Blog Archives

Geometry Optimizations for Excited States


Electronic excitations are calculated vertically according to the Frank—Condon principle, this means that the geometry does not change upon the excitation and we merely calculate the energy required to reach the next electronic state. But for some instances, say calculating not only the absorption spectra but also the emission, it is important to know what the geometry minimum of this final state looks like, or if it even exists at all (Figure 1). Optimizing the geometry of a given excited state requires the prior calculation of the vertical excitations whether via a multireference method, quantum Monte Carlo, or the Time Dependent Density Functional Theory, TD-DFT, which due to its lower computational cost is the most widespread method.

Most single-reference treatments, ab initio or density based, yield good agreement with experiments for lower states, but not so much for the higher excitations or process that involve the excitation of two electrons. Of course, an appropriate selection of the method ensures the accuracy of the obtained results, and the more states are considered, the better their description although it becomes more computationally demanding in turn.

Figure 1. The vertical excitation does not match the minimum on the excited state

In Gaussian 09 and 16, the argument to the ROOT keyword selects a given excited state to be optimized. In the following example, five excited states are calculated and the optimization is requested upon the second excited state. If no ROOT is specified, then the optimization would be carried out by default on the first excited state (Where L.O.T. stands for Level of Theory).

#p opt TD=(nstates=5,root=2) L.O.T.

Gaussian16 includes now the calculation of analytic second derivatives which allows for the calculation of vibrational frequencies for IR and Raman spectra, as well as transition state optimization and IRC calculations in excited states opening thus an entire avenue for the computation of photochemistry.

If you already computed the excited states and just want to optimize one of them from a previous calculation, you can read the previous results with the following input :

#p opt TD=(Read,Root=N) L.O.T. Density=Current Guess=Read Geom=AllCheck

Common problems. The following error message is commonly observed in excited state calculations whether in TD-DFT, CIS or other methods:

No map to state XX, you need to solve for more vectors in order to follow this state.

This message usually means you need to increase the number of excited states to be calculated for a proper description of the one you’re interested in. Increase the number N for nstates=N in the route section at higher computational cost. A rule of thumb is to request at least 2 more states than the state of interest. This message can also reflect the fact that during the optimization the energy ordering changes between states, and can also mean that the ground state wave function is unstable, i.e., the energy of the excited state becomes lower than that of the ground state, in this case a single determinant approach is unviable and CAS should be used if the size of the molecule allows it. Excited state optimizations are tricky this way, in some cases the optimization may cross from one PES to another making it hard to know if the resulting geometry corresponds to the state of interest or another. Gaussian recommends changing the step size of the optimization from the default 0.3 Bohr radius to 0.1, but obviously this will make the calculation take longer.

Opt=(MaxStep=10)

If the minimum on the excited state potential energy surface (PES) doesn’t exist, then the excited state is not bound; take for example the first excited state of the H2 molecule which doesn’t show a minimum, and therefore the optimized geometry would correspond to both H atoms moving away from each other indefinitely (Figure 2). Nevertheless, a failed optimization doesn’t necessarily means the minimum does not exist and further analysis is required, for instance, checking the gradient is converging to zero while the forces do not.

Figure 2. An unbound excited state with no minima ensures the dissociation of the system along the reaction coordinate

Grimme’s Dispersion DFT-D3 in Gaussian #CompChem


I was just asked if it is possible to perform DFT-D3 calculations in Gaussian and my first answer was to use the following  keyword:

EmpiricalDispersion=GD3

which is available in G16 and G09 only in revision D, apparently.

There are also some overlays that can be used to invoke the use dispersion in various scenarios:

IOp(3/74=x) Exchange and Correlation Potentials

-77

-76

-60

-59

DSD-PBEP86 (double hybrid, DFT-D3).

PW6B95-D3.

B2PLYP-D3 (double hybrid, DFT-D3).

B97-D (DFT-D3).

IOp(3/76=x) Mixing of HF and DFT.

-33 PW6B95 and PW6B95-D3 coefficients.

IOp(3/124=x) Empirical dispersion term.

30

40

50

Force dispersion type 3 (Grimme DFT-D3).

Force dispersion type 4 (Grimme DFT-D3(BJ)).

Force dispersion type 5 (Grimme D3, PM7 version).

 

The D3 correction method of Grimme defines the van der Waals energy like:

