G09

go back to Main Page, Computational Resources, Chemistry & More, Computational Codes, GAUSSIAN

---

You MUST use the "p" keyword for extra-print in all your calculations.

This is an essential requirement to use afterwards the Chemical Repository, SCIPIO, where you will need to save your data.

---

How to submit g09 jobs on kimik2

-- For Feliu Maseras Group --

for the 12 processors queue:

qs g09.c12m24 nameinput.in

for the 8 processors queue:

qs g09.c8m24 nameinput.in

for the 4 processors queue:

qs g09.c4m8 nameinput.in

or

qs g09.cq4m4 nameinput.in

Revision Issues

• Revision A1 and A2;
  • do NOT perform the Dispersion calculation properly. You need the revision B1 or onwards.
  • implementation of the QZVP basis set is incorrect, it has not the polarization functions. You need the revision B1 or onwards.

• Revision C1;
  • problems with the modredundant section. To freeze a distance you can not use "B atom1 atom2 value F", as this does not do anything. Instead, you will have to move the geometry to the value and use "B atom1 atom2 F".
  • for some geometries the calculation stops after reading the initial structure.

• Revision D1;
  • several issues from previous versions have been sorted, and some new interesting features added, such as DFT-D3 for most functionals, ...
  • It still does not behave properly with the modredundant section: you can not use "B atom1 atom2 value F" (see above)
  • A severe error concerning minima optimizations and frozen geometries of specific systems has been observed. The error consist on a call as if the chk file was incorrect/corrupted "Operation on file out of range...dumping /fiocom/..a large list of numbers...Error termination in NtrErr: NtrErr Called from FileIO." We don't know yet exactly what is special in the systems that crash, although we suspect that it might have something to do with having to connected aromatic/planar rings. The error does not occur in TS searches. There are several solutions to overcome this depending on the project needs: i) to use the version g09.a02 if the the method of choice is available on that version. ii) to use a guess=always, although it is expensive. iii) the use of opt=cartesian also helps out to diminish the probability to get the error although it is does not always work and it is also expensive.

How to

• to view frequencies in gaussview3

• how to create .fchk files

• how to compute BSSE

Basis sets

