Running (Only) A Single-Point Energy Calculation In Crystal06/09, Proper Input Format For Long-Range Dispersion Contributions In Crystal09, And Removing The MPICH2 Content From The Output File In Pcrystal

Now enjoying the benefits of dispersion-corrected solid-state density functional theory (and a proper MPICH2 implementation for infrared intensity calculations, although this now a problem for reasons to be addressed in an upcoming post) in Crystal09, three issues in recent calculations caused me to think hard enough about keyword formats and job runs that I have opted to post briefly about what to do in case google and bing are your preferred methods of manual searching.

1. How To Run Only A Single-Point Energy Calculation In Crystal06/Crystal09

This had never come up before and, by the time I needed to find an input file to see what do to, the first google search provided Civalleri's Total Energy Calculation page that currently has broken links to .zip files. There is quite a bit about the different geometry optimization approaches in the manual, but a search for "single-point" provides no information about what to do for only single-point energy calculations.

The solution, it should be obvious after, is simply to not include the geometry optimization section in the input file. What would otherwise be the following (with arbitrary geometry optimization-like info between [COORDINATES] and [BASIS SETS]…

[COORDINATES]
OPTGEOM
TOLDEG
0.000005
TOLDEX
0.000020
END
END
[BASIS SETS]

becomes…

[COORDINATES]
[BASIS SETS]

One problem solved by simply not having any optimization parameters (again, makes sense and is now google-able).

2. Proper GRIMME Input Format For Long-Range Dispersion Contributions In Crystal09

This is another example where one's first efforts in translating the manual into calculations may lead to considerable confusion until the proper format is finally identified (by which time you've run many pruned-down input tests).

GRIMME
1.05 20. 25.
1.05 20. 25. s6 (scaling factor) d (steepness) Rcut (cutoff radius)
5
1  0.14 1.001 Hydrogen Conventional Atomic number , C6 , Rvdw
6  1.75 1.452 Carbon Conventional Atomic number , C6 , Rvdw
7  1.23 1.397 Nitrogen Conventional Atomic number , C6 , Rvdw
8  0.70 1.342 Oxygen Conventional Atomic number , C6 , Rvdw
17 5.07 1.639 Chlorine Conventional Atomic number , C6 ,'Rvdw

I'm not even sure where the final ,'Rvdw comes from. Your .out file may terminate with the following error (or something similar)…

rank 7 in job 8  korterquad_51438   caused collective abort of all ranks
  exit status of rank 7: killed by signal 9

And the ERROR.peN file with any content will show the following, clearly pointing to a GRIMME-specific error…

 ERROR **** GRIMME_INPUT **** ELEMENT NOT DEFINED:           1

The problem is the additional content within the manual pages for the GRIMME keyword that require pruning (or, at least, some identifier to show what is and what is not needed). The proper GRIMME section above is properly provided in the INPUT file as…

GRIMME
1.05 20. 25.
5
1  0.14 1.001
6  1.75 1.452
7  1.23 1.397
8  0.70 1.342
17 5.07 1.639

Where (see page 88 of the Crystal09 manual)…

GRIMME <- keyword is called
1.05 20. 25. <- scaling factor, steepness, cutoff distance
5 <- number of elements in the list (not the total number of atoms)
1  0.14 1.001 <- atomic number, dispersion coefficient, van der Waals radius
...

When all is properly run, the bottom of your output file will look something like the following:

 CYC  43 ETOT(AU) -5.784662098123E+03 DETOT  1.18E-11 tst  8.17E-15 PX  6.73E-08

 == SCF ENDED - CONVERGENCE ON ENERGY      E(AU) -5.7846620981229E+03 CYCLES  43

 ENERGY EXPRESSION=HARTREE+FOCK EXCH*0.20000+(BECKE  EXCH)*0.80000+LYP    CORR

 TOTAL ENERGY(DFT)(AU)( 43) -5.7846620981229E+03 DE 1.2E-11 tester 8.2E-15
 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT EDFT        TELAPSE     4705.82 TCPU     4651.41

 *******************************************************************************

 GRIMME DISPERSION ENERGY CORRECTION

 SCALE FACTOR (S6):     1.0500

 GRIMME DISPERSION ENERGY (AU) -1.9723347118951E-01
 TOTAL ENERGY + DISP (AU) -5.7848593315941E+03

 *******************************************************************************

