Gaussian How-tos

(updated )

Comments

Licensing

XSEDE was updated to ACCESS in Fall 2022, see this “cheatsheet” page for help transitioning. To use Gaussian on XSEDE machines:

Single-point energy calculations / geometry optimizations

Below is the sequence of steps for running a Gaussian calculation to compute a single-point energy or to optimize a structure.

  1. Use GaussView or any other suitable software to prepare the structure to be optimized, export it in the XYZ / MOL / MOL2 format, where you can extract the Cartesian coordinates.
  2. Create an appropriate folder hierarchy and construct a Gaussian input file called gaussian.gjf: 
    cd /path/to/projects/ABC/moleculeX
vim gaussian.gjf
  3. A basic Gaussian input file may look like the following. You can read about the different keywords on its website. Important notes and common changes include:
    • Gaussian has a very strict formatting: the empty lines between different sections, the four stars and two empty lines after the basis set definition.
    • Lines starting with a exclamation mark are comment lines and ignored.
    • DFT functional: M062X, M06L.
    • Basis set and Gen keyword: AUG-cc-pVDZ, AUG-cc-pVTZ; 3-21G, 6-311+G(3df,3pd); etc.
    • Charge and spin states: 0 1 means a neutral singlet.
    • Fill in the Cartesian coordinates you prepared earlier.
    • nprocshared and mem will be adjusted automatically by the submission script below.
    %nprocshared=128
%mem=253000MB
%chk=gaussian.chk
# M062X/Gen/Auto
# Opt
!# Opt=ReadFC
!# Geom=(Check, NewDefinition) Guess=Read
# Integral(Grid=UltraFine) SCF=(Tight,IntRep,XQC) NoSymm
!#p Integral(Grid=UltraFine) SCF=(Conver=4,VarAcc,IntRep,XQC,MaxCycle=512,VShift=400) NoSymm

Molecule

0   1
 H               -118.70001270  -91.91738745 -127.75908914
... More lines of Cartesian coordinates
 O               -120.07157063  -94.57259547 -129.34664170

H O 0
AUG-cc-pVTZ
****
  4. Prepare submission scripts. The best practice is to write a generic submission script and put it in the local asset directory (e.g., /path/to/projects/ABC/bin/). Then create symbolic links: 
    ln -sf ../bin/sub-gaussian.script .  # Run in current directory /path/to/projects/ABC/moleculeX

    NOTE: these symlinks must use relative paths instead of absolute paths, because the local asset directory will be synced to the remote target under your home directory, which may have a different absolute path.

  5. Submit the job using shotgun 
    shotgun submit -c 8 -m 2000 -n 1 -t 24 -r bridges2 -q RM-shared -j ../../bin/sub-gaussian-bridges2.script  # submitting the job

    NOTE: the ./ part in ./sub-gaussian-bridges2.script is required. See shotgun --help submit for a list of options that can be used.

  6. If the job runs out of time and you want to continue the calculation, you have two choices:
    • Copy the Cartesian coordinates of the last configuration from output file gaussian.log to gaussian.gjf. There is a script that may help with this task:
    /home/common/sample-input/gaussian/bin/g_copy_last_geom.sh gaussian.log gaussian.gjf
    • Read frequencies and coordinates from the checkpoint file gaussian.chk by commenting and un-commenting the relevant lines.
    !# Opt
# Opt=ReadFC
# Geom=(Check, NewDefinition) Guess=Read

    You also need to comment out or delete the coordinates lines. When all is done, rename the output files and submit the job again:

    mv gaussian.log equil1.log
shotgun submit -c 8 -m 2000 -n 1 -t 24 -r bridges2 -q RM-shared -j ../../bin/sub-gaussian-bridges2.script  # submitting the job

    If you keep the structure of the sample input file, you can try using a script to make this change:

    /home/common/sample-input/gaussian/bin/g_change.sh -c gaussian.gjf  # -c for continuation, absence of the flag will switch back to the normal calculation

