Modeling the Behavior of Siloxane Polymers as Emulsifiers using Molecular Dynamics Simulations:

And Development of AMBER Parameters

Justin M. Shorb

Dr. Brent P. Krueger

Dr. Michael E. Silver

Hope College, Holland, MI 49423

February, 2004


Due to the fact that silicon-containing molecules are being used more frequently in organic chemistry, it is necessary to extend current organic modeling tools to incorporate silicon. Specifically, this research aims to model a new class of siloxane polymers with a silicon-oxygen backbone that behave as emulsifiers in aqueous solution. Current work is being done using molecular dynamics packages that include MMFF94 and MM2, two force fields that are already parameterized for silicon. These results provide initial insight into the behavior of these novel materials as well as comparison for the development of parameters specific to modeling silicon-oxygen polymers in the AMBER molecular dynamics program. AMBER was chosen due to its widespread use and success in modeling the aqueous behavior of biopolymers.; We have obtained accurate starting values via ab initio ;calculations using the Gaussian98 quantum mechanics program. Using these initial values, molecular dynamics parameters have been optimized to match experimental IR spectral measurements using a genetic algorithm program: PARMSCAN.



The nature of a polymer containing a silicon-oxygen backbone (siloxanes) 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 Field1 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.


Modeling Techniques

Currently, work is being done to model the siloxanes in various silicon parameter containing force fields for two initial goals: primarily to begin modeling the novel siloxane polymer of Dr. Michael Silver and secondly to help verify the work done to develop AMBER's parameter set. At the present time little work has been done with the other force fields, but below can be seen a sample of a molecule minimized in vacuum using the MOE modeling package6; note the position of the sidechains.Work is being done to choose an appropriate force field for investigation as well as determining suitable environment for modeling.

Here, the darker gray atoms represent silicon and the lighter gray atoms are carbon (solely in sidechains). The smaller picture is a representation of the molecule prior to minimization where the sidechains are more clearly visible. For those with a keen eye, it is apparent that the two pictures actually do not contain the same sidechains. This is due to the project still being at its early stages. Looking at the denatured strand at the right, it is clear,however, that there is a correlation to the structure of a polypeptide. The current efforts focus on determining the best way to integrate these new 'psuedo-peptides,' or siloxane monomers, into current protein building software in order to create an easy to build and manipulate structure compatible with various other biomodeling software.


AMBER Force Field and Parameters

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

An example of a silicone-based monomer is shown below. With a slightly more complex molecule shown from VMD2 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 nanoparticle3. 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


The Spring Parameters: Bond Length and Angle

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 Gaussian984 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.


Sample Molecules

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.

  • Si-O

  • Si-C

  • O-Si-C

  • C-Si-C

  • Si-H

  • O-Si-H

  • C-Si-H

  • O-Si-O

  • Si-O-Si

  • Torsional

  • van der Waals


Initial Parameters from ab initio Calculations

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 θ0 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θ 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)

Comparison of Computational Models


B3LYP 6-311*G

B3LYP 6-31G



MP2 6-311*G

MP2 6-31G

MP2 cc-pVDZ

MP2 cc-pVTZ

Si-H Bond Length









H-Si-H Angle









Energy (Hartrees)









Energy (Kcal/mol)










Matching Reality with a Genetic Algorithm: PARMSCAN

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, parmscan5, the initial values can be entered into a script. Currently, the PARMSCAN program is being used to match AMBER derived normal modes to the IR spectrum calculated by Gaussian with B3LYP/6-311**.

Move mouse over spectra.

The red line indicates parmscan derived; the blue line indicates B3LYP, scaled from literature. Note that these are preliminary results. PARMSCAN weights error linearly, therefore the largest errors are found in the higher wavenumber values.


Development of the van der Waals Force Parameters


Development of van der Waals (vdW) parameters is the greatest obstacle in parameter development. This is the reasoning behind using alternate force fields for comparison as well as the genetic algorithm.



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.

Current work is being done to investigate the viability of MOE versus Macro-Model for their ability to incorporate new residues. Then, the protein building functions can be used to build expanded polymers of various lengths and composition.

Initial Parameters

Bond Lengths



Si - O



Si - C



Si - H



Si - O (1)



Bond Angles



C - Si - C



O - Si - O



C - Si - O



C - Si - H



O - Si - H



Si - O - Si




Acknowledgements and References

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.

  6. Chemical Computing Group