• About 5d/6d basis sets. The keywords 5D, 6D, 7F, and 10F are used to specify use of Cartesian or pure d and f (and higher) functions. The defaults are different depending if writing the basis sets in the gen section (at the end of the input file) or in the command line (#p opt belyp/6-31+g*). When using the Gen keyword the defaults are 5D and 7F (pure functions). When specifying the basis in the command line the defaults are 6D and 10F functions (Cartesian functions). The defaults can be changed typing 5d or 6d in the command line accordingly. It's very important to keep the use of one or another number of functions in all calculations, as if not there will be different basis sets and that will induce an error.
• About the QZVP basis set. The implementation is different between the A.02 (and A.01) and the B.01 versions of g09. The basis set in A.01 and A.02 miss the polarization functions, which means that these versions stored a QZV basis set. If you want to use the true QZVP basis, you should use the B.01 version. DO NOT RUN GAUSSIAN QZVP with REVISIONS A1 and A2 = KIMIK2

Dispersion

In revisions before G09.D1: One of the new features in Gaussian G09, is the possibility to include dispersion. This can be done in two different ways.

• By using the B97D functional, which includes dispersion. With this you can optimize or do single point calculations.

• With the iop: iop(3/124=3),which allows to compute the Dispersion on a given structure optimized at your favourite level of theory (to add Dispersion as single point). The good point of this is that can be done using any set of basis and method (using SDDall you can save lot's of computer time).

The issue is that the first and second revision of Gaussian A1 and A2 do NOT perform the Dispersion calculation properly. You need the revision B1 (not in kimik2 look at cesca).

For revision G09.D1: There are different ways of including dispersion:

• By using B97D or B97D3

• By using the keyword EmpiricalDispersion = and choosing a dispersion model

Large systems frequency calculations

For some reason, it seems a numerical problem, frequency calculations with large systems take too long. To reduce the computer time you can try to use the keyword Int=Acc2E=11 According to the manual it increases the accuracy from ^10 to ^11, but probably due to some bug, decreases the computational cost.

ONIOM in Gaussian09

• New force field implantation for ONIOM(QM:MM)

• ONIOM tips

• Convergence problems with ONIOM in G09

• Transition state searches with ONIOM

Issues concerning the coordinate system

Some extra checks of the completeness of the coordinate system were added in Gaussian09. Thus in some cases you might get the error message "Error in internal coordinate system" where it did not appear in Gaussian03.

Typical cases where this error message appear, are for linear system. E.g. when three covalently bond atoms are arranged on a line. For organometallic systems a typical case is for carbonyl ligands coordinated to a metal centre, M-CO. But you can also find it in many other cases.

• If your system contains any linearly arranged atoms and you get this error message, a solution might be to add linear bend coordinates to these atoms. In your modredundant section you write: L atom1 atom2 atom3 -1 A

• In the error message you will get the dihedral/angle affected, you can try to kill the problematic angle/dihedral atom1 atom2 atom3 K

• Modify a bit the geometry, by changing a bit the angles/dihedrals or even just deleting from the third decimal of the xyz coordinates.

• To optimize in cartesian coordinates could be also a solution to this problem although your calculations are going to be slower and you are not going to use modredun (not possible to freeze and scan coordinated).

• Another solution is to use a Z-matrix and add ghost atoms, that will allow you to define the system without angles and dihedrals close to 0 and 180.

• Revision C1 of Gaussian09 might have improved the way in which the internal coordinates are defined, as for some systems exactly the same input stops for RevA1 and runs without problems with RevC1. Be aware that RevC1 has issues with the modredundant section (see above)

• The Error does of course depend on the geometry, but also on the method/basis, so the same geometry might give the error with, let's say, M06 and not with B3LYP. So you can optimize the geometry "a bit" with another method and then try again with the method you whish.

Solvent Effects

As with Gaussian03, the SCRF keyword requests that a calculation be performed in the presence of a solvent by placing the solute in a cavity within the solvent reaction field.

• The integral equation formalism variant IEFPCM, is the default SCRF method. It has not changed from Gaussian03, BUT the formalism used and its implementation has changed. That is: you will get different results in G03 and G09, using the same method.

• The default RADII used in G03 was UAO , whilst now is UFF with spheres placed by default on all Hydrogen atoms. No need to use the SPHEREONH= keyword, except if you use UAO or another radii model that does not have them explicitly.

• It seems to give no "convergence failure" problems, but maybe its too early to say so.

• It is able to perform frequency calculations in solvent, giving enthalpies, free energies, ZPE corrections...

• The G09 output for PCM gives only one energy which corresponds, by default, to the electrostic energy. The non-electroestatic components of the energy of solvation are not included. If you wish to include them you will need to include "scrf=(read,PCM,...)" in the command line, and Cav Dis Rep (in separate lines) at the end of the input.

• The SMD method includes automatically the non-electrostatic terms.

The RI-DFT method

From the Gaussian manual:

Gaussian 09 provides the density fitting approximation for pure DFT calculations. This approach expands the density in a set of atom-centered functions when computing the Coulomb interaction instead of computing all of the two-electron integrals. It provides significant performance gains for pure DFT calculations on medium sized systems too small to take advantage of the linear scaling algorithms without a significant degradation in the accuracy of predicted structures, relative energies and molecular properties. Gaussian 09 can generate an appropriate fitting basis automatically from the AO basis, or you may select one of the built-in fitting sets.

What all this means is that for all calculations that have more than about 1000 basis functions, it will speed up the calculation greatly if you use the RI-DFT method. Generally, the greater the number of basis functions, the better it compares to normal DFT. The accuracy of this approximation is very good, and generally you can afford to use a much higher basis set than you would with a normal DFT calcualtion (so you get rid of BSSE errors etc.).

How to set up an RI-DFT calculation.

You have to use a pure DFT functional... no hybrids allowed...

You also have to specify TWO basis sets... The first one is the normal basis set specification (GENECP is supported for complicated setups with transition metals and ECPs), then the second one is the "density fitting basis", which you can specify yourself or Gaussian can generate automatically.

For example:

# PBEPBE/TZVP/TZVPFit

would specify a calculation with the PBEPBE functional, TZVP "normal" basis set, and TZVPFit "density fitting" basis set. The / sign between each section is essential.

For an example of an input file, and more information on how to set up the calculation, see the following page... RI-DFT input

MECP calculations (Old Prof. Harvey version)

To find Minimum Energy Crossing Points in Gaussian you can use the program developed by Prof. Jeremy Noel Harvey. To use the program for first time you will need a bit of patience.

1.- Download the program from the repository:

• Connect to the repository: http://aliga.iciq.es/fsn/ ask Martín for the username and password.
• Go to the directory /Gaussian/MECP/ and copy all files and directories into a new directory in your home.

2.- There are two directories the PROGRAM and the SHELL. There is also a "README" file that can help if you follow the instructions carefully, and another file called "envia" which is a submitting script.

2.1.- In the PROGRAM directory you will find the fortran program, there you have to compile the program every time that the number of atoms of your system changes. That is done by:

• Delete all *.o files and the MECP.x file,
• Edit the file MainProgram.f: Change the value of the variable Natom, that is change number 5 in the line "parameter (Natom=5)" by the number of atoms of your system.
• Re-compile the program by typing "make". This will generate a new version of the executable file MECP.x.
• Copy the MECP.x generated in the directory you want to work in.

2.2.- In the SHELL directory there is a number of files you will need to mount your system. Copy all of them to your working directory. Go to your working directory and start editing the files:

• Change the number of atoms in files first_progfile and master_script.
• Add your starting geometry in the file geom. The geometry should be in the same format as in the model file: number of the element and xyz.
• Read the file geom in the file first_progfile.
• Copy the first_progfile to ProgFile. The first_progfile will not be changed along the calculation but the ProgFile will.
• Edit the 'Input_Header_A' and 'Input_Header_B' according to your methodology requirements. One Input_Header for each spin state. They should have different names for the .chk files.
• Use the 'footer' file to include the basis sets if necessary.
• Once the Input_Header files, the geom and the footer are ready type "csh -f Make_First_Inputs" to generate automatically the first input files: Job0_A.gjf and Job0_B.gjf. If you don't have the .chk files available, open the two input files Job0_A.gjf and Job0_B.gjf and delete "guess(read)".
• Submit the calculations with a submitting script. I use the "envia" file you will find in the main directory. Please edit the "envia" file according to the machines you want to use. To submit it to the queue just type "qsub envia".

3.- If you do not use DFT or you have special needs you might have to change the extract_energy or extract_gradient files by the ones available in the SHELL directory, or the main_script file.

MECP calculations (Python version)

Jaime Rodríguez-Guerra and Ignacio Funes have adapted the old-FORTRAN script to find Minimum Energy Crossing Points in Gaussian (by Prof. Jeremy Noel Harvey). Now, a little script governs all the process and you do not need to recompile and copy the files when your number of atoms change. In addition, some issues such as the generation of initial input files, the selection of the Gaussian version (g09 or g16), the thresholds, maximum number of steps, have been automatized.

1.- Download the program from the easyMECP's GitHub repository:

• Go to the following web-page and download the program. You just need to copy the easymecp.py script to your working directory. The script has many arguments that can be easily specified (Gaussian version, thresholds...), you can find all the information by typing python easymecp.py -h. Remember that if you want to execute python scripts in your kimik2 home directory, you should activate the corresponding module (module load python/3.5.2).

2.- Prepare your input files. It should contain the following:

• Input_header_A and Input_header_B: this represents your upper part of the calculation (up to charge/multiplicity) for each multiplicity.
• geom: this file contains your geometry. You do not need to format it with atomic numbers: if your input has the atomic symbols, now the program works fine!
• footer: contains additional information (normally basis set).

3. Run the program: you need to load the Python module in the machine, so the submit script should be modified. Here you have examples that you can modify according to the machine you are using. In addition, this script calls directly the program and it should have the corresponding flags to define thresholds, gaussian version. One example for g09 and other for g16 can be download here: File:Mecpy pruebas.tgz

Also, if you need to restart your calculation, just review your file geom and send the script again, it automatically save your old results in JOBS1, JOBS2... directories. You can find more useful information in the GitHub webpage!

BOMD calculations

This keyword requests a classical trajectory calculation using a Born-Oppenheimer molecular dynamics model.

http://www.gaussian.com/g_tech/g_ur/k_bomd.htm

• Giving initial velocities: If you want to give initial velocities, Gaussian is able to generate them providing the energy you want to add into a given normal mode,

by using the keyword nsample and adding after the coordinates:Mode-num, VibEng(Mode-num), ... for as many modes you wish.

• • As example: this will provide an energy of -16.4 kcal/mol to mode 1

BOMD(Phase=(1,14),RTemp=300,nsample=1,update=12,maxpoints=2,ntraj=10)

title & charge-multiplicity & coordinates

0
1 -16.4

• • Please note that in the current version, the sign of the energy is not changing the direction of the velocities generated. So, to change the direction, you will have to run a calculation just to generate the number of velocities you wish ntraj = XX (with maxpoints=2 is very fast). Then take all XX velocities MW Cartesian velocity: from the output (remove all text and leave only the three columns of velocities). If the direction is not the desired one change the sign of your velocities (to know it make a test with a small number of steps to each direction). Then read the velocities with the sign changed in a new input after the coordinates (with the NPath on the top and a blank line separating velocities), you will have to use the keyword ReadMWVelocity in the BOMD specifications. You need a set of velocities for each trajectory if not all trajectories will be equal.

• Visualizing BOMD trajectories: You can visualize each trajectory separately with GaussView.
  • The geometry of your last structure in GaussView does not correspond to the last geometry of your trajectory but to the last geometry found in the output file.
  • If the number of steps of your trajectory is very large you might want to kill your GaussView and open it again than wait until GaussView closes your file.

Gaussian09 vs Gaussian03

• Scaling and New Functional Testing

• Solvation models, Polarizable Continuum Model (PCM)

• Tests on solvation energies, PCM in G03 and G09 and SMD