Pseudopotentials for heavy atoms

You may want to use effective core potentials starting from the first-row transition metal elements. You can use the [GenECP keyword and specify the flavor of pseudopotential. A basic example is given below. Again note the empty lines and four-star separators.

# M06L/GenECP/Auto
...skipping the remaining input until the last Cartesian coordinates line
Cr 0 -33.587176 33.370038 -45.351159

H C O F 0
AUG-cc-pVDZ
****
Cr 0
SDD
****

Cr 0
SDD

Constrained optimizations

  • On your workstation
  1. In Materials Studio, you can constrain atoms by clicking the menu item Modify -> Constraints and tick Fix Cartesian position in the dialogue.
  2. File -> Export the structure, choose the .xsd format, and click Options to select version 6.1.
  3. Use rsync / scp / WinSCP (Windows only) to upload the file to group server.
  • On the group server
  1. Convert the .xsd file to Gaussian format. 
    cd /path/to/projects/ABC/moleculeX  # go to the directory where the .xsd file is uploaded.
module add vasp-tools/linux materialsstudio  # in order to use the vasp-tools Materials Studio integration
msvasp ms2gaussian moleculeX.xsd gaussian.gjf  # converts the .xsd file to Gaussian format
  2. Edit gaussian.gjf to include other input parameters.

Stretching/bending/torsion scans

See the tutorial in the linked section header above to perform Gaussian scans.

Free energy calculations using implicit solvent models

In order to find the free energy of your system in the solvent, follow the following steps:

  1. Add the following keyword to the Route section of the input file to calculate the harmonic frequencies at the end of a geometry optimization job. Thermochemistry will be calculated at 298K unless specified otherwise.

    NOTE Please read and understand the document written by Joseph W. Ochterski from Gaussian, Inc. on our group’s Google Drive or on the Gaussian website

    # Freq
# Temperature=523
  2. Add the following keyword to the Route section of the input file with an appropriate solvent model. The SMD solvation model developed by Cramer and Truhlar at Minnesota is often a good choice. 
    # SCRF=(SMD,Solvent=n-Octanol)
  3. In certain cases (see discussions here and here), depending on the system and the solvation model, it may be more appropriate to perform a gas-phase optimization and a single-point energy calculation at the gas phase optimized structure. The solvent phase free energy is then calculated as Gsolv = Esolv + ΔGcorr_gas, where, Gsolv is the free energy of the system in the solvent, Esolv is the electronic energy of the system in the solvent, and ΔGcorr_gas is the gas phase thermal correction term to the free energy.
  4. If you keep the structure of the sample input file, you can use a script to switch on and off solvation and frequency calculations: 
    /home/common/sample-input/gaussian/bin/g_change.sh -f -s gaussian.gjf  # -f for frequency, -g for solvation, absence of a flag will turn off the corresponding option

Calculations of redox potentials

On Zotero, the Electrochemistry/Redox folder includes a few useful papers. Of particular interest is Ho, J.; Coote, M. L.; Cramer, C. J.; Truhlar, D. G. Chapter 6 Theoretical Calculation of Reduction Potentials. In Organic Electrochemistry; CRC Press, 2016; pp 229–259.

Here we use an example from J. Am. Chem. Soc. 2019, 141, 6392, specifically one of the reactions illustrated in the top left panel of the figure: L2Ni(+1)Br + e-→ L2Ni(0) + Br- (all species in DMF solvent). This is written as a reduction reaction. The reaction free energy can therefore be converted to the reduction potential via the Nernst equation.

To calculate the reaction free energy, three geometry optimization jobs need to be performed for the three chemical species using implicit DMF solvent. The following excerpt can be found near the end of the output files for L2Ni(+1)Br:

 SCF Done:  E(UM06L) =  -1332.49281641
...
...
 Zero-point correction=                           0.428097 (Hartree/Particle)
 Thermal correction to Energy=                    0.458133
 Thermal correction to Enthalpy=                  0.459077
 Thermal correction to Gibbs Free Energy=         0.362588
 Sum of electronic and zero-point Energies=          -1332.064719
 Sum of electronic and thermal Energies=             -1332.034684
 Sum of electronic and thermal Enthalpies=           -1332.033739
 Sum of electronic and thermal Free Energies=        -1332.130229

There is a script to extract these numbers, but please make sure to double check the first few times you use the script for a new type of calculations.

/home/common/sample-input/gaussian/bin/g_grep_result.sh -g gaussian.log

Collecting the results into the table below, where all species were computed in DMF solution, except the last column, we have the electronic energy, enthalpy, and free energy at T = 298 K and ρ0(1atm) = p/RT = 101325/8.3144621/298 = 40.89 mol/m3. From the last two columns, we can in principle compute the solvation free energy for Br- [i.e., the reaction free energy for Br- (gas, ρ0) → Br- (DMF, ρ0) ]:  − 13.554610 − (−13.468305) =  − 0.086305 Ha =  − 0.086305 * 2625.4997356709 kJ/mol =  − 226.59 kJ/mol. However, in standard definition, solution-phase species should use a reference state with a density of ρ1M = 1 mol/L = 1000 mol/m3 and we should compute the free energy for Br- (gas, ρ0) → Br- (DMF, 1 M)]. To do this, we need to add the conversion factor G[Br-(DMF, ρ1M)] - G[Br-(DMF, ρ0)] = RTln (ρ1M/ρ0) = 8.3144621 * 298 * ln (1000/40.89) J/mol = 7.92 kJ/mol, to give  − 226.59 + 7.92 =  − 218.7 kJ/mol.

