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.crystal.Crystal;
41 import ffx.openmm.Force;
42 import ffx.openmm.IntArray;
43 import ffx.openmm.amoeba.VdwForce;
44 import ffx.potential.bonded.Atom;
45 import ffx.potential.extended.ExtendedSystem;
46 import ffx.potential.nonbonded.NonbondedCutoff;
47 import ffx.potential.nonbonded.VanDerWaals;
48 import ffx.potential.nonbonded.VanDerWaalsForm;
49 import ffx.potential.parameters.ForceField;
50 import ffx.potential.parameters.VDWPairType;
51 import ffx.potential.parameters.VDWType;
52
53 import java.util.HashMap;
54 import java.util.Map;
55 import java.util.logging.Level;
56 import java.util.logging.Logger;
57
58 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_AlchemicalMethod.OpenMM_AmoebaVdwForce_Annihilate;
59 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_AlchemicalMethod.OpenMM_AmoebaVdwForce_Decouple;
60 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_CutoffPeriodic;
61 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_NoCutoff;
62 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
63 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
64 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_False;
65 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
66 import static java.lang.String.format;
67
68
69
70
71 public class AmoebaVdwForce extends VdwForce {
72
73 private static final Logger logger = Logger.getLogger(AmoebaVdwForce.class.getName());
74
75
76
77
78
79
80
81
82 private int vdWClassForNoInteraction = 0;
83
84
85
86
87 private final Map<Integer, Integer> vdwClassToOpenMMType = new HashMap<>();
88
89
90
91
92
93
94 public AmoebaVdwForce(OpenMMEnergy openMMEnergy) {
95 VanDerWaals vdW = openMMEnergy.getVdwNode();
96 if (vdW == null) {
97 destroy();
98 return;
99 }
100
101 VanDerWaalsForm vdwForm = vdW.getVDWForm();
102 NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
103 Crystal crystal = openMMEnergy.getCrystal();
104
105 double radScale = 1.0;
106 if (vdwForm.radiusSize == VDWType.RADIUS_SIZE.DIAMETER) {
107 radScale = 0.5;
108 }
109
110 ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
111
112 Map<String, VDWType> vdwTypes = forceField.getVDWTypes();
113
114 for (VDWType vdwType : vdwTypes.values()) {
115 int atomClass = vdwType.atomClass;
116 if (!vdwClassToOpenMMType.containsKey(atomClass)) {
117 double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
118 double rad = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
119
120 if (rad == 0) {
121 rad = OpenMM_NmPerAngstrom * radScale;
122 }
123 int type = addParticleType(rad, eps);
124 vdwClassToOpenMMType.put(atomClass, type);
125 if (atomClass <= vdWClassForNoInteraction) {
126 vdWClassForNoInteraction = atomClass - 1;
127 }
128 }
129 }
130
131
132 int type = addParticleType(OpenMM_NmPerAngstrom, 0.0);
133 vdwClassToOpenMMType.put(vdWClassForNoInteraction, type);
134
135 Map<String, VDWPairType> vdwPairTypeMap = forceField.getVDWPairTypes();
136 for (VDWPairType vdwPairType : vdwPairTypeMap.values()) {
137 int c1 = vdwPairType.atomClasses[0];
138 int c2 = vdwPairType.atomClasses[1];
139 int type1 = vdwClassToOpenMMType.get(c1);
140 int type2 = vdwClassToOpenMMType.get(c2);
141 double rMin = vdwPairType.radius * OpenMM_NmPerAngstrom;
142 double eps = vdwPairType.wellDepth * OpenMM_KJPerKcal;
143 addTypePair(type1, type2, rMin, eps);
144 addTypePair(type2, type1, rMin, eps);
145 }
146
147 ExtendedSystem extendedSystem = vdW.getExtendedSystem();
148 double[] vdwPrefactorAndDerivs = new double[3];
149
150 int[] ired = vdW.getReductionIndex();
151 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
152 int nAtoms = atoms.length;
153 for (int i = 0; i < nAtoms; i++) {
154 Atom atom = atoms[i];
155 VDWType vdwType = atom.getVDWType();
156 int atomClass = vdwType.atomClass;
157 type = vdwClassToOpenMMType.get(atomClass);
158 int isAlchemical = atom.applyLambda() ? 1 : 0;
159 double scaleFactor = 1.0;
160 if (extendedSystem != null) {
161 extendedSystem.getVdwPrefactor(i, vdwPrefactorAndDerivs);
162 scaleFactor = vdwPrefactorAndDerivs[0];
163 }
164 addParticle_1(ired[i], type, vdwType.reductionFactor, isAlchemical, scaleFactor);
165 }
166
167 setCutoffDistance(nonbondedCutoff.off * OpenMM_NmPerAngstrom);
168 if (vdW.getDoLongRangeCorrection()) {
169 setUseDispersionCorrection(OpenMM_True);
170 } else {
171 setUseDispersionCorrection(OpenMM_False);
172 }
173
174 if (crystal.aperiodic()) {
175 setNonbondedMethod(OpenMM_AmoebaVdwForce_NoCutoff);
176 } else {
177 setNonbondedMethod(OpenMM_AmoebaVdwForce_CutoffPeriodic);
178 }
179
180 if (vdW.getLambdaTerm()) {
181 boolean annihilate = vdW.getIntramolecularSoftcore();
182 if (annihilate) {
183 setAlchemicalMethod(OpenMM_AmoebaVdwForce_Annihilate);
184 } else {
185 setAlchemicalMethod(OpenMM_AmoebaVdwForce_Decouple);
186 }
187 setSoftcoreAlpha(vdW.getAlpha());
188 setSoftcorePower((int) vdW.getBeta());
189 }
190
191 int[][] bondMask = vdW.getMask12();
192 int[][] angleMask = vdW.getMask13();
193
194
195 IntArray exclusions = new IntArray(0);
196 for (int i = 0; i < nAtoms; i++) {
197 exclusions.append(i);
198 final int[] bondMaski = bondMask[i];
199 for (int value : bondMaski) {
200 exclusions.append(value);
201 }
202 final int[] angleMaski = angleMask[i];
203 for (int value : angleMaski) {
204 exclusions.append(value);
205 }
206 setParticleExclusions(i, exclusions);
207 exclusions.resize(0);
208 }
209 exclusions.destroy();
210
211 int forceGroup = forceField.getInteger("VDW_FORCE_GROUP", 1);
212 setForceGroup(forceGroup);
213
214 logger.log(Level.INFO, vdW.toString());
215 logger.log(Level.FINE, format(" Force group:\t\t%d\n", forceGroup));
216 }
217
218
219
220
221
222
223
224 public static Force constructForce(OpenMMEnergy openMMEnergy) {
225 VanDerWaals vdW = openMMEnergy.getVdwNode();
226 if (vdW == null) {
227 return null;
228 }
229 return new AmoebaVdwForce(openMMEnergy);
230 }
231
232
233
234
235
236
237
238 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
239 VanDerWaals vdW = openMMEnergy.getVdwNode();
240 VanDerWaalsForm vdwForm = vdW.getVDWForm();
241 double radScale = 1.0;
242 if (vdwForm.radiusSize == VDWType.RADIUS_SIZE.DIAMETER) {
243 radScale = 0.5;
244 }
245
246 ExtendedSystem extendedSystem = vdW.getExtendedSystem();
247 double[] vdwPrefactorAndDerivs = new double[3];
248
249 int[] ired = vdW.getReductionIndex();
250 for (Atom atom : atoms) {
251 int index = atom.getXyzIndex() - 1;
252 VDWType vdwType = atom.getVDWType();
253
254
255 int type = vdwClassToOpenMMType.get(vdwType.atomClass);
256 if (!atom.getUse()) {
257
258 type = vdwClassToOpenMMType.get(vdWClassForNoInteraction);
259 }
260 int isAlchemical = atom.applyLambda() ? 1 : 0;
261 double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
262 double rad = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
263
264 double scaleFactor = 1.0;
265 if (extendedSystem != null) {
266 extendedSystem.getVdwPrefactor(index, vdwPrefactorAndDerivs);
267 scaleFactor = vdwPrefactorAndDerivs[0];
268 }
269
270 setParticleParameters(index, ired[index], rad, eps, vdwType.reductionFactor, isAlchemical, type, scaleFactor);
271 }
272 updateParametersInContext(openMMEnergy.getContext());
273 }
274
275 }