# Blog Archives

## Percentage of Molecular Orbital Composition – G09,G16

Canonical Molecular Orbitals are–by construction–delocalized over the various atoms making up a molecule. In some contexts it is important to know how much of any given orbital is made up by a particular atom or group of atoms, and while you could calculate it by hand given the coefficients of each MO in terms of every AO (or basis set function) centered on each atom there is a straightforward way to do it in Gaussian.

If we’re talking about ‘dividing’ a molecular orbital into atomic components, we’re most definitely talking about population analysis calculations, so we’ll resort to the pop keyword and the orbitals option in the standard syntax:

#p M052x/cc-pVDZ pop=orbitals

This will produce the following output right after the Mulliken population analysis section:

Atomic contributions to Alpha molecular orbitals:
Alpha occ 140 OE=-0.314 is Pt1-d=0.23 C38-p=0.16 C31-p=0.16 C36-p=0.16 C33-p=0.15
Alpha occ 141 OE=-0.313 is Pt1-d=0.41
Alpha occ 142 OE=-0.308 is Cl2-p=0.25
Alpha occ 143 OE=-0.302 is Cl2-p=0.72 Pt1-d=0.18
Alpha occ 144 OE=-0.299 is Cl2-p=0.11
Alpha occ 145 OE=-0.298 is C65-p=0.11 C58-p=0.11 C35-p=0.11 C30-p=0.11
Alpha occ 146 OE=-0.293 is C58-p=0.10
Alpha occ 147 OE=-0.291 is C22-p=0.09
Alpha occ 148 OE=-0.273 is Pt1-d=0.18 C11-p=0.12 C7-p=0.11
Alpha occ 149 OE=-0.273 is Pt1-d=0.18
Alpha vir 150 OE=-0.042 is C9-p=0.18 C13-p=0.18
Alpha vir 151 OE=-0.028 is C7-p=0.25 C16-p=0.11 C44-p=0.11
Alpha vir 152 OE=0.017 is Pt1-p=0.10
Alpha vir 153 OE=0.021 is C36-p=0.15 C31-p=0.14 C63-p=0.12 C59-p=0.12 C38-p=0.11 C33-p=0.11
Alpha vir 154 OE=0.023 is C36-p=0.13 C31-p=0.13 C63-p=0.11 C59-p=0.11
Alpha vir 155 OE=0.027 is C65-p=0.11 C58-p=0.10
Alpha vir 156 OE=0.029 is C35-p=0.14 C30-p=0.14 C65-p=0.12 C58-p=0.11
Alpha vir 157 OE=0.032 is C52-p=0.09
Alpha vir 158 OE=0.040 is C50-p=0.14 C22-p=0.13 C45-p=0.12 C17-p=0.11
Alpha vir 159 OE=0.044 is C20-p=0.15 C48-p=0.14 C26-p=0.12 C54-p=0.11


Alpha and Beta densities are listed separately only in unrestricted calculations, otherwise only the first is printed. Each orbital is listed sequentially (occ = occupied; vir = virtual) with their energy value (OE = orbital energy) in atomic units following and then the fraction with which each atom contributes to each MO.

By default only the ten highest occupied orbitals and ten lowest virtual orbitals will be assessed, but the number of MOs to be analyzed can be modified with orbitals=N, if you want to have all orbitals analyzed then use the option AllOrbitals instead of just orbitals. Also, the threshold used for printing the composition is set to 10% but it can be modified with the option ThreshOrbitals=N, for the same compound as before here’s the output lines for HOMO and LUMO (MOs 149, 150) with ThreshOrbitals set to N=1, i.e. 1% as occupation threshold (ThreshOrbitals=1):

Alpha occ 149 OE=-0.273 is Pt1-d=0.18 N4-p=0.08 N6-p=0.08 C20-p=0.06 C13-p=0.06 C48-p=0.06 C9-p=0.06 C24-p=0.05 C52-p=0.05 C16-p=0.04 C44-p=0.04 C8-p=0.03 C15-p=0.03 C17-p=0.03 C45-p=0.02 C46-p=0.02 C18-p=0.02 C26-p=0.02 C54-p=0.02 N5-p=0.01 N3-p=0.01
Alpha vir 150 OE=-0.042 is C9-p=0.18 C13-p=0.18 C44-p=0.08 C16-p=0.08 C15-p=0.06 C8-p=0.06 N6-p=0.04 N4-p=0.04 C52-p=0.04 C24-p=0.04 N5-p=0.03 N3-p=0.03 C46-p=0.03 C18-p=0.03 C48-p=0.02 C20-p=0.02