Compared to the experimental value, -241.8 kJ/mol (Marcus, Y. Ions in solutions and their solvation, 2015 John Wiley & Sons, Inc.), the implicit solvation calculation is underestimated by 23 kJ/mol. It is common for these small ions to use the experimental solvation free energy instead. Therefore, we take G[Br-(DMF, ρ1M)] =  − 13.468305 − 241.8/2625.4997356709 =  − 13.5604018 Ha. Finally, the reaction free energy is G[Br-(DMF, ρ1M)] + G[L2Ni(0)(DMF, ρ1M)] - G[L2Ni(+1)Br(DMF, ρ1M)] = G[Br-(DMF, ρ1M)] + G[L2Ni(0)(DMF, ρ0)] - G[L2Ni(+1)Br(DMF, ρ0)] =  − 13.5604018 − 1318.665092 − (−1332.130229) =  − 0.0952648 Ha =  − 250.12 kJ/mol.

Species

L2Ni(+1)Br

L2Ni(0)

Br-

Br- (gas)

E

-1332.492816

-1319.034577

-13.538435

-13.452129

H

-1332.033739

-1318.579726

-13.536074

-13.449769

G

-1332.130229

-1318.665092

-13.554610

-13.468305

Using the Nernst equation, we can convert this standard-state reaction free energy to standard reduction potential:  − (−250.12) * 1000/96485.33289 = 2.59 V. Finally, if we want to compute reduction potential relative to the standard hydrogen electrode, we need to know the standard-state reaction free energy for H+ (DMF, 1 M) + e-→ 1/2H2 (gas, 1 atm), and its corresponding standard reduction potential. This is a completely analogous calculation. We use the experimental solvation free energy for H+ in DMF, -1118.98 kJ/mol. Reaction free energy:  − 1.17094/2 * 2625.4997356709 − (−0.01*2625.4997356709−1118.98) =  − 391.92 kJ/mol, which is  − (−391.92) * 1000/96485.33289 = 4.06 V. Therefore, the reduction potential for L2Ni(+1)Br + e-→ L2Ni(0) + Br- is 2.59 − 4.06 =  − 1.47 V. If we want to express this to other reference such as Ag/AgNO3 (10 mM), we just need to further subtract the potential of Ag/AgNO3 (10 mM) relative to SHE, which is  − 1.47 − 0.5 =  − 2.0 V.