The Crystal09 manual refers you to Table 1 of the Stefan Grimme paper, "Semiempirical GGA-type density functional constructed with a long-range dispersion correction" (Journal of Computational Chemistry, Volume 27, Issue 15, Pages 1787 – 1799), which I've put together into the proper format below. Be sure to (1) delete the elements in parentheses ( -> get rid of the (H) <- ), (2) remove those atoms you do not need, (3) be sure to change the "number of elements" number for your structure, and (4) get and reference the Grimme paper so you have the proper C6 parameters and van der Waals radii accounted for (you'll be the right nitwit if I mis-copied something and you ran with it (although I trust my input), and you should have the reference regardless).

( H)   1   0.14 1.001
(Li)   3   1.61 0.825
(Na)  11   5.71 1.144
( K)  19  10.80 1.485
(Rb)  37  24.67 1.628
(Be)   4   1.61 1.408
(Mg)  12   5.71 1.364
(Ca)  20  10.80 1.474
(Sr)  38  24.67 1.606
( B)   5   3.13 1.485
(Al)  13  10.79 1.639
(Ga)  31  16.99 1.650
(In)  49  37.32 1.672
( C)   6   1.75 1.452
(Si)  14   9.23 1.716
(Ge)  32  17.10 1.727
(Sn)  50  38.71 1.804
( N)   7   1.23 1.397
( P)  15   7.84 1.705
(As)  33  16.37 1.760
(Sb)  51  38.44 1.881
( O)   8   0.70 1.342
( S)  16   5.57 1.683
(Se)  34  12.64 1.771
(Te)  52  31.74 1.892
( F)   9   0.75 1.287
(Cl)  17   5.07 1.639
(Br)  35  12.47 1.749
( I)  53  31.50 1.892
(He)   2   0.08 1.012
(Ne)  10   0.63 1.243
(Ar)  18   4.61 1.595
(Kr)  36  12.01 1.727
(Xe)  54  29.99 1.881
Y-Cd      24.67 1.639
Sc-Zn     10.80 1.562

Note that the d-block is identical for each row (so no atom numbers provided).

3. Removing The MPICH2 Content From The Output File In Pcrystal(/09)

This final issue does not occur in Pcrystal(/06) but does in Pcrystal(/09), with the reason being (I assume) the new use of MPICH2 in Pcrystal(/09) instead of MPICH in Pcrystal(/06).  The problem comes from running the following set of commands at the terminal window in MPICH2:

mpiexec -machinefile machine -np N /path/to/Pcrystal &>FILENAME.out &

Embedded within the FILENAME.out file will be all flavors of MPI-specific output, perhaps such as the following (in this case errors, but it happens in proper output as well):

application called MPI_Abort(MPI_COMM_WORLD, 1) - process 4
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 7
rank 7 in job 9  korterquad_51438   caused collective abort of all ranks
 exit status of rank 7: return code 1 
rank 4 in job 9  korterquad_51438   caused collective abort of all ranks
 exit status of rank 4: killed by signal 9 

or…

mpiexec_machine (handle_stdin_input 1089): stdin problem; if pgm is run in background...
mpiexec_machine (handle_stdin_input 1090):     e.g.: mpiexec -n 4 a.out < /dev/null &

The solution is to break up the mpiexec output from the Pcrystal output, performed by directing the mpiexec-specific content to, in this case, /dev/null (because it is not necessary except for diagnostic purposes).

mpiexec -machinefile machine -np N /path/to/Pcrystal < /dev/null &>FILENAME.out &

Which removes all traces of mpi-specific output from FILENAME.out.

Terahertz Spectroscopic Investigation Of S-(+)-Ketamine Hydrochloride And Vibrational Assignment By Density Functional Theory, "Function Follows Functional Follows Formalism"