$\displaystyle E_{\rm disp} = -\frac{1}{2} \sum_{i=1}^{N_{at}} \sum_{j=1}^{N_{at...
...{6ij}} {r_{ij,{L}}^6} +f_{d,8}(r_{ij,L})\,\frac{C_{8ij}} {r_{ij,L}^8} \right ),$

where coefficients $ C_{6ij}$ are adjusted depending on the geometry of atoms i and j. The damping D3 function for is:

$\displaystyle f_{d,n}(r_{ij}) = \frac{s_n}{1+6(r_{ij}/(s_{R,n}R_{0ij}))^{-\alpha_{n}}},$

where the values of s are adjustable parameters fit for the exchange-correlation functionals used in each calculation.

Partial Optimizations with Gaussian09


Sometimes you just need to optimize some fragment or moiety of your molecule for a number of reasons -whether because of its size, your current interest, or to skew the progress of a previous optimization- or maybe you want just some kind of atoms to have their positions optimized. I usually optimize hydrogen atoms when working with crystallographic files but that for some reason I want to preserve the rest of the molecule as refined, in order to keep it under a crystalline field of sorts.
Asking Gaussian to optimize some of the atoms in your molecule requires you to make a list albeit the logic behind it is not quite straightforward to me. This list is invoked by the ReadOptimize keyword in the route section and it includes all atoms by default, you can then further tell G09 which atoms are to be included or excluded from the optimization.
So, for example you want to optimize all atoms EXCEPT hydrogens, then your input should bear the ReadOptimize keyword in the route section and then, at the end of the molecule specification, the following line:

atoms notatoms=H

If you wish to selectively add some atoms to the list while excluding others, here’s an example:

atoms=C H S notatoms=5-8

This list adds, and therefore optimizes, all carbon, hydrogen and sulfur atoms, except atoms 5, 6, 7 and 8, should they be any of the previous elements in the C H S list.
The way I selectively optimize hydrogen atoms is by erasing all atoms from the list -using the noatoms instruction- and then selecting which are to be included in the list -with atoms=H-, but I haven’t tried it with only selecting hydrogen atoms from the start, as in atoms=H

noatoms atoms=H

I probably get very confused because I learned to do this with the now obsolete ReadFreeze keyword; now it sometimes may seem to me like I’m using double negatives or something – please do not optimize all atoms except if they are hydrogen atoms. You can include numbers, ranks or symbols in this list as a final line of your input file.

Common errors (by common I mean I’ve got them):

Lets look at the end of an input I just was working with:

>  AtmSel:  Line=”P  0″
>  Maximum list size exceeded in AddBin.
>  Error termination via Lnk1e in…

AtmSel is the routine which reads the atoms list and I was using a pseudopotential on phosphorous atoms, I placed the atoms list at the end of the file but it should be placed right after the coordinates and the connectivity matrix, should there be one, and thus before any external basis set or pseudopotential or any other specification to be read by Gaussian.

As a sort of test you can use the instruction:

%kjob l103
%chk=myfile.chk
...

at the Link0 section (where your checkpoint is defined). This will kill the job after the link 103 is finished, thus you will only get a list of what parameters were frozen and which were active. Then, if things look ok, you can run the job without the %kjob l103 instruction and get it done.

As usual I hope this helps. Thanks for reading except to those who didn’t read it except for the parts they did read.

Transition State Search (QST2 & QST3) and IRC with Gaussian09


Theoretical evaluation of a reaction mechanism is all about finding the right transition states (TS) but there are no guarantees within the available methods to actually find the one we need. Chemical intuition in the proposal of a mechanism is paramount. Let’s remember that a TS is a critical point on a Potential Energy Surface (PES) that is a minimum in every dimension but one. For a PES with more than two degrees of freedom, a hyper-surface, envisioning the location of a TS is a bit tricky, in the case of a three dimensional PES (two degrees of freedom) the saddle point constitutes the location of the TS as depicted in figure 1 by a section of a revolution hyperboloid.

400px-Saddle_point

Fig1. Saddle point on a surface (min in one direction; max in the other)

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

The following procedure considers gas phase calculations. Nevertheless, the use of the SCRF keyword activates the implicit solvent calculation of choice in order to evaluate to some degree the solvent influence on the reaction energetics at different temperatures with the use of the temperature keyword.

