A few weeks back we wrote about using WFN(X) files with MultiWFN in order to find σ-holes in halogen atoms by calculating the maximum potential on a given surface. We later found out that using a chk file to generate a wfn(x) file using the guess=(read,only) keyword didn’t retrieve the MP2 wavefunction but rather the HF wavefunction! Luckily we realized this problem very quickly and were able to fix it. We tried to generate the wfn(x) file with the following keywords at the route section
#p guess=(read,only) density=current
but we kept retrieving the HF values, which we noticed by running the corresponding HF calculation and noticing that every value extracted from the WFN file was exactly the same.
So, if you want a WFN(X) file for post processing an MP2 (or any other post-HartreFock calculation for that matter) ask for it from the beginning of your calculation in the same job. I still don’t know how to work around this or but will be happy to report it whenever I do.
PS. A sincere apology to all subscribers for getting a notification to this post when it wasn’t still finished.
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:
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
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.
Editing large molecules on a seemingly simple visualizer as GaussView can be a bit daunting. I’m working on a follow up of that project we recently published in JACS but now we require to attach two macrocycles to the organometallic moiety; the only caveat is that this time we don’t have any crystallographic data with which to start. Generating a 3D model of this structure is already hard enough and even when you managed to do it there are many degrees of freedom that in some cases can lead to unrealistic geometries after optimization.
I recently came across a simple way to edit a large complicated molecule by optimizing the fragments separately and then joining them in a new molecule by using the clipboard. This rather simple method, that I for one had never exploited has just saved me a few good hours.
Copy a molecule (CTRL+C) and it will go to the clipboard as a molecular fragment for which you can define a new hot atom and thus bind it to the other fragment as you would with the regular builder. I strongly suggest to use a “New Molecule Group” instead of editing over an existing molecule. Also, if you are using the “paste” button, observe that it has three different options; you may want to use the last one “append to existing molecule” or you will have your fragments in different windows.
And remember, dihedral angles are your best friends.
As far as population analysis methods goes, the Quantum Theory of Atoms in Molecules (QTAIM) a.k.a Atoms in Molecules (AIM) has become a popular option for defining atomic properties in molecular systems, however, its calculation is a bit tricky and maybe not as straightforward as Mulliken’s or NBO.
Personally I find AIM a philosophical question since, after the introduction of the molecule concept by Stanislao Cannizzaro in 1860 (although previously developed by Amadeo Avogadro who was dead at the time of the Karlsruhe congress), the questions of whether or not an atom retains its identity when bound to others? where does an atom end and the next begins? What are the connections between atoms in a molecule? are truly interesting and far deeper than we usually consider because it takes a big mental leap to think about how matter is organized to give rise to substances. Particularly I’m very interested with the concept of a Molecular Graph which in turn is concerned with the way we “draw lines” to form conceptual molecules. Perhaps in a different post we can go into the detail of the method, which is based in the Laplacian operator of the electron density, but today, I just want to collect the basic steps in getting the most basic AIM answers for any given molecule. Recently, my good friend Pezhman Zarabadi-Poor and I have used rather extensively the following procedure. We hope to have a couple of manuscripts published later on. Therefore, I’ve asked Pezhman to write a sort of guest post on how to run AIMALL, which is our selected program for the integration algorithm.
The first thing we need is a WFN or WFX file, which contains the wavefunction in a Fortran unformatted file on which the Laplacian integration is to be performed. This is achieved in Gaussian09 by incluiding the keyword output=wfn or output=wfx in the route section and adding a name for this file at the bottom line of the input file, e.g.
(NOTE: WFX is an eXtended version of WFN; particularly necessary when using pseudopotentials or ECP’s)
Analyzing this file requires the use of a third party software such as AIMALL suite of programs, of which the standard version is free of charge upon registration to their website.
OpenAIMStudio (the accompanying graphical interface) and select the AIMQB program from the run menu as shown in figure 1.
Select your WFN/WFX file on which the calculation is to be run. (Figure 2)
You can control several options for the integration of the Laplacian of the electron density as well as other features. If your molecules are simple enough, you may go through with a successful and meaningful calculation using the default settings. After the calculation is finished, several result files are obtained. We’ll work in this tutorial only with *.mpgviz (which contains information about the molecular graph, MG) and *.sum (which contains all of needed numerical data).
Visualization of the MG yields different kinds of critical points, such as: 1) Nuclear Attractor Critical Points (NACP); 2) Bond Critical Points (BCP); 3) Ring CP’s (RCP); and 4) Cage CP’s (CCP).
Of the above, BCP are the ones that indicate the presence of a chemical bond between two atoms, although this conclusion is not without controversy as pointed out by Foroutan-Njead in his paper: C. Foroutan-Nejad, S. Shahbazian and R. Marek, Chemistry – A European Journal, 2014, 20, 10140-10152. However, at a first approximation, BCP’s can help us to explore chemical interactions.
Now, let’s go back to visualizing those MGs (in our examples we’ve used methane and ethylene and acetylene). We open the corresponding *.mpgviz file in AIMStudio and export the image from the file menu and using the save as picture option (figure 3).
The labeled atoms are NACP’s while the green dots correspond to BCP’s. Multiplicity of a bond cannot be discerned within the MG; in order to find out whether a bond is a single, double or triple bond we have to look into the *.sum file, in which we’ll take a look at the bond orders between pairs of atoms in the section labeled “Diatomic Electron Pair Contributions and Delocalization Data” (Figure 4).
Delocalization indexes, DI’s, show the approximate number of electrons shared between two atoms. From the above examples we get the following DI(C,C) values: 1.93 for C2H4 and 2.87 for C2H2; on the other hand, DI(C,H) values are 0.98 for CH4, 0.97 in C2H4 and 0.96 in C2H2. These are our usual bond orders.
This is the first part of a crash tutorial on AIM, in my opinion this is the very basics anyone needs to get started with this interesting and widespread method. Thanks to all who asked about QTAIM, now you have your long answer.
Thanks a lot to my good friend Dr Pezhman Zarabadi-Poor for providing this contribution to the blog, we hope you all find it helpful. Please share and comment.
A couple of weeks ago I posted a solution for a common error regarding .fchk files that will display the error below when opened with GaussView5.0. As I expected, this error has to do with the use of diffuse functions in the basis set and is related to a change of format between Gaussian versions.
CConnectionGFCHK::Parse_GFCHK() Missing or bad data: Alpha Orbital Energies Line Number 1234
Although the method described in the previous post works just fine, the following update is a better approach. Due to a change of spelling between G03 and G09 (which has been corrected for G09 but not available for GV versions prior to 5.0.9) one must change “independent” for “independant”
To make the change directly from the terminal the following command is needed:
sed -i 's/independent/independant/g' file.fchk
Alternatively you can redirect the output to a new file
sed -e 's/independent/independant/g' file.fchk > newfile.fchk
if you want to keep the old version and work with a new one.
Of course this edition can be performed manually with any text editor available (for example if you work in Windows) but solutions from the terminal always seem easier and a lot more fun to me.
Thanks to Dr. Fernando Cortés for sharing his insight into this issue.
I’ve found the following error regarding the opening of .fchk files in GaussView5.0.
CConnectionGFCHK::Parse_GFCHK() Missing or bad data: Alpha Orbital Energies Line Number 1234
The error is prevented to a first approximation (i.e. it at least will allow GV to open and visualize the file but other issues may arise) by opening the file and modifying the number of basis functions to equal the number of independent functions (which is lower)
FILE HEADER FOpt RM062X 6-311++G(d,p) Number of atoms I 75 Info1-9 I N= 9 163 163 0 0 0 110 2 18 -502 Charge I 0 Multiplicity I 1 Number of electrons I 314 Number of alpha electrons I 157 Number of beta electrons I 157 Number of basis functions I 1199 Number of independent functions I 1199 Number of point charges in /Mol/ I 0 Number of translation vectors I 0 Atomic numbers I N= 75 ... ... ... ...
Once both numbers match you can open the file normally and work with it. My guess is this will continue to happen with highly polarized basis sets but I need to run some tests.
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.
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  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.
#p opt=(qst2,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)
Title card for reagents
Cartesian Coordinates for reagents
Title card for products
Cartesian Coordinates for products
#p opt=(qst3,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)
Title Card for reagents
Cartesian Coordinates for reagents
Title card for products
Cartesian Coordinates for products
Title card for TS
Cartesian Coordinates for TS
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.
#p m062x/6-31++G(d,p) IRC=(Maxpoints=50,RCFC,phase=(2,1))Temperature=373.15 SCRF=(Solvent=Water) geom=allcheck
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.
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 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.
Calculating both Polarizability and the Hyperpolarizability in Gaussian is actually very easy and straightforward. However, interpreting the results requires a deeper understanding of the underlying physics of such phenomena. Herein I will try to describe the most common procedures for calculating both quantities in Gaussian09 and the way to interpret the results; if possible I will also try to address some of the most usual problems associated with their calculation.
The dipole moment of a molecule changes when is placed under a static electric field, and this change can be calculated as
pe = pe,0 + α:E + (1/2) β:EE + … (1)
where pe,0 is the dipole moment in the absence of an electric field; α is a second rank tensor called the polarizability tensor and β is the first in an infinite series of dipole hiperpolarizabilities. The molecular potential energy changes as well with the influence of an external field in the following way
U = U0 – pe.E – (1/2) α:EE – (1/6) β:EEE – … (2)
Route Section Keyword: Polar
This keyword requests calculation of the polarizability and, if available, hyperpolarizability for the molecule under study. This keyword is both available for DFT and HF methods. Hyperpolarizabilities are NOT available for methods that lack analytic derivatives, for example CCSD(T), QCISD, MP4 and other post Hartree-Fock methods.
Frequency dependent polarizabilities may be calculated by including CPHF=RdFreq in the route section and then specifying the frequency (expressed in Hartrees!!!) to which the calculation should be performed, after the molecule specification preceded by a blank line. Example:
#HF/6-31G(d) Polar CPHF=RdFreq Title Section Charge Multiplicity Molecular coordinates ==blank line== 0.15
In this example 0.15 is the frequency in Hartrees to which the calculation is to be performed. By default the output file will also include the static calculation, that is, ω = 0.0. Below you can find an example of the output when the CPHF=RdFreq is employed (taken from Gaussian’s website) Notice that the second section is performed at ω = 0.1 Ha
SCF Polarizability for W= 0.000000: 1 2 3 1 0.482729D+01 2 0.000000D+00 0.112001D+02 3 0.000000D+00 0.000000D+00 0.165696D+02 Isotropic polarizability for W= 0.000000 10.87 Bohr**3. SCF Polarizability for W= 0.100000: 1 2 3 1 0.491893D+01 2 0.000000D+00 0.115663D+02 3 0.000000D+00 0.000000D+00 0.171826D+02 Isotropic polarizability for W= 0.100000 11.22 Bohr**3.
You may have noticed now that the polarizabilities are expressed in volume units (Bohr^3) and the reason is the following:
Consider the simplest case of an atom with nuclear charge Q, radius r, and subjected to an electric field, E, which creates a force QE, and displaces the nucleus by a distance d. According to Gauss’ law this latter force is given by:
(dQ^2)/(4πεr^3) = QE (Hey! WordPress! I could really use an equation editor in here!)
if the polarizability is defined by Qd/E then we can rearrange the previous equation and yield
α = 4πεr^3 which in atomic units yields volume units, r^3, since 4πε = 1. This is why polarizabilities are usually referred to as ‘polarizability volumes’.
****THIS POST IS STILL IN PROGRESS. WILL COMPLETE IT IN SHORT. SORRY FOR ANY INCONVENIENCE****
The use of double zeta quality basis sets is paramount but it also makes these calculations more time consuming. Polarization functions on the basis set functions are a requirement for good results.
As usual, please rate/comment/share this post if you found it useful or if you think someone else might find it useful. Thanks for reading!
I am frequently asked how to include an extra set of basis functions in a calculation or how to use an entirely external basis set. Sometimes this question also implies the explicit declaration of an external pseudopotential or Effective Core Potential (ECP).
New basis sets and ECPs are published continuously in specialized journals all the time. The same happens with functionals for DFT calculations. The format in which they are published is free and usually only a list of coefficients and exponents are shown and one has to figure out how to introduce it in ones calculation. The EMSL Basis Set Exchange site helps you get it right! It has a clickable periodic table and a list of many (not all) different basis sets at the left side. Below the periodic table there is a menu from which one can select which program we want our basis set for; finally we click on “get basis set” and a pop-up window shows the result in the selected format along with the corresponding references for citation. A multiple query can be performed by selecting more than one element on the table, which generates a list that almost sure can be used as input without further manipulations. Dr. David Feller is to be thanked for leading the creation of this repository. More on the history and mission of the EMSL can be found on their About page. Because of my experience, the rest of the post addresses the inclusion of external basis sets in Gaussian, other programs such as NwChem will be addressed in a different post soon.
The correct format for inclusion of an external basis set is exemplified below with the inclusion of the 3-21G basis set for Carbon as obtained from the EMSL Basis Set Exchange site (blank lines are marked explicitly just to emphasize their location:
spin multiplicity Molecular coordinates - blank line - C 0 S 3 1.00 172.2560000 0.0617669 25.9109000 0.3587940 5.5333500 0.7007130 SP 2 1.00 3.6649800 -0.3958970 0.2364600 0.7705450 1.2158400 0.8606190 SP 1 1.00 0.1958570 1.0000000 1.0000000 **** - blank line -
The use of four stars ‘****’ is mandatory to indicate the end of the basis set specification for any given atom. If a basis set is to be declared for a second atom, it should be included after the **** line without any blank line in between.
WARNING! Sometimes we can find more than one basis set in a single file this is due to different representations, spherical or cartesian basis sets. Gaussian by default uses cartesian (5D,7F) functions. Pure gaussian use 6 functions for d-type orbitals and 10 for f-type orbitals (6D, 10F). Calculations must be consistent throughout, hence all basis functions should be either cartesian or pure.
Inclusion of a pseudopotential allows for more computational resources to be used for calculation of the electronic structure of the valence shell by replacing the inner electrons for a set of functions which simulate the presence of these and their effect (such as shielding) on the valence electrons. There are full core pseudopotentialas, which replace the entire core (kernel). There are also medium core pseudopotentials which only replace the previous kernel to the full one, allowing for the outermost core electrons to be explicitly calculated. The correct inclusion of a pseudopotential is shown below exemplified by the LANL2DZ ECP by Hay and Wadt for the Chlorine atom.
spin multiplicity Molecular coordinates - blank line - basis set for atom1 **** basis set for atom2 (if there is any) **** - blank line - CL 0 CL-ECP 2 10 d potential 5 1 94.8130000 -10.0000000 2 165.6440000 66.2729170 2 30.8317000 -28.9685950 2 10.5841000 -12.8663370 2 3.7704000 -1.7102170 s-d potential 5 0 128.8391000 3.0000000 1 120.3786000 12.8528510 2 63.5622000 275.6723980 2 18.0695000 115.6777120 2 3.8142000 35.0606090 p-d potential 6 0 216.5263000 5.0000000 1 46.5723000 7.4794860 2 147.4685000 613.0320000 2 48.9869000 280.8006850 2 13.2096000 107.8788240 2 3.1831000 15.3439560
If a second ECP is to be introduced, it should be placed right after the first one without any blank line! If a blank line is detected then the program will assume it’s done reading all ECPs and Basis Sets.
Finally, here is an example of a combination of both keywords. If a second ECP was needed then we’d place it at the end of the first one without a blank line. The molecule is any given chlorinated hydrocarbon (H, C and Cl atoms exclusively)
#P B3LYP/gen pseudo=read ADDITIONAL-KEYWORDS - blank line - 0 1 Molecular Coordinates - blank line - H 0 S 3 1.00 19.2384000 0.0328280 2.8987000 0.2312040 0.6535000 0.8172260 S 1 1.00 0.1776000 1.0000000 **** C 0 S 7 1.00 4233.0000000 0.0012200 634.9000000 0.0093420 146.1000000 0.0454520 42.5000000 0.1546570 14.1900000 0.3588660 5.1480000 0.4386320 1.9670000 0.1459180 S 2 1.00 5.1480000 -0.1683670 0.4962000 1.0600910 S 1 1.00 0.1533000 1.0000000 P 4 1.00 18.1600000 0.0185390 3.9860000 0.1154360 1.1430000 0.3861880 0.3594000 0.6401140 P 1 1.00 0.1146000 1.0000000 **** Cl 0 S 2 1.00 2.2310000 -0.4900589 0.4720000 1.2542684 S 1 1.00 0.1631000 1.0000000 P 2 1.00 6.2960000 -0.0635641 0.6333000 1.0141355 P 1 1.00 0.1819000 1.0000000 **** - blank line - CL 0 CL-ECP 2 10 d potential 5 1 94.8130000 -10.0000000 2 165.6440000 66.2729170 2 30.8317000 -28.9685950 2 10.5841000 -12.8663370 2 3.7704000 -1.7102170 s-d potential 5 0 128.8391000 3.0000000 1 120.3786000 12.8528510 2 63.5622000 275.6723980 2 18.0695000 115.6777120 2 3.8142000 35.0606090 p-d potential 6 0 216.5263000 5.0000000 1 46.5723000 7.4794860 2 147.4685000 613.0320000 2 48.9869000 280.8006850 2 13.2096000 107.8788240 2 3.1831000 15.3439560 - blank line -
If you like this post or found it useful please leave a comment, share it or just give it a like. It is as much fun to find out people is reading as it is finding the answer to ones questions in someone else’s blog 🙂