Species

H+ (gas)

H2 (gas)

E

0

-1.16958

H

0.00236

-1.15615

G

-0.01

-1.17094

You can find this and other computed results in the figure on pg. S27 of the Supporting Information.

Adsorbed Molecular Free Energy Correction

Unlike gas molecules, adsorbed molecules form chemical bonds with substrate, which limits the translational and rotational freedom of the adsorbed molecules. So the contribution of translation and rotation to entropy and enthalpy is significantly reduced (so called hindered translator / hindered rotor model). This does not mean no translational or rotational contribution.

One common method is to attribute the translational or rotational part of the contribution to vibration, that is, the 3N vibrations of the surface-adsorbing molecules (except the virtual frequency) are all used to calculate the correction of the thermo energy.

The small the vibration frequencies have large contribution to entropy. It is very likely that a small vibration frequency will lead to abnormal entropy and free energy correction. So, it suggests that when the free energy of the surface adsorption molecule is corrected, the contribution of the frequency below 50 cm-1 is calculated as 50 cm-1, and VASPKIT also does this. And VASPKIT neglects the electron motion because of its small contribution.

And the PV term is also delete. At last VASPKIT will output:

1. Zero-point energy εZPE;

2. Thermal correction to U(T) = H(T) = εZPE + ΔU0 → T

3. Thermal correction to G(T) = εZPE + ΔU0 → T + TS = εZPE + ΔG0 → T

Example: Thermo energy correction for adsorbed O on Au(111)

(1) Fix all slab metal atoms, do frequency calculation for the adsorbed molecule. Check the frequencies:

   $ grep cm OUTCAR
      1 f  =   10.904836 THz    68.517103 2PiTHz  363.746154 cm-1    45.098792 meV
      2 f  =   10.827330 THz    68.030119 2PiTHz  361.160834 cm-1    44.778253 meV
      3 f  =   10.751668 THz    67.554721 2PiTHz  358.637024 cm-1    44.465340 meV

(2) In the same folder, run VASPKIT 501. Input temperature(K).

   501
    +-------------------------- Warm Tips --------------------------+
    This Feature Was Contributed by Nan XU, Qiang LI and Jincheng LIU.
        See An Example in vaspkit/examples/thermo_correction/ORR.
    Only vibrations! No Translation & Rotation & Electron contributions.
    +---------------------------------------------------------------+

    Please Enter The Temperature (K):

    ------------>>
   298.15
     -->> (01) Reading OUTCAR File...
    +-------------------------- Summary ----------------------------+
    Neglect PV contribution to translation for adsorbed molecules.
    To avoid abnormal entropy contribution,
    frequencies less than 50 cm-1 are set to 50 cm-1.
    U(T) = H(T) = ZPE + Delat_U(0->T)
    G(T) = H(T) - TS = ZPE + Delat_U(0->T) - TS

    Temperature (K): 298.1
    Zero-point energy E_ZPE   :       1.549 kcal/mol    0.067171 eV
    Thermal correction to U(T):       2.206 kcal/mol    0.095670 eV
    Thermal correction to H(T):       2.206 kcal/mol    0.095670 eV
    Thermal correction to G(T):       1.207 kcal/mol    0.052342 eV
    Entropy S                 :      14.021 J/(mol*K)   0.000145 eV

So, ZPE is 0.067171 eV. Thermal correction to U(T), H(T), and G(T) are 0.095670 eV, 0.187167 eV, and  − 0.446723 eV, respectively.

Calculation of CM5 charges (Gaussian/cm5pac)

See the tutorial in the linked section header above to perform CM5 charge calculations.

NEB calculations

See the tutorial in the linked section header above to perform NEB calculations.

Comments

Comments powered by Disqus