The generation of isotopically-substituted molecular crystal spectra has become a point of interest, which means blog post. To be clear, this is for cases where isotopic substitution does not affect the crystal geometry – the crystal cell does not change significantly upon deuteration (and for those who believe isotopic substitution never leads to significant changes in the solid, I refer you Zhou, Kye, and Harbison's article on Isotopomeric Polymprphism and their work on 4-methylpyridine pentachlorophenol, which changes dramatically upon deuteration. I beat on this point because blindly assuming of the crystal cell geometry in such cases will produce spectra noticeably different than measured. It's NOT the calculation's fault!).
The generation of isotopically-substituted spectra and intensities in Crystal09 is trivial provided that you KEEP THE FREQINFO.DAT FILE. In fact, you need keep ONLY the FREQINFO.DAT to generate these spectra, which greatly reduces file transfer loads and allows for the scripted calculation of new vibrational spectra and thermodynamic data post-frequency calculation.
As my example system, I'm using the dispersion-corrected crystal cell of alpha-HMX (I have it handy, it's a small system, and having anything about HMX on your website is proven to increase traffic) at the B3LYP/6-31G(d,p) level of theory. Original input file (the one where the original normal mode analysis is performed) is below:
Test - alpha-HMX 6-31Gdp set DFT/B3LYP FREQ
CRYSTAL
0 0 0
43
15.14 23.89 5.913 124.3
14
6 1.016493675797E-01 -4.109909899348E-02 -3.351438244488E-03
6 -6.539109813231E-02 -6.180633576707E-02 -1.110575784790E-02
1 9.149797846691E-02 -4.382919469310E-02 -1.860042940246E-01
1 1.558888705857E-01 -6.829708099502E-02 4.595161229829E-02
1 -5.138242817334E-02 -5.844587273099E-02 -1.920922064181E-01
1 -9.781600273101E-02 -1.015710562102E-01 2.063738273292E-02
7 1.992579327285E-02 -5.951921578598E-02 1.040704228546E-01
7 1.232154652110E-01 1.634305404407E-02 5.951841980010E-02
7 2.220759010770E-02 -7.142100857312E-02 3.299259852838E-01
7 2.054067942916E-01 2.817244373261E-02 1.473285310628E-01
8 -4.761487685316E-02 -8.656669456613E-02 4.192568497756E-01
8 9.327421157186E-02 -6.479426971916E-02 4.286363161888E-01
8 2.563441491059E-01 -1.128705054032E-02 1.760581823035E-01
8 2.225071782791E-01 7.736574474011E-02 1.903699942346E-01
FREQCALC
INTENS
END
END
8 4
0 0 6 2.0 1.0
5484.671700 0.1831100000E-02
825.2349500 0.1395010000E-01
188.0469600 0.6844510000E-01
52.96450000 0.2327143000
16.89757000 0.4701930000
5.799635300 0.3585209000
0 1 3 6.0 1.0
15.53961600 -0.1107775000 0.7087430000E-01
3.599933600 -0.1480263000 0.3397528000
1.013761800 1.130767000 0.7271586000
0 1 1 0.0 1.0
0.2700058000 1.000000000 1.000000000
0 3 1 0.0 1.0
0.800000000 1.00000000
7 4
0 0 6 2.0 1.0
4173.51100 0.183480000E-02
627.457900 0.139950000E-01
142.902100 0.685870000E-01
40.2343300 0.232241000
12.8202100 0.469070000
4.39043700 0.360455000
0 1 3 5.0 1.0
11.6263580 -0.114961000 0.675800000E-01
2.71628000 -0.169118000 0.323907000
0.772218000 1.14585200 0.740895000
0 1 1 0.0 1.0
0.212031300 1.00000000 1.00000000
0 3 1 0.0 1.0
0.800000000 1.00000000
6 4
0 0 6 2.0 1.0
.3047524880D+04 .1834737130D-02
.4573695180D+03 .1403732280D-01
.1039486850D+03 .6884262220D-01
.2921015530D+02 .2321844430D+00
.9286662960D+01 .4679413480D+00
.3163926960D+01 .3623119850D+00
0 1 3 4.0 1.0
.7868272350D+01 -.1193324200D+00 .6899906660D-01
.1881288540D+01 -.1608541520D+00 .3164239610D+00
.5442492580D+00 .1143456440D+01 .7443082910D+00
0 1 1 0.0 1.0
.1687144782D+00 .1000000000D+01 .1000000000D+01
0 3 1 0.0 1.0
.8000000000D+00 .1000000000D+01
1 3
0 0 3 1.0 1.0
.1873113696D+02 .3349460434D-01
.2825394365D+01 .2347269535D+00
.6401216923D+00 .8137573262D+00
0 0 1 0.0 1.0
.1612777588D+00 .1000000000D+01
0 2 1 0.0 1.0
.1100000000D+01 .1000000000D+01
99 0
END
DFT
B3LYP
XLGRID
END
EXCHSIZE
10654700
BIPOSIZE
10654700
TOLINTEG
8 8 8 8 16
SCFDIR
MAXCYCLE
100
TOLDEE
11
GRIMME
1.05 20. 25.
4
1 0.14 1.001
6 1.75 1.452
7 1.23 1.397
8 0.70 1.342
SHRINK
8 8
LEVSHIFT
5 0
FMIXING
50
END
END
Upon completion of this run, you need only the FREQINFO.DAT file, the last set of coordinates from the .OUT file (for atom counting purposes) and an input file which is modified from the original only in the specification of the ISOTOPES section and which includes a RESTART.
Question – how does one deal with isotopically-labeling atoms when it breaks the space group symmetry? If I isotopically label Atom 1 in the asymmetric unit, what happens to the other N symmetry-related atoms?
Answer – Crystal09, in its infinite wisdom, does not consider the asymmetric unit in the isotopic substitution scheme. If you've 14 atoms in the asymmetric unit (the symmetry-unique atoms you provide in the input file)…
14 6 1.016493675797E-01 -4.109909899348E-02 -3.351438244488E-03 6 -6.539109813231E-02 -6.180633576707E-02 -1.110575784790E-02 ... 8 2.563441491059E-01 -1.128705054032E-02 1.760581823035E-01 8 2.225071782791E-01 7.736574474011E-02 1.903699942346E-01
and 56 atoms in the full unit cell…
ATOMS IN THE ASYMMETRIC UNIT 14 - ATOMS IN THE UNIT CELL: 56
ATOM X/A Y/B Z/C
*******************************************************************************
1 T 6 C -1.460999048177E-01 1.393970283287E-01 6.390170683069E-02
2 F 6 C 1.393970283287E-01 -1.460999048177E-01 -5.719883034171E-02
3 F 6 C 3.071988303417E-01 1.860982931693E-01 1.106029716713E-01
4 F 6 C 1.860982931693E-01 3.071988303417E-01 3.960999048177E-01
...
53 T 8 O 4.522856069554E-02 3.355114277736E-01 1.095029287847E-01
54 F 8 O 3.355114277736E-01 4.522856069554E-02 -4.902429172538E-01
55 F 8 O -2.597570827462E-01 1.404970712153E-01 -8.551142777356E-02
56 F 8 O 1.404970712153E-01 -2.597570827462E-01 2.047714393045E-01
your ISOTOPES section relies on the numbering of the atoms in the "56 atom" list.
The input file below will calculate an isotopically-labeled vibrational spectrum for 8 of the hydrogen atoms that ends up breaking the unit cell symmetry (which will be more obvious from the produced mode energies). Again, the atom numbers come from the "ATOMS IN THE ASYMMETRIC UNIT" part of the original optimization by which you performed the original normal mode analysis (hopefully).
Test - alpha-HMX 6-31Gdp set DFT/B3LYP FREQ - Isotopic Substitution
CRYSTAL
0 0 0
43
15.14 23.89 5.913 124.3
14
6 1.016493675797E-01 -4.109909899348E-02 -3.351438244488E-03
6 -6.539109813231E-02 -6.180633576707E-02 -1.110575784790E-02
1 9.149797846691E-02 -4.382919469310E-02 -1.860042940246E-01
1 1.558888705857E-01 -6.829708099502E-02 4.595161229829E-02
1 -5.138242817334E-02 -5.844587273099E-02 -1.920922064181E-01
1 -9.781600273101E-02 -1.015710562102E-01 2.063738273292E-02
7 1.992579327285E-02 -5.951921578598E-02 1.040704228546E-01
7 1.232154652110E-01 1.634305404407E-02 5.951841980010E-02
7 2.220759010770E-02 -7.142100857312E-02 3.299259852838E-01
7 2.054067942916E-01 2.817244373261E-02 1.473285310628E-01
8 -4.761487685316E-02 -8.656669456613E-02 4.192568497756E-01
8 9.327421157186E-02 -6.479426971916E-02 4.286363161888E-01
8 2.563441491059E-01 -1.128705054032E-02 1.760581823035E-01
8 2.225071782791E-01 7.736574474011E-02 1.903699942346E-01
FREQCALC
RESTART
ISOTOPES
8
9 2
10 2
11 2
13 2
14 2
15 2
16 2
18 2
INTENS
END
END
8 4
0 0 6 2.0 1.0
5484.671700 0.1831100000E-02
825.2349500 0.1395010000E-01
188.0469600 0.6844510000E-01
52.96450000 0.2327143000
16.89757000 0.4701930000
5.799635300 0.3585209000
0 1 3 6.0 1.0
15.53961600 -0.1107775000 0.7087430000E-01
3.599933600 -0.1480263000 0.3397528000
1.013761800 1.130767000 0.7271586000
0 1 1 0.0 1.0
0.2700058000 1.000000000 1.000000000
0 3 1 0.0 1.0
0.800000000 1.00000000
7 4
0 0 6 2.0 1.0
4173.51100 0.183480000E-02
627.457900 0.139950000E-01
142.902100 0.685870000E-01
40.2343300 0.232241000
12.8202100 0.469070000
4.39043700 0.360455000
0 1 3 5.0 1.0
11.6263580 -0.114961000 0.675800000E-01
2.71628000 -0.169118000 0.323907000
0.772218000 1.14585200 0.740895000
0 1 1 0.0 1.0
0.212031300 1.00000000 1.00000000
0 3 1 0.0 1.0
0.800000000 1.00000000
6 4
0 0 6 2.0 1.0
.3047524880D+04 .1834737130D-02
.4573695180D+03 .1403732280D-01
.1039486850D+03 .6884262220D-01
.2921015530D+02 .2321844430D+00
.9286662960D+01 .4679413480D+00
.3163926960D+01 .3623119850D+00
0 1 3 4.0 1.0
.7868272350D+01 -.1193324200D+00 .6899906660D-01
.1881288540D+01 -.1608541520D+00 .3164239610D+00
.5442492580D+00 .1143456440D+01 .7443082910D+00
0 1 1 0.0 1.0
.1687144782D+00 .1000000000D+01 .1000000000D+01
0 3 1 0.0 1.0
.8000000000D+00 .1000000000D+01
1 3
0 0 3 1.0 1.0
.1873113696D+02 .3349460434D-01
.2825394365D+01 .2347269535D+00
.6401216923D+00 .8137573262D+00
0 0 1 0.0 1.0
.1612777588D+00 .1000000000D+01
0 2 1 0.0 1.0
.1100000000D+01 .1000000000D+01
99 0
END
DFT
B3LYP
XLGRID
END
EXCHSIZE
10654700
BIPOSIZE
10654700
TOLINTEG
8 8 8 8 16
SCFDIR
MAXCYCLE
100
TOLDEE
11
GRIMME
1.05 20. 25.
4
1 0.14 1.001
6 1.75 1.452
7 1.23 1.397
8 0.70 1.342
SHRINK
8 8
LEVSHIFT
5 0
FMIXING
50
END
END
The difference is in the FREQCALC section, which calls RESTART (to use the FREQINFO.DAT file), ISOTOPES (obvious), the total number of atoms that are having their isotopes changed (8), then the list, containing the atom number and the new mass (here, 2 for deuterium).
The proof is in the high-frequency region, where the last 16 modes (H-atom motion) in the non-deuterated form…
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
MODES EIGV FREQUENCIES IRREP IR INTENS RAMAN
(HARTREE**2) (CM**-1) (THZ) (KM/MOL)
...
153- 153 0.2003E-03 3106.1384 93.1197 (A2 ) I ( 0.00) A
154- 154 0.2003E-03 3106.5054 93.1307 (B1 ) A ( 0.02) A
155- 155 0.2004E-03 3106.5586 93.1323 (A1 ) A ( 0.23) A
156- 156 0.2004E-03 3106.8420 93.1408 (B2 ) A ( 0.48) A
157- 157 0.2017E-03 3117.1664 93.4503 (B2 ) A ( 1.13) A
158- 158 0.2018E-03 3117.4901 93.4600 (B1 ) A ( 2.33) A
159- 159 0.2021E-03 3120.2876 93.5439 (A1 ) A ( 115.24) A
160- 160 0.2022E-03 3120.7805 93.5586 (A2 ) I ( 0.00) A
161- 161 0.2131E-03 3203.6552 96.0432 (A1 ) A ( 44.59) A
162- 162 0.2131E-03 3203.6581 96.0433 (B2 ) A ( 115.98) A
163- 163 0.2132E-03 3204.6505 96.0730 (B1 ) A ( 15.30) A
164- 164 0.2132E-03 3204.8874 96.0801 (A2 ) I ( 0.00) A
165- 165 0.2157E-03 3223.4669 96.6371 (A1 ) A ( 44.98) A
166- 166 0.2157E-03 3223.5803 96.6405 (B2 ) A ( 27.02) A
167- 167 0.2158E-03 3223.8536 96.6487 (B1 ) A ( 35.26) A
168- 168 0.2158E-03 3224.3355 96.6631 (A2 ) I ( 0.00) A
change to the following last 16 modes (H/D-atom motion) upon deuteration. Note the mode energies split and the mode symmetries go from (A1,A2,B1,B2) to (A). Also note your IR mode intensities change, giving you the complete picture upon isotopic substitution.
HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH
MODES EIGV FREQUENCIES IRREP IR INTENS RAMAN
(HARTREE**2) (CM**-1) (THZ) (KM/MOL)
...
153- 153 0.1074E-03 2274.8942 68.1996 (A ) A ( 1.07) A
154- 154 0.1075E-03 2275.5949 68.2206 (A ) A ( 3.75) A
155- 155 0.1075E-03 2275.7008 68.2238 (A ) A ( 2.93) A
156- 156 0.1099E-03 2300.7446 68.9746 (A ) A ( 4.68) A
157- 157 0.1148E-03 2351.7846 70.5047 (A ) A ( 11.32) A
158- 158 0.1183E-03 2387.0269 71.5613 (A ) A ( 36.17) A
159- 159 0.1183E-03 2387.2610 71.5683 (A ) A ( 16.04) A
160- 160 0.1184E-03 2387.6687 71.5805 (A ) A ( 3.73) A
161- 161 0.2006E-03 3108.6223 93.1942 (A ) A ( 0.93) A
162- 162 0.2009E-03 3110.5061 93.2506 (A ) A ( 12.43) A
163- 163 0.2009E-03 3110.7567 93.2581 (A ) A ( 13.67) A
164- 164 0.2039E-03 3134.0133 93.9554 (A ) A ( 40.48) A
165- 165 0.2147E-03 3215.5160 96.3987 (A ) A ( 19.38) A
166- 166 0.2157E-03 3223.4291 96.6360 (A ) A ( 35.29) A
167- 167 0.2157E-03 3223.5925 96.6409 (A ) A ( 29.50) A
168- 168 0.2158E-03 3223.8729 96.6493 (A ) A ( 8.37) A