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