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 com.sun.jna.ptr.DoubleByReference;
41 import com.sun.jna.ptr.IntByReference;
42 import ffx.crystal.Crystal;
43 import ffx.openmm.CustomBondForce;
44 import ffx.openmm.CustomNonbondedForce;
45 import ffx.openmm.DoubleArray;
46 import ffx.openmm.IntSet;
47 import ffx.potential.bonded.Atom;
48 import ffx.potential.nonbonded.NonbondedCutoff;
49 import ffx.potential.nonbonded.VanDerWaals;
50 import ffx.potential.nonbonded.VanDerWaalsForm;
51 import ffx.potential.parameters.ForceField;
52
53 import java.util.logging.Level;
54 import java.util.logging.Logger;
55
56 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
57 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
58 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_CutoffPeriodic;
59 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_NoCutoff;
60 import static ffx.potential.nonbonded.VanDerWaalsForm.EPSILON_RULE.GEOMETRIC;
61 import static ffx.potential.nonbonded.VanDerWaalsForm.RADIUS_RULE.ARITHMETIC;
62 import static ffx.potential.nonbonded.VanDerWaalsForm.VDW_TYPE.LENNARD_JONES;
63 import static java.lang.String.format;
64
65
66
67
68
69
70
71
72
73
74 public class FixedChargeAlchemicalForces {
75
76 private static final Logger logger = Logger.getLogger(RestrainPositionsForce.class.getName());
77
78 private CustomNonbondedForce fixedChargeSoftcoreForce;
79
80 private CustomBondForce alchemicalAlchemicalStericsForce;
81
82 private CustomBondForce nonAlchemicalAlchemicalStericsForce;
83
84 public FixedChargeAlchemicalForces(OpenMMEnergy openMMEnergy,
85 FixedChargeNonbondedForce fixedChargeNonBondedForce) {
86
87 VanDerWaals vdW = openMMEnergy.getVdwNode();
88 if (vdW == null) {
89 return;
90 }
91
92
93
94
95
96 VanDerWaalsForm vdwForm = vdW.getVDWForm();
97 if (vdwForm.vdwType != LENNARD_JONES || vdwForm.radiusRule != ARITHMETIC
98 || vdwForm.epsilonRule != GEOMETRIC) {
99 logger.info(format(" VDW Type: %s", vdwForm.vdwType));
100 logger.info(format(" VDW Radius Rule: %s", vdwForm.radiusRule));
101 logger.info(format(" VDW Epsilon Rule: %s", vdwForm.epsilonRule));
102 logger.log(Level.SEVERE, " Unsupported van der Waals functional form.");
103 return;
104 }
105
106
107 String stericsMixingRules = " epsilon = sqrt(epsilon1*epsilon2);";
108 stericsMixingRules += " rmin = 0.5 * (sigma1 + sigma2) * 1.122462048309372981;";
109
110
111 String stericsEnergyExpression = "(vdw_lambda^beta)*epsilon*x*(x-2.0);";
112
113 stericsEnergyExpression += " x = 1.0 / (alpha*(1.0-vdw_lambda)^2.0 + (r/rmin)^6.0);";
114
115 String energyExpression = stericsEnergyExpression + stericsMixingRules;
116
117 fixedChargeSoftcoreForce = new CustomNonbondedForce(energyExpression);
118
119
120 OpenMMSystem openMMSystem = openMMEnergy.getSystem();
121
122 double alpha = openMMSystem.getVdWSoftcoreAlpha();
123 double beta = openMMSystem.getVdwSoftcorePower();
124
125 fixedChargeSoftcoreForce.addGlobalParameter("vdw_lambda", 1.0);
126 fixedChargeSoftcoreForce.addGlobalParameter("alpha", alpha);
127 fixedChargeSoftcoreForce.addGlobalParameter("beta", beta);
128 fixedChargeSoftcoreForce.addPerParticleParameter("sigma");
129 fixedChargeSoftcoreForce.addPerParticleParameter("epsilon");
130
131
132 IntSet alchemicalGroup = new IntSet();
133 IntSet nonAlchemicalGroup = new IntSet();
134 DoubleByReference charge = new DoubleByReference();
135 DoubleByReference sigma = new DoubleByReference();
136 DoubleByReference eps = new DoubleByReference();
137
138 int index = 0;
139 DoubleArray parameters = new DoubleArray(0);
140 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
141 for (Atom atom : atoms) {
142 if (atom.applyLambda()) {
143 alchemicalGroup.insert(index);
144 } else {
145 nonAlchemicalGroup.insert(index);
146 }
147
148 fixedChargeNonBondedForce.getParticleParameters(index, charge, sigma, eps);
149 double sigmaValue = sigma.getValue();
150 double epsValue = eps.getValue();
151
152
153 if (sigmaValue == 0.0) {
154 sigmaValue = 1.0;
155 epsValue = 0.0;
156 }
157
158 parameters.append(sigmaValue);
159 parameters.append(epsValue);
160 fixedChargeSoftcoreForce.addParticle(parameters);
161 parameters.resize(0);
162 index++;
163 }
164 parameters.destroy();
165 fixedChargeSoftcoreForce.addInteractionGroup(alchemicalGroup, alchemicalGroup);
166 fixedChargeSoftcoreForce.addInteractionGroup(alchemicalGroup, nonAlchemicalGroup);
167 alchemicalGroup.destroy();
168 nonAlchemicalGroup.destroy();
169
170 Crystal crystal = openMMEnergy.getCrystal();
171 if (crystal.aperiodic()) {
172 fixedChargeSoftcoreForce.setNonbondedMethod(OpenMM_CustomNonbondedForce_NoCutoff);
173 } else {
174 fixedChargeSoftcoreForce.setNonbondedMethod(OpenMM_CustomNonbondedForce_CutoffPeriodic);
175 }
176
177 NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
178 double off = nonbondedCutoff.off;
179 double cut = nonbondedCutoff.cut;
180 if (cut == off) {
181 logger.warning(" OpenMM does not properly handle cutoffs where cut == off!");
182 if (cut == Double.MAX_VALUE || cut == Double.POSITIVE_INFINITY) {
183 logger.info(" Detected infinite or max-value cutoff; setting cut to 1E+40 for OpenMM.");
184 cut = 1E40;
185 } else {
186 logger.info(format(" Detected cut %8.4g == off %8.4g; scaling cut to 0.99 of off for OpenMM.", cut, off));
187 cut *= 0.99;
188 }
189 }
190 fixedChargeSoftcoreForce.setCutoffDistance(OpenMM_NmPerAngstrom * off);
191 fixedChargeSoftcoreForce.setUseSwitchingFunction(OpenMM_True);
192 fixedChargeSoftcoreForce.setSwitchingDistance(OpenMM_NmPerAngstrom * cut);
193
194 ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
195 int forceGroup = forceField.getInteger("VDW_FORCE_GROUP", 1);
196 fixedChargeSoftcoreForce.setForceGroup(forceGroup);
197
198
199 alchemicalAlchemicalStericsForce = new CustomBondForce(stericsEnergyExpression);
200
201
202 nonAlchemicalAlchemicalStericsForce = new CustomBondForce(stericsEnergyExpression);
203
204
205 alchemicalAlchemicalStericsForce.addPerBondParameter("rmin");
206 alchemicalAlchemicalStericsForce.addPerBondParameter("epsilon");
207 alchemicalAlchemicalStericsForce.addGlobalParameter("vdw_lambda", 1.0);
208 alchemicalAlchemicalStericsForce.addGlobalParameter("alpha", alpha);
209 alchemicalAlchemicalStericsForce.addGlobalParameter("beta", beta);
210
211 nonAlchemicalAlchemicalStericsForce.addPerBondParameter("rmin");
212 nonAlchemicalAlchemicalStericsForce.addPerBondParameter("epsilon");
213 nonAlchemicalAlchemicalStericsForce.addGlobalParameter("vdw_lambda", 1.0);
214 nonAlchemicalAlchemicalStericsForce.addGlobalParameter("alpha", alpha);
215 nonAlchemicalAlchemicalStericsForce.addGlobalParameter("beta", beta);
216
217 int range = fixedChargeNonBondedForce.getNumExceptions();
218
219 IntByReference atomi = new IntByReference();
220 IntByReference atomj = new IntByReference();
221 int[][] torsionMask = vdW.getMask14();
222
223 for (int i = 0; i < range; i++) {
224 fixedChargeNonBondedForce.getExceptionParameters(i, atomi, atomj, charge, sigma, eps);
225
226
227
228 fixedChargeSoftcoreForce.addExclusion(atomi.getValue(), atomj.getValue());
229
230
231 int[] maskI = torsionMask[atomi.getValue()];
232 int jID = atomj.getValue();
233 boolean epsException = false;
234
235 for (int mask : maskI) {
236 if (mask == jID) {
237 epsException = true;
238 break;
239 }
240 }
241
242 if (epsException) {
243 Atom atom1 = atoms[atomi.getValue()];
244 Atom atom2 = atoms[atomj.getValue()];
245 boolean bothAlchemical = false;
246 boolean oneAlchemical = false;
247 if (atom1.applyLambda() && atom2.applyLambda()) {
248 bothAlchemical = true;
249 } else if ((atom1.applyLambda() && !atom2.applyLambda()) || (!atom1.applyLambda() && atom2.applyLambda())) {
250 oneAlchemical = true;
251 }
252 if (bothAlchemical || oneAlchemical) {
253 DoubleArray bondParameters = new DoubleArray(0);
254 bondParameters.append(sigma.getValue() * 1.122462048309372981);
255 bondParameters.append(eps.getValue());
256 if (bothAlchemical) {
257 alchemicalAlchemicalStericsForce.addBond(atomi.getValue(), atomj.getValue(), bondParameters);
258 } else {
259 nonAlchemicalAlchemicalStericsForce.addBond(atomi.getValue(), atomj.getValue(), bondParameters);
260 }
261 bondParameters.destroy();
262 }
263 }
264 }
265 alchemicalAlchemicalStericsForce.setForceGroup(forceGroup);
266 nonAlchemicalAlchemicalStericsForce.setForceGroup(forceGroup);
267 logger.log(Level.INFO, format(" Added fixed charge softcore force \t%d", forceGroup));
268 logger.log(Level.INFO, format(" Alpha = %8.6f and beta = %8.6f", alpha, beta));
269 }
270
271 public CustomNonbondedForce getFixedChargeSoftcoreForce() {
272 return fixedChargeSoftcoreForce;
273 }
274
275 public CustomBondForce getAlchemicalAlchemicalStericsForce() {
276 return alchemicalAlchemicalStericsForce;
277 }
278
279 public CustomBondForce getNonAlchemicalAlchemicalStericsForce() {
280 return nonAlchemicalAlchemicalStericsForce;
281 }
282
283 }