Using Fragments in Gaussian In this video, we will demonstrate how to
set up fragments in Gaussian.
In doing so, we will model antiferromagnetic coupling using a ferrodoxin compound.
Additionally, we will model phenol dimer by setting up fragments and running a counterpoise
correction calculation.
Additionally, we will examine counterpoise corrections for weakly-bound systems by modeling
phenol dimer.
Ferredoxins are a class of proteins that facilitate electron transfer in organisms.
Our molecule of interest is FeS(SCH3)2 – dimer.
This molecule has an overall charge of -2 and a spin multiplicity of 1.
Although it is a singlet, the charge and spin multiplicities must be localized in order
to generate a reasonable initial wavefunction for the calculations.
Here is our molecule.
Its two central iron atoms are weakly bonded to their respective methylsulfanyl group.
The central iron atoms are also bonded to the central sulfur atoms.
We begin by running an optimization.
Since our molecule is a singlet, this will be a restricted calculation, which means the
molecule will be modeled as a closed-shell species.
When we use the APFD/6-311+G(2d,p) model chemistry to run an Opt+Freq calculation, the calculation
succeeds.
However, although Gaussian located a stationary point, the Results Summary tab indicates that
the structure has two imaginary frequencies, which means that it is a second order saddle
point rather than a minimum.
Running a follow-on stability calculation indicates the molecule has an RHF-to-UHF instability.
Thus, the structure is distorted from the minimum and has an unstable wavefunction.
For a singlet, merely requesting an unrestricted method will still run a restricted, closed-shell
calculation, producing the exact same results.
Instead, we need to specify an initial guess which will result in an accurate representation
of the molecule’s electron configuration in the final results.
While the molecule has an overall charge of -2 and is a singlet, each Fe(III) atom has
five unpaired electrons.
One of them Alpha electrons, and the other Beta electrons.
To represent the localization of charges and spin distribution across the molecule, we
will use fragments to define separate regions with different charge and spin characteristics.
Here is our original molecule.
In order to set up fragments, we navigate to the Tools menu and select "Atom Groups."
Once the Atom Group Editor Tab opens, we select Gaussian Fragment from the Atom Group Class
drop-down menu in the upper left.
We will be focusing on the six highlighted columns.
We will define four fragments: two containing an iron atom and its attached methlysulfanyl
group, and two consisting of a single sulfur atom each.
Currently, all of the atoms are in Gaussian Fragment (1).
To add atoms to a new fragment, we select the atom marquee selection tool and highlight
one iron atom and its respective methlysulfanyl group.
We then click plus sign in the Gaussian Fragment (2) row to add these selected atoms to this
new fragment.
We do this again for the other iron and methlysulfanyl group, adding these atoms to Fragment (3).
To make a fourth fragment group, we navigate to Group Actions in the top right of the Atom
Group Editor panel and select "New Group.”
We keep the default name of Gaussian Fragment (4) and press "Ok."
We then use the single atom selector tool to select a sulfur atom and add it to Gaussian
Fragment (4).
Now that we have created our four fragment groups, we select the highlight boxes so that
we can clearly see each defined group.
We go on to specify the charge, spin multiplicities, and unpaired spins type for each fragment
in the initial guess.
In fragments 1 and 4 ,which contain the sulfur atoms, we assign a charge of -2 and select
Singlet (spin multiplicity of 1).
For fragments 2 and 3, we change the charge to 1 and set the spin multiplicity to Sextet.
Finally, we change fragment 3's Unpaired Spins to beta.
To continue to see the different fragments highlighted after we close out of the Atom
Group Editor Tab, we check the Persistent Visuals box in the bottom left of the panel.
As shown, the four fragments are now clearly highlighted and defined.
Next, we navigate to the Gaussian Calculation Setup Dialog.
In the Charge and Multiplicity section at the top left of the window, the values reflecting
each fragments' individual charge and spin multiplicity have been automatically entered
by GaussView to reflect our fragment definitions.
In the Job Type Panel, we specify the job type as "Energy."
Next, we move to the method panel and change the method to "Unrestricted."
Notice that the APFD keyword at the top of the window now has a "U" at the beginning
to represent the unrestricted method selected.
In the General Panel, we verify that the "Write Gaussian Fragment Data" box is selected.
Moving to the Guess panel, we check "Only do guess (No SCF)" as well as "Use fragments
(atom groups) for generating guess."
When selecting these options, the Keywords have changed to include "guess= only" and
"fragments = 4" to specify that our four defined fragments will be used in generating our wavefunction.
The job is ready to be run.
We run the Guess=Only job as an unrestricted calculation and the job succeeds.
The follow-on stability calculation confirms that the antiferromagnetic singlet wavefunction
is stable and is an appropriate initial guess for this molecule.
We can now run an optimization starting from this wavefunction.
This is the input file for that job.
It reads the wavefunction, molecule specification and initial guess from the fragment guess
job’s checkpoint file, and then runs the optimization.
The optimization successfully finds a minimum.
This table summarizes the results of our various calculations.
The minimum has an energy that is .2 hartrees lower than the second order saddle point,
and the iron sulfur bond lengths are significantly longer.
We will next visualizing the spin density.
We open the Results menu and select "Surfaces/Contours."
We create a New Cube of Type "Spin Density."
Once our new cube is available, we set the isovalue in the Density field below the surface
list to 0.004.
Then we go to the Surface Actions menu and choose "New Surface."
The spin density surface will be displayed.
The spin density is mainly on the iron atoms, with a lesser amount on the methylated sulfur
atoms.
For a more quantitative view of the spin density, we examine the Mulliken spin density results
using the Charge Distribution visualization feature.
We check Show Numbers to include the numeric value in the View window.
We see that the predicted spin density on the iron atoms is about plus-or-minus 3.7
and that there is a small amount of spin density of the same polarity on the bonded methysulfanyl
sulfur atoms.
The spin density on the central sulfur atoms is 0.
Another use of fragments is for evaluating the importance of basis set superposition
error (if any) for weakly-bound complexes.
In a weakly-bound complex, basis set superposition error occurs when basis functions on one molecule
overlap with ones on another molecule.
For example, in a dimer where the monomers are at close intermolecular distances, basis
functions from one monomer can inadvertently describe the other.
This overlap results in an artificial--and therefore erroneous--stabilization of the
molecular structure.
A counterpoise calculation computes a correction for the basis set superposition error.
We will use phenol dimer as our example molecule.
We will run a counterpoise correction calculation for the APFD method using the 3-21G basis
set.
We choose this basis set, which is of course far smaller than what would ever be used in
a production calculation, in order to illustrate an extreme case of this phenomenon.
In order to do so, we must again create fragments.
This time, we only need two fragments, one for each phenol group.
We open the Atom Group Editor.
All atoms are automatically placed in Gaussian Fragment (1).
We again use the atom marquee selector tool to select all of the atoms in one phenol group,
which we add to Gaussian Fragment (2).
In the Gaussian Calculation Setup Window, we specify our job type as Opt+Freq.
In the Method panel, we select the APFD functional and select the 3-21G Basis Set.
In the General tab, we verify that "Use Counterpoise" is checked.
Notice that in the Keywords section at the top of the window, the counterpoise parameter
is set as two to because there are two fragments.
When the job completes, The counterpoise-corrected energy can be found in the output file.
The first highlighted line contains the predicted total energy including the counterpoise correction.
The following line contains the counterpoise correction amount itself.
The final 2 lines in this excerpt report the predicted complexation energy with and without
the counterpoise correction.
We can see that the difference is over 6 kcal/mol.
The predicted value for the angle between the two rings is 67.6 degrees, which we determined
by examining the C-O-O-C dihedral angle.
When optimized with the same model chemistry without a counterpoise correction, the angle
in the resulting structure was 95.6 degrees.
The predicted binding energy including the counterpoise correction is about -12 kcal/mol.
The counterpoise correction makes up 55% of that value.
In general, the larger the basis set, the less need for counterpoise correction.
When using the larger 6-31+G(d,p) basis set, the energy correction decreases significantly,
and it is a much smaller percentage of the value.
The difference in ring angles with and without counterpoise correction also narrows.
Finally, our usual 6-311+G(2d,p) basis set has an even smaller counterpoise correction
and ring angles that are very close with and without it..
Basis set superposition error is not a significant factor at this level of theory.
In summary, use fragments to construct an appropriate initial guess for molecular systems
with non-standard electronic configurations.
You will need to run a sequence of jobs: a fragment guess calculation, a stability calculation
to verify the stability of the wavefunction, and then finally the jobs you need to predict
the items of interest.
You should compute counterpoise corrections for weakly-bound complexes when you cannot
afford to model them with a large basis set.
You can use the results as a diagnostic to determine whether basis set superposition
error is a significant factor for this system and model chemistry.