Development of Parameters for Modeling Silicon-Oxygen Polymers in AMBER Molecular Dynamics Simulations

Due to the fact that silcon-containing molecules are being used more frequently in materials chemistry, it is necessary to extend current organic modeling tools to incorporate silicon. Specifically, this research aims to develop accurate parameters for modeling silicon in the AMBER molecular dynamics program. This was accomplished by first obtaining accurate initial values for force constants and equilibrium values via ab initio calculations using the Gaussian quantum mechanics program. We investigated various methods including PM3,DFT/B3-LYP,MP2,and MP4, also basis sets including 6-31G(d),6-311+G(d,p), cc-pVDZ, and cc-pVTZ. After comparison to experimental vibrational spectra, the MP2/cc-pVTZ model chemistry was chosen for the majority of calcuations. Then these initial force-field values were optimized within the AMBER program to match experimental IR spectrum measurements. Van der Waal's parameters will be optimized in OPLS-fashion by matching the AMBER MD simulation liquid properties to experimental data.

The nature of a polymer containing a silicon-oxygen backbone (silicones) is unique in that it contains the ability to act in many ways like a carbon backbone, yet in instances to behave entirely differently. One property that silicone-based polymers have is their ability to be in a liquid state for any length of polymer. This property allows liquid-based dynamic simulations to be performed on polymers of great lengths. These long chained silicone-based polymers, with hydrophilic and hydrophobic organic sidechains, can fold and bend fluidly to form structures whose properties are based largely on the properties of the sidechains. Therefore, finding a molecular dynamics modeling tool which can accurately model organic sidechains in the liquid state becomes necessary. Also, in order for the polymer model to be created efficiently, a program which allows creation of residues (repeat groups for the polymer) and easy linking of a large number of residues would be beneficial. The AMBER Force Field¹ and associated molecular dynamics modeling programs were chosen for their accuracy and efficiency, yet in order to use AMBER, the Force Field parameters for the silicon-oxygen backbone must be derived.

An example of a silicone-based monomer is shown below. With a slightly more complex molecule shown from VMD² the complexity of the folding of the monomer can be seen. The Yellow atoms represent silicon. The addition of a hydrophilic group and hydrophobic group to the molecule creates a folding in on itself. It is theorized that this polymer will, in aqueous solution, form a nanoparticle³. This is just one example of a use for molecular dynamics involving covalently bonded silicon.

Necessary Parameters Bond Parameters Si-O Si-C Si-H Angle Parameters C-Si-C C-Si-O C-Si-H O-Si-H O-Si-O Si-O-Si Torsional Angles Van der Waals Parameters

For easy reference to the parameter variables to be developed. Below is the the generic form of the AMBER Potential Energy equation¹.

The parameters used to model the potential for bond length and angle can be seen above to be of the form of a spring equation. Since the spring equation is simply a parabolic function, this is an ideal first step in parameterization.

Strategy: In order to derive a good estimate of the parameters for the spring equation, first Gaussian98⁴ is utilized to obtain ab initio calculation estimates. Then, a genetic algorithm will be incorporated in order to refine those values to match reality in AMBER, where possible.

The Gaussian98 program has the ability to do large Quantum Mechanical (QM) calculations, but at a great temperal cost. Obviously then, whole polymers cannot be used to calculate the initial parameters. The molecules used, and the initial parameters derived from them, are shown below.

Using the Gaussian98 program, it is possible to equilibriate the a structure using the geometry optimization calculation. This equilibriated structure determined the initial equilibrium bond lengths and angles (r0, or &theta0 from the AMBER potential function). Then, using the 'scan' keyword, it is possible to vary the parameter in question and get each new molecule's respective energies. These energies can then be plotted vs. the bond lenght/angle to reveal a very parabolic set of points, as shown to the right. This can then be fit using a least-squares regression technique found in the program 'gnuplot' to determine the force constant for the spring equation (kr or k&theta from the AMBER potential function).

The plot above shows the energies (hartrees) vs. various angles of H-Si-H in Silane.

