1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.openmm;
39
40 import ffx.openmm.DoubleArray;
41 import ffx.openmm.Force;
42 import ffx.openmm.CustomGBForce;
43 import ffx.potential.MolecularAssembly;
44 import ffx.potential.bonded.Atom;
45 import ffx.potential.nonbonded.GeneralizedKirkwood;
46 import ffx.potential.parameters.MultipoleType;
47
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
52 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
53 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_ComputationType.OpenMM_CustomGBForce_ParticlePair;
54 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_ComputationType.OpenMM_CustomGBForce_ParticlePairNoExclusions;
55 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_ComputationType.OpenMM_CustomGBForce_SingleParticle;
56 import static java.lang.String.format;
57
58
59
60
61 public class FixedChargeGBForce extends CustomGBForce {
62
63 private static final Logger logger = Logger.getLogger(FixedChargeGBForce.class.getName());
64
65
66
67
68
69
70 public FixedChargeGBForce(OpenMMEnergy openMMEnergy) {
71 GeneralizedKirkwood gk = openMMEnergy.getGK();
72 if (gk == null) {
73 destroy();
74 return;
75 }
76
77 double sTens = 0.0;
78 if (gk.getNonPolarModel() == GeneralizedKirkwood.NonPolarModel.BORN_SOLV
79 || gk.getNonPolarModel() == GeneralizedKirkwood.NonPolarModel.BORN_CAV_DISP) {
80 sTens = gk.getSurfaceTension();
81 sTens *= OpenMM_KJPerKcal;
82 sTens *= 100.0;
83
84
85 }
86
87 addPerParticleParameter("q");
88 addPerParticleParameter("radius");
89 addPerParticleParameter("scale");
90 addPerParticleParameter("surfaceTension");
91 addGlobalParameter("solventDielectric", gk.getSolventPermittivity());
92 addGlobalParameter("soluteDielectric", 1.0);
93 addGlobalParameter("dOffset", gk.getDielecOffset() * OpenMM_NmPerAngstrom);
94 addGlobalParameter("probeRadius", gk.getProbeRadius() * OpenMM_NmPerAngstrom);
95
96 addComputedValue("I", """
97 0.5*((1/L^3-1/U^3)/3.0+(1/U^4-1/L^4)/8.0*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);
98 U=r+sr2;
99 C=2/3*(1/or1^3-1/L^3)*step(sr2-r-or1);
100 L = step(sr2 - r1r)*sr2mr + (1 - step(sr2 - r1r))*L;
101 sr2mr = sr2 - r;
102 r1r = radius1 + r;
103 L = step(r1sr2 - r)*radius1 + (1 - step(r1sr2 - r))*L;
104 r1sr2 = radius1 + sr2;
105 L = r - sr2;
106 sr2 = scale2 * radius2;
107 or1 = radius1;
108 or2 = radius2;
109 """, OpenMM_CustomGBForce_ParticlePairNoExclusions);
110
111 addComputedValue("B", """
112 step(BB-radius)*BB + (1 - step(BB-radius))*radius;
113 BB = 1 / ( (3.0*III)^(1.0/3.0) );
114 III = step(II)*II + (1 - step(II))*1.0e-9/3.0;
115 II = maxI - I;
116 maxI = 1/(3.0*radius^3);
117 """, OpenMM_CustomGBForce_SingleParticle);
118
119 addEnergyTerm(
120 "surfaceTension*(radius+probeRadius+dOffset)^2*((radius+dOffset)/B)^6/6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B",
121 OpenMM_CustomGBForce_SingleParticle);
122
123
124 addEnergyTerm("""
125 -138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;
126 f=sqrt(r^2+B1*B2*exp(-r^2/(2.455*B1*B2)));
127 """, OpenMM_CustomGBForce_ParticlePair);
128
129 double[] baseRadii = gk.getBaseRadii();
130 double[] overlapScale = gk.getOverlapScale();
131 DoubleArray doubleArray = new DoubleArray(0);
132 MolecularAssembly molecularAssembly = openMMEnergy.getMolecularAssembly();
133 Atom[] atoms = molecularAssembly.getAtomArray();
134 int nAtoms = atoms.length;
135 for (int i = 0; i < nAtoms; i++) {
136 MultipoleType multipoleType = atoms[i].getMultipoleType();
137 doubleArray.append(multipoleType.charge);
138 doubleArray.append(OpenMM_NmPerAngstrom * baseRadii[i]);
139 doubleArray.append(overlapScale[i]);
140 doubleArray.append(sTens);
141 addParticle(doubleArray);
142 doubleArray.resize(0);
143 }
144 doubleArray.destroy();
145
146 double cut = gk.getCutoff();
147 setCutoffDistance(OpenMM_NmPerAngstrom * cut);
148
149 int forceGroup = molecularAssembly.getForceField().getInteger("GK_FORCE_GROUP", 1);
150 setForceGroup(forceGroup);
151 logger.log(Level.INFO, format(" Custom generalized Born force \t%d", forceGroup));
152 }
153
154
155
156
157
158
159
160 public static Force constructForce(OpenMMEnergy openMMEnergy) {
161 GeneralizedKirkwood gk = openMMEnergy.getGK();
162 if (gk == null) {
163 return null;
164 }
165 return new FixedChargeGBForce(openMMEnergy);
166 }
167
168
169
170
171
172
173
174 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
175 GeneralizedKirkwood gk = openMMEnergy.getGK();
176 double[] baseRadii = gk.getBaseRadii();
177 double[] overlapScale = gk.getOverlapScale();
178 boolean nea = gk.getNativeEnvironmentApproximation();
179
180 double sTens = 0.0;
181 if (gk.getNonPolarModel() == GeneralizedKirkwood.NonPolarModel.BORN_SOLV
182 || gk.getNonPolarModel() == GeneralizedKirkwood.NonPolarModel.BORN_CAV_DISP) {
183 sTens = gk.getSurfaceTension();
184 sTens *= OpenMM_KJPerKcal;
185 sTens *= 100.0;
186 }
187
188 DoubleArray parameters = new DoubleArray(0);
189 double lambdaElec = openMMEnergy.getSystem().getLambdaElec();
190 for (Atom atom : atoms) {
191 int index = atom.getXyzIndex() - 1;
192 double chargeUseFactor = 1.0;
193 if (!atom.getUse() || !atom.getElectrostatics()) {
194 chargeUseFactor = 0.0;
195 }
196 double lambdaScale = lambdaElec;
197 if (!atom.applyLambda()) {
198 lambdaScale = 1.0;
199 }
200
201 chargeUseFactor *= lambdaScale;
202 MultipoleType multipoleType = atom.getMultipoleType();
203 double charge = multipoleType.charge * chargeUseFactor;
204 double surfaceTension = sTens * chargeUseFactor;
205
206 double overlapScaleUseFactor = nea ? 1.0 : chargeUseFactor;
207 double oScale = overlapScale[index] * overlapScaleUseFactor;
208 double baseRadius = baseRadii[index];
209
210 parameters.append(charge);
211 parameters.append(OpenMM_NmPerAngstrom * baseRadius);
212 parameters.append(oScale);
213 parameters.append(surfaceTension);
214 setParticleParameters(index, parameters);
215 parameters.resize(0);
216 }
217 parameters.destroy();
218 updateParametersInContext(openMMEnergy.getContext());
219 }
220
221 }