Physical Interaction plugins calculate the physical interactions between atomic and molecular entities, such as attraction and repulsion.
1.1. AIREBO
The Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential is an extension of the REBO potential described in the REBO_MBM plugin description. The extensions that the AIREBO potential provides are:
Non-bonded, intermolecular interactions, modelled with a Lennard-Jones (LJ) 12-6 potential
Four-body torsional interaction, modelled with a torsional potential with a single minimum
The resulting system energy is then giving by
Equation 5.1. The AIREBO Potential
The reference and citation to use for this AIREBO plugin and any publications that result from simulations with this potential is
S. J. Stuart, A. B. Tutein and J. A. Harrison, "A Reactive Potential for Hydrocarbons with Intermolecular Interactions", J. Chem. Phys.112 6472-6486 (2000).
1.1.1. Configuration
The following is an example of how to configure the AIREBO plugin in the simulation specification file.
This refers to the Δt/τT term in the van Gunstern-Berendsen Thermostat used by the plugin. The range is [0.01 < Δt/τT < 1.2] where 1.2 greatly dampens any thermal activity in the system, and 0.01 barely dampens thermal activity. A value of 1.0 seems to work well.
1.1.2. Operation Notes
Parameterized for hydrocarbons only. Additional parameterizations will be made available in future releases.
1.2. BondCalculator
This plugin discovers molecules in the simulation space and specifies bonds to describe them. It uses a simple scheme based on the distance between atoms to decide bonds. The bond length data comes from here [chemviz.ncsa.uiuc.edu].
1.3. MPQC_SClib
This quantum mechanics plugin encapsulates the functionality of the Scientific Computing Toolkit (SC) used by the Massively Parallel Quantum Chemistry Program (MPQC), hence MPQC_SClib. The following levels of theory are supported:
Closed shell, unrestricted, and general restricted open shell Hartree-Fock energies and gradients
Closed shell, unrestricted, and general restricted open shell Kohn-Sham density functional theory (DFT) energies and gradients
Second order closed shell Møller-Plesset perturbation theory energies and gradients.
Currently, the plugin can be used to calculate system energies, gradient dynamics, and electrostatic potential (ESP) planes.
1.3.1. Common Configuration
The following is an example of the common elements of each MPQC_SClib configuration in the simulation specification file.
The SClib component of the plugin writes status and result information as it runs. This specifies the directory where SClib should write its output file (SClib_out.txt).
Note
To watch detailed output from the MPQC_SClib plugin as it is running, execute (from a Unix command line) tail -f <logDirectory>/SClib_out.txt
dataDirectory
Where to find the MPQC_SClib atomic information and basis set data files. All distributions come with a set under the data/MPQC_SClib directory as shown in the example configuration.
multiplicity
Multiplicity is determined by the number of unpaired electrons in the system you're calculating. Here's a table showing the multiplicity depending on the number of unpaired electrons:
Unpaired Electrons
Multiplicity
0
1 (default)
1
2
2
3
3
4
Note, however, that you need not calculate with the maximum multiplicity as shown in the table, you can get away with multiplicities of either 1 or 2 for every system. Choose multiplicity=1 for systems with an even number of unpaired electrons, choose multiplicity=2 for systems with an odd number of unpaired electrons.
You can further reduce things if you don't happen to know how many unpaired electrons there are, as long as you know how many electrons there are. So choose multiplicity=1 for systems with an even number of electrons, choose multiplicity=2 for systems with an odd number of electrons.
method
The method to use, note that the effective multiplicity will ultimately determine the specific method. Possible values are
HF - general Hartree-Fock. With multiplicity=1, closed-shell Hartree-Fock is used, with multiplicity>1, unrestricted Hartree-Fock is used.
RHF - restricted Hartree-Fock. With multiplicity=1, closed-shell Hartree-Fock is used, with multiplicity>1, high-spin open-shell Hartree-Fock is used.
UHF - unrestricted Hartree-Fock regardless of multiplicity.
KS - Kohn-Sham (DFT). With multiplicity=1, closed-shell Kohn-Sham is used, with multiplicity>1, unrestricted Kohn-Sham is used. Note that a functional must be specified with this method.
RKS - restricted Kohn-Sham (DFT). With multiplicity=1, closed-shell Kohn-Sham is used, with multiplicity>1, high-spin open-shell Kohn-Sham is used. Note that a functional must be specified with this method.
UKS - unrestricted Kohn-Sham (DFT) regardless of multiplicity. Note that a functional must be specified with this method.
MP2 - second order closed shell Møller-Plesset perturbation theory. multiplicity must be 1.
basis
The Gaussian basis set to use. The following is a table of available basis sets (from the SClib documentation) listing the supported elements for each basis and the number of basis functions for H (n0), first row (n1), and second row (n2) atoms.
Basis Set
Elements
n0
n1
n2
STO-2G
H-Ca
1
5
9
STO-3G
H-Kr
1
5
9
STO-3G*
H-Ar
1
5
14
STO-6G
H-Be, C-Kr
1
5
9
MINI (Huzinaga)
H-Ca
1
5
9
MINI (Scaled)
H-Ca
1
5
9
MIDI (Huzinaga)
H-Na
2
9
DZ (Dunning)
H, Li, B-Ne, Al-Cl
2
10
18
DZP (Dunning)
H, Li, B-Ne, Al-Cl
5
16
24
DZP + Diffuse (Dunning)
H, B-Ne
6
19
3-21G
H-Kr
2
9
13
3-21G*
H-Ar
2
9
19
3-21++G
H-Ar
3
13
17
3-21++G*
H-Ar
3
13
23
4-31G
H-Ne, P-Cl
2
9
13
4-31G*
H-Ne, P-Cl
2
15
19
4-31G**
H-Ne, P-Cl
5
15
19
6-31G
H-Ar
2
9
13
6-31G*
H-Ar
2
15
19
6-31G**
H-Ar
5
15
19
6-31+G*
H-Ar
2
19
23
6-31++G
H-Ar
3
13
17
6-31++G*
H-Ar
3
19
23
6-31++G**
H-Ar
6
19
23
6-311G
H-Ar, Ga-Kr
3
13
21
6-311G*
H-Ar, Ga-Kr
3
18
26
6-311G**
H-Ar, Ga-Kr
6
18
26
6-311G(2df,2pd)
H-Ne
14
30
6-311++G**
H-Ne
7
23
6-311++G(2d,2p)
H-Ne
10
29
6-311++G(3df,3pd)
H-Ar
19
42
50
cc-pVDZ
H, He, B-Ne, Al-Ar, Ga-Kr
5
14
18
cc-pVTZ
H, He, B-Ne, Al-Ar, Ga-Se
14
30
34
cc-pVQZ
H, He, B-Ne, Al-Ar
30
55
59
cc-pV5Z
H-Ne, Al-Ar
55
91
95
aug-cc-pVDZ
H, He, B-Ne, Al-Ar
9
23
27
aug-cc-pVTZ
H, He, B-Ne, Al-Ar
23
46
50
aug-cc-pVQZ
H, He, B-Ne, Al-Ar
46
80
84
aug-cc-pV5Z
H, He, B-Ne, Al-Ar
80
127
131
cc-pCVDZ
B-Ne
18
cc-pCVTZ
B-Ne
43
cc-pCVQZ
B-Ne
84
cc-pCV5Z
B-Ne
145
aug-cc-pCVDZ
B-F
27
aug-cc-pCVTZ
B-Ne
59
aug-cc-pCVQZ
B-Ne
109
aug-cc-pCV5Z
B-F
181
NASA Ames ANO
H, B-Ne, Al, P, Ti, Fe, Ni
30
55
59
functional
The exchange-correlation functional to use for DFT methods. The table below (from the SClib documentation) lists the functional names and the equivalent functionals in other packages. The Description column gives the sub-functionals used to build up the functional and its coefficient if it is other than one. The G98 column lists the equivalent functional in Gaussian 98 A.6. The NWChem column lists the equivalent functional in NWChem 3.3.1.
Specifies the grid to be used for numerical integrations. The following values can be given: xcoarse, coarse, medium, fine, xfine, and ultrafine. Default: fine
maxEnergyConvergenceAttempts
Specifies the maximum number of attempts that should be made to achieve the desired energy calculation accuracy (desiredEnergyAccuracy.) Default: -1 indicating no limit to the number of attempts.
desiredEnergyAccuracy
Specifies the desired energy calculation accuracy in Hartrees (1 Hartree = 4.359 744 17(75) x 10-18 Joules.) Default: 1.0e-6
1.3.2. What Level of Theory Should I Use?
Quantum chemistry methods is a large body of science, but here's a fast break-down, from least to most accurate, with respect to the MPQC_SClib plugin.
The most basic method is the Hartree-Fock theory - each electron sees the average effect of all the other electrons.
Density Functional Theory (DFT) - the exact electron density can be determined by finding the density that minimizes a density functional. The catch is that the exact form of the functional isn't known, but there's several functionals available to approximate it.
Correlated methods - electrons respond to the dynamic motion of each of the other electrons, not just the average effect. Møller-Plesset perturbation theory is an example of this method.
The accuracy trade-off is computational intensity, ie, time.
Each of the methods listed above utilize a set of basis functions called a basis set. The basis functions describe orbitals, which are descriptions of the behavior of an electron in an atom according to quantum mechanics. The more basis functions there are in a basis set, the more accurate the description of the orbitals, and the more costly the calculation.
Note
Each basis set only has functionals for a sub-set of atoms.
Below is a brief comparison of various levels of theory used to calculate the energy of an H2O molecule, including some notes
[1]
about the basis set used.
Table 5.3. MPQC_SClib Comparison of Various Levels of Theory
Method
Basis Set (# basis functions)
Basis Set Notes
Functional
Calculation Time
Energy (Hartrees)
HF
3-21G* (30)
A small basis set.
n/a
< 1 s
-75.5788429667
HF
6-311G** (50)
Usually acceptable for HF and DFT.
n/a
2 s
-76.0314305762
KS
3-21G* (30)
A small basis set.
B3LYP
32 s
-75.9719803768
KS
6-311G** (50)
Usually acceptable for HF and DFT.
B3LYP
41 s
-76.4390839802
MP2
cc-pVTZ (78)
Designed for correlated methods such as MP2.
n/a
37 s
-76.309362958967
MP2
cc-pVQZ (144)
n/a
9 min 40 s
-76.338134762910
1.3.3. MPQC_SClib Gradient Dynamics Configuration
Gradient dynamics uses the calculated energy gradient to determine atom velocities. It's very similar to geometry optimization, but allows for pathing, ie, the controlled motion of molecular tools, to act on the system. The Cytosine-Guanine Example below uses gradient dynamics where the guanine molecule is "pushed" towards the cytosine molecule and subsequently hydrogen bonding takes place.
The following fragment is an example of how to configure the MPQC_SClib plugin for gradient dynamics in the simulation specification file.
Whether or not to perform gradient dynamics on the system. Either "yes" or "no".
deltaTbyTau
This refers to the Δt/τT term in the van Gunstern-Berendsen Thermostat used by the plugin. The range is [0.01 < Δt/τT < 1.2] where 1.2 greatly dampens any thermal activity in the system, and 0.01 barely dampens thermal activity. A value of 1.0 seems to work well.
1.3.4. Cytosine-Guanine Example
In this example we use the MPQC_SClib gradient dynamics mechanism to show the hydrogen bonding of a DNA base-pair: cytosine and guanine. The cytosine molecule is briefly pushed towards the guanine molecule to give it some momentum in the right direction.
All the files for this simulation are included in the installed cytosine-guanine example, including the resultant HiveKeeper visualization file. You can open the simulation with HiveKeeper and play the visualization (without having to re-run the simulation, which would take a few days.) While playing the visualization in HiveKeeper, you can rotate the system and see the characteristic "twisting" of the two molecules as they bond.
A video of the simulation is also available here: [MPEG] 688 Kb
The potential energy associated with the charge distribution of a molecule is the electrostatic potential (ESP). The ESP at some point due to a point charge is given by
Equation 5.2. Electrostatic Potential
where
E
is the electrostatic potential
k
is (4πεo)-1
q
is the value of the point charge
r
is the distance between q and the point where E is being calculated
In our molecules and systems we have many "point" charges comprising protons and electron densities. The total electrostatic potential at a point may be found by measuring the strengths of the individual electrostatic potentials due to each of the individual charges and then summing them. To form a visualization of ESP we assign a color to each point depending on how positive or negative the ESP is at that point.
The image to the right shows the ESP of a water molecule, which has two very polar covalent bonds. The ESP coloring shows how the majority of electron density (blue color) resides on the oxygen atom, with a lack of electron density (red color) on the hydrogen atoms. The white color in this visualization is used to show high potentials - the white inside the red is interpreted as very positive potential.
The ESP of a molecule can help give an understanding of such properties as its melting/boiling point, solubility, acidity, and current flow characteristics
[2].
The ESP Plane and the ESP Window
The MPQC_SClib plugin can calculate the ESP for a bounded window lying on a plane in space. While the bounded window is not a plane, if the window is big enough, everything outside the window, along the plane, would effectively have insignificant ESP. And so we refer to the window as a plane in some sense.
The following fragment is an example of how to configure the MPQC_SClib plugin for ESP plane output in the simulation specification file.
This specifies the resolution of the window. The total number of ESP data points in the window will number resolution2.
centerPoint
Specifies the point where the ESP window shall be centered at.
normalPoint
Combined with the centerPoint, specifies the orientation of the ESP window.
outputLength
The length of one side of the ESP window. The area of the window will be outputLength2
cutoffHeight
The distance above and below the window in which to include atoms for ESP calculation.
cutoffWidth
The distance around the window in which to include atoms for ESP calculation.
All measurements are in meters.
This illustration shows the parameters and how they specify the ESP window.
1.3.6. Electrostatic Potential of a 5,5 Carbon Nanotube Example
Using the MPQC_SClib plugin, we specified a flat square that visually encodes the electrostatic potential (ESP) of all atoms, within a specified range of the square, as color - red for positive ESP and blue for negative. We then pulled a 5,5 carbon nanotube through the square to observe its ESP.
ESP Window and Nanotube Path
In our example, we pull a nanotube through a window that shows the ESP of everything passing through it. The image to the left shows the window and how the nanotube will be pulled through it.
The dynamics of the nanotube itself will be simulated with the REBO_MBM plugin.
Here is an MPEG showing the path of the nanotube only (no ESP visualization.)
Nanotube ESP
The figure to the right shows the ESP of a slice of the nanotube as it begins to travel through the ESP window. We can see positive ESP around the pairs of carbon atoms at level with the window, and negative ESP radiating out from in between carbon pairs. The ESP is not symmetric because the nanotube is at an angle to the ESP window, as is shown in the above diagram. Note that there is no negative ESP inside the nanotube.
This plugin is being deprecated - use the AIREBO plugin instead.
This plugin employs D. Brenner et al's second-generation Reactive Empirical Bond Order (REBO) potential energy function for solid carbon and hydrocarbon molecules
[3]
. This potential allows for covalent bond breaking and forming with associated changes in atomic hybridization. It has been used to model many processes including reactive surface and cluster dynamics
[4]
, carbon nanotube properties
[5]
, polycrystalline diamond structure
[6]
, brittle fracture dynamics
[7]
, and processes associated with nanoindentation
[8]
.
The general form of this potential, originally derived by Abell from chemical pseudopotential theory
[9]
is:
Equation 5.3. REBO Potential
where
Eb
is the chemical binding energy
i, j
are nearest neighbor atoms
rij
is the distance between atoms i and j
V R
is the interatomic, pair-additive repulsive function (core-core, etc)
V A
is the interatomic, pair-additive attractive function (from valence electrons, etc)
bij
is a bond order between atoms i and j
Tersoff introduced analytic parameterized forms for the bond order that were able to describe bonding for a number elements in various phases
[10]
. The value of the bond order was assumed to depend both on the local coordination and on bond angles.
Brenner, et al, established that Tersoff's solid-state scheme can be applied to both solid-state and molecular systems
[11]
, and developed the first-generation form of the potential
[12]
. The second-generation form contains improved analytic functions and an extended database relative to the earlier version.
This plugin was developed with code from P. McCluskey's C port of Brenner et al's FORTRAN version, which incorporated C code from J. Michelsen, hence the name REBO_MBM. The license for some of the code used for the plugin asks me to give the following acknowledgement:
"This product includes software developed by Zyvex LLC."
1.4.1. Configuration
The following is an example of how to configure the REBO_MBM plugin in the simulation specification file.
Where to find the REBO spline data files. All distributions come with a set under the data/REBO_MBM directory as shown in the example configuration.
1.4.2. Operation Notes
This plugin primarily supports molecules made up totally of either Carbon atoms, or Carbon and Hydrogen atoms. The minimal Silicon atom support is best for positioning Carbon atoms in a place where they will form stronger bonds with other Carbon atoms. Non-C/H/Si atoms are ignored completely.
The original code this plugin uses has a hard-coded value of 15,000 for the number of atoms it can handle. Work will be done to remove this limit for the next version.
Works with RC_Traverser
Updates the totalEnergy, kineticEnergy, potentialEnergy, and idealTemp entries in the MeasurementSet.
1.4.3. Mechanosynthesis of Diamond Simulation Example
Using the REBO_MBM Interaction Plugin described above, D. Huang et al's proposed reaction sequence of epitaxial diamond growth
[13]
, Sinnott et al's hydrogen abstraction
[14]
, and adaptations of Ralph Merkle's proposed hydrocarbon assembly tools
[15]
, we created an example of the mechanosynthesis of diamond with the simulator.
Diamond Seed
We start with a small molecule of diamond.
Carbon atoms are yellow; hydrogen atoms are white
The diamond seed is anchored by six of its hydrogens, indicated by green stars. The ability to anchor atoms in space is a feature of NanoHive-1.
The brighter carbon atoms indicate the primary initial diamond structure that will be extended
The first phase of the mechanosynthesis is to abstract the two hydrogens, indicated as 1 and 2, in preparation for the new carbons.
Hydrogen Abstraction
Here we see hydrogen 1 has just been abstracted from the diamond seed.
The upper molecule (less the bottom-most hydrogen, which has just been abstracted) is the hydrogen abstractor
The abstractor was lowered down at 0.2 nm/ps directly over the target hydrogen, until the bottom carbon was 0.16 nm above the target hydrogen. The abstractor was moved by pathing, or controlling the motion of, three of the hydrogens at the top of the abstractor. The ability to path atoms is a feature NanoHive-1.
The bond indicated in green is a double bond between the two carbons
The target hydrogen is abstracted because it forms a stronger bond with the double-bonded carbon, at the end of the abstractor tool, than with the single-bonded carbon in the diamond seed.
Hydrogen 2 is abstracted in the same manner.
Dimer Deposition
This dimer deposition tool is used to deposit the first of two carbon dimers onto the diamond seed.
The blue atoms are silicon
The molecule is really two tools that act as one until the dimer is deposited, at which point the left and right silicon structures act independently to complete the deposition
The silicon-carbon/silicon-hydrogen model used in the REBO_MBM plugin is only an approximation to facilitate the dimer deposition - for this reason, the tool is not hydrogenated like the hydrogen abstractor was
The hydrogen acts as a catalyst to the deposition reaction
Here we see the first dimer has been deposited onto the diamond seed and the two tools are being removed.
The second dimer is deposited in the same manner.
Resulting Diamond Molecule
The second dimer has just been deposited (it did not have an additional hydrogen), and the tool removed, leaving the resulting (larger) diamond molecule.
[1] Janssen, C. L., "Quantum Chemistry on Clusters", ClusterWorld2 18 (August 2004).
[2] Liang, G. C., A. W. Ghosh, M. Paulsson, and S. Datta (2004), Phys. Rev. B69 No.3, 0353XX.
[3] Brenner, D. W., O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott (2002), J. Phys.: Condens. Matter14 783-802.
[4] Sinnot, S. B., L. F. Qi, O. A. Shenderova, and D. W. Brenner (1999), Advances in Classical Trajectory Methods: Molecular Dynamics of Clusters, Surfaces, Liquids and interfaces vol 4, ed W Hase (Stamford, CT: JAI Press) pp 1-26 ch 1
[5] Wijesundara, M. B. J., L. Hanley, B. Ni, S. B. Sinnott (2000), Proc. Natl Acad. Sci. USA97 23.
[6] Shenderova, O. A., D. W. Brenner, and L. H. Yang (1999), Phys. Rev. B60 7043.
[7] Shenderova, O. A., D. W. Brenner, A. Omeltchenko, X. Su, and L. H. Yang (2000), Phys. Rev. B61 3877.
[8] Sinnott, S. B., R. J. Colton, C. T. White, O. A. Shenderova, D. W. Brenner, and J. A. Harrison (1997), J. Vac. Sci. Technol. A15 936.