For this method of determining the initial parameters, many different QM methods were tried including PM3, DFT/B3-LYP, MP2, and MP4, also basis sets including 6-31G(d), 6-311+G(d,p), cc-pVDZ, and cc-pVTZ. After comparison to vibrational spectra, B3-LYP and MP2 with higher level basis sets, such as 6-311+G(d,p), cc-pVDZ, and cc-pVTZ proved to be fairly accurate based off of Aldrich and NIST values. The comparison of various models is shown below. The values differed slightly. For the majority of calculations, the chemical model MP2/cc-pVTZ was used for the hight level accuracy and the experienced stability with the molucules in question. (results shown in conclusion)

Now that each of the spring parameters have been given an initial value, it is important to match these values to reality. Using a genetic algorithm program, parmscan⁵, the initial values can be entered into a script. This genetic algorithm uses the initial values, restraints, and allowable ranges for parameters to alter the parameters to best match experimental data. For the purposes of determining the bond length and angle parameters, using vibrational frequencies was deemed most appropriate. Therefore, by entering in IR-spectrum values, parmscan will run variations of the molecules in question through nmode calculations with the AMBER force field to creat a best fit of parameters. This program can also be set to match a number of other physical constants as well. These options are discussed under van der Waals: Matching Reality.

Development of van der Waals (vdW) parameters is the greatest obstacle in parameter development. The vdW initial parameters were derived by using rough estimates based off of other silicon vdW functions in the MMX, MM3, and MMFF. By fitting these functions with the AMBER Lennard-Jones 6-12 equation, a rough estimate can be determined. Since the vdW parameters are frequently considered a "catch-all" for adjusting errors within the force field parameter set, then a certain amount of disagreement is to be expected. In order to account for systematic differences between force fields, this same method was applied to carbon, sulfur, phosphorous and oxygen vdW parameters and a mean value for silicon was extrapolated with a certain amount of error removed. The table below reveals the data used.

Just as mentioned under the Spring Parameters section, parmscan can be used to optimize the vdW parameters as well. Parmscan has the option of matching vibrational frequencies but also can be used to optimize parameters to match heat capacity, enthalpy of vaporization, entropy or density. As some of these properties are highly sensitive to vdW forces, parmscan can optimize vdW via the same genetic algorithm for optimizing the spring-like parameters. This step in the research has just begun and therefore no results can be revealed currently.

Another method attempted thus far to determine the vdW parameters was based solely off of the initial values from the table above. By using the molecule silane (SiH₄) and the initial parameters for the vdW potential function, along with derived bond and angle parameters from the previously described methods, a periodic boundary system of 4089 silane molecules were then constructed in the LEaP module of the AMBER package. This 'box' of silanes was then attempted to be equilibriated using the AMBER force field at various temperatures. The attempted equilibriation can be seen below.

After trials and various temperatures and various initial allowable box-sizes, it was concluded that this method must be reserved until after better vdW parameters can be found to use. This will be another test to determine the validity of the parameters parmscan derives.

So far, initial values for nearly all of the spring-like parameters have been determined and are shown below. The O-Si-O and the Si-O-Si angles require molecules too large to be optimized in Gaussian98, and therefore will be determined in parmscan, as will the torsions. Current work is now being done to create an efficient way to optimize these parameters with parmscan, as there are many parameters to be optimized. It has also been determined that unlike most protonic bound hydrogens, the silicon is less electronegative than hydrogen and thus forms a hydridic hydrogen. This change in polarity will also be investigated for changes in vdW parameters.

Many thanks to Sherman-Fairchild Foundation and the National Science Foundation for funding as well as the Hope College Department of Chemistry. Also, to Junmei Wang for correspondence and the use of his parmscan algorithm. Many of the computations were performed at the National Center for Supercomputing Applications (NCSA) at UIUC. Thanks are also due to them for computer time. Thanks to Matthew Zwier for help in the creation of this site.

References and Resources:

1. AMBER Molecular Dynamics Package
2. VMD - Visualizing Molecular Dynamics
3. Research of Dr. Michael Silver of the Hope College Chemistry Department - contact info upon request
4. Gaussian98 Quantum Mechanics Program
5. Wang, Junmei, and Peter A. Kollman. J Comput Chem 2001, 22, 12, 1228.