The fragment=n label in the coordinates can be used as in BSSE Counterpoise calculations and the output will show the orbital composition by fragments with the label "Fr", grouping all contributions to the MO by the AOs centered on the atoms in that fragment.

As always, thanks for reading, sharing, and rating. I hope someone finds this useful.

## Fixing the error: Bad data into FinFrg

I found this error in the calculation of two interacting fragments, both with unpaired electrons. So, two radicals interact at a certain distance and the full system is deemed as a singlet, therefore the unpaired electron on each fragment have opposite spins. The problem came when trying to calculate the Basis Set Superposition Error (BSSE) because in the Counterpoise method you need to assign a charge and multiplicity to each fragment, however it’s not obvious how to assign opposite spins.

The core of the problem is related to the guess construction; normally a Counterpoise calculation would look like the following example:

#p B3LYP/6-31G(d,p) counterpoise=2

-2,1 -1,2 -1,2
C(Fragment=1)        0.00   0.00   0.00
O(Fragment=2)        1.00   1.00   1.00
...

In which the first pair of charge-multiplicity numbers correspond to the whole molecule and the following to those of each fragment in increasing order of N (in this case, N = 2). So for this hypothetical example we have two anions (but could easily be two cations) each with an unpaired electron, yielding a complex of charge = -2 and a singlet multiplicity which implies those two unpaired electrons have opposite spin. But if the guess (the initial trial wavefunction from which the SCF will begin) has a problem understanding this then the title error shows up:

Bad data into FinFrg
Error termination via Lnk1e ...

The solution to this problem is as simple as it may be obscure: Create a convenient guess wavefunction by placing a negative sign to the multiplicity of one of the fragments in the following example. You may then use the guess as the starting point of other calculations since it will be stored in the checkpoint file. By using this negative sign we’re not requesting a negative multiplicity, but a given multiplicity of opposite spin to the other fragment.

#p B3LYP/6-31G(d,p) guess=(only,fragment=2)

-2,1 -1,2 -1,-2
C(Fragment=1)        0.00   0.00   0.00
O(Fragment=2)        1.00   1.00   1.00
...

This way, the second fragment will have the opposite spin (but the same multiplicity) as the first fragment. The only keyword tells gaussian to only calculate the guess wave function and then exit the program. You may then use that guess as the starting point for other calculations such as my failed Counterpoise one.

## 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

## Error for Gaussian16 .log files and GaussView5

There’s an error message when opening some Gaussian16 output files in GaussView5 for which the message displayed is the following:

ConnectionGLOG::Parse_Gauss_Coord().
Failure reading oriented atomic coordinates. Line Number

We have shared some solutions to the GaussView handling of *chk and *.fchk files in teh past but never for *.log files, and this time Dr. Davor Šakić from the University of Zagreb in Croatia has brought to my attention a fix for this error. If “Dipole orientation” with subsequent orientation is removed, the file becomes again readable by GaussView5.

Here you can download a script to fix the file without any hassle. The usage from the command line is simply:

˜$chmod 777 Fg16TOgv5 ˜$ ./Fg16TOgv5 name.log

The first line is to change and grant all permissions to the script (use at your discretion/own risk), which in turn will take the output file name.log and yield two more files: gv5_name.log and and name.arch; the latter archive allows for easy generation of SI files while the former is formatted for GaussView5.x.

Thanks to Dr. Šakić for his script and insight, we hope you find it useful and if indeed you do please credit him whenever its due, also, if you find this or other posts in the blog useful, please let us know by sharing, staring and commenting in all of them, your feedback is incredibly helpful in justifying to my bosses the time I spent curating this blog.

## 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: where coefficients are adjusted depending on the geometry of atoms i and j. The damping D3 function for is: where the values of s are adjustable parameters fit for the exchange-correlation functionals used in each calculation.