The first step consists of a high level optimization of all minimums involved, such as reagents, products and intermediates, with a subsequent frequency analysis that includes no imaginary eigenvalues.

In order to find the structures of the transition states we use in Gaussian the Synchronous Transit-guided Quasi-Newton method [1] through the keywords QST2 or QST3. In the former case, coordinates for the reagents and products are needed as input; for the latter keyword, coordinates for the TS structure guess is needed also.

QST2)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst2,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line

QST3)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst3,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title Card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line–
Title card for TS
Q M
Cartesian Coordinates for TS
blank line

NOTE: It is fundamental that the numbering order is kept constant throughout the molecular specifications of all two, or three, input structures. Hence, I recommend to build a set of molecules, save their structure, and then modified the coordinates on the same file to produce the following structure, that way the number for every atom will remain the same for every step.

As I wrote above, there are no guarantees of finding the right TS so many attempts are probably needed. Once we have the optimized structures for all the species involved in our mechanistic proposal we can plot their energies very simply with MS Excel the way we’ve previously described in this blog (reblogged from eutactic.wordpress.com)

Once we’ve succeeded in finding the structure of our TS we may run an Internal Reaction Coordinate (IRC) calculation. This calculation will connect the TS structure to those of the products and the reagents. Initial constant forces are required and these are commonly retrieved from the TS calculation checkpoint file through the RCFC keyword.

%chk=QST3_2p.chk
%nprocshared=8

#p m062x/6-31++G(d,p) IRC=(Maxpoints=50,RCFC,phase=(2,1))Temperature=373.15 SCRF=(Solvent=Water) geom=allcheck

Title Card

Q M
blank line

Finally, the IRC path can be visualized with GaussView from the Results menu. A successful IRC will link both structures along a single reaction coordinate proving that both reagents and products are linked by the obtained TS.

Hat tip to Howard Diaz who has become quite skillful in calculating these mechanisms as proven by his recent poster at the XII RMFQT a couple of weeks back. And as usual thanks to everyone who reads, comments, likes, recommends, rates and shares my silly little posts.

Delta G of solvation in Gaussian09


How to calculate the Delta G of solvation? This is a question that I get a lot in this blog, so it is about time I wrote a (mini)post on it, and at the same time put an end to this posting drought which has lasted for quite a few months due to a lot of pending work with which I’ve had to catch up. Therefore, this is another post in the series of SCRF calculations that are so popular in this blog. For the other posts on this subjects remember to click here and here.

SMD

SMD is the keyword you want to use when performing a Self Consistent Reaction Field (SCRF) calculation with G09. This keyword was only made available in this last version of the program and it corresponds to Truhlar’s and coworkers solvation model which is recommended by Gaussian itself as the preferred model to calculate Delta G of solvation. The syntax used is the standard way used in any other Gaussian input files as follows:

# 'route section keywords' SCRF=SMD

Separately, we must either perform a gas phase calculation or use the DoVacuum keyword within the same SCRF input, and then take the energy difference between gas phase and solvated models.

# 'route section keywords' SCRF=(SMD,DoVacuum)

No solvation or cavity model should be defined since, by definition, SMD will use the IEFPCM model which is a synonym for PCM.

As opposed to the previous versions of Gaussian, the output energy already contains all corrections, this is why we must take the difference between both values (remember to calculate them both at the same level of theory if calculated separately!). Nevertheless, when using the SMD keyword we get a separate line, just below the energy, stating the SMD-CDS non electrostatic value in kCal/mol.

The radii were also defined in the original paper by Truhlar; I’m not sure if using the keyword RADII with any of its options yields a different result or if it even ends in an error. Its worth the try!

Some calculation variations are not available when using SMD, such as Dis (calculation of the solute-solvent dispersion interaction energy), Rep (solute-solvent repulsion interaction energy) and Cav (inclusion of the solute cavitation energy in the total energy). I guess the reason for this might be that the SMD model is highly parametrized.

Have you found any issue with any item listed above? Pleases share your thoughts in the comments section below. As usual I hope this post was useful and that you all rate it, like it and comment.

References

A. V. Marenich, C. J. Cramer, and D. G. Truhlar, “Universal solvation model based on solute electron density and a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions,” J. Phys. Chem. B, 113 (2009) 6378-96.
%d bloggers like this: