Contents
Resources
OpenMM Documentation (of these links, the User Guide is the most important)
OpenMM Script Builder (Useful for Running Simple Simulations)
AMOEBA Model Overview (useful for understanding multipole conventions used by OpenMM)
Overview
This tutorial will teach you to, first, implement custom force fields in OpenMM and, second, to run simple simulations using these force fields. Here we will use a model for CO2 as an example; however, the content of this tutorial should be applicable to many different force fields and systems.
For this tutorial, it's expected that you've read through the OpenMM User's Guide and are familiar with its contents. Furthermore, you should have already installed and tested OpenMM. Finally, for the examples in this tutorial you'll need the following files:
Untar the file and make sure you can run the first example by executing the following commands:
$ tar -xvf openmm_tutorial.tgz $ cd 01_flexible_isotropic_model $ python run_co2.py
The sample output you get should be similar (but likely not identical) to the sample output code located in each directory.
Implementing New Potentials
Electrostatics
He we will be using the files found in the first subdirectory, 01_flexible_isotropic_model. Looking inside the file sapt_co2.xml (which contains our force field), the following block of xml code represents the model for electrostatics:
code xml
<AmoebaMultipoleForce direct11Scale="1.0" direct12Scale="1.0" direct13Scale="1.0" direct14Scale="1.0"
mpole12Scale="0.5" mpole13Scale="0.5" mpole14Scale="0.8" mpole15Scale="0.8" mutual11Scale="0.0"
mutual12Scale="0.0" mutual13Scale="1.0" mutual14Scale="1.0" polar12Scale="1.0" polar13Scale="0.0"
polar14Intra="0.5" polar14Scale="1.0" polar15Scale="1.0" >
<Multipole type="1" kz="2" kx="0" c0="0.884358" d1="0" d2="0" d3="2.64589e-07" q11="6.41139e-05" q21="0" q22="6.41139e-05" q31="0" q32="0" q33="-0.000128228" />
<Multipole type="2" kz="1" kx="0" c0="-0.442179" d1="0" d2="0" d3="-0.00621566" q11="-3.68839e-05" q21="0" q22="-3.68839e-05" q31="0" q32="0" q33="7.37679e-05" />
<Multipole type="4" kz="0" kx="0" c0="0.0" d1="0" d2="0" d3="0" q11="0" q21="0" q22="0" q31="0" q32="0" q33="0" />
<Polarize type="1" polarizability="1.333661e-04" thole="0.33" pgrp1="2" />
<Polarize type="2" polarizability="1.227116e-03" thole="0.33" pgrp1="1" />
<Polarize type="4" polarizability="1e-9" thole="0.33" pgrp1="1" />
</AmoebaMultipoleForce>This OpenMM Force encodes AMOEBA-style electrostatics and polarization, which is very similar to the functional forms used by ORIENT. See Ponder2010 for a good reference on the functional forms and conventions used by AMOEBA.
Several points are worth mention. First, each multipole tag represents the permanent electrostatics for a particular atomtype. Each tag represents the following:
Type: Atomtype defined by the multipole tag. This should correspond to the atomtype name given in the <AtomTypes> section. Note that, while in general OpenMM allows the user to give arbitrary names to atomtypes, AmoebaMultipoleForce requires that all atom types be integers. In the example above, the first multipole tag is for carbon, atomtype 1.
kz Atomtype defining the z-axis of the local reference frame. In the example above, kz="2" implies that the local reference frame for carbon (atomtype 1) has a z axis pointing from the central carbon to atomtype 2 (oxygen). Note that, in this case, there are two matching atoms that could define the z-axis, and internally OpenMM might choose either one. It's up to to ensure that your electrostatic energies don't depend on which atom is chosen as the z-axis, but there should be no problem so long as your multipoles obey the local symmetry of the system.
kx Atomtype defining the xz-plane of the local reference frame. In the example above, neither atomtype 1 nor atomtype 2 should be sensitive to the choice of z-axis, so we have simply written kx="0". There isn't (and shouldn't be!) an atomtype 0; instead, this is a flag to tell OpenMM to use the global coordinate frame for the x-axis.
c0, d1, d2, etc. Cartesian multipole moments for the atomtype, expressed in the local reference frame.
Important
Depending on the axis convention used, you may need more flexibility in defining the z and x axes. See the following resources for details on how to do this:
Important
At the time of this writing, the scale factors represented in the AmoebaMultipoleForce tag are hard coded in, and can't actually be modified in the xml file. See the OpenMM forum for details.
To obtain parameters for the AmoebaMultipoleForce, it is necessary to convert the multipole moments output by ORIENT (which are usually given in atomic units, and in the global reference frame) into the input format for OpenMM (SI units, local reference frame). The following script will perform the necessary conversion, such that usually only the atomtype and kz/kx fields need manual editing:
generate_electrostatic_parametrs.tgz
This script is run as follows (again using co2 as an example):
code bash $ python convert_mom_to_xml.py co2_ISA-GRID_L2.mom co2.axes <Multipole type="C0" kz="fill" kx="fill" c0="0.884358" d1="0" d2="0" d3="2.64589e-07" q11="6.41139e-05" q21="0" q22="6.41139e-05" q31="0" q32="0" q33="-0.000128228" /> <Multipole type="O1" kz="fill" kx="fill" c0="-0.44221" d1="0" d2="0" d3="-0.00621566" q11="-3.68839e-05" q21="0" q22="-3.68839e-05" q31="0" q32="0" q33="7.37679e-05" /> <Multipole type="O2" kz="fill" kx="fill" c0="-0.442204" d1="0" d2="0" d3="-0.00621518" q11="-3.68877e-05" q21="0" q22="-3.68877e-05" q31="0" q32="0" q33="7.37753e-05" />
Important
My .axes file follows a slightly different format than ORIENT's. Also note that, for this system, the partial charges don't add up precisely to zero. Before simulating this system in OpenMM, you may need to manually edit the c0 fields such that you end up with a neutral system.
Polarization
Polarization is also described using the AmoebaMultipoleForce, and is given by the Polarize tags. Each polarize tag contains the following information:
type: same as above
polarizability: Isotropic induced dipole polarizability, given in units of nm^3
thole: Dimensionless Thole parameter, following the damping function used by AMOEBA.
pgrp1, pgrp2, etc.: Polarization groups for that atomtype. See Shi2013 for details. In brief, the permanent multipoles of atoms within the same polarization group do not polarize each other. In this example, C and O both belong in the same polarization group.
Warning
As of this writing, sometimes using a polarizability of exactly zero leads to divergences in the PME calculation. If you use a polarizability of exactly zero for an atom, and see infinite energies in the AmoebaMultipoleForce, using a very small polarizability (say 1e-9) will fix this problem.
Non-electrostatic Pairwise Energies
Isotropic Models
Aside from electrostatic and polarization energies (described above), custom pair potentials in OpenMM can be implemented using CustomNonbondedForce. The OpenMM documentation does a good job of describing how to use this force and generate your own xml files, so you should read over this documentation before proceeding.
Below I've copied part of the sapt_co2.xml file, which shows the isotropic pair potential using a functional form given in VanVleet2016:
code xml
<CustomNonbondedForce bondCutoff="4"
energy="A*K2*exBr - Adi*(f6*C6/(r^6) - f8*C8/(r^8) - f10*C10/(r^10) - f12*C12/(r^12));
A=Aex-Ael-Ain-Adh;
Aex=(Aexch1*Aexch2);
Ael=(Aelec1*Aelec2);
Ain=(Aind1*Aind2);
Adh=(Adhf1*Adhf2);
Adi=(Adisp1*Adisp2);
K2=(Br^2)/3 + Br + 1;
f12 = f10 - exX*( (1/39916800)*(X^11)*(1 + X/12) );
f10 = f8 - exX*( (1/362880)*(X^9)*(1 + X/10 ) );
f8 = f6 - exX*( (1/5040)*(X^7)*(1 + X/8 ) );
f6 = 1 - exX*(1 + X * (1 + (1/2)*X*(1 + (1/3)*X*(1 + (1/4)*X*(1 + (1/5)*X*(1 + (1/6)*X ) ) ) ) ) );
exX = exp(-X);
X = Br - r * ( 2*(B^2)*r + 3*B )/(Br^2 + 3*Br + 3) ;
exBr = exp(-Br);
Br = B*r;
B=sqrt(Bexp1*Bexp2);
C6=sqrt(C61*C62); C8=sqrt(C81*C82); C10=sqrt(C101*C102); C12=sqrt(C121*C122)"
>
<PerParticleParameter name="Aexch"/>
<PerParticleParameter name="Aelec"/>
<PerParticleParameter name="Aind"/>
<PerParticleParameter name="Adhf"/>
<PerParticleParameter name="Adisp"/>
<PerParticleParameter name="Bexp"/>
<PerParticleParameter name="C6"/>
<PerParticleParameter name="C8"/>
<PerParticleParameter name="C10"/>
<PerParticleParameter name="C12"/>
<Atom class="C_CO2" Aexch="7.073403e+01" Aelec="3.871424e+01" Aind="3.020657e+01" Adhf="4.204415e+00" Adisp="7.114610e-01" Bexp="3.820882e+01" C6="1.262157e-03" C8="6.962318e-05" C10="3.262527e-06" C12="2.019001e-07"/>
<Atom class="O_CO2" Aexch="2.215524e+02" Aelec="1.257688e+02" Aind="1.276994e+01" Adhf="5.232102e+01" Adisp="1.072543e+00" Bexp="4.621758e+01" C6="9.540046e-04" C8="4.693183e-05" C10="3.163703e-06" C12="1.458385e-07"/>
</CustomNonbondedForce>In particular, this section of the xml file describes the exponentially-decaying charge penetration terms and the damped dispersion expansion. The OpenMM documentation is quite good for CustomNonbondedForce, so consult this guide for further questions.
Anisotropic Models
Anisotropic pair potentials are also (albeit in a limited sense) available in OpenMM. The Schmidt group is actively working on a developing a new force type (CustomAnisotropicNonbondedForce) which will allow anisotropic force fields to be developed in a general sense. In the meantime, a workaround (which uses CustomHbondForce) is available, and can be used on small molecules. Note that CustomHbondForce has not been optimized very well for the CPU platform, and the GPU code for running the anisotropic potential (see subdirectory 02_flexible_anisotropic_potential in the example files given above) runs at least 100x faster than the CPU code. In short, contact the Schmidt group before implementing your own anisotropic potentials, as we soon hope to have better implementations for these models.
Using CustomHbondForce, the anisotropic pair potential for CO2 looks as follows:
code xml
<CustomHbondForce particlesPerDonor="2" particlesPerAcceptor="2" bondCutoff="4"
energy="scale*(A*K2*exBr - (Adi-Adi_lr)*(f6*C6/(r^6) + f8*C8/(r^8) + f10*C10/(r^10) + f12*C12/(r^12)));
A=Aex-Ael-Ain-Adh;
;
Aex=(Aexch1*Aexch2*Aexch1_sph*Aexch2_sph);
Aexch1_sph= 1 + aexch_y101*y101 + aexch_y201*y201 + aexch_y22c1*y22c1;
Aexch2_sph= 1 + aexch_y102*y102 + aexch_y202*y202 + aexch_y22c2*y22c2;
;
Ael=(Aelec1*Aelec2*Aelec1_sph*Aelec2_sph);
Aelec1_sph= 1 + aelec_y101*y101 + aelec_y201*y201 + aelec_y22c1*y22c1;
Aelec2_sph= 1 + aelec_y102*y102 + aelec_y202*y202 + aelec_y22c2*y22c2;
;
Ain=(Aind1*Aind2*Aind1_sph*Aind2_sph);
Aind1_sph= 1 + aind_y101*y101 + aind_y201*y201 + aind_y22c1*y22c1;
Aind2_sph= 1 + aind_y102*y102 + aind_y202*y202 + aind_y22c2*y22c2;
;
Adh=(Adhf1*Adhf2*Adhf1_sph*Adhf2_sph);
Adhf1_sph= 1 + adhf_y101*y101 + adhf_y201*y201 + adhf_y22c1*y22c1;
Adhf2_sph= 1 + adhf_y102*y102 + adhf_y202*y202 + adhf_y22c2*y22c2;
;
Adi_lr=Adisp1*Adisp2;
Adi=(Adisp1*Adisp2*Adisp1_sph*Adisp2_sph);
Adisp1_sph= 1 + adisp_y101*y101 + adisp_y201*y201 + adisp_y22c1*y22c1;
Adisp2_sph= 1 + adisp_y102*y102 + adisp_y202*y202 + adisp_y22c2*y22c2;
;
K2=(Br^2)/3 + Br + 1;
f12 = f10 - exX*( (1/39916800)*(X^11)*(1 + X/12) );
f10 = f8 - exX*( (1/362880)*(X^9)*(1 + X/10 ) );
f8 = f6 - exX*( (1/5040)*(X^7)*(1 + X/8 ) );
f6 = 1 - exX*(1 + X * (1 + (1/2)*X*(1 + (1/3)*X*(1 + (1/4)*X*(1 + (1/5)*X*(1 + (1/6)*X ) ) ) ) ) );
exX = exp(-X);
X = Br - r * ( 2*(B^2)*r + 3*B )/(Br^2 + 3*Br + 3) ;
exBr = exp(-Br);
Br = B*r;
y101 = cos(phi1);
y102 = cos(phi2);
y201 = 0.5*(3*cos(phi1)^2 - 1);
y202 = 0.5*(3*cos(phi2)^2 - 1);
y22c1 = sqrt(0.75)*sin(phi1)^2*cos(2*theta1);
y22c2 = sqrt(0.75)*sin(phi2)^2*cos(2*theta2);
phi1=angle(d2,d1,a1);
phi2=angle(a2,a1,d1);
theta1=0;
theta2=0;
r = distance(a1,d1);
B=sqrt(Bexp1*Bexp2);
C6=sqrt(C61*C62); C8=sqrt(C81*C82); C10=sqrt(C101*C102); C12=sqrt(C121*C122)">
Scale factor to avoid double counting:
<GlobalParameter name="scale" defaultValue="0.5"/>
<PerDonorParameter name="Aexch1"/>
<PerDonorParameter name="aexch_y101"/>
<PerDonorParameter name="aexch_y201"/>
<PerDonorParameter name="aexch_y22c1"/>
<PerDonorParameter name="Aelec1"/>
<PerDonorParameter name="aelec_y101"/>
<PerDonorParameter name="aelec_y201"/>
<PerDonorParameter name="aelec_y22c1"/>
<PerDonorParameter name="Aind1"/>
<PerDonorParameter name="aind_y101"/>
<PerDonorParameter name="aind_y201"/>
<PerDonorParameter name="aind_y22c1"/>
<PerDonorParameter name="Adhf1"/>
<PerDonorParameter name="adhf_y101"/>
<PerDonorParameter name="adhf_y201"/>
<PerDonorParameter name="adhf_y22c1"/>
<PerDonorParameter name="Adisp1"/>
<PerDonorParameter name="adisp_y101"/>
<PerDonorParameter name="adisp_y201"/>
<PerDonorParameter name="adisp_y22c1"/>
<PerDonorParameter name="Bexp1"/>
<PerDonorParameter name="C61"/>
<PerDonorParameter name="C81"/>
<PerDonorParameter name="C101"/>
<PerDonorParameter name="C121"/>
<PerAcceptorParameter name="Aexch2"/>
<PerAcceptorParameter name="aexch_y102"/>
<PerAcceptorParameter name="aexch_y202"/>
<PerAcceptorParameter name="aexch_y22c2"/>
<PerAcceptorParameter name="Aelec2"/>
<PerAcceptorParameter name="aelec_y102"/>
<PerAcceptorParameter name="aelec_y202"/>
<PerAcceptorParameter name="aelec_y22c2"/>
<PerAcceptorParameter name="Aind2"/>
<PerAcceptorParameter name="aind_y102"/>
<PerAcceptorParameter name="aind_y202"/>
<PerAcceptorParameter name="aind_y22c2"/>
<PerAcceptorParameter name="Adhf2"/>
<PerAcceptorParameter name="adhf_y102"/>
<PerAcceptorParameter name="adhf_y202"/>
<PerAcceptorParameter name="adhf_y22c2"/>
For both Donor's and acceptors, class1 should be the particle, class2 should be the z-axis atom, and class3 should be the x-axis atom.
<Donor class1="C_CO2" class2="O1_CO2"
Aexch1="1.414656e+02" aexch_y101="0.000000e+00" aexch_y201="4.032290e-01" aexch_y22c1="0.000000e+00"
Aelec1="1.421757e+02" aelec_y101="0.000000e+00" aelec_y201="9.823030e-01" aelec_y22c1="0.000000e+00"
Aind1="1.919441e+01" aind_y101="0.000000e+00" aind_y201="-1.000000e+00" aind_y22c1="0.000000e+00"
Adhf1="1.680096e+01" adhf_y101="0.000000e+00" adhf_y201="-1.000000e+00" adhf_y22c1="0.000000e+00"
Adisp1="1.017705e+00" adisp_y101="0.000000e+00" adisp_y201="-5.366200e-01" adisp_y22c1="0.000000e+00"
Bexp1="4.348231e+01" C61="1.262157e-03" C81="6.962318e-05" C101="3.262527e-06" C121="2.019001e-07"/>
<Acceptor class1="C_CO2" class2="O1_CO2"
Aexch2="1.414656e+02" aexch_y102="0.000000e+00" aexch_y202="4.032290e-01" aexch_y22c2="0.000000e+00"
Aelec2="1.421757e+02" aelec_y102="0.000000e+00" aelec_y202="9.823030e-01" aelec_y22c2="0.000000e+00"
Aind2="1.919441e+01" aind_y102="0.000000e+00" aind_y202="-1.000000e+00" aind_y22c2="0.000000e+00"
Adhf2="1.680096e+01" adhf_y102="0.000000e+00" adhf_y202="-1.000000e+00" adhf_y22c2="0.000000e+00"
Adisp2="1.017705e+00" adisp_y102="0.000000e+00" adisp_y202="-5.366200e-01" adisp_y22c2="0.000000e+00"
Bexp2="4.348231e+01" C62="1.262157e-03" C82="6.962318e-05" C102="3.262527e-06" C122="2.019001e-07"/>
<Donor class1="O1_CO2" class2="C_CO2"
Aexch1="1.868046e+02" aexch_y101="2.045300e-02" aexch_y201="-8.284400e-02" aexch_y22c1="0.000000e+00"
Aelec1="8.714577e+01" aelec_y101="-2.134580e-01" aelec_y201="-2.281160e-01" aelec_y22c1="0.000000e+00"
Aind1="3.459401e+01" aind_y101="1.000000e+00" aind_y201="4.099520e-01" aind_y22c1="0.000000e+00"
Adhf1="2.856014e+01" adhf_y101="-1.000000e+00" adhf_y201="-5.779180e-01" adhf_y22c1="0.000000e+00"
Adisp1="7.304340e-01" adisp_y101="-1.000000e+00" adisp_y201="-3.965300e-01" adisp_y22c1="0.000000e+00"
Bexp1="4.450987e+01" C61="9.540046e-04" C81="4.693183e-05" C101="3.163703e-06" C121="1.458385e-07"/>
<Acceptor class1="O1_CO2" class2="C_CO2"
Aexch2="1.868046e+02" aexch_y102="2.045300e-02" aexch_y202="-8.284400e-02" aexch_y22c2="0.000000e+00"
Aelec2="8.714577e+01" aelec_y102="-2.134580e-01" aelec_y202="-2.281160e-01" aelec_y22c2="0.000000e+00"
Aind2="3.459401e+01" aind_y102="1.000000e+00" aind_y202="4.099520e-01" aind_y22c2="0.000000e+00"
Adhf2="2.856014e+01" adhf_y102="-1.000000e+00" adhf_y202="-5.779180e-01" adhf_y22c2="0.000000e+00"
Adisp2="7.304340e-01" adisp_y102="-1.000000e+00" adisp_y202="-3.965300e-01" adisp_y22c2="0.000000e+00"
Bexp2="4.450987e+01" C62="9.540046e-04" C82="4.693183e-05" C102="3.163703e-06" C122="1.458385e-07"/>
<Donor class1="O2_CO2" class2="C_CO2"
Aexch1="1.868046e+02" aexch_y101="2.045300e-02" aexch_y201="-8.284400e-02" aexch_y22c1="0.000000e+00"
Aelec1="8.714577e+01" aelec_y101="-2.134580e-01" aelec_y201="-2.281160e-01" aelec_y22c1="0.000000e+00"
Aind1="3.459401e+01" aind_y101="1.000000e+00" aind_y201="4.099520e-01" aind_y22c1="0.000000e+00"
Adhf1="2.856014e+01" adhf_y101="-1.000000e+00" adhf_y201="-5.779180e-01" adhf_y22c1="0.000000e+00"
Adisp1="7.304340e-01" adisp_y101="-1.000000e+00" adisp_y201="-3.965300e-01" adisp_y22c1="0.000000e+00"
Bexp1="4.450987e+01" C61="9.540046e-04" C81="4.693183e-05" C101="3.163703e-06" C121="1.458385e-07"/>
<Acceptor class1="O2_CO2" class2="C_CO2"
Aexch2="1.868046e+02" aexch_y102="2.045300e-02" aexch_y202="-8.284400e-02" aexch_y22c2="0.000000e+00"
Aelec2="8.714577e+01" aelec_y102="-2.134580e-01" aelec_y202="-2.281160e-01" aelec_y22c2="0.000000e+00"
Aind2="3.459401e+01" aind_y102="1.000000e+00" aind_y202="4.099520e-01" aind_y22c2="0.000000e+00"
Adhf2="2.856014e+01" adhf_y102="-1.000000e+00" adhf_y202="-5.779180e-01" adhf_y22c2="0.000000e+00"
Adisp2="7.304340e-01" adisp_y102="-1.000000e+00" adisp_y202="-3.965300e-01" adisp_y22c2="0.000000e+00"
Bexp2="4.450987e+01" C62="9.540046e-04" C82="4.693183e-05" C102="3.163703e-06" C122="1.458385e-07"/>
</CustomHbondForce>Several points are of note (most of which showcase why CustomAnisotropicNonbonedForce is being developed in the first place). All interactions between Donor Particles (here C and each oxygen) are Acceptor Particles (again C and O1,O2) are counted, so a scale factor of 1/2 is necessary to avoid double counting. Second, CustomHbondForce allows angles and dihedrals between atoms as input variables, and the local coordinate system can be calculated in this manner. Otherwise, this is a standard use of CustomHbondForce, and the OpenMM documentation for this class does a good job of describing how to use it.
Long-range Corrections
As a final implementation note, it is generally important in simulations to use long-range corrections for both the electrostatic and dispersion terms. Available long-range corrections differ based on which force is being used, but I've put some good defaults in the run_co2.py script provided in the example directory above:
code python
# Set distance cutoffs, constraints, and other force-specific options for each
# force we might encounter
for force in system.getForces():
if isinstance(force, mm.CustomHbondForce):
force.setNonbondedMethod(mm.CustomHbondForce.CutoffPeriodic)
force.setCutoffDistance(two_body_cutoff)
elif isinstance(force, mm.CustomNonbondedForce):
force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)
force.setCutoffDistance(two_body_cutoff)
force.setUseLongRangeCorrection(True)
elif isinstance(force, mm.AmoebaMultipoleForce):
force.setNonbondedMethod(mm.AmoebaMultipoleForce.PME)
elif isinstance(force, mm.NonbondedForce):
force.setNonbondedMethod(mm.NonbondedForce.LJPME)
force.setCutoffDistance(two_body_cutoff)
elif isinstance(force, mm.CustomManyParticleForce):
force.setNonbondedMethod(mm.CustomManyParticleForce.CutoffPeriodic)
force.setCutoffDistance(three_body_cutoff)
else:
passHere you can see that we're using PME for electrostatics. For dispersion, this script uses a cutoff with a set cutoff distance, and additionally corrects for energies beyond the cutoff using a long range correction term (called via the command "force.setUseLongRangeCorrection(True)" and as described in the OpenMM documentation.
Finally, note that, for anisotropic force fields, CustomHbondForce has no option for a long-range correction. In the example directory 02_flexible_anisotropic_potential, I've subtracted off the isotropic dispersion energy from CustomHbondForce and added it to CustomNonbondedForce, which then allows me to add a long-range correction.
Combination Rules
In general, OpenMM assumes that you will use combination rules to describe cross-interactions in the pair potential. I have not tried using explicit cross-terms, but this is in principle possible. See this example xml file for an example where this is done.