Accepted in the Journal of Physical Chemistry A, with my fingers crossed for pulling off the rare double-header in an upcoming print edition of the journal (having missed it by three intermediate articles with the Cs2B12H12 and HMX papers back in 2006 (you'd keep track, too). A fortuitous overlap of scheduled defense dates between P. Hakey, Ph.D. and M. Hudson, A.B.D.). A brief summary of interesting points from this study is provided below, including what I think is a useful point about how to most easily interpret AND represent solid-state vibrational spectra for publications.

1. AS USUAL, YOU CANNOT USE GAS-PHASE CALCULATIONS TO ASSIGN SOLID-STATE TERAHERTZ SPECTRA. It will take a phenomenal piece of data and one helluvan interpretation to convince me otherwise. As a more subtle point (for those attempting an even worse job of vibrational mode assignment), if the molecule exists in its protonated form in the solid-state, do not use the neutral form for your gas-phase calculation (this is a point that came up as part of an MDMA re-assignment published (and posted here) previously).

2. It is very difficult to find what I would consider to be "complete data sets" for molecules and solids being studied by spectroscopic and computational methods. For many molecular solids, the influences of thermal motion are not important to providing a proper vibrational analysis by solid-state density functional theory methods. Heating a crystal may make spectral lines broader, but phase changes and unusual spectral features do not often result when heating a sample from cryogenic (say, liquid nitrogen) to room temperature. Yes, there are thousands of cases where this is not true, but several fold more cases where it is. We are fortunate to live in a temperature regime where characterization is reasonably straightforward and yet we can modify a system to observe its subtle changes under standard laboratory conditions. The THz spectrum of S-(+)-Ketamine Hydrochloride gets a bit cleaner upon cooling, which makes the assignment easier. As the ultimate goal is to be able to characterize these systems in a person's pocket instead of their liquid nitrogen thermos, the limited observed change to the spectrum upon cooling is important to note.

3. Crystal06 vs. DMol3 – This paper contains what is hoped to be a level, pragmatic discussion about the strengths and weaknesses of computational tools available to terahertz spectroscopists for use in their efforts to assign spectra. This type of discussion is, as a computational chemist using tools and not developing tools, a touchy subject to present on not because of the finger-pointing of limitations with software, but because the Crystal06 team and Accelrys (through Delley's initial DMol3 code) clearly are doing things that the vast majority of their users (myself included) could in no way do by themselves. The analysis for the theory-minded terahertz spectroscopist is presented comparing two metrics – speed and functionality (specifically, infra-red intensity prediction). What is observed as the baseline is that both DMol3 and Crystal06 make available density functionals and basis sets that, when used at high levels of theory and rigorous convergence criteria, produce simulated terahertz spectra with vibrational mode energies that are in good (if not very good) agreement with each other. For the terahertz spectroscopist, Crystal06 provides as output (although this is system size- and basis set size-dependent) rigorous infrared intensity predictions for vibrational modes, inseparable from mode energy as "the most important" pieces of information for mode assignments. While DMol3 does not produce infrared intensities (the many previous terahertz papers I've worked on employed difference-dipole calculations that are, at best, a guesstimate), DMol3 produces very good mode energy predictions in 1/6th to (I've seen it happen) 1/10th the time of a comparable Crystal06 calculation. This is the reason DMol3 has been the go-to program for all of the neutron scattering spectroscopy papers cited on this blog (where intensity is determined by normal mode eigenvectors, which are provided by both (and any self-respecting quantum chemical code) programs).

Now, it should be noted that this difference in functionality has NOTHING to do with formalism. Both codes are excellent for what they are intended to do. To the general assignment-minded spectroscopist (the target audience of the Discussion in the paper), any major problem with Crystal06 likely originates with the time to run calculations (and, quite frankly, the time it takes to run a calculation is the worst possible reason for not running a calculation if you need that data. Don't blame the theory, blame the deadline). In my past exchanges with George Fitzgerald of Accelrys, the issue of DMol3 infrared intensities came up as a feature request that would greatly improve the (this) user experience and Dr. Fitzgerald is very interested (of course) in making a great code that much better. Neither code will be disappearing from my toolbox anytime soon.

4. The Periodicity Of The Molecular Solid Doesn't Care What The Space Group Is – One of the more significant problems facing the assignment-minded spectroscopist is the physical description of molecular motion in a vibrational mode. In the simplest motions involving the most weakly interacting molecules, translational and rotational motions are often quite easy to pick out and state as such. When the molecules are very weakly interacting, often the intramolecular vibrational modes are easy to identify as well, as they are largely unchanged from their gas-phase descriptions. In ionic solids or strongly hydrogen-bonded systems, it is often much harder to separate out individual molecular motions from "group modes" involving the in- and out-of-phase motions of multiple molecules. In the unit cells of molecular solids, it can be the case that these group modes appear, by inspection, to be extremely complicated, sometimes too involved to easily describe in the confines of a table in a journal article.

S-(+)-Ketamine Hydrochloride is one such example where a great simplification in vibrational mode description comes from thinking, well, "outside the box." The image below shows two cells and the surrounding molecules of S-(+)-Ketamine Hydrochloride. As it is difficult to see why the mode descriptions are complex from just an image, assume that I am right in this statement of complexity. Part of this complexity comes from the fact that the two molecules in the unit cell are not strongly interacting, instead packed together by van der Waals and dispersion forces more than anything else. The key to a greatly simplified assignment comes from the realization that the most polar fragments of these molecules are aligned on the edges of the unit cell.

An alternate view of molecular vibrational motion comes from considering not the contents of the defined unit cell but the hydrogen-bonding and ionic bonding arrangement that exists between pairs of molecules between unit cells. The colorized image below shows two distinct chains (red and blue) that, when the predicted vibrational modes are animated, become trivial to characterize as the relative motions of a hydrogen/ionic-bonded chain. Rotational motions appear as spinning motions of the chains, translational motions as either chain sliding motions or chain breathing modes. It appears as a larger macromolecule undergoing very "molecular" vibrations. In optical vibrational spectroscopy, selection rules and the unit cell arrangement do not produce in- and out-of-phase motions of the red and blue chains, as only one "chain" exists in the periodicity of the unit cell. In neutron scattering spectroscopy, these relative motions between red and blue would appear in the phonon region. This same discussion was had, in part, in a previous post on the solid-state terahertz assignment of ephedrine (with a nicer picture).

So, look at the cell contents, then see if there's more structure than crystal packing would indicate. It greatly simplifies the assignment (which, in turn. greatly simplifies the reader's digestion of the vibrational motions).

Patrick M. Hakey, Damian G. Allis, Matthew R. Hudson, Wayne Ouellette, and Timothy M. Korter

Department of Chemistry, Syracuse University, Syracuse, New York 13244-4100

Abstract: The terahertz (THz) spectrum of (S)-(+)-ketamine hydrochloride has been investigated from 10 to 100 cm-1 (0.3-3.0 THz) at both liquid-nitrogen (78 K) and room (294 K) temperatures. Complete solid-state density functional theory structural analyses and normal-mode analyses are performed using a single hybrid density functional (B3LYP) and three generalized gradient approximation density functionals (BLYP, PBE, PW91). An assignment of the eight features present in the well-resolved cryogenic spectrum is provided based upon solid-state predictions at a PW91/6-31G(d,p) level of theory. The simulations predict that a total of 13 infrared- active vibrational modes contribute to the THz spectrum with 26.4% of the spectral intensity originating from external lattice vibrations.

pubs.acs.org/journal/jpcafh
www.somewhereville.com/?p=29
www.somewhereville.com/?p=26
www.somewhereville.com/?p=126
en.wikipedia.org/wiki/Density_functional_theory
en.wikipedia.org/wiki/Ketamine
www.crystal.unito.it
accelrys.com/products/materials-studio/quantum-and-catalysis-software.html
en.wikipedia.org/wiki/Time_domain_terahertz_spectroscopy
en.wikipedia.org/wiki/Computational_chemistry
accelrys.com
en.wikipedia.org/wiki/Inelastic_neutron_scattering
en.wikipedia.org/wiki/Vibrational_spectroscopy
www.somewhereville.com/?p=680