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;
39
40 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_3D_DoubleArray_create;
41 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_3D_DoubleArray_destroy;
42 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_3D_DoubleArray_set;
43 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGKCavitationForce_NonbondedMethod.OpenMM_AmoebaGKCavitationForce_NoCutoff;
44 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGKCavitationForce_addParticle;
45 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGKCavitationForce_create;
46 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGKCavitationForce_setNonbondedMethod;
47 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGKCavitationForce_setParticleParameters;
48 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGKCavitationForce_updateParametersInContext;
49 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_addParticle_1;
50 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_create;
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setDielectricOffset;
52 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setIncludeCavityTerm;
53 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setParticleParameters_1;
54 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setProbeRadius;
55 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setSoluteDielectric;
56 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setSolventDielectric;
57 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setSurfaceAreaFactor;
58 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setTanhParameters;
59 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_setTanhRescaling;
60 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaGeneralizedKirkwoodForce_updateParametersInContext;
61 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent12;
62 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent13;
63 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent14;
64 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_Covalent15;
65 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_CovalentType.OpenMM_AmoebaMultipoleForce_PolarizationCovalent11;
66 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_Bisector;
67 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_NoAxisType;
68 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ThreeFold;
69 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZBisect;
70 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZOnly;
71 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_MultipoleAxisTypes.OpenMM_AmoebaMultipoleForce_ZThenX;
72 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_NonbondedMethod.OpenMM_AmoebaMultipoleForce_NoCutoff;
73 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_NonbondedMethod.OpenMM_AmoebaMultipoleForce_PME;
74 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Direct;
75 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Extrapolated;
76 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_PolarizationType.OpenMM_AmoebaMultipoleForce_Mutual;
77 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_addMultipole;
78 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_create;
79 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setAEwald;
80 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setCovalentMap;
81 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setCutoffDistance;
82 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setEwaldErrorTolerance;
83 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setExtrapolationCoefficients;
84 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setMultipoleParameters;
85 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setMutualInducedMaxIterations;
86 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setMutualInducedTargetEpsilon;
87 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setNonbondedMethod;
88 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setPmeGridDimensions;
89 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_setPolarizationType;
90 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaMultipoleForce_updateParametersInContext;
91 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaTorsionTorsionForce_addTorsionTorsion;
92 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaTorsionTorsionForce_create;
93 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaTorsionTorsionForce_setTorsionTorsionGrid;
94 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_AlchemicalMethod.OpenMM_AmoebaVdwForce_Decouple;
95 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_CutoffPeriodic;
96 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_NonbondedMethod.OpenMM_AmoebaVdwForce_NoCutoff;
97 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_addParticleType;
98 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_addParticle_1;
99 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_addTypePair;
100 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_create;
101 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setAlchemicalMethod;
102 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setCutoffDistance;
103 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setNonbondedMethod;
104 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setParticleExclusions;
105 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setParticleParameters;
106 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setSoftcoreAlpha;
107 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setSoftcorePower;
108 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_setUseDispersionCorrection;
109 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaVdwForce_updateParametersInContext;
110 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_addParticle;
111 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_create;
112 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setAwater;
113 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setDispoff;
114 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setEpsh;
115 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setEpso;
116 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setParticleParameters;
117 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setRminh;
118 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setRmino;
119 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setShctd;
120 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_setSlevy;
121 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AmoebaWcaDispersionForce_updateParametersInContext;
122 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AngstromsPerNm;
123 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
124 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KcalPerKJ;
125 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
126 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_RadiansPerDegree;
127 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_AndersenThermostat_create;
128 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_AndersenThermostat_setDefaultCollisionFrequency;
129 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_AndersenThermostat_setDefaultTemperature;
130 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_BondArray_append;
131 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_BondArray_create;
132 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_BondArray_destroy;
133 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_False;
134 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
135 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CMMotionRemover_create;
136 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_applyConstraints;
137 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_create_2;
138 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_destroy;
139 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_getState;
140 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_reinitialize;
141 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_setParameter;
142 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_setPeriodicBoxVectors;
143 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_setPositions;
144 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Context_setVelocities;
145 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomAngleForce_addAngle;
146 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomAngleForce_addPerAngleParameter;
147 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomAngleForce_create;
148 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomAngleForce_setAngleParameters;
149 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomAngleForce_updateParametersInContext;
150 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomBondForce_addBond;
151 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomBondForce_addGlobalParameter;
152 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomBondForce_addPerBondParameter;
153 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomBondForce_create;
154 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomBondForce_setBondParameters;
155 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomBondForce_updateParametersInContext;
156 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCentroidBondForce_addBond;
157 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCentroidBondForce_addGroup;
158 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCentroidBondForce_addPerBondParameter;
159 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCentroidBondForce_create;
160 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCentroidBondForce_setUsesPeriodicBoundaryConditions;
161 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCompoundBondForce_addBond;
162 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCompoundBondForce_addGlobalParameter;
163 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCompoundBondForce_addPerBondParameter;
164 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCompoundBondForce_create;
165 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCompoundBondForce_setBondParameters;
166 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomCompoundBondForce_updateParametersInContext;
167 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomExternalForce_addGlobalParameter;
168 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomExternalForce_addParticle;
169 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomExternalForce_addPerParticleParameter;
170 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomExternalForce_create;
171 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_ComputationType.OpenMM_CustomGBForce_ParticlePair;
172 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_ComputationType.OpenMM_CustomGBForce_ParticlePairNoExclusions;
173 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_ComputationType.OpenMM_CustomGBForce_SingleParticle;
174 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_addComputedValue;
175 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_addEnergyTerm;
176 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_addGlobalParameter;
177 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_addParticle;
178 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_addPerParticleParameter;
179 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_create;
180 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_setCutoffDistance;
181 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_setParticleParameters;
182 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomGBForce_updateParametersInContext;
183 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomIntegrator_addComputePerDof;
184 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomIntegrator_addConstrainPositions;
185 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomIntegrator_addConstrainVelocities;
186 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomIntegrator_addGlobalVariable;
187 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomIntegrator_addPerDofVariable;
188 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomIntegrator_addUpdateContextState;
189 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomIntegrator_create;
190 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_addExclusion;
191 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_addGlobalParameter;
192 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_addInteractionGroup;
193 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_addParticle;
194 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_addPerParticleParameter;
195 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_create;
196 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_setCutoffDistance;
197 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_setNonbondedMethod;
198 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_setSwitchingDistance;
199 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_setUseSwitchingFunction;
200 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_DoubleArray_append;
201 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_DoubleArray_create;
202 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_DoubleArray_destroy;
203 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_DoubleArray_resize;
204 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_DoubleArray_set;
205 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Force_setForceGroup;
206 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Force_setName;
207 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_HarmonicBondForce_addBond;
208 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_HarmonicBondForce_create;
209 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_HarmonicBondForce_setBondParameters;
210 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_HarmonicBondForce_updateParametersInContext;
211 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntArray_append;
212 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntArray_create;
213 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntArray_destroy;
214 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntArray_resize;
215 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntArray_set;
216 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntSet_create;
217 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntSet_destroy;
218 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_IntSet_insert;
219 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Integrator_destroy;
220 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Integrator_setConstraintTolerance;
221 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Integrator_step;
222 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_LangevinIntegrator_create;
223 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_LangevinIntegrator_setRandomNumberSeed;
224 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_LocalEnergyMinimizer_minimize;
225 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_MonteCarloBarostat_create;
226 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_MonteCarloBarostat_setRandomNumberSeed;
227 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_addParticle;
228 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_create;
229 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_createExceptionsFromBonds;
230 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_getExceptionParameters;
231 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_getNumExceptions;
232 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_getParticleParameters;
233 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setCutoffDistance;
234 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setExceptionParameters;
235 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setNonbondedMethod;
236 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setPMEParameters;
237 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setParticleParameters;
238 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setSwitchingDistance;
239 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setUseDispersionCorrection;
240 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_setUseSwitchingFunction;
241 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_updateParametersInContext;
242 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_PeriodicTorsionForce_addTorsion;
243 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_PeriodicTorsionForce_create;
244 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_PeriodicTorsionForce_setTorsionParameters;
245 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_PeriodicTorsionForce_updateParametersInContext;
246 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Platform_getNumPlatforms;
247 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Platform_getOpenMMVersion;
248 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Platform_getPlatformByName;
249 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Platform_getPluginLoadFailures;
250 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Platform_loadPluginsFromDirectory;
251 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Platform_setPropertyDefaultValue;
252 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Energy;
253 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Forces;
254 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Positions;
255 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_DataType.OpenMM_State_Velocities;
256 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_destroy;
257 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_getForces;
258 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_getKineticEnergy;
259 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_getPeriodicBoxVectors;
260 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_getPositions;
261 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_getPotentialEnergy;
262 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_State_getVelocities;
263 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_StringArray_destroy;
264 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_StringArray_get;
265 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_StringArray_getSize;
266 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_addConstraint;
267 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_addForce;
268 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_addParticle;
269 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_create;
270 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_destroy;
271 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_getNumConstraints;
272 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_setDefaultPeriodicBoxVectors;
273 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_System_setParticleMass;
274 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Vec3Array_append;
275 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Vec3Array_create;
276 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Vec3Array_destroy;
277 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Vec3Array_get;
278 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_VerletIntegrator_create;
279 import static ffx.potential.nonbonded.GeneralizedKirkwood.NonPolarModel.GAUSS_DISP;
280 import static ffx.potential.nonbonded.VanDerWaalsForm.EPSILON_RULE.GEOMETRIC;
281 import static ffx.potential.nonbonded.VanDerWaalsForm.RADIUS_RULE.ARITHMETIC;
282 import static ffx.potential.nonbonded.VanDerWaalsForm.RADIUS_SIZE.RADIUS;
283 import static ffx.potential.nonbonded.VanDerWaalsForm.RADIUS_TYPE.R_MIN;
284 import static ffx.potential.nonbonded.VanDerWaalsForm.VDW_TYPE.LENNARD_JONES;
285 import static ffx.utilities.Constants.KCAL_TO_GRAM_ANG2_PER_PS2;
286 import static ffx.utilities.Constants.KCAL_TO_KJ;
287 import static ffx.utilities.Constants.R;
288 import static ffx.utilities.Constants.kB;
289 import static java.lang.Double.isFinite;
290 import static java.lang.Double.isInfinite;
291 import static java.lang.Double.isNaN;
292 import static java.lang.Math.exp;
293 import static java.lang.String.format;
294 import static java.util.Arrays.copyOfRange;
295 import static org.apache.commons.math3.util.FastMath.PI;
296 import static org.apache.commons.math3.util.FastMath.abs;
297 import static org.apache.commons.math3.util.FastMath.cos;
298 import static org.apache.commons.math3.util.FastMath.pow;
299 import static org.apache.commons.math3.util.FastMath.sqrt;
300 import static org.apache.commons.math3.util.FastMath.toRadians;
301
302 import com.sun.jna.Memory;
303 import com.sun.jna.Pointer;
304 import com.sun.jna.ptr.DoubleByReference;
305 import com.sun.jna.ptr.IntByReference;
306 import com.sun.jna.ptr.PointerByReference;
307 import edu.rit.mp.CharacterBuf;
308 import edu.rit.pj.Comm;
309 import edu.uiowa.jopenmm.OpenMMLibrary;
310 import edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_CustomNonbondedForce_NonbondedMethod;
311 import edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_NonbondedForce_NonbondedMethod;
312 import edu.uiowa.jopenmm.OpenMMUtils;
313 import edu.uiowa.jopenmm.OpenMM_Vec3;
314 import ffx.crystal.Crystal;
315 import ffx.potential.bonded.Angle;
316 import ffx.potential.bonded.AngleTorsion;
317 import ffx.potential.bonded.Atom;
318 import ffx.potential.bonded.Bond;
319 import ffx.potential.bonded.ImproperTorsion;
320 import ffx.potential.bonded.OutOfPlaneBend;
321 import ffx.potential.bonded.PiOrbitalTorsion;
322 import ffx.potential.bonded.RestraintBond;
323 import ffx.potential.bonded.RestraintTorsion;
324 import ffx.potential.bonded.StretchBend;
325 import ffx.potential.bonded.StretchTorsion;
326 import ffx.potential.bonded.Torsion;
327 import ffx.potential.bonded.TorsionTorsion;
328 import ffx.potential.bonded.UreyBradley;
329 import ffx.potential.extended.ExtendedSystem;
330 import ffx.potential.nonbonded.RestrainPosition;
331 import ffx.potential.nonbonded.GeneralizedKirkwood;
332 import ffx.potential.nonbonded.GeneralizedKirkwood.NonPolarModel;
333 import ffx.potential.nonbonded.NonbondedCutoff;
334 import ffx.potential.nonbonded.ParticleMeshEwald;
335 import ffx.potential.nonbonded.pme.Polarization;
336 import ffx.potential.nonbonded.pme.SCFAlgorithm;
337 import ffx.potential.nonbonded.ReciprocalSpace;
338 import ffx.potential.nonbonded.RestrainGroups;
339 import ffx.potential.nonbonded.VanDerWaals;
340 import ffx.potential.nonbonded.VanDerWaalsForm;
341 import ffx.potential.nonbonded.implicit.ChandlerCavitation;
342 import ffx.potential.nonbonded.implicit.DispersionRegion;
343 import ffx.potential.nonbonded.implicit.GaussVol;
344 import ffx.potential.parameters.AngleType;
345 import ffx.potential.parameters.AngleType.AngleFunction;
346 import ffx.potential.parameters.AngleType.AngleMode;
347 import ffx.potential.parameters.BondType;
348 import ffx.potential.parameters.BondType.BondFunction;
349 import ffx.potential.parameters.ForceField;
350 import ffx.potential.parameters.ImproperTorsionType;
351 import ffx.potential.parameters.MultipoleType;
352 import ffx.potential.parameters.OutOfPlaneBendType;
353 import ffx.potential.parameters.PiOrbitalTorsionType;
354 import ffx.potential.parameters.PolarizeType;
355 import ffx.potential.parameters.TorsionTorsionType;
356 import ffx.potential.parameters.TorsionType;
357 import ffx.potential.parameters.UreyBradleyType;
358 import ffx.potential.parameters.VDWPairType;
359 import ffx.potential.parameters.VDWType;
360 import ffx.potential.utils.EnergyException;
361 import ffx.potential.utils.PotentialsFunctions;
362 import ffx.potential.utils.PotentialsUtils;
363 import ffx.utilities.Constants;
364 import java.io.File;
365 import java.io.IOException;
366 import java.time.LocalDateTime;
367 import java.time.format.DateTimeFormatter;
368 import java.util.ArrayList;
369 import java.util.Arrays;
370 import java.util.HashMap;
371 import java.util.LinkedHashMap;
372 import java.util.List;
373 import java.util.Map;
374 import java.util.logging.Level;
375 import java.util.logging.Logger;
376 import java.util.stream.Collectors;
377 import java.util.stream.IntStream;
378 import org.apache.commons.configuration2.CompositeConfiguration;
379 import org.apache.commons.io.FilenameUtils;
380
381
382
383
384
385
386
387 @SuppressWarnings("deprecation")
388 public class ForceFieldEnergyOpenMM extends ForceFieldEnergy {
389
390 private static final Logger logger = Logger.getLogger(ForceFieldEnergyOpenMM.class.getName());
391
392 public final int enforcePBC;
393
394 private final Atom[] atoms;
395
396 private final int nAtoms;
397
398 private final double finiteDifferenceStepSize;
399
400 private Context context;
401
402 private System system;
403
404
405
406
407
408 private double lambdaStart = 0.0;
409
410 private boolean twoSidedFiniteDifference = true;
411
412 private final boolean freeOpenMM;
413
414
415
416
417
418
419
420
421
422 protected ForceFieldEnergyOpenMM(MolecularAssembly molecularAssembly, Platform requestedPlatform, int nThreads) {
423 super(molecularAssembly, nThreads);
424
425 Crystal crystal = getCrystal();
426 int symOps = crystal.spaceGroup.getNumberOfSymOps();
427 if (symOps > 1) {
428 logger.info("");
429 logger.severe(" OpenMM does not support symmetry operators.");
430 }
431
432 logger.info("\n Initializing OpenMM");
433
434 ForceField forceField = molecularAssembly.getForceField();
435 context = new Context(forceField, requestedPlatform);
436
437 atoms = molecularAssembly.getAtomArray();
438 nAtoms = atoms.length;
439 system = new System(molecularAssembly);
440
441 boolean aperiodic = super.getCrystal().aperiodic();
442 boolean pbcEnforced = forceField.getBoolean("ENFORCE_PBC", !aperiodic);
443 enforcePBC = pbcEnforced ? OpenMM_True : OpenMM_False;
444
445 finiteDifferenceStepSize = forceField.getDouble("FD_DLAMBDA", 0.001);
446 twoSidedFiniteDifference = forceField.getBoolean("FD_TWO_SIDED", twoSidedFiniteDifference);
447
448 freeOpenMM = forceField.getBoolean("FREE_OPENMM", true);
449 }
450
451
452
453
454
455
456
457
458
459 private static int getDefaultDevice(CompositeConfiguration props) {
460 String availDeviceProp = props.getString("availableDevices", props.getString("CUDA_DEVICES"));
461 if (availDeviceProp == null) {
462 int nDevs = props.getInt("numCudaDevices", 1);
463 availDeviceProp = IntStream.range(0, nDevs).mapToObj(Integer::toString)
464 .collect(Collectors.joining(" "));
465 }
466 availDeviceProp = availDeviceProp.trim();
467
468 String[] availDevices = availDeviceProp.split("\\s+");
469 int nDevs = availDevices.length;
470 int[] devs = new int[nDevs];
471 for (int i = 0; i < nDevs; i++) {
472 devs[i] = Integer.parseInt(availDevices[i]);
473 }
474
475 logger.info(format(" Number of CUDA devices: %d.", nDevs));
476
477 int index = 0;
478 try {
479 Comm world = Comm.world();
480 if (world != null) {
481 int size = world.size();
482
483
484 int messageLen = 100;
485 String host = world.host();
486
487 host = host.substring(0, Math.min(messageLen, host.length()));
488
489 host = String.format("%-100s", host);
490 char[] messageOut = host.toCharArray();
491 CharacterBuf out = CharacterBuf.buffer(messageOut);
492
493
494 char[][] incoming = new char[size][messageLen];
495 CharacterBuf[] in = new CharacterBuf[size];
496 for (int i = 0; i < size; i++) {
497 in[i] = CharacterBuf.buffer(incoming[i]);
498 }
499
500 try {
501 world.allGather(out, in);
502 } catch (IOException ex) {
503 logger.severe(
504 String.format(" Failure at the allGather step for determining rank: %s\n%s", ex,
505 Utilities.stackTraceToString(ex)));
506 }
507 int ownIndex = -1;
508 int rank = world.rank();
509 boolean selfFound = false;
510
511 for (int i = 0; i < size; i++) {
512 String hostI = new String(incoming[i]);
513 if (hostI.equalsIgnoreCase(host)) {
514 ++ownIndex;
515 if (i == rank) {
516 selfFound = true;
517 break;
518 }
519 }
520 }
521 if (!selfFound) {
522 logger.severe(
523 String.format(" Rank %d: Could not find any incoming host messages matching self %s!",
524 rank, host.trim()));
525 } else {
526 index = ownIndex % nDevs;
527 }
528 }
529 } catch (IllegalStateException ise) {
530
531 }
532 return devs[index];
533 }
534
535
536
537
538
539
540
541 private static Pointer pointerForString(String string) {
542 Pointer pointer = new Memory(string.length() + 1);
543 pointer.setString(0, string);
544 return pointer;
545 }
546
547
548
549
550
551
552
553
554 private static String stringFromArray(PointerByReference stringArray, int i) {
555 Pointer platformPtr = OpenMM_StringArray_get(stringArray, i);
556 if (platformPtr == null) {
557 return null;
558 }
559 return platformPtr.getString(0);
560 }
561
562
563
564
565
566
567
568
569
570
571
572
573 public void createContext(String integratorString, double timeStep, double temperature,
574 boolean forceCreation) {
575 context.create(integratorString, timeStep, temperature, forceCreation);
576 }
577
578
579
580
581
582
583
584
585
586
587
588
589 public State createState(boolean positions, boolean energies, boolean forces, boolean velocities) {
590 return new State(positions, energies, forces, velocities);
591 }
592
593
594 @Override
595 public boolean destroy() {
596 boolean ffxFFEDestroy = super.destroy();
597 if (freeOpenMM) {
598 free();
599 logger.fine(" Destroyed the Context, Integrator, and OpenMMSystem.");
600 }
601 return ffxFFEDestroy;
602 }
603
604
605 @Override
606 public double energy(double[] x) {
607 return energy(x, false);
608 }
609
610
611 @Override
612 public double energy(double[] x, boolean verbose) {
613
614 if (lambdaBondedTerms) {
615 return 0.0;
616 }
617
618
619 context.getContextPointer();
620
621 updateParameters(atoms);
622
623
624 unscaleCoordinates(x);
625
626 setCoordinates(x);
627
628 State state = new State(false, true, false, false);
629 double e = state.potentialEnergy;
630 state.free();
631
632 if (!isFinite(e)) {
633 String message = String.format(" Energy from OpenMM was a non-finite %8g", e);
634 logger.warning(message);
635 if (lambdaTerm) {
636 system.printLambdaValues();
637 }
638 throw new EnergyException(message);
639 }
640
641 if (verbose) {
642 logger.log(Level.INFO, String.format("\n OpenMM Energy: %14.10g", e));
643 }
644
645
646 scaleCoordinates(x);
647
648 return e;
649 }
650
651
652 @Override
653 public double energyAndGradient(double[] x, double[] g) {
654 return energyAndGradient(x, g, false);
655 }
656
657
658 @Override
659 public double energyAndGradient(double[] x, double[] g, boolean verbose) {
660 if (lambdaBondedTerms) {
661 return 0.0;
662 }
663
664
665
666
667 unscaleCoordinates(x);
668
669
670 context.getContextPointer();
671
672 setCoordinates(x);
673
674 State state = new State(false, true, true, false);
675 double e = state.potentialEnergy;
676 g = state.getGradient(g);
677 state.free();
678
679 if (!isFinite(e)) {
680 String message = format(" Energy from OpenMM was a non-finite %8g", e);
681 logger.warning(message);
682 if (lambdaTerm) {
683 system.printLambdaValues();
684 }
685 throw new EnergyException(message);
686 }
687
688
689
690
691
692
693
694
695
696
697 if (maxDebugGradient < Double.MAX_VALUE) {
698 boolean extremeGrad = Arrays.stream(g)
699 .anyMatch((double gi) -> (gi > maxDebugGradient || gi < -maxDebugGradient));
700 if (extremeGrad) {
701 File origFile = molecularAssembly.getFile();
702 String timeString = LocalDateTime.now()
703 .format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
704
705 String filename = format("%s-LARGEGRAD-%s.pdb",
706 FilenameUtils.removeExtension(molecularAssembly.getFile().getName()), timeString);
707 PotentialsFunctions ef = new PotentialsUtils();
708 filename = ef.versionFile(filename);
709
710 logger.warning(
711 format(" Excessively large gradients detected; printing snapshot to file %s", filename));
712 ef.saveAsPDB(molecularAssembly, new File(filename));
713 molecularAssembly.setFile(origFile);
714 }
715 }
716
717 if (verbose) {
718 logger.log(Level.INFO, format("\n OpenMM Energy: %14.10g", e));
719 }
720
721
722 scaleCoordinatesAndGradient(x, g);
723
724 return e;
725 }
726
727
728
729
730
731
732
733
734 public double energyAndGradientFFX(double[] x, double[] g) {
735 return super.energyAndGradient(x, g, false);
736 }
737
738
739
740
741
742
743
744
745
746 public double energyAndGradientFFX(double[] x, double[] g, boolean verbose) {
747 return super.energyAndGradient(x, g, verbose);
748 }
749
750
751
752
753
754
755
756 public double energyFFX(double[] x) {
757 return super.energy(x, false);
758 }
759
760
761
762
763
764
765
766
767 public double energyFFX(double[] x, boolean verbose) {
768 return super.energy(x, verbose);
769 }
770
771
772
773
774
775
776 public Context getContext() {
777 return context;
778 }
779
780
781
782
783
784
785 @Override
786 public double[] getGradient(double[] g) {
787 State state = new State(false, false, true, false);
788 g = state.getGradient(g);
789 state.free();
790 return g;
791 }
792
793
794 @Override
795 public Platform getPlatform() {
796 return context.platform;
797 }
798
799
800
801
802
803
804 public System getSystem() {
805 return system;
806 }
807
808
809 @Override
810 public double getd2EdL2() {
811 return 0.0;
812 }
813
814
815 @Override
816 public double getdEdL() {
817
818 if (!lambdaTerm) {
819 return 0.0;
820 }
821
822
823 double[] x = new double[getNumberOfVariables()];
824 getCoordinates(x);
825
826 double currentLambda = getLambda();
827 double width = finiteDifferenceStepSize;
828 double ePlus;
829 double eMinus;
830
831 if (twoSidedFiniteDifference) {
832 if (currentLambda + finiteDifferenceStepSize > 1.0) {
833 setLambda(currentLambda - finiteDifferenceStepSize);
834 eMinus = energy(x);
835 setLambda(currentLambda);
836 ePlus = energy(x);
837 } else if (currentLambda - finiteDifferenceStepSize < 0.0) {
838 setLambda(currentLambda + finiteDifferenceStepSize);
839 ePlus = energy(x);
840 setLambda(currentLambda);
841 eMinus = energy(x);
842 } else {
843
844 setLambda(currentLambda + finiteDifferenceStepSize);
845 ePlus = energy(x);
846 setLambda(currentLambda - finiteDifferenceStepSize);
847 eMinus = energy(x);
848 width *= 2.0;
849 setLambda(currentLambda);
850 }
851 } else {
852
853 if (currentLambda + finiteDifferenceStepSize > 1.0) {
854 setLambda(currentLambda - finiteDifferenceStepSize);
855 eMinus = energy(x);
856 setLambda(currentLambda);
857 ePlus = energy(x);
858 } else {
859 setLambda(currentLambda + finiteDifferenceStepSize);
860 ePlus = energy(x);
861 setLambda(currentLambda);
862 eMinus = energy(x);
863 }
864 }
865
866
867 double dEdL = (ePlus - eMinus) / width;
868
869
870
871 return dEdL;
872 }
873
874
875 @Override
876 public void getdEdXdL(double[] gradients) {
877
878 }
879
880
881 public void setActiveAtoms() {
882 system.updateAtomMass();
883
884
885 }
886
887
888
889
890
891
892 @Override
893 public void setCoordinates(double[] x) {
894
895 context.setOpenMMPositions(x);
896 }
897
898
899 @Override
900 public void setCrystal(Crystal crystal) {
901 super.setCrystal(crystal);
902 context.setPeriodicBoxVectors();
903 }
904
905
906 @Override
907 public void setLambda(double lambda) {
908
909 if (!lambdaTerm) {
910 logger.fine(" Attempting to set lambda for a ForceFieldEnergyOpenMM with lambdaterm false.");
911 return;
912 }
913
914
915 if (lambda < 0.0 || lambda > 1.0) {
916 String message = format(" Lambda value %8.3f is not in the range [0..1].", lambda);
917 logger.warning(message);
918 return;
919 }
920
921 super.setLambda(lambda);
922
923
924 double mappedLambda = lambda;
925 if (lambdaStart > 0) {
926 double windowSize = 1.0 - lambdaStart;
927 mappedLambda = lambdaStart + lambda * windowSize;
928 }
929
930 if (system != null) {
931 system.setLambda(mappedLambda);
932 if (atoms != null) {
933 List<Atom> atomList = new ArrayList<>();
934 for (Atom atom : atoms) {
935 if (atom.applyLambda()) {
936 atomList.add(atom);
937 }
938 }
939
940 updateParameters(atomList.toArray(new Atom[0]));
941 } else {
942 updateParameters(null);
943 }
944 }
945 }
946
947
948
949
950
951
952 public void updateParameters(Atom[] atoms) {
953 if (atoms == null) {
954 atoms = this.atoms;
955 }
956 system.updateParameters(atoms);
957 }
958
959
960 private void free() {
961 if (context != null) {
962 context.free();
963 context = null;
964 }
965 if (system != null) {
966 system.free();
967 system = null;
968 }
969 }
970
971
972
973
974
975
976
977
978
979
980
981
982
983 public class Context {
984
985
986 private final Platform platform;
987
988 private final Integrator integrator;
989
990 private final double constraintTolerance = ForceFieldEnergy.DEFAULT_CONSTRAINT_TOLERANCE;
991
992 private PointerByReference contextPointer = null;
993
994 private String integratorString = "VERLET";
995
996 private double timeStep = 0.001;
997
998 private PointerByReference platformPointer = null;
999
1000 private double temperature = 298.15;
1001
1002
1003
1004
1005
1006
1007
1008 Context(ForceField forceField, Platform requestedPlatform) {
1009 loadPlatform(requestedPlatform);
1010 platform = requestedPlatform;
1011 integrator = new Integrator(forceField, constraintTolerance);
1012 }
1013
1014
1015
1016
1017
1018
1019 public PointerByReference getContextPointer() {
1020 if (contextPointer == null) {
1021 create(integratorString, timeStep, temperature, true);
1022 }
1023 return contextPointer;
1024 }
1025
1026
1027
1028
1029
1030
1031 public PointerByReference getIntegrator() {
1032 return integrator.getIntegratorPointer();
1033 }
1034
1035
1036
1037
1038
1039
1040 public void integrate(int numSteps) {
1041 OpenMM_Integrator_step(integrator.getIntegratorPointer(), numSteps);
1042 }
1043
1044
1045
1046
1047
1048
1049
1050 public void optimize(double eps, int maxIterations) {
1051 OpenMM_LocalEnergyMinimizer_minimize(contextPointer,
1052 eps / (OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ), maxIterations);
1053 }
1054
1055
1056
1057
1058
1059
1060 public void setOpenMMPositions(double[] x) {
1061 PointerByReference positions = OpenMM_Vec3Array_create(0);
1062 OpenMM_Vec3.ByValue coords = new OpenMM_Vec3.ByValue();
1063 double[] d = new double[3];
1064 int index = 0;
1065 for (Atom a : atoms) {
1066 if (a.isActive()) {
1067 a.moveTo(x[index++], x[index++], x[index++]);
1068 a.getXYZ(d);
1069 coords.x = d[0] * OpenMM_NmPerAngstrom;
1070 coords.y = d[1] * OpenMM_NmPerAngstrom;
1071 coords.z = d[2] * OpenMM_NmPerAngstrom;
1072 OpenMM_Vec3Array_append(positions, coords);
1073 } else {
1074
1075 coords.x = a.getX() * OpenMM_NmPerAngstrom;
1076 coords.y = a.getY() * OpenMM_NmPerAngstrom;
1077 coords.z = a.getZ() * OpenMM_NmPerAngstrom;
1078 OpenMM_Vec3Array_append(positions, coords);
1079 }
1080 }
1081 OpenMM_Context_setPositions(contextPointer, positions);
1082 if (freeOpenMM) {
1083 logger.finer(" Free OpenMM positions.");
1084 OpenMM_Vec3Array_destroy(positions);
1085 logger.finer(" Free OpenMM positions completed.");
1086 }
1087 }
1088
1089
1090
1091
1092
1093
1094 public void setOpenMMVelocities(double[] v) {
1095 PointerByReference velocities = OpenMM_Vec3Array_create(0);
1096 OpenMM_Vec3.ByValue vel = new OpenMM_Vec3.ByValue();
1097 int index = 0;
1098 double[] velocity = new double[3];
1099 for (Atom a : atoms) {
1100 if (a.isActive()) {
1101 a.setVelocity(v[index++], v[index++], v[index++]);
1102 a.getVelocity(velocity);
1103 vel.x = velocity[0] * OpenMM_NmPerAngstrom;
1104 vel.y = velocity[1] * OpenMM_NmPerAngstrom;
1105 vel.z = velocity[2] * OpenMM_NmPerAngstrom;
1106 OpenMM_Vec3Array_append(velocities, vel);
1107 } else {
1108
1109 a.setVelocity(0.0, 0.0, 0.0);
1110 vel.x = 0.0;
1111 vel.y = 0.0;
1112 vel.z = 0.0;
1113 OpenMM_Vec3Array_append(velocities, vel);
1114 }
1115 }
1116 OpenMM_Context_setVelocities(contextPointer, velocities);
1117
1118 if (freeOpenMM) {
1119 logger.finer(" Free OpenMM velocities.");
1120 OpenMM_Vec3Array_destroy(velocities);
1121 logger.finer(" Free OpenMM velocities completed.");
1122 }
1123 }
1124
1125
1126 public void setPeriodicBoxVectors() {
1127 Crystal crystal = getCrystal();
1128 if (!crystal.aperiodic()) {
1129 OpenMM_Vec3 a = new OpenMM_Vec3();
1130 OpenMM_Vec3 b = new OpenMM_Vec3();
1131 OpenMM_Vec3 c = new OpenMM_Vec3();
1132 double[][] Ai = crystal.Ai;
1133 a.x = Ai[0][0] * OpenMM_NmPerAngstrom;
1134 a.y = Ai[0][1] * OpenMM_NmPerAngstrom;
1135 a.z = Ai[0][2] * OpenMM_NmPerAngstrom;
1136 b.x = Ai[1][0] * OpenMM_NmPerAngstrom;
1137 b.y = Ai[1][1] * OpenMM_NmPerAngstrom;
1138 b.z = Ai[1][2] * OpenMM_NmPerAngstrom;
1139 c.x = Ai[2][0] * OpenMM_NmPerAngstrom;
1140 c.y = Ai[2][1] * OpenMM_NmPerAngstrom;
1141 c.z = Ai[2][2] * OpenMM_NmPerAngstrom;
1142 OpenMM_Context_setPeriodicBoxVectors(contextPointer, a, b, c);
1143 }
1144 }
1145
1146
1147 @Override
1148 public String toString() {
1149 return String.format(
1150 " OpenMM context with integrator %s, timestep %9.3g fsec, temperature %9.3g K, constraintTolerance %9.3g",
1151 integratorString, timeStep, temperature, constraintTolerance);
1152 }
1153
1154
1155 void free() {
1156 if (integrator != null) {
1157 integrator.free();
1158 }
1159 if (contextPointer != null && freeOpenMM) {
1160 logger.fine(" Free OpenMM Context.");
1161 OpenMM_Context_destroy(contextPointer);
1162 logger.fine(" Free OpenMM Context completed.");
1163 contextPointer = null;
1164 }
1165 }
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176 Context create(String integratorString, double timeStep, double temperature,
1177 boolean forceCreation) {
1178
1179 if (contextPointer != null && !forceCreation) {
1180 if (this.temperature == temperature && this.timeStep == timeStep
1181 && this.integratorString.equalsIgnoreCase(integratorString)) {
1182
1183 return this;
1184 }
1185 }
1186
1187 this.integratorString = integratorString;
1188 this.timeStep = timeStep;
1189 this.temperature = temperature;
1190
1191 if (contextPointer != null) {
1192 logger.fine(" Free OpenMM Context.");
1193 OpenMM_Context_destroy(contextPointer);
1194 logger.fine(" Free OpenMM Context completed.");
1195 contextPointer = null;
1196 }
1197
1198 logger.info("\n Creating OpenMM Context");
1199
1200 PointerByReference integratorPointer = integrator.createIntegrator(integratorString,
1201 this.timeStep, temperature);
1202
1203
1204 double currentLambda = getLambda();
1205
1206 if (lambdaTerm) {
1207 ForceFieldEnergyOpenMM.this.setLambda(1.0);
1208 }
1209
1210
1211 contextPointer = OpenMM_Context_create_2(system.getSystem(), integratorPointer,
1212 platformPointer);
1213
1214
1215 if (lambdaTerm) {
1216 ForceFieldEnergyOpenMM.this.setLambda(currentLambda);
1217 }
1218
1219
1220 int nVar = ForceFieldEnergyOpenMM.super.getNumberOfVariables();
1221 double[] x = new double[nVar];
1222 double[] v = new double[nVar];
1223 double[] vel3 = new double[3];
1224 int index = 0;
1225 for (Atom a : atoms) {
1226 if (a.isActive()) {
1227 a.getVelocity(vel3);
1228
1229 x[index] = a.getX();
1230 v[index++] = vel3[0];
1231
1232 x[index] = a.getY();
1233 v[index++] = vel3[1];
1234
1235 x[index] = a.getZ();
1236 v[index++] = vel3[2];
1237 }
1238 }
1239
1240
1241 setPeriodicBoxVectors();
1242
1243
1244 setOpenMMPositions(x);
1245
1246
1247 setOpenMMVelocities(v);
1248
1249
1250 OpenMM_Context_applyConstraints(contextPointer, constraintTolerance);
1251
1252
1253
1254 State state = new State(true, false, false, true);
1255 state.getPositions(x);
1256 state.getVelocities(v);
1257 state.free();
1258
1259 return this;
1260 }
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272 void reinitContext() {
1273 if (contextPointer != null) {
1274 int preserveState = 1;
1275 OpenMM_Context_reinitialize(contextPointer, preserveState);
1276 }
1277 }
1278
1279
1280
1281
1282
1283
1284
1285 void setParameter(String name, double value) {
1286 if (contextPointer != null) {
1287 OpenMM_Context_setParameter(contextPointer, name, value);
1288 }
1289 }
1290
1291
1292 private void loadPlatform(Platform requestedPlatform) {
1293
1294 OpenMMUtils.init();
1295
1296 logger.log(Level.INFO, " Loaded from:\n {0}", OpenMMLibrary.JNA_NATIVE_LIB.toString());
1297
1298
1299 Pointer version = OpenMM_Platform_getOpenMMVersion();
1300 logger.log(Level.INFO, " Version: {0}", version.getString(0));
1301
1302
1303 logger.log(Level.FINE, " Lib Directory: {0}", OpenMMUtils.getLibDirectory());
1304
1305 PointerByReference libs = OpenMM_Platform_loadPluginsFromDirectory(
1306 OpenMMUtils.getLibDirectory());
1307 int numLibs = OpenMM_StringArray_getSize(libs);
1308 logger.log(Level.FINE, " Number of libraries: {0}", numLibs);
1309 for (int i = 0; i < numLibs; i++) {
1310 String libString = stringFromArray(libs, i);
1311 logger.log(Level.FINE, " Library: {0}", libString);
1312 }
1313 OpenMM_StringArray_destroy(libs);
1314
1315
1316 logger.log(Level.INFO, "\n Plugin Directory: {0}", OpenMMUtils.getPluginDirectory());
1317
1318 PointerByReference plugins = OpenMM_Platform_loadPluginsFromDirectory(
1319 OpenMMUtils.getPluginDirectory());
1320 int numPlugins = OpenMM_StringArray_getSize(plugins);
1321 logger.log(Level.INFO, " Number of Plugins: {0}", numPlugins);
1322 boolean cuda = false;
1323 for (int i = 0; i < numPlugins; i++) {
1324 String pluginString = stringFromArray(plugins, i);
1325 logger.log(Level.INFO, " Plugin: {0}", pluginString);
1326 if (pluginString != null) {
1327 pluginString = pluginString.toUpperCase();
1328 boolean amoebaCudaAvailable = pluginString.contains("AMOEBACUDA");
1329 if (amoebaCudaAvailable) {
1330 cuda = true;
1331 }
1332 }
1333 }
1334 OpenMM_StringArray_destroy(plugins);
1335
1336 int numPlatforms = OpenMM_Platform_getNumPlatforms();
1337 logger.log(Level.INFO, " Number of Platforms: {0}", numPlatforms);
1338
1339 if (requestedPlatform == Platform.OMM_CUDA && !cuda) {
1340 logger.severe(" The OMM_CUDA platform was requested, but is not available.");
1341 }
1342
1343
1344 if (logger.isLoggable(Level.FINE)) {
1345 PointerByReference pluginFailers = OpenMM_Platform_getPluginLoadFailures();
1346 int numFailures = OpenMM_StringArray_getSize(pluginFailers);
1347 for (int i = 0; i < numFailures; i++) {
1348 String pluginString = stringFromArray(pluginFailers, i);
1349 logger.log(Level.FINE, " Plugin load failure: {0}", pluginString);
1350 }
1351 OpenMM_StringArray_destroy(pluginFailers);
1352 }
1353
1354 String defaultPrecision = "mixed";
1355 String precision = molecularAssembly.getForceField().getString("PRECISION", defaultPrecision)
1356 .toLowerCase();
1357 precision = precision.replace("-precision", "");
1358 switch (precision) {
1359 case "double", "mixed", "single" ->
1360 logger.info(String.format(" Precision level: %s", precision));
1361 default -> {
1362 logger.info(
1363 String.format(" Could not interpret precision level %s, defaulting to %s", precision,
1364 defaultPrecision));
1365 precision = defaultPrecision;
1366 }
1367 }
1368
1369 if (cuda && requestedPlatform != Platform.OMM_REF) {
1370 int defaultDevice = getDefaultDevice(molecularAssembly.getProperties());
1371 platformPointer = OpenMM_Platform_getPlatformByName("CUDA");
1372 int deviceID = molecularAssembly.getForceField().getInteger("CUDA_DEVICE", defaultDevice);
1373 String deviceIDString = Integer.toString(deviceID);
1374
1375 OpenMM_Platform_setPropertyDefaultValue(platformPointer, pointerForString("CudaDeviceIndex"),
1376 pointerForString(deviceIDString));
1377 OpenMM_Platform_setPropertyDefaultValue(platformPointer, pointerForString("Precision"),
1378 pointerForString(precision));
1379 logger.info(String.format(" Platform: AMOEBA CUDA (Device ID %d)", deviceID));
1380 try {
1381 Comm world = Comm.world();
1382 if (world != null) {
1383 logger.info(String.format(" Running on host %s, rank %d", world.host(), world.rank()));
1384 }
1385 } catch (IllegalStateException ise) {
1386 logger.fine(" Could not find the world communicator!");
1387 }
1388 } else {
1389 platformPointer = OpenMM_Platform_getPlatformByName("Reference");
1390 logger.info(" Platform: AMOEBA CPU Reference");
1391 }
1392 }
1393 }
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404 private class Integrator {
1405
1406
1407 private final double constraintTolerance;
1408
1409 private final double frictionCoeff;
1410
1411 private PointerByReference integratorPointer = null;
1412
1413
1414
1415
1416
1417
1418
1419 Integrator(ForceField forceField, double constraintTolerance) {
1420 this.constraintTolerance = constraintTolerance;
1421 frictionCoeff = forceField.getDouble("FRICTION_COEFF", 91.0);
1422 }
1423
1424
1425
1426
1427
1428
1429 public PointerByReference getIntegratorPointer() {
1430 return integratorPointer;
1431 }
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441 PointerByReference createIntegrator(String integratorString, double timeStep,
1442 double temperature) {
1443 switch (integratorString) {
1444 default -> createVerletIntegrator(timeStep);
1445 case "LANGEVIN" -> createLangevinIntegrator(timeStep, temperature, frictionCoeff);
1446 case "MTS" -> createCustomMTSIntegrator(timeStep);
1447 case "LANGEVIN-MTS" ->
1448 createCustomMTSLangevinIntegrator(timeStep, temperature, frictionCoeff);
1449 }
1450
1451 return integratorPointer;
1452 }
1453
1454
1455
1456
1457
1458
1459
1460
1461 private void createLangevinIntegrator(double dt, double temperature, double frictionCoeff) {
1462 free();
1463 integratorPointer = OpenMM_LangevinIntegrator_create(temperature, frictionCoeff, dt);
1464 CompositeConfiguration properties = molecularAssembly.getProperties();
1465 if (properties.containsKey("integrator-seed")) {
1466 int randomSeed = properties.getInt("integrator-seed", 0);
1467 OpenMM_LangevinIntegrator_setRandomNumberSeed(integratorPointer, randomSeed);
1468 }
1469
1470 OpenMM_Integrator_setConstraintTolerance(integratorPointer, constraintTolerance);
1471 logger.info(" Langevin Integrator");
1472 logger.info(format(" Target Temperature: %6.2f (K)", temperature));
1473 logger.info(format(" Friction Coefficient: %6.2f (1/psec)", frictionCoeff));
1474 logger.info(format(" Time step: %6.2f (fsec)", dt * 1000));
1475 logger.info(format(" Degrees of Freedom: %6d", system.calculateDegreesOfFreedom()));
1476 }
1477
1478
1479
1480
1481
1482
1483 private void createCustomMTSIntegrator(double dt) {
1484 createCustomIntegrator(dt);
1485
1486 int n = 4;
1487
1488
1489 int[] forceGroups = {1, 0};
1490
1491 int[] subSteps = {1, 4};
1492 if (system.amoebaCavitationForce != null) {
1493 n = 8;
1494
1495
1496
1497 forceGroups = new int[] {2, 1, 0};
1498
1499
1500
1501 subSteps = new int[] {1, 2, 8};
1502 }
1503
1504 OpenMM_CustomIntegrator_addPerDofVariable(integratorPointer, "x1", 0.0);
1505 OpenMM_CustomIntegrator_addUpdateContextState(integratorPointer);
1506 createMTSSubStep(1, forceGroups, subSteps);
1507 OpenMM_CustomIntegrator_addConstrainVelocities(integratorPointer);
1508 logger.info(" Custom MTS Integrator");
1509 logger.info(format(" Time step: %6.2f (fsec)", dt * 1000));
1510 logger.info(format(" Inner Time step: %6.2f (fsec)", dt / n * 1000));
1511 logger.info(format(" Friction Coefficient: %6.2f", frictionCoeff));
1512 logger.info(format(" Degrees of Freedom: %6d", system.calculateDegreesOfFreedom()));
1513 }
1514
1515
1516
1517
1518
1519
1520
1521
1522 private void createMTSSubStep(int parentSubsteps, int[] forceGroups, int[] subSteps) {
1523 int forceGroup = forceGroups[0];
1524 int steps = subSteps[0];
1525 int stepsPerParentStep = steps / parentSubsteps;
1526 if (stepsPerParentStep < 1 || steps % parentSubsteps != 0) {
1527 throw new IllegalArgumentException(
1528 "The number for substeps for each group must be a multiple of the number for the previous group");
1529 }
1530 if (forceGroup < 0 || forceGroup > 31) {
1531 throw new IllegalArgumentException("Force group must be between 0 and 31");
1532 }
1533 for (int i = 0; i < stepsPerParentStep; i++) {
1534 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "v",
1535 "v+0.5*(dt/" + steps + ")*f" + forceGroup + "/m");
1536 if (forceGroups.length == 1) {
1537 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "x", "x+(dt/" + steps + ")*v");
1538 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "x1", "x");
1539 OpenMM_CustomIntegrator_addConstrainPositions(integratorPointer);
1540 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "v",
1541 "v+(x-x1)/(dt/" + steps + ")");
1542 OpenMM_CustomIntegrator_addConstrainVelocities(integratorPointer);
1543 } else {
1544 createMTSSubStep(steps, copyOfRange(forceGroups, 1, forceGroups.length),
1545 copyOfRange(subSteps, 1, subSteps.length));
1546 }
1547 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "v",
1548 "v+0.5*(dt/" + steps + ")*f" + forceGroup + "/m");
1549 }
1550 }
1551
1552
1553
1554
1555
1556
1557
1558
1559 private void createCustomMTSLangevinIntegrator(double dt, double temperature,
1560 double frictionCoeff) {
1561 createCustomIntegrator(dt);
1562
1563 int n = 4;
1564
1565
1566 int[] forceGroups = {1, 0};
1567
1568 int[] subSteps = {1, 4};
1569 if (system.amoebaCavitationForce != null) {
1570 n = 8;
1571
1572
1573
1574 forceGroups = new int[] {2, 1, 0};
1575
1576
1577
1578 subSteps = new int[] {1, 2, 8};
1579 }
1580
1581 OpenMM_CustomIntegrator_addGlobalVariable(integratorPointer, "a",
1582 exp(-frictionCoeff * dt / n));
1583 OpenMM_CustomIntegrator_addGlobalVariable(integratorPointer, "b",
1584 sqrt(1.0 - exp(-2.0 * frictionCoeff * dt / n)));
1585 OpenMM_CustomIntegrator_addGlobalVariable(integratorPointer, "kT",
1586 R * temperature * KCAL_TO_KJ);
1587 OpenMM_CustomIntegrator_addPerDofVariable(integratorPointer, "x1", 0.0);
1588 StringBuilder sb = new StringBuilder(" Update Context State\n");
1589 OpenMM_CustomIntegrator_addUpdateContextState(integratorPointer);
1590
1591 createMTSLangevinSubStep(1, forceGroups, subSteps, sb);
1592
1593 logger.finest(" Langevin-MTS steps:" + sb);
1594 OpenMM_CustomIntegrator_addConstrainVelocities(integratorPointer);
1595 logger.info(" Custom MTS Langevin Integrator");
1596 logger.info(format(" Time step: %6.2f (fsec)", dt * 1000));
1597 logger.info(format(" Inner Time step: %6.2f (fsec)", dt / n * 1000));
1598 logger.info(format(" Friction Coefficient: %6.2f (1/psec)", frictionCoeff));
1599 logger.info(format(" Degrees of Freedom: %6d", system.calculateDegreesOfFreedom()));
1600 }
1601
1602
1603
1604
1605
1606
1607
1608
1609 private void createMTSLangevinSubStep(int parentSubsteps, int[] forceGroups, int[] subSteps,
1610 StringBuilder sb) {
1611 int forceGroup = forceGroups[0];
1612 int steps = subSteps[0];
1613 int stepsPerParentStep = steps / parentSubsteps;
1614 if (stepsPerParentStep < 1 || steps % parentSubsteps != 0) {
1615 throw new IllegalArgumentException(
1616 "The number for substeps for each group must be a multiple of the number for the previous group");
1617 }
1618 if (forceGroup < 0 || forceGroup > 31) {
1619 throw new IllegalArgumentException("Force group must be between 0 and 31");
1620 }
1621
1622 sb.append(" Force Group: ").append(forceGroup).append(" ForceGroup length: ")
1623 .append(forceGroups.length).append(" Steps: ").append(steps)
1624 .append(" Step Per Parent Step: ").append(stepsPerParentStep).append(" Parent Sub Steps: ")
1625 .append(parentSubsteps).append("\n");
1626
1627 for (int i = 0; i < stepsPerParentStep; i++) {
1628 String step = "v+0.5*(dt/" + steps + ")*f" + forceGroup + "/m";
1629 sb.append(" v = ").append(step).append("\n");
1630 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "v", step);
1631
1632 if (forceGroups.length == 1) {
1633 step = "x+(dt/" + 2 * steps + ")*v";
1634 sb.append(" x = ").append(step).append("\n");
1635 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "x", step);
1636 step = "a*v + b*sqrt(kT/m)*gaussian";
1637 sb.append(" v = ").append(step).append("\n");
1638 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "v", step);
1639 step = "x+(dt/" + 2 * steps + ")*v";
1640 sb.append(" x = ").append(step).append("\n");
1641 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "x", step);
1642 step = "x";
1643 sb.append(" x1 = ").append(step).append("\n");
1644 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "x1", step);
1645 sb.append(" Constrain Positions\n");
1646 OpenMM_CustomIntegrator_addConstrainPositions(integratorPointer);
1647 step = "v+(x-x1)/(dt/" + steps + ")";
1648 sb.append(" v = ").append(step).append("\n");
1649 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "v", step);
1650 sb.append(" Constrain Velocities\n");
1651 OpenMM_CustomIntegrator_addConstrainVelocities(integratorPointer);
1652 } else {
1653 createMTSLangevinSubStep(steps, copyOfRange(forceGroups, 1, forceGroups.length),
1654 copyOfRange(subSteps, 1, subSteps.length), sb);
1655 }
1656 step = "v+0.5*(dt/" + steps + ")*f" + forceGroup + "/m";
1657 sb.append(" v = ").append(step).append("\n");
1658 OpenMM_CustomIntegrator_addComputePerDof(integratorPointer, "v", step);
1659 }
1660 }
1661
1662
1663
1664
1665
1666
1667 private void createVerletIntegrator(double dt) {
1668 free();
1669 integratorPointer = OpenMM_VerletIntegrator_create(dt);
1670 OpenMM_Integrator_setConstraintTolerance(integratorPointer, constraintTolerance);
1671 logger.info(" Verlet Integrator");
1672 logger.info(format(" Time step: %6.2f (fsec)", dt * 1000));
1673 logger.info(format(" Degrees of Freedom: %6d", system.calculateDegreesOfFreedom()));
1674 }
1675
1676
1677
1678
1679
1680
1681 private void createCustomIntegrator(double dt) {
1682 free();
1683 integratorPointer = OpenMM_CustomIntegrator_create(dt);
1684 OpenMM_Integrator_setConstraintTolerance(integratorPointer, constraintTolerance);
1685 }
1686
1687
1688 private void free() {
1689 if (integratorPointer != null && freeOpenMM) {
1690 logger.fine(" Free OpenMM Integrator.");
1691 OpenMM_Integrator_destroy(integratorPointer);
1692 logger.fine(" Free OpenMM Integrator completed.");
1693 integratorPointer = null;
1694 }
1695 }
1696 }
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712 public class System {
1713
1714 private static final double DEFAULT_MELD_SCALE_FACTOR = -1.0;
1715 private final double meldScaleFactor;
1716
1717 private final double collisionFreq;
1718
1719
1720
1721
1722 private final boolean useMeld;
1723
1724 ForceField forceField;
1725
1726 Atom[] atoms;
1727
1728 private PointerByReference system;
1729
1730 private PointerByReference ommBarostat = null;
1731
1732 private PointerByReference commRemover = null;
1733
1734
1735
1736
1737 private boolean updateBondedTerms = false;
1738
1739
1740
1741
1742 private final boolean manyBodyTitration;
1743
1744 private PointerByReference bondForce = null;
1745
1746 private PointerByReference angleForce = null;
1747
1748 private PointerByReference stretchBendForce = null;
1749
1750 private PointerByReference inPlaneAngleForce = null;
1751
1752 private PointerByReference ureyBradleyForce = null;
1753
1754 private PointerByReference outOfPlaneBendForce = null;
1755
1756 private PointerByReference piTorsionForce = null;
1757
1758 private PointerByReference torsionForce = null;
1759 private PointerByReference[] restraintTorsions = null;
1760
1761 private PointerByReference improperTorsionForce = null;
1762
1763 private PointerByReference amoebaVDWForce = null;
1764
1765 private PointerByReference amoebaMultipoleForce = null;
1766
1767 private PointerByReference amoebaGeneralizedKirkwoodForce = null;
1768
1769 private PointerByReference amoebaWcaDispersionForce = null;
1770
1771 private PointerByReference amoebaCavitationForce = null;
1772
1773 private PointerByReference customGBForce = null;
1774
1775 private PointerByReference fixedChargeNonBondedForce = null;
1776
1777 private boolean softcoreCreated = false;
1778
1779 private boolean[] chargeExclusion;
1780
1781 private boolean[] vdWExclusion;
1782
1783 private double[] exceptionChargeProd;
1784
1785 private double[] exceptionEps;
1786
1787
1788
1789 private Map<Integer, Integer> vdwClassToOpenMMType;
1790
1791
1792
1793
1794
1795 private int vdWClassForNoInteraction;
1796
1797
1798
1799
1800 private boolean elecLambdaTerm;
1801
1802
1803
1804
1805 private boolean vdwLambdaTerm;
1806
1807
1808
1809
1810 private boolean torsionLambdaTerm;
1811
1812 private double lambdaVDW = 1.0;
1813
1814 private double lambdaElec = 1.0;
1815
1816 private double lambdaTorsion = 1.0;
1817
1818
1819
1820
1821
1822
1823
1824 private double electrostaticStart = 0.6;
1825
1826 private double electrostaticLambdaPower;
1827
1828 private double vdWSoftcoreAlpha = 0.25;
1829
1830 private PointerByReference ommThermostat = null;
1831
1832 private double vdwSoftcorePower = 3.0;
1833
1834 private double torsionalLambdaPower = 2.0;
1835
1836
1837
1838
1839
1840
1841 System(MolecularAssembly molecularAssembly) {
1842
1843 system = OpenMM_System_create();
1844 logger.info("\n System created");
1845
1846 forceField = molecularAssembly.getForceField();
1847 atoms = molecularAssembly.getAtomArray();
1848
1849
1850 try {
1851 addAtoms();
1852 } catch (Exception e) {
1853 logger.severe(" Atom without mass encountered.");
1854 }
1855
1856
1857 meldScaleFactor = forceField.getDouble("MELD_SCALE_FACTOR", DEFAULT_MELD_SCALE_FACTOR);
1858 if (meldScaleFactor <= 1.0 && meldScaleFactor > 0.0) {
1859 useMeld = true;
1860 elecLambdaTerm = true;
1861 vdwLambdaTerm = true;
1862 torsionLambdaTerm = true;
1863 } else {
1864 useMeld = false;
1865 elecLambdaTerm = false;
1866 vdwLambdaTerm = false;
1867 torsionLambdaTerm = false;
1868 }
1869
1870
1871 elecLambdaTerm = forceField.getBoolean("ELEC_LAMBDATERM", elecLambdaTerm);
1872 vdwLambdaTerm = forceField.getBoolean("VDW_LAMBDATERM", vdwLambdaTerm);
1873 torsionLambdaTerm = forceField.getBoolean("TORSION_LAMBDATERM", torsionLambdaTerm);
1874
1875 manyBodyTitration = forceField.getBoolean("MANYBODY_TITRATION", false);
1876
1877 if (!forceField.getBoolean("LAMBDATERM", false)) {
1878 lambdaTerm = (elecLambdaTerm || vdwLambdaTerm || torsionLambdaTerm);
1879 } else {
1880 lambdaTerm = true;
1881 }
1882
1883 VanDerWaals vdW = ForceFieldEnergyOpenMM.super.getVdwNode();
1884 if (vdW != null) {
1885 vdWSoftcoreAlpha = vdW.getAlpha();
1886 vdwSoftcorePower = (int) vdW.getBeta();
1887 }
1888
1889
1890 lambdaStart = forceField.getDouble("LAMBDA_START", 0.0);
1891 if (lambdaStart > 1.0) {
1892 lambdaStart = 1.0;
1893 } else if (lambdaStart < 0.0) {
1894 lambdaStart = 0.0;
1895 }
1896
1897 electrostaticStart = forceField.getDouble("PERMANENT_LAMBDA_START", electrostaticStart);
1898 if (electrostaticStart > 1.0) {
1899 electrostaticStart = 1.0;
1900 } else if (electrostaticStart < 0.0) {
1901 electrostaticStart = 0.0;
1902 }
1903 electrostaticLambdaPower = forceField.getDouble("PERMANENT_LAMBDA_EXPONENT", 2.0);
1904
1905 if (useMeld) {
1906
1907 lambdaStart = 0.0;
1908
1909 electrostaticStart = 0.0;
1910
1911 electrostaticLambdaPower = 1.0;
1912
1913 vdwSoftcorePower = 1;
1914
1915 vdWSoftcoreAlpha = 0.0;
1916
1917 torsionalLambdaPower = 1.0;
1918
1919 twoSidedFiniteDifference = false;
1920 }
1921
1922 collisionFreq = forceField.getDouble("COLLISION_FREQ", 0.1);
1923
1924
1925
1926 boolean rigidHydrogen = forceField.getBoolean("RIGID_HYDROGEN", false);
1927 boolean rigidBonds = forceField.getBoolean("RIGID_BONDS", false);
1928 boolean rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
1929 if (rigidHydrogen) {
1930 addHydrogenConstraints();
1931 }
1932 if (rigidBonds) {
1933 addUpBondConstraints();
1934 }
1935 if (rigidHydrogenAngles) {
1936 setUpHydrogenAngleConstraints();
1937 }
1938
1939 logger.info("\n Bonded Terms\n");
1940
1941
1942 if (rigidBonds) {
1943 logger.info(" Not creating AmoebaBondForce because bonds are constrained.");
1944 } else {
1945 addBondForce();
1946 }
1947
1948
1949 addAngleForce();
1950 addInPlaneAngleForce();
1951
1952
1953 addStretchBendForce();
1954
1955
1956 addUreyBradleyForce();
1957
1958
1959 addOutOfPlaneBendForce();
1960
1961
1962 addTorsionForce();
1963
1964
1965 addImproperTorsionForce();
1966
1967
1968 addPiTorsionForce();
1969
1970
1971 addTorsionTorsionForce();
1972
1973
1974 addHarmonicRestraintForce();
1975
1976
1977 addRestraintBondForce();
1978
1979
1980 addRestrainGroupsForce();
1981
1982
1983 addStretchTorsionForce();
1984
1985
1986 addAngleTorsionForce();
1987
1988 setDefaultPeriodicBoxVectors();
1989
1990 addRestraintTorsions();
1991
1992 if (vdW != null) {
1993 logger.info("\n Non-Bonded Terms\n");
1994 VanDerWaalsForm vdwForm = vdW.getVDWForm();
1995 if (vdwForm.vdwType == LENNARD_JONES) {
1996 addFixedChargeNonBondedForce();
1997 } else {
1998
1999 addAmoebaVDWForce();
2000
2001
2002 addAmoebaMultipoleForce();
2003 }
2004 }
2005
2006 if (lambdaTerm) {
2007 logger.info(format("\n Lambda path start: %6.3f", lambdaStart));
2008 logger.info(format(" Lambda scales torsions: %s", torsionLambdaTerm));
2009 if (torsionLambdaTerm) {
2010 logger.info(format(" torsion lambda power: %6.3f", torsionalLambdaPower));
2011 }
2012 logger.info(format(" Lambda scales vdW interactions: %s", vdwLambdaTerm));
2013 if (vdwLambdaTerm) {
2014 logger.info(format(" van Der Waals alpha: %6.3f", vdWSoftcoreAlpha));
2015 logger.info(format(" van Der Waals lambda power: %6.3f", vdwSoftcorePower));
2016 }
2017 logger.info(format(" Lambda scales electrostatics: %s", elecLambdaTerm));
2018
2019 if (elecLambdaTerm) {
2020 logger.info(format(" Electrostatics start: %6.3f", electrostaticStart));
2021 logger.info(format(" Electrostatics lambda power: %6.3f", electrostaticLambdaPower));
2022 }
2023 logger.info(format(" Using Meld: %s", useMeld));
2024 if (useMeld) {
2025 logger.info(format(" Meld scale factor: %6.3f", meldScaleFactor));
2026 }
2027 }
2028 }
2029
2030
2031
2032
2033
2034
2035 public void addAndersenThermostatForce(double targetTemp) {
2036 addAndersenThermostatForce(targetTemp, collisionFreq);
2037 }
2038
2039
2040
2041
2042
2043
2044
2045 public void addAndersenThermostatForce(double targetTemp, double collisionFreq) {
2046 if (ommThermostat == null) {
2047 ommThermostat = OpenMM_AndersenThermostat_create(targetTemp, collisionFreq);
2048 OpenMM_System_addForce(system, ommThermostat);
2049 logger.info("\n Adding an Andersen thermostat");
2050 logger.info(format(" Target Temperature: %6.2f (K)", targetTemp));
2051 logger.info(format(" Collision Frequency: %6.2f (1/psec)", collisionFreq));
2052 } else {
2053 OpenMM_AndersenThermostat_setDefaultTemperature(ommThermostat, targetTemp);
2054 OpenMM_AndersenThermostat_setDefaultCollisionFrequency(ommThermostat, collisionFreq);
2055 logger.fine(" Updated the Andersen thermostat");
2056 logger.fine(format(" Target Temperature: %6.2f (K)", targetTemp));
2057 logger.fine(format(" Collision Frequency: %6.2f (1/psec)", collisionFreq));
2058 }
2059 }
2060
2061
2062 public void addCOMMRemoverForce() {
2063 int frequency = 100;
2064 if (commRemover == null) {
2065 commRemover = OpenMM_CMMotionRemover_create(frequency);
2066 OpenMM_System_addForce(system, commRemover);
2067 logger.info("\n Adding a center of mass motion remover");
2068 logger.info(format(" Frequency: %6d", frequency));
2069 }
2070 }
2071
2072
2073
2074
2075
2076
2077
2078
2079 public void addMonteCarloBarostatForce(double targetPressure, double targetTemp, int frequency) {
2080 if (ommBarostat == null) {
2081 double pressureInBar = targetPressure * Constants.ATM_TO_BAR;
2082 ommBarostat = OpenMM_MonteCarloBarostat_create(pressureInBar, targetTemp, frequency);
2083 CompositeConfiguration properties = molecularAssembly.getProperties();
2084 if (properties.containsKey("barostat-seed")) {
2085 int randomSeed = properties.getInt("barostat-seed", 0);
2086 logger.info(format(" Setting random seed %d for Monte Carlo Barostat", randomSeed));
2087 OpenMM_MonteCarloBarostat_setRandomNumberSeed(ommBarostat, randomSeed);
2088 }
2089 OpenMM_System_addForce(system, ommBarostat);
2090 logger.info("\n Adding a Monte Carlo barostat");
2091 logger.info(format(" Target Pressure: %6.2f (atm)", targetPressure));
2092 logger.info(format(" Target Temperature: %6.2f (K)", targetTemp));
2093 logger.info(format(" MC Move Frequency: %6d", frequency));
2094 } else {
2095 logger.fine("\n Updating the Monte Carlo barostat");
2096 logger.fine(format(" Target Pressure: %6.2f (atm)", targetPressure));
2097 logger.fine(format(" Target Temperature: %6.2f (K)", targetTemp));
2098 logger.fine(format(" MC Move Frequency: %6d", frequency));
2099 }
2100 }
2101
2102
2103
2104
2105
2106
2107 public int calculateDegreesOfFreedom() {
2108
2109 int dof = getNumberOfVariables();
2110
2111 dof = dof - OpenMM_System_getNumConstraints(system);
2112
2113 if (commRemover != null) {
2114 dof -= 3;
2115 }
2116 return dof;
2117 }
2118
2119
2120 public void free() {
2121 if (system != null && freeOpenMM) {
2122 logger.fine(" Free OpenMM system.");
2123 OpenMM_System_destroy(system);
2124 logger.fine(" Free OpenMM system completed.");
2125 system = null;
2126 }
2127 }
2128
2129
2130 public void printLambdaValues() {
2131 logger.info(format("\n Lambda Values\n Torsion: %6.3f vdW: %6.3f Elec: %6.3f ", lambdaTorsion,
2132 lambdaVDW, lambdaElec));
2133 }
2134
2135
2136
2137
2138
2139
2140 public void setLambda(double lambda) {
2141
2142
2143 lambdaTorsion = 1.0;
2144
2145
2146 lambdaVDW = 1.0;
2147
2148
2149 lambdaElec = 1.0;
2150
2151 if (torsionLambdaTerm) {
2152
2153 lambdaTorsion = pow(lambda, torsionalLambdaPower);
2154 if (useMeld) {
2155 lambdaTorsion = meldScaleFactor + lambda * (1.0 - meldScaleFactor);
2156 }
2157 }
2158
2159 if (elecLambdaTerm && vdwLambdaTerm) {
2160
2161 if (lambda < electrostaticStart) {
2162
2163 lambdaElec = 0.0;
2164 } else {
2165
2166 double elecWindow = 1.0 - electrostaticStart;
2167 lambdaElec = (lambda - electrostaticStart) / elecWindow;
2168 lambdaElec = pow(lambdaElec, electrostaticLambdaPower);
2169 }
2170 lambdaVDW = lambda;
2171 if (useMeld) {
2172 lambdaElec = sqrt(meldScaleFactor + lambda * (1.0 - meldScaleFactor));
2173 lambdaVDW = meldScaleFactor + lambda * (1.0 - meldScaleFactor);
2174 }
2175 } else if (vdwLambdaTerm) {
2176
2177 lambdaElec = 0.0;
2178 lambdaVDW = lambda;
2179 if (useMeld) {
2180 lambdaVDW = meldScaleFactor + lambda * (1.0 - meldScaleFactor);
2181 }
2182
2183 } else if (elecLambdaTerm) {
2184
2185 lambdaElec = lambda;
2186 if (useMeld) {
2187 lambdaElec = sqrt(meldScaleFactor + lambda * (1.0 - meldScaleFactor));
2188 }
2189 }
2190 }
2191
2192 public void setUpdateBondedTerms(boolean updateBondedTerms) {
2193 this.updateBondedTerms = updateBondedTerms;
2194 }
2195
2196
2197
2198
2199
2200
2201 PointerByReference getSystem() {
2202 return system;
2203 }
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215 private void setDefaultPeriodicBoxVectors() {
2216 Crystal crystal = getCrystal();
2217 if (!crystal.aperiodic()) {
2218 OpenMM_Vec3 a = new OpenMM_Vec3();
2219 OpenMM_Vec3 b = new OpenMM_Vec3();
2220 OpenMM_Vec3 c = new OpenMM_Vec3();
2221 double[][] Ai = crystal.Ai;
2222 a.x = Ai[0][0] * OpenMM_NmPerAngstrom;
2223 a.y = Ai[0][1] * OpenMM_NmPerAngstrom;
2224 a.z = Ai[0][2] * OpenMM_NmPerAngstrom;
2225 b.x = Ai[1][0] * OpenMM_NmPerAngstrom;
2226 b.y = Ai[1][1] * OpenMM_NmPerAngstrom;
2227 b.z = Ai[1][2] * OpenMM_NmPerAngstrom;
2228 c.x = Ai[2][0] * OpenMM_NmPerAngstrom;
2229 c.y = Ai[2][1] * OpenMM_NmPerAngstrom;
2230 c.z = Ai[2][2] * OpenMM_NmPerAngstrom;
2231 OpenMM_System_setDefaultPeriodicBoxVectors(system, a, b, c);
2232 }
2233 }
2234
2235
2236
2237
2238
2239
2240 void updateParameters(Atom[] atoms) {
2241
2242 if (vdwLambdaTerm) {
2243 if (fixedChargeNonBondedForce != null) {
2244 if (!softcoreCreated) {
2245 addCustomNonbondedSoftcoreForce();
2246
2247 context.reinitContext();
2248 softcoreCreated = true;
2249 }
2250 context.setParameter("vdw_lambda", lambdaVDW);
2251 } else if (amoebaVDWForce != null) {
2252 context.setParameter("AmoebaVdwLambda", lambdaVDW);
2253 if (softcoreCreated) {
2254
2255
2256 if (!torsionLambdaTerm && !elecLambdaTerm) {
2257 return;
2258 }
2259 } else {
2260 softcoreCreated = true;
2261 }
2262 }
2263 }
2264
2265
2266
2267
2268 if (updateBondedTerms) {
2269 if (bondForce != null) {
2270 updateBondForce();
2271 }
2272 if (angleForce != null) {
2273 updateAngleForce();
2274 }
2275 if (stretchBendForce != null) {
2276 updateStretchBendForce();
2277 }
2278 if (inPlaneAngleForce != null) {
2279 updateInPlaneAngleForce();
2280 }
2281 if (ureyBradleyForce != null) {
2282 updateUreyBradleyForce();
2283 }
2284 if (outOfPlaneBendForce != null) {
2285 updateOutOfPlaneBendForce();
2286 }
2287 if (piTorsionForce != null) {
2288 updatePiTorsionForce();
2289 }
2290 }
2291
2292 if (torsionLambdaTerm || updateBondedTerms) {
2293 if (torsionForce != null) {
2294 updateTorsionForce();
2295 }
2296 if (improperTorsionForce != null) {
2297 updateImproperTorsionForce();
2298 }
2299 }
2300
2301 if (restraintTorsions != null && restraintTorsions.length > 0) {
2302 updateRestraintTorsions();
2303 }
2304
2305 if (atoms == null || atoms.length == 0) {
2306 return;
2307 }
2308
2309
2310 if (fixedChargeNonBondedForce != null) {
2311 updateFixedChargeNonBondedForce(atoms);
2312 }
2313
2314
2315 if (customGBForce != null) {
2316 updateCustomGBForce(atoms);
2317 }
2318
2319
2320 if (amoebaVDWForce != null) {
2321 updateAmoebaVDWForce(atoms);
2322 }
2323
2324
2325 if (amoebaMultipoleForce != null) {
2326 updateAmoebaMultipoleForce(atoms);
2327 }
2328
2329
2330 if (amoebaGeneralizedKirkwoodForce != null) {
2331 updateGeneralizedKirkwoodForce(atoms);
2332 }
2333
2334
2335 if (amoebaWcaDispersionForce != null) {
2336 updateWCAForce(atoms);
2337 }
2338
2339
2340 if (amoebaCavitationForce != null) {
2341 updateCavitationForce(atoms);
2342 }
2343 }
2344
2345
2346
2347
2348
2349 private void addAtoms() throws Exception {
2350 double totalMass = 0.0;
2351 for (Atom atom : atoms) {
2352 double mass = atom.getMass();
2353 totalMass += mass;
2354 if (mass < 0.0) {
2355 throw new Exception(" Atom with mass less than 0.");
2356 }
2357 if (mass == 0.0) {
2358 logger.info(format(" Atom %s has zero mass.", atom));
2359 }
2360 OpenMM_System_addParticle(system, mass);
2361 }
2362 logger.log(Level.INFO, format(" Atoms \t\t%6d", nAtoms));
2363 logger.log(Level.INFO, format(" Mass \t\t%12.3f", totalMass));
2364 }
2365
2366
2367 private void updateAtomMass() {
2368 int index = 0;
2369 for (Atom atom : atoms) {
2370 double mass = 0.0;
2371 if (atom.isActive()) {
2372 mass = atom.getMass();
2373 }
2374 OpenMM_System_setParticleMass(system, index++, mass);
2375 }
2376 }
2377
2378
2379 private void addBondForce() {
2380 Bond[] bonds = getBonds();
2381 if (bonds == null || bonds.length < 1) {
2382 return;
2383 }
2384
2385 String energy;
2386 BondType bondType = bonds[0].getBondType();
2387 if (bondType.bondFunction == BondFunction.QUARTIC) {
2388 energy = format("k*(d^2 + %.15g*d^3 + %.15g*d^4); d=r-r0",
2389 bondType.cubic / OpenMM_NmPerAngstrom,
2390 bondType.quartic / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
2391 } else {
2392 energy = "k*(d^2); d=r-r0";
2393 }
2394 bondForce = OpenMM_CustomBondForce_create(energy);
2395 OpenMM_CustomBondForce_addPerBondParameter(bondForce, "r0");
2396 OpenMM_CustomBondForce_addPerBondParameter(bondForce, "k");
2397 OpenMM_Force_setName(bondForce, "AmoebaBond");
2398
2399 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
2400
2401 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2402 for (Bond bond : bonds) {
2403 int i1 = bond.getAtom(0).getXyzIndex() - 1;
2404 int i2 = bond.getAtom(1).getXyzIndex() - 1;
2405 bondType = bond.bondType;
2406 double r0 = bondType.distance * OpenMM_NmPerAngstrom;
2407 double k = kParameterConversion * bondType.forceConstant * bond.bondType.bondUnit;
2408 OpenMM_DoubleArray_append(parameters, r0);
2409 OpenMM_DoubleArray_append(parameters, k);
2410 OpenMM_CustomBondForce_addBond(bondForce, i1, i2, parameters);
2411 OpenMM_DoubleArray_resize(parameters, 0);
2412 }
2413 OpenMM_DoubleArray_destroy(parameters);
2414
2415 int forceGroup = forceField.getInteger("BOND_FORCE_GROUP", 0);
2416 OpenMM_Force_setForceGroup(bondForce, forceGroup);
2417 OpenMM_System_addForce(system, bondForce);
2418 logger.log(Level.INFO, format(" Bonds \t\t%6d\t\t%1d", bonds.length, forceGroup));
2419 }
2420
2421
2422 private void updateBondForce() {
2423 Bond[] bonds = getBonds();
2424 if (bonds == null || bonds.length < 1) {
2425 return;
2426 }
2427
2428 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
2429
2430 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2431 int index = 0;
2432 for (Bond bond : bonds) {
2433 int i1 = bond.getAtom(0).getXyzIndex() - 1;
2434 int i2 = bond.getAtom(1).getXyzIndex() - 1;
2435 BondType bondType = bond.bondType;
2436 double r0 = bondType.distance * OpenMM_NmPerAngstrom;
2437 double k = kParameterConversion * bondType.forceConstant * bondType.bondUnit;
2438 OpenMM_DoubleArray_append(parameters, r0);
2439 OpenMM_DoubleArray_append(parameters, k);
2440 OpenMM_CustomBondForce_setBondParameters(bondForce, index++, i1, i2, parameters);
2441 OpenMM_DoubleArray_resize(parameters, 0);
2442 }
2443 OpenMM_DoubleArray_destroy(parameters);
2444
2445 if (context.contextPointer != null) {
2446 OpenMM_CustomBondForce_updateParametersInContext(bondForce, context.contextPointer);
2447 }
2448 }
2449
2450
2451 private void addAngleForce() {
2452 Angle[] angles = getAngles();
2453 if (angles == null || angles.length < 1) {
2454 return;
2455 }
2456 boolean rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
2457 String energy;
2458 AngleType angleType = angles[0].angleType;
2459 if (angleType.angleFunction == AngleFunction.SEXTIC) {
2460 energy = format(
2461 "k*(d^2 + %.15g*d^3 + %.15g*d^4 + %.15g*d^5 + %.15g*d^6); d=%.15g*theta-theta0",
2462 angleType.cubic, angleType.quartic, angleType.pentic, angleType.sextic, 180.0 / PI);
2463 } else {
2464 energy = format("k*(d^2); d=%.15g*theta-theta0", 180.0 / PI);
2465 }
2466 angleForce = OpenMM_CustomAngleForce_create(energy);
2467 OpenMM_CustomAngleForce_addPerAngleParameter(angleForce, "theta0");
2468 OpenMM_CustomAngleForce_addPerAngleParameter(angleForce, "k");
2469 OpenMM_Force_setName(angleForce, "Angle");
2470
2471 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2472 int angleCount = 0;
2473 for (Angle angle : angles) {
2474 AngleMode angleMode = angle.angleType.angleMode;
2475
2476 if (!manyBodyTitration && angleMode == AngleMode.IN_PLANE) {
2477
2478 } else if (isHydrogenAngle(angle) && rigidHydrogenAngles) {
2479 logger.log(Level.INFO, " Constrained angle %s was not added the AngleForce.", angle);
2480 } else {
2481 int i1 = angle.getAtom(0).getXyzIndex() - 1;
2482 int i2 = angle.getAtom(1).getXyzIndex() - 1;
2483 int i3 = angle.getAtom(2).getXyzIndex() - 1;
2484 double theta0 = angle.angleType.angle[angle.nh];
2485 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
2486 if (angleMode == AngleMode.IN_PLANE) {
2487
2488
2489 k = 0.0;
2490 }
2491 OpenMM_DoubleArray_append(parameters, theta0);
2492 OpenMM_DoubleArray_append(parameters, k);
2493 OpenMM_CustomAngleForce_addAngle(angleForce, i1, i2, i3, parameters);
2494 angleCount++;
2495 OpenMM_DoubleArray_resize(parameters, 0);
2496 }
2497 }
2498 OpenMM_DoubleArray_destroy(parameters);
2499
2500 int forceGroup = forceField.getInteger("ANGLE_FORCE_GROUP", 0);
2501 OpenMM_Force_setForceGroup(angleForce, forceGroup);
2502 OpenMM_System_addForce(system, angleForce);
2503 logger.log(Level.INFO, format(" Angles \t\t%6d\t\t%1d", angleCount, forceGroup));
2504 }
2505
2506
2507 private void updateAngleForce() {
2508 Angle[] angles = getAngles();
2509 if (angles == null || angles.length < 1) {
2510 return;
2511 }
2512 boolean rigidHydrogenAngles = forceField.getBoolean("RIGID_HYDROGEN_ANGLES", false);
2513 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2514 int index = 0;
2515 for (Angle angle : angles) {
2516 AngleMode angleMode = angle.angleType.angleMode;
2517 if (!manyBodyTitration && angleMode == AngleMode.IN_PLANE) {
2518
2519 }
2520
2521 else if (!rigidHydrogenAngles || !isHydrogenAngle(angle)) {
2522 int i1 = angle.getAtom(0).getXyzIndex() - 1;
2523 int i2 = angle.getAtom(1).getXyzIndex() - 1;
2524 int i3 = angle.getAtom(2).getXyzIndex() - 1;
2525 double theta0 = angle.angleType.angle[angle.nh];
2526 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
2527 if (angleMode == AngleMode.IN_PLANE) {
2528
2529 k = 0.0;
2530 }
2531 OpenMM_DoubleArray_append(parameters, theta0);
2532 OpenMM_DoubleArray_append(parameters, k);
2533 OpenMM_CustomAngleForce_setAngleParameters(angleForce, index++, i1, i2, i3, parameters);
2534 OpenMM_DoubleArray_resize(parameters, 0);
2535 }
2536 }
2537 OpenMM_DoubleArray_destroy(parameters);
2538
2539 if (context.contextPointer != null) {
2540 OpenMM_CustomAngleForce_updateParametersInContext(angleForce, context.contextPointer);
2541 }
2542 }
2543
2544
2545 private void addInPlaneAngleForce() {
2546 Angle[] angles = getAngles();
2547 if (angles == null || angles.length < 1) {
2548 return;
2549 }
2550
2551 AngleType angleType = angles[0].angleType;
2552 String energy = format("""
2553 k*(d^2 + %.15g*d^3 + %.15g*d^4 + %.15g*d^5 + %.15g*d^6);
2554 d=theta-theta0;
2555 theta = %.15g*pointangle(x1, y1, z1, projx, projy, projz, x3, y3, z3);
2556 projx = x2-nx*dot;
2557 projy = y2-ny*dot;
2558 projz = z2-nz*dot;
2559 dot = nx*(x2-x3) + ny*(y2-y3) + nz*(z2-z3);
2560 nx = px/norm;
2561 ny = py/norm;
2562 nz = pz/norm;
2563 norm = sqrt(px*px + py*py + pz*pz);
2564 px = (d1y*d2z-d1z*d2y);
2565 py = (d1z*d2x-d1x*d2z);
2566 pz = (d1x*d2y-d1y*d2x);
2567 d1x = x1-x4;
2568 d1y = y1-y4;
2569 d1z = z1-z4;
2570 d2x = x3-x4;
2571 d2y = y3-y4;
2572 d2z = z3-z4;
2573 """, angleType.cubic, angleType.quartic, angleType.pentic, angleType.sextic, 180.0 / PI);
2574 inPlaneAngleForce = OpenMM_CustomCompoundBondForce_create(4, energy);
2575 OpenMM_CustomCompoundBondForce_addPerBondParameter(inPlaneAngleForce, "theta0");
2576 OpenMM_CustomCompoundBondForce_addPerBondParameter(inPlaneAngleForce, "k");
2577 OpenMM_Force_setName(inPlaneAngleForce, "InPlaneAngle");
2578
2579 PointerByReference particles = OpenMM_IntArray_create(0);
2580 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2581 for (Angle angle : angles) {
2582 AngleMode angleMode = angle.angleType.angleMode;
2583
2584 if (!manyBodyTitration && angleMode == AngleMode.NORMAL) {
2585
2586 } else {
2587 double theta0 = angle.angleType.angle[angle.nh];
2588 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
2589 int i1 = angle.getAtom(0).getXyzIndex() - 1;
2590 int i2 = angle.getAtom(1).getXyzIndex() - 1;
2591 int i3 = angle.getAtom(2).getXyzIndex() - 1;
2592 int i4 = 0;
2593 if (angleMode == AngleMode.NORMAL) {
2594
2595
2596 k = 0.0;
2597 Atom fourthAtom = angle.getFourthAtomOfTrigonalCenter();
2598 if (fourthAtom != null) {
2599 i4 = fourthAtom.getXyzIndex() - 1;
2600 } else {
2601 while (i1 == i4 || i2 == i4 || i3 == i4) {
2602 i4++;
2603 }
2604 }
2605 } else {
2606 i4 = angle.getAtom4().getXyzIndex() - 1;
2607 }
2608 OpenMM_IntArray_append(particles, i1);
2609 OpenMM_IntArray_append(particles, i2);
2610 OpenMM_IntArray_append(particles, i3);
2611 OpenMM_IntArray_append(particles, i4);
2612 OpenMM_DoubleArray_append(parameters, theta0);
2613 OpenMM_DoubleArray_append(parameters, k);
2614 OpenMM_CustomCompoundBondForce_addBond(inPlaneAngleForce, particles, parameters);
2615 OpenMM_IntArray_resize(particles, 0);
2616 OpenMM_DoubleArray_resize(parameters, 0);
2617 }
2618 }
2619 OpenMM_IntArray_destroy(particles);
2620 OpenMM_DoubleArray_destroy(parameters);
2621
2622 int forceGroup = forceField.getInteger("IN_PLANE_ANGLE_FORCE_GROUP", 0);
2623 OpenMM_Force_setForceGroup(inPlaneAngleForce, forceGroup);
2624 OpenMM_System_addForce(system, inPlaneAngleForce);
2625 logger.log(Level.INFO, format(" In-Plane Angles \t%6d\t\t%1d", angles.length, forceGroup));
2626 }
2627
2628
2629 private void updateInPlaneAngleForce() {
2630 Angle[] angles = getAngles();
2631 if (angles == null || angles.length < 1) {
2632 return;
2633 }
2634 PointerByReference particles = OpenMM_IntArray_create(0);
2635 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2636 int index = 0;
2637 for (Angle angle : angles) {
2638 AngleMode angleMode = angle.angleType.angleMode;
2639 if (!manyBodyTitration && angleMode == AngleMode.NORMAL) {
2640
2641 } else {
2642 double theta0 = angle.angleType.angle[angle.nh];
2643 double k = OpenMM_KJPerKcal * angle.angleType.angleUnit * angle.angleType.forceConstant;
2644 int i1 = angle.getAtom(0).getXyzIndex() - 1;
2645 int i2 = angle.getAtom(1).getXyzIndex() - 1;
2646 int i3 = angle.getAtom(2).getXyzIndex() - 1;
2647
2648 int i4 = 0;
2649 if (angleMode == AngleMode.NORMAL) {
2650
2651 k = 0.0;
2652 Atom fourthAtom = angle.getFourthAtomOfTrigonalCenter();
2653 if (fourthAtom != null) {
2654 i4 = fourthAtom.getXyzIndex() - 1;
2655 } else {
2656 while (i1 == i4 || i2 == i4 || i3 == i4) {
2657 i4++;
2658 }
2659 }
2660 } else {
2661 i4 = angle.getAtom4().getXyzIndex() - 1;
2662 }
2663 OpenMM_IntArray_append(particles, i1);
2664 OpenMM_IntArray_append(particles, i2);
2665 OpenMM_IntArray_append(particles, i3);
2666 OpenMM_IntArray_append(particles, i4);
2667 OpenMM_DoubleArray_append(parameters, theta0);
2668 OpenMM_DoubleArray_append(parameters, k);
2669 OpenMM_CustomCompoundBondForce_setBondParameters(inPlaneAngleForce, index++, particles,
2670 parameters);
2671 OpenMM_IntArray_resize(particles, 0);
2672 OpenMM_DoubleArray_resize(parameters, 0);
2673 }
2674 }
2675 OpenMM_IntArray_destroy(particles);
2676 OpenMM_DoubleArray_destroy(parameters);
2677
2678 if (context.contextPointer != null) {
2679 OpenMM_CustomCompoundBondForce_updateParametersInContext(inPlaneAngleForce,
2680 context.contextPointer);
2681 }
2682 }
2683
2684
2685 private void addUreyBradleyForce() {
2686 UreyBradley[] ureyBradleys = getUreyBradleys();
2687 if (ureyBradleys == null || ureyBradleys.length < 1) {
2688 return;
2689 }
2690
2691 ureyBradleyForce = OpenMM_HarmonicBondForce_create();
2692 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
2693
2694 for (UreyBradley ureyBradley : ureyBradleys) {
2695 int i1 = ureyBradley.getAtom(0).getXyzIndex() - 1;
2696 int i2 = ureyBradley.getAtom(2).getXyzIndex() - 1;
2697 UreyBradleyType ureyBradleyType = ureyBradley.ureyBradleyType;
2698 double length = ureyBradleyType.distance * OpenMM_NmPerAngstrom;
2699
2700
2701 double k =
2702 2.0 * ureyBradleyType.forceConstant * ureyBradleyType.ureyUnit * kParameterConversion;
2703 OpenMM_HarmonicBondForce_addBond(ureyBradleyForce, i1, i2, length, k);
2704 }
2705
2706 int forceGroup = forceField.getInteger("UREY_BRADLEY_FORCE", 0);
2707 OpenMM_Force_setForceGroup(ureyBradleyForce, forceGroup);
2708 OpenMM_System_addForce(system, ureyBradleyForce);
2709 logger.log(Level.INFO,
2710 format(" Urey-Bradleys \t%6d\t\t%1d", ureyBradleys.length, forceGroup));
2711 }
2712
2713
2714 private void updateUreyBradleyForce() {
2715 UreyBradley[] ureyBradleys = getUreyBradleys();
2716 if (ureyBradleys == null || ureyBradleys.length < 1) {
2717 return;
2718 }
2719
2720 double kParameterConversion = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
2721
2722 int index = 0;
2723 for (UreyBradley ureyBradley : ureyBradleys) {
2724 int i1 = ureyBradley.getAtom(0).getXyzIndex() - 1;
2725 int i2 = ureyBradley.getAtom(2).getXyzIndex() - 1;
2726 UreyBradleyType ureyBradleyType = ureyBradley.ureyBradleyType;
2727 double length = ureyBradleyType.distance * OpenMM_NmPerAngstrom;
2728
2729
2730 double k =
2731 2.0 * ureyBradleyType.forceConstant * ureyBradleyType.ureyUnit * kParameterConversion;
2732 OpenMM_HarmonicBondForce_setBondParameters(ureyBradleyForce, index++, i1, i2, length, k);
2733 }
2734
2735 if (context.contextPointer != null) {
2736 OpenMM_HarmonicBondForce_updateParametersInContext(ureyBradleyForce, context.contextPointer);
2737 }
2738 }
2739
2740
2741 private void addOutOfPlaneBendForce() {
2742 OutOfPlaneBend[] outOfPlaneBends = getOutOfPlaneBends();
2743 if (outOfPlaneBends == null || outOfPlaneBends.length < 1) {
2744 return;
2745 }
2746
2747 OutOfPlaneBendType outOfPlaneBendType = outOfPlaneBends[0].outOfPlaneBendType;
2748 String energy = format(
2749 "k*(theta^2 + %.15g*theta^3 + %.15g*theta^4 + %.15g*theta^5 + %.15g*theta^6); "
2750 + "theta = %.15g*pointangle(x2, y2, z2, x4, y4, z4, projx, projy, projz); "
2751 + "projx = x2-nx*dot; projy = y2-ny*dot; projz = z2-nz*dot; "
2752 + "dot = nx*(x2-x3) + ny*(y2-y3) + nz*(z2-z3); "
2753 + "nx = px/norm; ny = py/norm; nz = pz/norm; " + "norm = sqrt(px*px + py*py + pz*pz); "
2754 + "px = (d1y*d2z-d1z*d2y); py = (d1z*d2x-d1x*d2z); pz = (d1x*d2y-d1y*d2x); "
2755 + "d1x = x1-x4; d1y = y1-y4; d1z = z1-z4; " + "d2x = x3-x4; d2y = y3-y4; d2z = z3-z4",
2756 outOfPlaneBendType.cubic, outOfPlaneBendType.quartic, outOfPlaneBendType.pentic,
2757 outOfPlaneBendType.sextic, 180.0 / PI);
2758 outOfPlaneBendForce = OpenMM_CustomCompoundBondForce_create(4, energy);
2759 OpenMM_CustomCompoundBondForce_addPerBondParameter(outOfPlaneBendForce, "k");
2760 OpenMM_Force_setName(outOfPlaneBendForce, "OutOfPlaneBend");
2761
2762 PointerByReference particles = OpenMM_IntArray_create(0);
2763 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2764 for (OutOfPlaneBend outOfPlaneBend : outOfPlaneBends) {
2765 outOfPlaneBendType = outOfPlaneBend.outOfPlaneBendType;
2766 int i1 = outOfPlaneBend.getAtom(0).getXyzIndex() - 1;
2767 int i2 = outOfPlaneBend.getAtom(1).getXyzIndex() - 1;
2768 int i3 = outOfPlaneBend.getAtom(2).getXyzIndex() - 1;
2769 int i4 = outOfPlaneBend.getAtom(3).getXyzIndex() - 1;
2770 double k =
2771 OpenMM_KJPerKcal * outOfPlaneBendType.forceConstant * outOfPlaneBendType.opBendUnit;
2772 OpenMM_IntArray_append(particles, i1);
2773 OpenMM_IntArray_append(particles, i2);
2774 OpenMM_IntArray_append(particles, i3);
2775 OpenMM_IntArray_append(particles, i4);
2776 OpenMM_DoubleArray_append(parameters, k);
2777 OpenMM_CustomCompoundBondForce_addBond(outOfPlaneBendForce, particles, parameters);
2778 OpenMM_IntArray_resize(particles, 0);
2779 OpenMM_DoubleArray_resize(parameters, 0);
2780 }
2781 OpenMM_IntArray_destroy(particles);
2782 OpenMM_DoubleArray_destroy(parameters);
2783 int forceGroup = forceField.getInteger("OUT_OF_PLANE_BEND_FORCE_GROUP", 0);
2784 OpenMM_Force_setForceGroup(outOfPlaneBendForce, forceGroup);
2785 OpenMM_System_addForce(system, outOfPlaneBendForce);
2786 logger.log(Level.INFO,
2787 format(" Out-of-Plane Bends \t%6d\t\t%1d", outOfPlaneBends.length, forceGroup));
2788 }
2789
2790
2791 private void updateOutOfPlaneBendForce() {
2792 OutOfPlaneBend[] outOfPlaneBends = getOutOfPlaneBends();
2793 if (outOfPlaneBends == null || outOfPlaneBends.length < 1) {
2794 return;
2795 }
2796
2797 PointerByReference particles = OpenMM_IntArray_create(0);
2798 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2799 int index = 0;
2800 for (OutOfPlaneBend outOfPlaneBend : outOfPlaneBends) {
2801 OutOfPlaneBendType outOfPlaneBendType = outOfPlaneBend.outOfPlaneBendType;
2802 int i1 = outOfPlaneBend.getAtom(0).getXyzIndex() - 1;
2803 int i2 = outOfPlaneBend.getAtom(1).getXyzIndex() - 1;
2804 int i3 = outOfPlaneBend.getAtom(2).getXyzIndex() - 1;
2805 int i4 = outOfPlaneBend.getAtom(3).getXyzIndex() - 1;
2806 double k =
2807 OpenMM_KJPerKcal * outOfPlaneBendType.forceConstant * outOfPlaneBendType.opBendUnit;
2808 OpenMM_IntArray_append(particles, i1);
2809 OpenMM_IntArray_append(particles, i2);
2810 OpenMM_IntArray_append(particles, i3);
2811 OpenMM_IntArray_append(particles, i4);
2812 OpenMM_DoubleArray_append(parameters, k);
2813 OpenMM_CustomCompoundBondForce_setBondParameters(outOfPlaneBendForce, index++, particles,
2814 parameters);
2815 OpenMM_IntArray_resize(particles, 0);
2816 OpenMM_DoubleArray_resize(parameters, 0);
2817 }
2818 OpenMM_IntArray_destroy(particles);
2819 OpenMM_DoubleArray_destroy(parameters);
2820
2821 if (context.contextPointer != null) {
2822 OpenMM_CustomCompoundBondForce_updateParametersInContext(outOfPlaneBendForce,
2823 context.contextPointer);
2824 }
2825
2826 }
2827
2828
2829 private void addStretchBendForce() {
2830 StretchBend[] stretchBends = getStretchBends();
2831 if (stretchBends == null || stretchBends.length < 1) {
2832 return;
2833 }
2834
2835 String energy = format(
2836 "(k1*(distance(p1,p2)-r12) + k2*(distance(p2,p3)-r23))*(%.15g*(angle(p1,p2,p3)-theta0))",
2837 180.0 / PI);
2838 stretchBendForce = OpenMM_CustomCompoundBondForce_create(3, energy);
2839 OpenMM_CustomCompoundBondForce_addPerBondParameter(stretchBendForce, "r12");
2840 OpenMM_CustomCompoundBondForce_addPerBondParameter(stretchBendForce, "r23");
2841 OpenMM_CustomCompoundBondForce_addPerBondParameter(stretchBendForce, "theta0");
2842 OpenMM_CustomCompoundBondForce_addPerBondParameter(stretchBendForce, "k1");
2843 OpenMM_CustomCompoundBondForce_addPerBondParameter(stretchBendForce, "k2");
2844 OpenMM_Force_setName(stretchBendForce, "AmoebaStretchBend");
2845
2846 PointerByReference particles = OpenMM_IntArray_create(0);
2847 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2848 for (StretchBend stretchBend : stretchBends) {
2849 int i1 = stretchBend.getAtom(0).getXyzIndex() - 1;
2850 int i2 = stretchBend.getAtom(1).getXyzIndex() - 1;
2851 int i3 = stretchBend.getAtom(2).getXyzIndex() - 1;
2852 double r12 = stretchBend.bond0Eq * OpenMM_NmPerAngstrom;
2853 double r23 = stretchBend.bond1Eq * OpenMM_NmPerAngstrom;
2854 double theta0 = stretchBend.angleEq * OpenMM_RadiansPerDegree;
2855 double k1 = stretchBend.force0 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
2856 double k2 = stretchBend.force1 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
2857 OpenMM_IntArray_append(particles, i1);
2858 OpenMM_IntArray_append(particles, i2);
2859 OpenMM_IntArray_append(particles, i3);
2860 OpenMM_DoubleArray_append(parameters, r12);
2861 OpenMM_DoubleArray_append(parameters, r23);
2862 OpenMM_DoubleArray_append(parameters, theta0);
2863 OpenMM_DoubleArray_append(parameters, k1);
2864 OpenMM_DoubleArray_append(parameters, k2);
2865 OpenMM_CustomCompoundBondForce_addBond(stretchBendForce, particles, parameters);
2866 OpenMM_IntArray_resize(particles, 0);
2867 OpenMM_DoubleArray_resize(parameters, 0);
2868 }
2869 OpenMM_IntArray_destroy(particles);
2870 OpenMM_DoubleArray_destroy(parameters);
2871
2872 int forceGroup = forceField.getInteger("STRETCH_BEND_FORCE_GROUP", 0);
2873 OpenMM_Force_setForceGroup(stretchBendForce, forceGroup);
2874 OpenMM_System_addForce(system, stretchBendForce);
2875 logger.log(Level.INFO,
2876 format(" Stretch-Bends \t%6d\t\t%1d", stretchBends.length, forceGroup));
2877 }
2878
2879
2880 private void updateStretchBendForce() {
2881 StretchBend[] stretchBends = getStretchBends();
2882 if (stretchBends == null || stretchBends.length < 1) {
2883 return;
2884 }
2885
2886 PointerByReference particles = OpenMM_IntArray_create(0);
2887 PointerByReference parameters = OpenMM_DoubleArray_create(0);
2888 int index = 0;
2889 for (StretchBend stretchBend : stretchBends) {
2890 int i1 = stretchBend.getAtom(0).getXyzIndex() - 1;
2891 int i2 = stretchBend.getAtom(1).getXyzIndex() - 1;
2892 int i3 = stretchBend.getAtom(2).getXyzIndex() - 1;
2893 double r12 = stretchBend.bond0Eq * OpenMM_NmPerAngstrom;
2894 double r23 = stretchBend.bond1Eq * OpenMM_NmPerAngstrom;
2895 double theta0 = stretchBend.angleEq * OpenMM_RadiansPerDegree;
2896 double k1 = stretchBend.force0 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
2897 double k2 = stretchBend.force1 * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
2898 OpenMM_IntArray_append(particles, i1);
2899 OpenMM_IntArray_append(particles, i2);
2900 OpenMM_IntArray_append(particles, i3);
2901 OpenMM_DoubleArray_append(parameters, r12);
2902 OpenMM_DoubleArray_append(parameters, r23);
2903 OpenMM_DoubleArray_append(parameters, theta0);
2904 OpenMM_DoubleArray_append(parameters, k1);
2905 OpenMM_DoubleArray_append(parameters, k2);
2906 OpenMM_CustomCompoundBondForce_setBondParameters(stretchBendForce, index++, particles,
2907 parameters);
2908 OpenMM_IntArray_resize(particles, 0);
2909 OpenMM_DoubleArray_resize(parameters, 0);
2910 }
2911 OpenMM_IntArray_destroy(particles);
2912 OpenMM_DoubleArray_destroy(parameters);
2913
2914 if (context.contextPointer != null) {
2915 OpenMM_CustomCompoundBondForce_updateParametersInContext(stretchBendForce,
2916 context.contextPointer);
2917 }
2918 }
2919
2920
2921 private void addTorsionForce() {
2922 Torsion[] torsions = getTorsions();
2923 if (torsions == null || torsions.length < 1) {
2924 return;
2925 }
2926
2927 torsionForce = OpenMM_PeriodicTorsionForce_create();
2928 for (Torsion torsion : torsions) {
2929 int a1 = torsion.getAtom(0).getXyzIndex() - 1;
2930 int a2 = torsion.getAtom(1).getXyzIndex() - 1;
2931 int a3 = torsion.getAtom(2).getXyzIndex() - 1;
2932 int a4 = torsion.getAtom(3).getXyzIndex() - 1;
2933 TorsionType torsionType = torsion.torsionType;
2934 int nTerms = torsionType.phase.length;
2935 for (int j = 0; j < nTerms; j++) {
2936 OpenMM_PeriodicTorsionForce_addTorsion(torsionForce, a1, a2, a3, a4, j + 1,
2937 torsionType.phase[j] * OpenMM_RadiansPerDegree,
2938 OpenMM_KJPerKcal * torsionType.torsionUnit * torsionType.amplitude[j]);
2939 }
2940
2941
2942 if (manyBodyTitration) {
2943 for (int j = nTerms; j < 6; j++) {
2944 OpenMM_PeriodicTorsionForce_addTorsion(torsionForce, a1, a2, a3, a4, j + 1, 0.0, 0.0);
2945 }
2946 }
2947 }
2948 int fGroup = forceField.getInteger("TORSION_FORCE_GROUP", 0);
2949
2950 OpenMM_Force_setForceGroup(torsionForce, fGroup);
2951 OpenMM_System_addForce(system, torsionForce);
2952
2953 logger.log(Level.INFO, format(" Torsions \t\t%6d\t\t%1d", torsions.length, fGroup));
2954 }
2955
2956
2957 private void updateTorsionForce() {
2958
2959 Torsion[] torsions = getTorsions();
2960 if (torsions == null || torsions.length < 1) {
2961 return;
2962 }
2963
2964 int index = 0;
2965 for (Torsion torsion : torsions) {
2966 TorsionType torsionType = torsion.torsionType;
2967 int nTerms = torsionType.phase.length;
2968 int a1 = torsion.getAtom(0).getXyzIndex() - 1;
2969 int a2 = torsion.getAtom(1).getXyzIndex() - 1;
2970 int a3 = torsion.getAtom(2).getXyzIndex() - 1;
2971 int a4 = torsion.getAtom(3).getXyzIndex() - 1;
2972 for (int j = 0; j < nTerms; j++) {
2973 double forceConstant =
2974 OpenMM_KJPerKcal * torsionType.torsionUnit * torsionType.amplitude[j] * lambdaTorsion;
2975 OpenMM_PeriodicTorsionForce_setTorsionParameters(torsionForce, index++, a1, a2, a3, a4,
2976 j + 1, torsionType.phase[j] * OpenMM_RadiansPerDegree, forceConstant);
2977 }
2978
2979
2980 if (manyBodyTitration) {
2981 for (int j = nTerms; j < 6; j++) {
2982 OpenMM_PeriodicTorsionForce_setTorsionParameters(torsionForce, index++, a1, a2, a3, a4,
2983 j + 1, 0.0, 0.0);
2984 }
2985 }
2986 }
2987
2988 if (context.contextPointer != null) {
2989 OpenMM_PeriodicTorsionForce_updateParametersInContext(torsionForce, context.contextPointer);
2990 }
2991 }
2992
2993
2994 private void addImproperTorsionForce() {
2995 ImproperTorsion[] improperTorsions = getImproperTorsions();
2996 if (improperTorsions == null || improperTorsions.length < 1) {
2997 return;
2998 }
2999
3000 improperTorsionForce = OpenMM_PeriodicTorsionForce_create();
3001 for (ImproperTorsion improperTorsion : improperTorsions) {
3002 int a1 = improperTorsion.getAtom(0).getXyzIndex() - 1;
3003 int a2 = improperTorsion.getAtom(1).getXyzIndex() - 1;
3004 int a3 = improperTorsion.getAtom(2).getXyzIndex() - 1;
3005 int a4 = improperTorsion.getAtom(3).getXyzIndex() - 1;
3006 ImproperTorsionType improperTorsionType = improperTorsion.improperType;
3007 OpenMM_PeriodicTorsionForce_addTorsion(improperTorsionForce, a1, a2, a3, a4,
3008 improperTorsionType.periodicity, improperTorsionType.phase * OpenMM_RadiansPerDegree,
3009 OpenMM_KJPerKcal * improperTorsion.improperType.impTorUnit * improperTorsion.scaleFactor
3010 * improperTorsionType.k);
3011 }
3012
3013 int forceGroup = forceField.getInteger("IMPROPER_TORSION_FORCE_GROUP", 0);
3014
3015 OpenMM_Force_setForceGroup(improperTorsionForce, forceGroup);
3016 OpenMM_System_addForce(system, improperTorsionForce);
3017
3018 logger.log(Level.INFO,
3019 format(" Improper Torsions \t%6d\t\t%1d", improperTorsions.length, forceGroup));
3020 }
3021
3022
3023 private void updateImproperTorsionForce() {
3024 ImproperTorsion[] improperTorsions = getImproperTorsions();
3025 if (improperTorsions == null || improperTorsions.length < 1) {
3026 return;
3027 }
3028
3029 int nImproperTorsions = improperTorsions.length;
3030 for (int i = 0; i < nImproperTorsions; i++) {
3031 ImproperTorsion improperTorsion = improperTorsions[i];
3032 int a1 = improperTorsion.getAtom(0).getXyzIndex() - 1;
3033 int a2 = improperTorsion.getAtom(1).getXyzIndex() - 1;
3034 int a3 = improperTorsion.getAtom(2).getXyzIndex() - 1;
3035 int a4 = improperTorsion.getAtom(3).getXyzIndex() - 1;
3036 ImproperTorsionType improperTorsionType = improperTorsion.improperType;
3037 double forceConstant =
3038 OpenMM_KJPerKcal * improperTorsion.improperType.impTorUnit * improperTorsion.scaleFactor
3039 * improperTorsionType.k * lambdaTorsion;
3040 OpenMM_PeriodicTorsionForce_setTorsionParameters(improperTorsionForce, i, a1, a2, a3, a4,
3041 improperTorsionType.periodicity, improperTorsionType.phase * OpenMM_RadiansPerDegree,
3042 forceConstant);
3043 }
3044
3045 if (context.contextPointer != null) {
3046 OpenMM_PeriodicTorsionForce_updateParametersInContext(improperTorsionForce,
3047 context.contextPointer);
3048 }
3049 }
3050
3051
3052 private void addPiTorsionForce() {
3053 PiOrbitalTorsion[] piOrbitalTorsions = getPiOrbitalTorsions();
3054 if (piOrbitalTorsions == null || piOrbitalTorsions.length < 1) {
3055 return;
3056 }
3057
3058 String energy = "2*k*sin(phi)^2;"
3059 + "phi = pointdihedral(x3+c1x, y3+c1y, z3+c1z, x3, y3, z3, x4, y4, z4, x4+c2x, y4+c2y, z4+c2z); "
3060 + "c1x = (d14y*d24z-d14z*d24y); c1y = (d14z*d24x-d14x*d24z); c1z = (d14x*d24y-d14y*d24x); "
3061 + "c2x = (d53y*d63z-d53z*d63y); c2y = (d53z*d63x-d53x*d63z); c2z = (d53x*d63y-d53y*d63x); "
3062 + "d14x = x1-x4; d14y = y1-y4; d14z = z1-z4; "
3063 + "d24x = x2-x4; d24y = y2-y4; d24z = z2-z4; "
3064 + "d53x = x5-x3; d53y = y5-y3; d53z = z5-z3; "
3065 + "d63x = x6-x3; d63y = y6-y3; d63z = z6-z3";
3066 piTorsionForce = OpenMM_CustomCompoundBondForce_create(6, energy);
3067 OpenMM_CustomCompoundBondForce_addPerBondParameter(piTorsionForce, "k");
3068 OpenMM_Force_setName(piTorsionForce, "PiTorsion");
3069
3070 PointerByReference particles = OpenMM_IntArray_create(0);
3071 PointerByReference parameters = OpenMM_DoubleArray_create(0);
3072 for (PiOrbitalTorsion piOrbitalTorsion : piOrbitalTorsions) {
3073 int a1 = piOrbitalTorsion.getAtom(0).getXyzIndex() - 1;
3074 int a2 = piOrbitalTorsion.getAtom(1).getXyzIndex() - 1;
3075 int a3 = piOrbitalTorsion.getAtom(2).getXyzIndex() - 1;
3076 int a4 = piOrbitalTorsion.getAtom(3).getXyzIndex() - 1;
3077 int a5 = piOrbitalTorsion.getAtom(4).getXyzIndex() - 1;
3078 int a6 = piOrbitalTorsion.getAtom(5).getXyzIndex() - 1;
3079 PiOrbitalTorsionType type = piOrbitalTorsion.piOrbitalTorsionType;
3080 double k =
3081 OpenMM_KJPerKcal * type.forceConstant * piOrbitalTorsion.piOrbitalTorsionType.piTorsUnit;
3082 OpenMM_IntArray_append(particles, a1);
3083 OpenMM_IntArray_append(particles, a2);
3084 OpenMM_IntArray_append(particles, a3);
3085 OpenMM_IntArray_append(particles, a4);
3086 OpenMM_IntArray_append(particles, a5);
3087 OpenMM_IntArray_append(particles, a6);
3088 OpenMM_DoubleArray_append(parameters, k);
3089 OpenMM_CustomCompoundBondForce_addBond(piTorsionForce, particles, parameters);
3090 OpenMM_IntArray_resize(particles, 0);
3091 OpenMM_DoubleArray_resize(parameters, 0);
3092 }
3093 OpenMM_IntArray_destroy(particles);
3094 OpenMM_DoubleArray_destroy(parameters);
3095
3096 int forceGroup = forceField.getInteger("PI_ORBITAL_TORSION_FORCE_GROUP", 0);
3097 OpenMM_Force_setForceGroup(piTorsionForce, forceGroup);
3098 OpenMM_System_addForce(system, piTorsionForce);
3099 logger.log(Level.INFO,
3100 format(" Pi-Orbital Torsions \t%6d\t\t%1d", piOrbitalTorsions.length, forceGroup));
3101 }
3102
3103
3104 private void updatePiTorsionForce() {
3105 PiOrbitalTorsion[] piOrbitalTorsions = getPiOrbitalTorsions();
3106 if (piOrbitalTorsions == null || piOrbitalTorsions.length < 1) {
3107 return;
3108 }
3109
3110 PointerByReference particles = OpenMM_IntArray_create(0);
3111 PointerByReference parameters = OpenMM_DoubleArray_create(0);
3112 int index = 0;
3113 for (PiOrbitalTorsion piOrbitalTorsion : piOrbitalTorsions) {
3114 int a1 = piOrbitalTorsion.getAtom(0).getXyzIndex() - 1;
3115 int a2 = piOrbitalTorsion.getAtom(1).getXyzIndex() - 1;
3116 int a3 = piOrbitalTorsion.getAtom(2).getXyzIndex() - 1;
3117 int a4 = piOrbitalTorsion.getAtom(3).getXyzIndex() - 1;
3118 int a5 = piOrbitalTorsion.getAtom(4).getXyzIndex() - 1;
3119 int a6 = piOrbitalTorsion.getAtom(5).getXyzIndex() - 1;
3120 PiOrbitalTorsionType type = piOrbitalTorsion.piOrbitalTorsionType;
3121 double k =
3122 OpenMM_KJPerKcal * type.forceConstant * piOrbitalTorsion.piOrbitalTorsionType.piTorsUnit;
3123 OpenMM_IntArray_append(particles, a1);
3124 OpenMM_IntArray_append(particles, a2);
3125 OpenMM_IntArray_append(particles, a3);
3126 OpenMM_IntArray_append(particles, a4);
3127 OpenMM_IntArray_append(particles, a5);
3128 OpenMM_IntArray_append(particles, a6);
3129 OpenMM_DoubleArray_append(parameters, k);
3130 OpenMM_CustomCompoundBondForce_setBondParameters(piTorsionForce, index++, particles,
3131 parameters);
3132 OpenMM_IntArray_resize(particles, 0);
3133 OpenMM_DoubleArray_resize(parameters, 0);
3134 }
3135 OpenMM_IntArray_destroy(particles);
3136 OpenMM_DoubleArray_destroy(parameters);
3137
3138 if (context.contextPointer != null) {
3139 OpenMM_CustomCompoundBondForce_updateParametersInContext(piTorsionForce,
3140 context.contextPointer);
3141 }
3142 }
3143
3144
3145 private void addTorsionTorsionForce() {
3146 TorsionTorsion[] torsionTorsions = getTorsionTorsions();
3147 if (torsionTorsions == null || torsionTorsions.length < 1) {
3148 return;
3149 }
3150
3151
3152 int nTypes = 0;
3153 LinkedHashMap<String, TorsionTorsionType> torTorTypes = new LinkedHashMap<>();
3154 PointerByReference amoebaTorsionTorsionForce = OpenMM_AmoebaTorsionTorsionForce_create();
3155 for (TorsionTorsion torsionTorsion : torsionTorsions) {
3156 int ia = torsionTorsion.getAtom(0).getXyzIndex() - 1;
3157 int ib = torsionTorsion.getAtom(1).getXyzIndex() - 1;
3158 int ic = torsionTorsion.getAtom(2).getXyzIndex() - 1;
3159 int id = torsionTorsion.getAtom(3).getXyzIndex() - 1;
3160 int ie = torsionTorsion.getAtom(4).getXyzIndex() - 1;
3161
3162 TorsionTorsionType torsionTorsionType = torsionTorsion.torsionTorsionType;
3163 String key = torsionTorsionType.getKey();
3164
3165
3166 int gridIndex = 0;
3167 if (torTorTypes.containsKey(key)) {
3168
3169
3170 int index = 0;
3171 for (String entry : torTorTypes.keySet()) {
3172 if (entry.equalsIgnoreCase(key)) {
3173 gridIndex = index;
3174 break;
3175 } else {
3176 index++;
3177 }
3178 }
3179 } else {
3180
3181 torTorTypes.put(key, torsionTorsionType);
3182 gridIndex = nTypes;
3183 nTypes++;
3184 }
3185
3186 Atom atom = torsionTorsion.getChiralAtom();
3187 int iChiral = -1;
3188 if (atom != null) {
3189 iChiral = atom.getXyzIndex() - 1;
3190 }
3191 OpenMM_AmoebaTorsionTorsionForce_addTorsionTorsion(amoebaTorsionTorsionForce, ia, ib, ic, id,
3192 ie, iChiral, gridIndex);
3193 }
3194
3195
3196 PointerByReference values = OpenMM_DoubleArray_create(6);
3197 int gridIndex = 0;
3198 for (String key : torTorTypes.keySet()) {
3199 TorsionTorsionType torTorType = torTorTypes.get(key);
3200 int nx = torTorType.nx;
3201 int ny = torTorType.ny;
3202 double[] tx = torTorType.tx;
3203 double[] ty = torTorType.ty;
3204 double[] f = torTorType.energy;
3205 double[] dx = torTorType.dx;
3206 double[] dy = torTorType.dy;
3207 double[] dxy = torTorType.dxy;
3208
3209
3210 PointerByReference grid3D = OpenMM_3D_DoubleArray_create(nx, ny, 6);
3211 int xIndex = 0;
3212 int yIndex = 0;
3213 for (int j = 0; j < nx * ny; j++) {
3214 int addIndex = 0;
3215 OpenMM_DoubleArray_set(values, addIndex++, tx[xIndex]);
3216 OpenMM_DoubleArray_set(values, addIndex++, ty[yIndex]);
3217 OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * f[j]);
3218 OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dx[j]);
3219 OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dy[j]);
3220 OpenMM_DoubleArray_set(values, addIndex, OpenMM_KJPerKcal * dxy[j]);
3221 OpenMM_3D_DoubleArray_set(grid3D, yIndex, xIndex, values);
3222 xIndex++;
3223 if (xIndex == nx) {
3224 xIndex = 0;
3225 yIndex++;
3226 }
3227 }
3228 OpenMM_AmoebaTorsionTorsionForce_setTorsionTorsionGrid(amoebaTorsionTorsionForce,
3229 gridIndex++, grid3D);
3230 OpenMM_3D_DoubleArray_destroy(grid3D);
3231 }
3232 OpenMM_DoubleArray_destroy(values);
3233
3234 int forceGroup = forceField.getInteger("TORSION_TORSION_FORCE_GROUP", 0);
3235
3236 OpenMM_Force_setForceGroup(amoebaTorsionTorsionForce, forceGroup);
3237 OpenMM_System_addForce(system, amoebaTorsionTorsionForce);
3238 logger.log(Level.INFO,
3239 format(" Torsion-Torsions \t%6d\t\t%1d", torsionTorsions.length, forceGroup));
3240 }
3241
3242
3243 private void addStretchTorsionForce() {
3244 StretchTorsion[] stretchTorsions = getStretchTorsions();
3245 if (stretchTorsions == null || stretchTorsions.length < 1) {
3246 return;
3247 }
3248
3249 PointerByReference stretchTorsionForce = OpenMM_CustomCompoundBondForce_create(4,
3250 StretchTorsion.stretchTorsionForm());
3251 OpenMM_CustomCompoundBondForce_addGlobalParameter(stretchTorsionForce, "phi1", 0);
3252 OpenMM_CustomCompoundBondForce_addGlobalParameter(stretchTorsionForce, "phi2", Math.PI);
3253 OpenMM_CustomCompoundBondForce_addGlobalParameter(stretchTorsionForce, "phi3", 0);
3254
3255 for (int m = 1; m < 4; m++) {
3256 for (int n = 1; n < 4; n++) {
3257 OpenMM_CustomCompoundBondForce_addPerBondParameter(stretchTorsionForce,
3258 String.format("k%d%d", m, n));
3259 }
3260 }
3261
3262 for (int m = 1; m < 4; m++) {
3263 OpenMM_CustomCompoundBondForce_addPerBondParameter(stretchTorsionForce,
3264 String.format("b%d", m));
3265 }
3266
3267 final double unitConv = OpenMM_KJPerKcal / OpenMM_NmPerAngstrom;
3268
3269 for (StretchTorsion strTors : stretchTorsions) {
3270 double[] constants = strTors.getConstants();
3271 PointerByReference strTorsParams = OpenMM_DoubleArray_create(0);
3272 for (int m = 0; m < 3; m++) {
3273 for (int n = 0; n < 3; n++) {
3274 int index = (3 * m) + n;
3275 double kmn = constants[index] * unitConv;
3276 OpenMM_DoubleArray_append(strTorsParams, kmn);
3277 }
3278 }
3279
3280 OpenMM_DoubleArray_append(strTorsParams, strTors.bondType1.distance * OpenMM_NmPerAngstrom);
3281 OpenMM_DoubleArray_append(strTorsParams, strTors.bondType2.distance * OpenMM_NmPerAngstrom);
3282 OpenMM_DoubleArray_append(strTorsParams, strTors.bondType3.distance * OpenMM_NmPerAngstrom);
3283
3284 PointerByReference strTorsParticles = OpenMM_IntArray_create(0);
3285 Atom[] atoms = strTors.getAtomArray(true);
3286 for (int i = 0; i < 4; i++) {
3287 OpenMM_IntArray_append(strTorsParticles, atoms[i].getXyzIndex() - 1);
3288 }
3289
3290 OpenMM_CustomCompoundBondForce_addBond(stretchTorsionForce, strTorsParticles, strTorsParams);
3291 OpenMM_DoubleArray_destroy(strTorsParams);
3292 OpenMM_IntArray_destroy(strTorsParticles);
3293 }
3294
3295 int forceGroup = forceField.getInteger("STRETCH_TORSION_FORCE_GROUP", 0);
3296
3297 OpenMM_Force_setForceGroup(stretchTorsionForce, forceGroup);
3298 OpenMM_System_addForce(system, stretchTorsionForce);
3299
3300 logger.log(Level.INFO,
3301 format(" Stretch-Torsions \t%6d\t\t%1d", stretchTorsions.length, forceGroup));
3302 }
3303
3304
3305 private void addAngleTorsionForce() {
3306 AngleTorsion[] angleTorsions = getAngleTorsions();
3307 if (angleTorsions == null || angleTorsions.length < 1) {
3308 return;
3309 }
3310
3311 PointerByReference angleTorsionForce = OpenMM_CustomCompoundBondForce_create(4,
3312 AngleTorsion.angleTorsionForm());
3313 OpenMM_CustomCompoundBondForce_addGlobalParameter(angleTorsionForce, "phi1", 0);
3314 OpenMM_CustomCompoundBondForce_addGlobalParameter(angleTorsionForce, "phi2", Math.PI);
3315 OpenMM_CustomCompoundBondForce_addGlobalParameter(angleTorsionForce, "phi3", 0);
3316
3317 for (int m = 1; m < 3; m++) {
3318 for (int n = 1; n < 4; n++) {
3319 OpenMM_CustomCompoundBondForce_addPerBondParameter(angleTorsionForce,
3320 format("k%d%d", m, n));
3321 }
3322 }
3323
3324 for (int m = 1; m < 3; m++) {
3325 OpenMM_CustomCompoundBondForce_addPerBondParameter(angleTorsionForce, format("a%d", m));
3326 }
3327
3328 for (AngleTorsion angleTorsion : angleTorsions) {
3329 double[] constants = angleTorsion.getConstants();
3330 PointerByReference atorsParams = OpenMM_DoubleArray_create(0);
3331 for (int m = 0; m < 2; m++) {
3332 for (int n = 0; n < 3; n++) {
3333 int index = (3 * m) + n;
3334 double kmn = constants[index] * OpenMM_KJPerKcal;
3335 OpenMM_DoubleArray_append(atorsParams, kmn);
3336 }
3337 }
3338
3339 Atom[] atoms = angleTorsion.getAtomArray(true);
3340
3341
3342
3343
3344
3345 OpenMM_DoubleArray_append(atorsParams,
3346 angleTorsion.angleType1.angle[0] * OpenMM_RadiansPerDegree);
3347 OpenMM_DoubleArray_append(atorsParams,
3348 angleTorsion.angleType2.angle[0] * OpenMM_RadiansPerDegree);
3349
3350 PointerByReference atorsParticles = OpenMM_IntArray_create(0);
3351 for (int i = 0; i < 4; i++) {
3352 OpenMM_IntArray_append(atorsParticles, atoms[i].getXyzIndex() - 1);
3353 }
3354
3355 OpenMM_CustomCompoundBondForce_addBond(angleTorsionForce, atorsParticles, atorsParams);
3356 OpenMM_DoubleArray_destroy(atorsParams);
3357 OpenMM_IntArray_destroy(atorsParticles);
3358 }
3359
3360 int forceGroup = forceField.getInteger("ANGLE_TORSION_FORCE_GROUP", 0);
3361
3362 OpenMM_Force_setForceGroup(angleTorsionForce, forceGroup);
3363 OpenMM_System_addForce(system, angleTorsionForce);
3364
3365 logger.log(Level.INFO,
3366 format(" Angle-Torsions \t%6d\t\t%1d", angleTorsions.length, forceGroup));
3367 }
3368
3369 private void addRestraintTorsions() {
3370 if (rTors != null && rTors.length > 0) {
3371 int nRT = rTors.length;
3372 restraintTorsions = new PointerByReference[nRT];
3373 for (int i = 0; i < nRT; i++) {
3374 PointerByReference rtOMM = OpenMM_PeriodicTorsionForce_create();
3375 RestraintTorsion rt = rTors[i];
3376 int a1 = rt.getAtom(0).getXyzIndex() - 1;
3377 int a2 = rt.getAtom(1).getXyzIndex() - 1;
3378 int a3 = rt.getAtom(2).getXyzIndex() - 1;
3379 int a4 = rt.getAtom(3).getXyzIndex() - 1;
3380 int nTerms = rt.torsionType.terms;
3381 for (int j = 0; j < nTerms; j++) {
3382 OpenMM_PeriodicTorsionForce_addTorsion(rtOMM, a1, a2, a3, a4, j + 1,
3383 rt.torsionType.phase[j] * OpenMM_RadiansPerDegree,
3384 OpenMM_KJPerKcal * rt.units * rt.torsionType.amplitude[j]);
3385 }
3386 int fGroup = forceField.getInteger("TORSION_FORCE_GROUP", 0);
3387
3388 OpenMM_Force_setForceGroup(rtOMM, fGroup);
3389 OpenMM_System_addForce(system, rtOMM);
3390 restraintTorsions[i] = rtOMM;
3391 }
3392 logger.info(format(" Added %d restraint torsions to OpenMM.", nRT));
3393 }
3394 }
3395
3396 private void updateRestraintTorsions() {
3397
3398
3399
3400
3401
3402 int nRT = restraintTorsions.length;
3403 for (int i = 0; i < nRT; i++) {
3404 RestraintTorsion rt = rTors[i];
3405 PointerByReference rtOMM = restraintTorsions[i];
3406 if (rt.applyLambda()) {
3407 int index = 0;
3408 TorsionType torsionType = rt.torsionType;
3409 int nTerms = torsionType.phase.length;
3410 int a1 = rt.getAtom(0).getXyzIndex() - 1;
3411 int a2 = rt.getAtom(1).getXyzIndex() - 1;
3412 int a3 = rt.getAtom(2).getXyzIndex() - 1;
3413 int a4 = rt.getAtom(3).getXyzIndex() - 1;
3414 for (int j = 0; j < nTerms; j++) {
3415 double forceConstant =
3416 OpenMM_KJPerKcal * rt.units * torsionType.amplitude[j] * rt.mapLambda(getLambda());
3417 OpenMM_PeriodicTorsionForce_setTorsionParameters(rtOMM, index++, a1, a2, a3, a4, j + 1,
3418 torsionType.phase[j] * OpenMM_RadiansPerDegree, forceConstant);
3419 }
3420 }
3421
3422 if (context.contextPointer != null) {
3423 OpenMM_PeriodicTorsionForce_updateParametersInContext(rtOMM, context.contextPointer);
3424 }
3425 }
3426 }
3427
3428
3429 private void addFixedChargeNonBondedForce() {
3430 VanDerWaals vdW = getVdwNode();
3431 if (vdW == null) {
3432 return;
3433 }
3434
3435
3436
3437
3438
3439 VanDerWaalsForm vdwForm = vdW.getVDWForm();
3440 if (vdwForm.vdwType != LENNARD_JONES || vdwForm.radiusRule != ARITHMETIC
3441 || vdwForm.epsilonRule != GEOMETRIC) {
3442 logger.info(format(" VDW Type: %s", vdwForm.vdwType));
3443 logger.info(format(" VDW Radius Rule: %s", vdwForm.radiusRule));
3444 logger.info(format(" VDW Epsilon Rule: %s", vdwForm.epsilonRule));
3445 logger.log(Level.SEVERE, " Unsupported van der Waals functional form.");
3446 return;
3447 }
3448
3449 fixedChargeNonBondedForce = OpenMM_NonbondedForce_create();
3450
3451
3452 double radScale = 1.0;
3453 if (vdwForm.radiusSize == RADIUS) {
3454 radScale = 2.0;
3455 }
3456
3457
3458 if (vdwForm.radiusType == R_MIN) {
3459 radScale /= 1.122462048309372981;
3460 }
3461
3462
3463 for (Atom atom : atoms) {
3464 VDWType vdwType = atom.getVDWType();
3465 double sigma = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
3466 double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
3467 double charge = 0.0;
3468 MultipoleType multipoleType = atom.getMultipoleType();
3469 if (multipoleType != null && atom.getElectrostatics()) {
3470 charge = multipoleType.charge;
3471 }
3472 OpenMM_NonbondedForce_addParticle(fixedChargeNonBondedForce, charge, sigma, eps);
3473 }
3474
3475
3476 double lj14Scale = vdwForm.getScale14();
3477 double coulomb14Scale = 1.0 / 1.2;
3478 ParticleMeshEwald pme = getPmeNode();
3479 if (pme != null) {
3480 coulomb14Scale = pme.getScale14();
3481 }
3482 Bond[] bonds = getBonds();
3483 if (bonds != null && bonds.length > 0) {
3484 PointerByReference bondArray = OpenMM_BondArray_create(0);
3485 for (Bond bond : bonds) {
3486 int i1 = bond.getAtom(0).getXyzIndex() - 1;
3487 int i2 = bond.getAtom(1).getXyzIndex() - 1;
3488 OpenMM_BondArray_append(bondArray, i1, i2);
3489 }
3490 OpenMM_NonbondedForce_createExceptionsFromBonds(fixedChargeNonBondedForce, bondArray,
3491 coulomb14Scale, lj14Scale);
3492 OpenMM_BondArray_destroy(bondArray);
3493
3494 int num = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
3495 chargeExclusion = new boolean[num];
3496 vdWExclusion = new boolean[num];
3497 exceptionChargeProd = new double[num];
3498 exceptionEps = new double[num];
3499
3500 IntByReference particle1 = new IntByReference();
3501 IntByReference particle2 = new IntByReference();
3502 DoubleByReference chargeProd = new DoubleByReference();
3503 DoubleByReference sigma = new DoubleByReference();
3504 DoubleByReference eps = new DoubleByReference();
3505
3506 for (int i = 0; i < num; i++) {
3507 OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, particle1,
3508 particle2, chargeProd, sigma, eps);
3509 if (abs(chargeProd.getValue()) > 0.0) {
3510 chargeExclusion[i] = false;
3511 exceptionChargeProd[i] = chargeProd.getValue();
3512 } else {
3513 exceptionChargeProd[i] = 0.0;
3514 chargeExclusion[i] = true;
3515 }
3516 if (abs(eps.getValue()) > 0.0) {
3517 vdWExclusion[i] = false;
3518 exceptionEps[i] = eps.getValue();
3519 } else {
3520 vdWExclusion[i] = true;
3521 exceptionEps[i] = 0.0;
3522 }
3523 }
3524 }
3525
3526 Crystal crystal = getCrystal();
3527 if (crystal.aperiodic()) {
3528 OpenMM_NonbondedForce_setNonbondedMethod(fixedChargeNonBondedForce,
3529 OpenMM_NonbondedForce_NonbondedMethod.OpenMM_NonbondedForce_NoCutoff);
3530 } else {
3531 OpenMM_NonbondedForce_setNonbondedMethod(fixedChargeNonBondedForce,
3532 OpenMM_NonbondedForce_NonbondedMethod.OpenMM_NonbondedForce_PME);
3533
3534 if (pme != null) {
3535
3536
3537 double aEwald = OpenMM_AngstromsPerNm * pme.getEwaldCoefficient();
3538 int nx = pme.getReciprocalSpace().getXDim();
3539 int ny = pme.getReciprocalSpace().getYDim();
3540 int nz = pme.getReciprocalSpace().getZDim();
3541 OpenMM_NonbondedForce_setPMEParameters(fixedChargeNonBondedForce, aEwald, nx, ny, nz);
3542 }
3543
3544 NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
3545 double off = nonbondedCutoff.off;
3546 double cut = nonbondedCutoff.cut;
3547 OpenMM_NonbondedForce_setCutoffDistance(fixedChargeNonBondedForce,
3548 OpenMM_NmPerAngstrom * off);
3549 OpenMM_NonbondedForce_setUseSwitchingFunction(fixedChargeNonBondedForce, OpenMM_True);
3550 if (cut == off) {
3551 logger.warning(" OpenMM does not properly handle cutoffs where cut == off!");
3552 if (cut == Double.MAX_VALUE || cut == Double.POSITIVE_INFINITY) {
3553 logger.info(" Detected infinite or max-value cutoff; setting cut to 1E+40 for OpenMM.");
3554 cut = 1E40;
3555 } else {
3556 logger.info(String.format(
3557 " Detected cut %8.4g == off %8.4g; scaling cut to 0.99 of off for OpenMM.", cut,
3558 off));
3559 cut *= 0.99;
3560 }
3561 }
3562 OpenMM_NonbondedForce_setSwitchingDistance(fixedChargeNonBondedForce,
3563 OpenMM_NmPerAngstrom * cut);
3564 }
3565
3566 OpenMM_NonbondedForce_setUseDispersionCorrection(fixedChargeNonBondedForce, OpenMM_False);
3567
3568 int forceGroup = forceField.getInteger("VDW_FORCE_GROUP", 1);
3569 int pmeGroup = forceField.getInteger("PME_FORCE_GROUP", 1);
3570 if (forceGroup != pmeGroup) {
3571 logger.severe(String.format(" ERROR: VDW-FORCE-GROUP is %d while PME-FORCE-GROUP is %d. "
3572 + "This is invalid for fixed-charge force fields with combined nonbonded forces.",
3573 forceGroup, pmeGroup));
3574 }
3575
3576 OpenMM_Force_setForceGroup(fixedChargeNonBondedForce, forceGroup);
3577 OpenMM_System_addForce(system, fixedChargeNonBondedForce);
3578
3579 logger.log(Level.INFO, format(" Fixed charge non-bonded force \t%1d", forceGroup));
3580
3581 GeneralizedKirkwood gk = getGK();
3582 if (gk != null) {
3583 addCustomGBForce();
3584 }
3585 }
3586
3587
3588
3589
3590
3591
3592 private void updateFixedChargeNonBondedForce(Atom[] atoms) {
3593 VanDerWaals vdW = getVdwNode();
3594
3595
3596 VanDerWaalsForm vdwForm = vdW.getVDWForm();
3597 if (vdwForm.vdwType != LENNARD_JONES || vdwForm.radiusRule != ARITHMETIC
3598 || vdwForm.epsilonRule != GEOMETRIC) {
3599 logger.log(Level.SEVERE, " Unsupported van der Waals functional form.");
3600 return;
3601 }
3602
3603
3604 double radScale = 1.0;
3605 if (vdwForm.radiusSize == RADIUS) {
3606 radScale = 2.0;
3607 }
3608
3609
3610 if (vdwForm.radiusType == R_MIN) {
3611 radScale /= 1.122462048309372981;
3612 }
3613
3614
3615 for (Atom atom : atoms) {
3616 int index = atom.getXyzIndex() - 1;
3617 boolean applyLambda = atom.applyLambda();
3618
3619 double charge = Double.MIN_VALUE;
3620 MultipoleType multipoleType = atom.getMultipoleType();
3621 if (multipoleType != null && atom.getElectrostatics()) {
3622 charge = multipoleType.charge;
3623 }
3624
3625 VDWType vdwType = atom.getVDWType();
3626 double sigma = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
3627 double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
3628
3629 if (applyLambda) {
3630
3631
3632 if (vdwLambdaTerm) {
3633 eps = 0.0;
3634 }
3635
3636 charge *= lambdaElec;
3637 }
3638
3639 if (!atom.getUse()) {
3640 eps = 0.0;
3641 charge = 0.0;
3642 }
3643
3644 OpenMM_NonbondedForce_setParticleParameters(fixedChargeNonBondedForce, index, charge, sigma,
3645 eps);
3646 }
3647
3648
3649 IntByReference particle1 = new IntByReference();
3650 IntByReference particle2 = new IntByReference();
3651 DoubleByReference chargeProd = new DoubleByReference();
3652 DoubleByReference sigma = new DoubleByReference();
3653 DoubleByReference eps = new DoubleByReference();
3654
3655 int numExceptions = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
3656 for (int i = 0; i < numExceptions; i++) {
3657
3658
3659 if (chargeExclusion[i] && vdWExclusion[i]) {
3660 continue;
3661 }
3662
3663 OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, particle1,
3664 particle2, chargeProd, sigma, eps);
3665
3666 int i1 = particle1.getValue();
3667 int i2 = particle2.getValue();
3668
3669 double qq = exceptionChargeProd[i];
3670 double epsilon = exceptionEps[i];
3671
3672 Atom atom1 = atoms[i1];
3673 Atom atom2 = atoms[i2];
3674
3675
3676
3677
3678
3679 double minEpsilon = 1.0e-12;
3680 double lambdaValue = lambdaElec;
3681
3682 if (lambdaValue < minEpsilon) {
3683 lambdaValue = minEpsilon;
3684 }
3685
3686 if (atom1.applyLambda()) {
3687 qq *= lambdaValue;
3688 if (vdwLambdaTerm) {
3689 epsilon = minEpsilon;
3690 }
3691 }
3692 if (atom2.applyLambda()) {
3693 qq *= lambdaValue;
3694 if (vdwLambdaTerm) {
3695 epsilon = minEpsilon;
3696 }
3697 }
3698 if (!atom1.getUse() || !atom2.getUse()) {
3699 qq = minEpsilon;
3700 epsilon = minEpsilon;
3701 }
3702 OpenMM_NonbondedForce_setExceptionParameters(fixedChargeNonBondedForce, i, i1, i2, qq,
3703 sigma.getValue(), epsilon);
3704 }
3705
3706 if (context.contextPointer != null) {
3707 OpenMM_NonbondedForce_updateParametersInContext(fixedChargeNonBondedForce,
3708 context.contextPointer);
3709 }
3710 }
3711
3712
3713
3714
3715
3716
3717
3718
3719 private void addCustomNonbondedSoftcoreForce() {
3720 VanDerWaals vdW = getVdwNode();
3721 if (vdW == null) {
3722 return;
3723 }
3724
3725
3726
3727
3728
3729 VanDerWaalsForm vdwForm = vdW.getVDWForm();
3730 if (vdwForm.vdwType != LENNARD_JONES || vdwForm.radiusRule != ARITHMETIC
3731 || vdwForm.epsilonRule != GEOMETRIC) {
3732 logger.info(format(" VDW Type: %s", vdwForm.vdwType));
3733 logger.info(format(" VDW Radius Rule: %s", vdwForm.radiusRule));
3734 logger.info(format(" VDW Epsilon Rule: %s", vdwForm.epsilonRule));
3735 logger.log(Level.SEVERE, " Unsupported van der Waals functional form.");
3736 return;
3737 }
3738
3739
3740 String stericsMixingRules = " epsilon = sqrt(epsilon1*epsilon2);";
3741 stericsMixingRules += " rmin = 0.5 * (sigma1 + sigma2) * 1.122462048309372981;";
3742
3743
3744 String stericsEnergyExpression = "(vdw_lambda^beta)*epsilon*x*(x-2.0);";
3745
3746 stericsEnergyExpression += " x = 1.0 / (alpha*(1.0-vdw_lambda)^2.0 + (r/rmin)^6.0);";
3747
3748 String energyExpression = stericsEnergyExpression + stericsMixingRules;
3749
3750 PointerByReference fixedChargeSoftcore = OpenMM_CustomNonbondedForce_create(energyExpression);
3751
3752
3753 double alpha = vdWSoftcoreAlpha;
3754 double beta = vdwSoftcorePower;
3755
3756 OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "vdw_lambda", 1.0);
3757 OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "alpha", alpha);
3758 OpenMM_CustomNonbondedForce_addGlobalParameter(fixedChargeSoftcore, "beta", beta);
3759 OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "sigma");
3760 OpenMM_CustomNonbondedForce_addPerParticleParameter(fixedChargeSoftcore, "epsilon");
3761
3762
3763 PointerByReference alchemicalGroup = OpenMM_IntSet_create();
3764 PointerByReference nonAlchemicalGroup = OpenMM_IntSet_create();
3765 DoubleByReference charge = new DoubleByReference();
3766 DoubleByReference sigma = new DoubleByReference();
3767 DoubleByReference eps = new DoubleByReference();
3768
3769 int index = 0;
3770 for (Atom atom : atoms) {
3771 if (atom.applyLambda()) {
3772 OpenMM_IntSet_insert(alchemicalGroup, index);
3773 } else {
3774 OpenMM_IntSet_insert(nonAlchemicalGroup, index);
3775 }
3776
3777 OpenMM_NonbondedForce_getParticleParameters(fixedChargeNonBondedForce, index, charge, sigma,
3778 eps);
3779 double sigmaValue = sigma.getValue();
3780 double epsValue = eps.getValue();
3781
3782
3783 if (sigmaValue == 0.0) {
3784 sigmaValue = 1.0;
3785 epsValue = 0.0;
3786 }
3787
3788 PointerByReference particleParameters = OpenMM_DoubleArray_create(0);
3789 OpenMM_DoubleArray_append(particleParameters, sigmaValue);
3790 OpenMM_DoubleArray_append(particleParameters, epsValue);
3791 OpenMM_CustomNonbondedForce_addParticle(fixedChargeSoftcore, particleParameters);
3792 OpenMM_DoubleArray_destroy(particleParameters);
3793
3794 index++;
3795 }
3796
3797 OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup,
3798 alchemicalGroup);
3799 OpenMM_CustomNonbondedForce_addInteractionGroup(fixedChargeSoftcore, alchemicalGroup,
3800 nonAlchemicalGroup);
3801 OpenMM_IntSet_destroy(alchemicalGroup);
3802 OpenMM_IntSet_destroy(nonAlchemicalGroup);
3803
3804 Crystal crystal = getCrystal();
3805 if (crystal.aperiodic()) {
3806 OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore,
3807 OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_NoCutoff);
3808 } else {
3809 OpenMM_CustomNonbondedForce_setNonbondedMethod(fixedChargeSoftcore,
3810 OpenMM_CustomNonbondedForce_NonbondedMethod.OpenMM_CustomNonbondedForce_CutoffPeriodic);
3811 }
3812
3813 NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
3814 double off = nonbondedCutoff.off;
3815 double cut = nonbondedCutoff.cut;
3816 if (cut == off) {
3817 logger.warning(" OpenMM does not properly handle cutoffs where cut == off!");
3818 if (cut == Double.MAX_VALUE || cut == Double.POSITIVE_INFINITY) {
3819 logger.info(" Detected infinite or max-value cutoff; setting cut to 1E+40 for OpenMM.");
3820 cut = 1E40;
3821 } else {
3822 logger.info(
3823 format(" Detected cut %8.4g == off %8.4g; scaling cut to 0.99 of off for OpenMM.", cut,
3824 off));
3825 cut *= 0.99;
3826 }
3827 }
3828
3829 OpenMM_CustomNonbondedForce_setCutoffDistance(fixedChargeSoftcore, OpenMM_NmPerAngstrom * off);
3830 OpenMM_CustomNonbondedForce_setUseSwitchingFunction(fixedChargeSoftcore, OpenMM_True);
3831 OpenMM_CustomNonbondedForce_setSwitchingDistance(fixedChargeSoftcore,
3832 OpenMM_NmPerAngstrom * cut);
3833
3834
3835
3836
3837
3838 int forceGroup = forceField.getInteger("VDW_FORCE_GROUP", 1);
3839
3840 OpenMM_Force_setForceGroup(fixedChargeSoftcore, forceGroup);
3841 OpenMM_System_addForce(system, fixedChargeSoftcore);
3842
3843
3844 PointerByReference alchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(
3845 stericsEnergyExpression);
3846
3847
3848 PointerByReference nonAlchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(
3849 stericsEnergyExpression);
3850
3851
3852 OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "rmin");
3853 OpenMM_CustomBondForce_addPerBondParameter(alchemicalAlchemicalStericsForce, "epsilon");
3854 OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "vdw_lambda", 1.0);
3855 OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "alpha", alpha);
3856 OpenMM_CustomBondForce_addGlobalParameter(alchemicalAlchemicalStericsForce, "beta", beta);
3857
3858 OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "rmin");
3859 OpenMM_CustomBondForce_addPerBondParameter(nonAlchemicalAlchemicalStericsForce, "epsilon");
3860 OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "vdw_lambda",
3861 1.0);
3862 OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "alpha", alpha);
3863 OpenMM_CustomBondForce_addGlobalParameter(nonAlchemicalAlchemicalStericsForce, "beta", beta);
3864
3865 int range = OpenMM_NonbondedForce_getNumExceptions(fixedChargeNonBondedForce);
3866
3867 IntByReference atomi = new IntByReference();
3868 IntByReference atomj = new IntByReference();
3869 int[][] torsionMask = vdW.getMask14();
3870
3871 for (int i = 0; i < range; i++) {
3872 OpenMM_NonbondedForce_getExceptionParameters(fixedChargeNonBondedForce, i, atomi, atomj,
3873 charge, sigma, eps);
3874
3875
3876
3877 OpenMM_CustomNonbondedForce_addExclusion(fixedChargeSoftcore, atomi.getValue(),
3878 atomj.getValue());
3879
3880
3881 int[] maskI = torsionMask[atomi.getValue()];
3882 int jID = atomj.getValue();
3883 boolean epsException = false;
3884
3885 for (int mask : maskI) {
3886 if (mask == jID) {
3887 epsException = true;
3888 break;
3889 }
3890 }
3891
3892 if (epsException) {
3893 Atom atom1 = atoms[atomi.getValue()];
3894 Atom atom2 = atoms[atomj.getValue()];
3895
3896 boolean bothAlchemical = false;
3897 boolean oneAlchemical = false;
3898
3899 if (atom1.applyLambda() && atom2.applyLambda()) {
3900 bothAlchemical = true;
3901 } else if ((atom1.applyLambda() && !atom2.applyLambda()) || (!atom1.applyLambda()
3902 && atom2.applyLambda())) {
3903 oneAlchemical = true;
3904 }
3905
3906 if (bothAlchemical) {
3907 PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
3908 OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
3909 OpenMM_DoubleArray_append(bondParameters, eps.getValue());
3910 OpenMM_CustomBondForce_addBond(alchemicalAlchemicalStericsForce, atomi.getValue(),
3911 atomj.getValue(), bondParameters);
3912 OpenMM_DoubleArray_destroy(bondParameters);
3913 } else if (oneAlchemical) {
3914 PointerByReference bondParameters = OpenMM_DoubleArray_create(0);
3915 OpenMM_DoubleArray_append(bondParameters, sigma.getValue() * 1.122462048309372981);
3916 OpenMM_DoubleArray_append(bondParameters, eps.getValue());
3917 OpenMM_CustomBondForce_addBond(nonAlchemicalAlchemicalStericsForce, atomi.getValue(),
3918 atomj.getValue(), bondParameters);
3919 OpenMM_DoubleArray_destroy(bondParameters);
3920 }
3921 }
3922 }
3923
3924 OpenMM_Force_setForceGroup(alchemicalAlchemicalStericsForce, forceGroup);
3925 OpenMM_Force_setForceGroup(nonAlchemicalAlchemicalStericsForce, forceGroup);
3926
3927
3928
3929
3930
3931
3932 OpenMM_System_addForce(system, alchemicalAlchemicalStericsForce);
3933 OpenMM_System_addForce(system, nonAlchemicalAlchemicalStericsForce);
3934
3935 logger.log(Level.INFO, format(" Added fixed charge softcore force \t%d", forceGroup));
3936 logger.log(Level.INFO, format(" Alpha = %8.6f and beta = %8.6f", alpha, beta));
3937 }
3938
3939
3940 private void addCustomGBForce() {
3941 GeneralizedKirkwood gk = getGK();
3942 if (gk == null) {
3943 return;
3944 }
3945
3946 double sTens = 0.0;
3947 if (gk.getNonPolarModel() == NonPolarModel.BORN_SOLV
3948 || gk.getNonPolarModel() == NonPolarModel.BORN_CAV_DISP) {
3949 sTens = gk.getSurfaceTension();
3950 sTens *= OpenMM_KJPerKcal;
3951 sTens *= 100.0;
3952
3953
3954 }
3955
3956 customGBForce = OpenMM_CustomGBForce_create();
3957 OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "q");
3958 OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "radius");
3959 OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "scale");
3960 OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "surfaceTension");
3961
3962 OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "solventDielectric",
3963 gk.getSolventPermittivity());
3964 OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "soluteDielectric", 1.0);
3965 OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "dOffset",
3966 gk.getDielecOffset() * OpenMM_NmPerAngstrom);
3967 OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "probeRadius",
3968 gk.getProbeRadius() * OpenMM_NmPerAngstrom);
3969
3970 OpenMM_CustomGBForce_addComputedValue(customGBForce, "I",
3971
3972
3973 "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);"
3974 + "U=r+sr2;"
3975
3976 + "C=2/3*(1/or1^3-1/L^3)*step(sr2-r-or1);"
3977
3978
3979 + "L = step(sr2 - r1r)*sr2mr + (1 - step(sr2 - r1r))*L;" + "sr2mr = sr2 - r;"
3980 + "r1r = radius1 + r;" + "L = step(r1sr2 - r)*radius1 + (1 - step(r1sr2 - r))*L;"
3981 + "r1sr2 = radius1 + sr2;" + "L = r - sr2;" + "sr2 = scale2 * radius2;"
3982 + "or1 = radius1; or2 = radius2", OpenMM_CustomGBForce_ParticlePairNoExclusions);
3983
3984 OpenMM_CustomGBForce_addComputedValue(customGBForce, "B",
3985
3986
3987 "step(BB-radius)*BB + (1 - step(BB-radius))*radius;" + "BB = 1 / ( (3.0*III)^(1.0/3.0) );"
3988 + "III = step(II)*II + (1 - step(II))*1.0e-9/3.0;" + "II = maxI - I;"
3989 + "maxI = 1/(3.0*radius^3)", OpenMM_CustomGBForce_SingleParticle);
3990
3991 OpenMM_CustomGBForce_addEnergyTerm(customGBForce,
3992 "surfaceTension*(radius+probeRadius+dOffset)^2*((radius+dOffset)/B)^6/6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B",
3993 OpenMM_CustomGBForce_SingleParticle);
3994
3995
3996 OpenMM_CustomGBForce_addEnergyTerm(customGBForce,
3997 "-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
3998 + "f=sqrt(r^2+B1*B2*exp(-r^2/(2.455*B1*B2)))", OpenMM_CustomGBForce_ParticlePair);
3999
4000 double[] baseRadii = gk.getBaseRadii();
4001 double[] overlapScale = gk.getOverlapScale();
4002 PointerByReference doubleArray = OpenMM_DoubleArray_create(0);
4003 for (int i = 0; i < nAtoms; i++) {
4004 MultipoleType multipoleType = atoms[i].getMultipoleType();
4005 OpenMM_DoubleArray_append(doubleArray, multipoleType.charge);
4006 OpenMM_DoubleArray_append(doubleArray, OpenMM_NmPerAngstrom * baseRadii[i]);
4007 OpenMM_DoubleArray_append(doubleArray, overlapScale[i]);
4008 OpenMM_DoubleArray_append(doubleArray, sTens);
4009
4010 OpenMM_CustomGBForce_addParticle(customGBForce, doubleArray);
4011 OpenMM_DoubleArray_resize(doubleArray, 0);
4012 }
4013 OpenMM_DoubleArray_destroy(doubleArray);
4014
4015 double cut = gk.getCutoff();
4016 OpenMM_CustomGBForce_setCutoffDistance(customGBForce, cut);
4017
4018 int forceGroup = forceField.getInteger("GK_FORCE_GROUP", 1);
4019 OpenMM_Force_setForceGroup(customGBForce, forceGroup);
4020 OpenMM_System_addForce(system, customGBForce);
4021
4022 logger.log(Level.INFO, format(" Custom generalized Born force \t%d", forceGroup));
4023 }
4024
4025
4026
4027
4028
4029
4030 private void updateCustomGBForce(Atom[] atoms) {
4031 GeneralizedKirkwood gk = getGK();
4032 double[] baseRadii = gk.getBaseRadii();
4033 double[] overlapScale = gk.getOverlapScale();
4034 boolean nea = gk.getNativeEnvironmentApproximation();
4035
4036 double sTens = 0.0;
4037 if (gk.getNonPolarModel() == NonPolarModel.BORN_SOLV
4038 || gk.getNonPolarModel() == NonPolarModel.BORN_CAV_DISP) {
4039 sTens = gk.getSurfaceTension();
4040 sTens *= OpenMM_KJPerKcal;
4041 sTens *= 100.0;
4042 }
4043
4044 PointerByReference doubleArray = OpenMM_DoubleArray_create(0);
4045 for (Atom atom : atoms) {
4046 int index = atom.getXyzIndex() - 1;
4047 double chargeUseFactor = 1.0;
4048 if (!atom.getUse() || !atom.getElectrostatics()) {
4049 chargeUseFactor = 0.0;
4050 }
4051 double lambdaScale = lambdaElec;
4052 if (!atom.applyLambda()) {
4053 lambdaScale = 1.0;
4054 }
4055
4056 chargeUseFactor *= lambdaScale;
4057 MultipoleType multipoleType = atom.getMultipoleType();
4058 double charge = multipoleType.charge * chargeUseFactor;
4059 double surfaceTension = sTens * chargeUseFactor;
4060
4061 double overlapScaleUseFactor = nea ? 1.0 : chargeUseFactor;
4062 double oScale = overlapScale[index] * overlapScaleUseFactor;
4063 double baseRadius = baseRadii[index];
4064
4065 OpenMM_DoubleArray_append(doubleArray, charge);
4066 OpenMM_DoubleArray_append(doubleArray, OpenMM_NmPerAngstrom * baseRadius);
4067 OpenMM_DoubleArray_append(doubleArray, oScale);
4068 OpenMM_DoubleArray_append(doubleArray, surfaceTension);
4069 OpenMM_CustomGBForce_setParticleParameters(customGBForce, index, doubleArray);
4070
4071
4072 OpenMM_DoubleArray_resize(doubleArray, 0);
4073 }
4074 OpenMM_DoubleArray_destroy(doubleArray);
4075
4076 if (context.contextPointer != null) {
4077 OpenMM_CustomGBForce_updateParametersInContext(customGBForce, context.contextPointer);
4078 }
4079 }
4080
4081
4082 private void addAmoebaVDWForce() {
4083 VanDerWaals vdW = getVdwNode();
4084 if (vdW == null) {
4085 return;
4086 }
4087
4088 amoebaVDWForce = OpenMM_AmoebaVdwForce_create();
4089
4090 VanDerWaalsForm vdwForm = vdW.getVDWForm();
4091 NonbondedCutoff nonbondedCutoff = vdW.getNonbondedCutoff();
4092 Crystal crystal = getCrystal();
4093
4094 double radScale = 1.0;
4095 if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
4096 radScale = 0.5;
4097 }
4098
4099
4100
4101
4102
4103
4104 vdWClassForNoInteraction = 0;
4105
4106 vdwClassToOpenMMType = new HashMap<>();
4107
4108 Map<String, VDWType> vdwTypes = forceField.getVDWTypes();
4109
4110 for (VDWType vdwType : vdwTypes.values()) {
4111 int atomClass = vdwType.atomClass;
4112 if (!vdwClassToOpenMMType.containsKey(atomClass)) {
4113 double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
4114 double rad = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
4115 int type = OpenMM_AmoebaVdwForce_addParticleType(amoebaVDWForce, rad, eps);
4116 vdwClassToOpenMMType.put(atomClass, type);
4117 if (atomClass <= vdWClassForNoInteraction) {
4118 vdWClassForNoInteraction = atomClass - 1;
4119 }
4120 }
4121 }
4122
4123
4124 int type = OpenMM_AmoebaVdwForce_addParticleType(amoebaVDWForce, OpenMM_NmPerAngstrom, 0.0);
4125 vdwClassToOpenMMType.put(vdWClassForNoInteraction, type);
4126
4127 Map<String, VDWPairType> vdwPairTypeMap = forceField.getVDWPairTypes();
4128 for (VDWPairType vdwPairType : vdwPairTypeMap.values()) {
4129 int c1 = vdwPairType.atomClasses[0];
4130 int c2 = vdwPairType.atomClasses[1];
4131 int type1 = vdwClassToOpenMMType.get(c1);
4132 int type2 = vdwClassToOpenMMType.get(c2);
4133 double rMin = vdwPairType.radius * OpenMM_NmPerAngstrom;
4134 double eps = vdwPairType.wellDepth * OpenMM_KJPerKcal;
4135 OpenMM_AmoebaVdwForce_addTypePair(amoebaVDWForce, type1, type2, rMin, eps);
4136 OpenMM_AmoebaVdwForce_addTypePair(amoebaVDWForce, type2, type1, rMin, eps);
4137 }
4138
4139 ExtendedSystem extendedSystem = vdW.getExtendedSystem();
4140 double[] vdwPrefactorAndDerivs = new double[3];
4141
4142 int[] ired = vdW.getReductionIndex();
4143 for (int i = 0; i < nAtoms; i++) {
4144 Atom atom = atoms[i];
4145 VDWType vdwType = atom.getVDWType();
4146 int atomClass = vdwType.atomClass;
4147 type = vdwClassToOpenMMType.get(atomClass);
4148 int isAlchemical = atom.applyLambda() ? 1 : 0;
4149 double scaleFactor = 1.0;
4150 if (extendedSystem != null) {
4151 extendedSystem.getVdwPrefactor(i, vdwPrefactorAndDerivs);
4152 scaleFactor = vdwPrefactorAndDerivs[0];
4153 }
4154
4155 OpenMM_AmoebaVdwForce_addParticle_1(amoebaVDWForce, ired[i], type, vdwType.reductionFactor,
4156 isAlchemical, scaleFactor);
4157 }
4158
4159 double cutoff = nonbondedCutoff.off * OpenMM_NmPerAngstrom;
4160 OpenMM_AmoebaVdwForce_setCutoffDistance(amoebaVDWForce, cutoff);
4161
4162 if (vdW.getDoLongRangeCorrection()) {
4163 OpenMM_AmoebaVdwForce_setUseDispersionCorrection(amoebaVDWForce, OpenMM_True);
4164 } else {
4165 OpenMM_AmoebaVdwForce_setUseDispersionCorrection(amoebaVDWForce, OpenMM_False);
4166 }
4167
4168 if (crystal.aperiodic()) {
4169 OpenMM_AmoebaVdwForce_setNonbondedMethod(amoebaVDWForce, OpenMM_AmoebaVdwForce_NoCutoff);
4170 } else {
4171 OpenMM_AmoebaVdwForce_setNonbondedMethod(amoebaVDWForce,
4172 OpenMM_AmoebaVdwForce_CutoffPeriodic);
4173 }
4174
4175 if (vdwLambdaTerm) {
4176 OpenMM_AmoebaVdwForce_setAlchemicalMethod(amoebaVDWForce, OpenMM_AmoebaVdwForce_Decouple);
4177 OpenMM_AmoebaVdwForce_setSoftcoreAlpha(amoebaVDWForce, vdWSoftcoreAlpha);
4178 OpenMM_AmoebaVdwForce_setSoftcorePower(amoebaVDWForce, (int) vdwSoftcorePower);
4179 }
4180
4181 int[][] bondMask = vdW.getMask12();
4182 int[][] angleMask = vdW.getMask13();
4183
4184
4185 PointerByReference exclusions = OpenMM_IntArray_create(0);
4186 for (int i = 0; i < nAtoms; i++) {
4187 OpenMM_IntArray_append(exclusions, i);
4188 final int[] bondMaski = bondMask[i];
4189 for (int value : bondMaski) {
4190 OpenMM_IntArray_append(exclusions, value);
4191 }
4192 final int[] angleMaski = angleMask[i];
4193 for (int value : angleMaski) {
4194 OpenMM_IntArray_append(exclusions, value);
4195 }
4196 OpenMM_AmoebaVdwForce_setParticleExclusions(amoebaVDWForce, i, exclusions);
4197 OpenMM_IntArray_resize(exclusions, 0);
4198 }
4199 OpenMM_IntArray_destroy(exclusions);
4200
4201 int forceGroup = forceField.getInteger("VDW_FORCE_GROUP", 1);
4202 OpenMM_Force_setForceGroup(amoebaVDWForce, forceGroup);
4203 OpenMM_System_addForce(system, amoebaVDWForce);
4204
4205 logger.log(Level.INFO, format(" AMOEBA van der Waals force \t\t%d", forceGroup));
4206 }
4207
4208
4209
4210
4211
4212
4213 private void updateAmoebaVDWForce(Atom[] atoms) {
4214 VanDerWaals vdW = getVdwNode();
4215 VanDerWaalsForm vdwForm = vdW.getVDWForm();
4216 double radScale = 1.0;
4217 if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
4218 radScale = 0.5;
4219 }
4220
4221 ExtendedSystem extendedSystem = vdW.getExtendedSystem();
4222 double[] vdwPrefactorAndDerivs = new double[3];
4223
4224 int[] ired = vdW.getReductionIndex();
4225 for (Atom atom : atoms) {
4226 int index = atom.getXyzIndex() - 1;
4227 VDWType vdwType = atom.getVDWType();
4228
4229
4230 int type = vdwClassToOpenMMType.get(vdwType.atomClass);
4231 if (!atom.getUse()) {
4232
4233 type = vdwClassToOpenMMType.get(vdWClassForNoInteraction);
4234 }
4235 int isAlchemical = atom.applyLambda() ? 1 : 0;
4236 double eps = OpenMM_KJPerKcal * vdwType.wellDepth;
4237 double rad = OpenMM_NmPerAngstrom * vdwType.radius * radScale;
4238
4239 double scaleFactor = 1.0;
4240 if (extendedSystem != null) {
4241 extendedSystem.getVdwPrefactor(index, vdwPrefactorAndDerivs);
4242 scaleFactor = vdwPrefactorAndDerivs[0];
4243 }
4244
4245 OpenMM_AmoebaVdwForce_setParticleParameters(amoebaVDWForce, index, ired[index], rad, eps,
4246 vdwType.reductionFactor, isAlchemical, type, scaleFactor);
4247 }
4248
4249 if (context.contextPointer != null) {
4250 OpenMM_AmoebaVdwForce_updateParametersInContext(amoebaVDWForce, context.contextPointer);
4251 }
4252 }
4253
4254
4255 private void addAmoebaMultipoleForce() {
4256 ParticleMeshEwald pme = getPmeNode();
4257 if (pme == null) {
4258 return;
4259 }
4260
4261 int[][] axisAtom = pme.getAxisAtoms();
4262 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
4263 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
4264 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
4265
4266 amoebaMultipoleForce = OpenMM_AmoebaMultipoleForce_create();
4267
4268 double polarScale = 1.0;
4269 SCFAlgorithm scfAlgorithm = null;
4270
4271 if (pme.getPolarizationType() != Polarization.MUTUAL) {
4272 OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce,
4273 OpenMM_AmoebaMultipoleForce_Direct);
4274 if (pme.getPolarizationType() == Polarization.NONE) {
4275 polarScale = 0.0;
4276 }
4277 } else {
4278 String algorithm = forceField.getString("SCF_ALGORITHM", "CG");
4279 try {
4280 algorithm = algorithm.replaceAll("-", "_").toUpperCase();
4281 scfAlgorithm = SCFAlgorithm.valueOf(algorithm);
4282 } catch (Exception e) {
4283 scfAlgorithm = SCFAlgorithm.CG;
4284 }
4285
4286 if (scfAlgorithm == SCFAlgorithm.EPT) {
4287
4288
4289
4290
4291 OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce,
4292 OpenMM_AmoebaMultipoleForce_Extrapolated);
4293 PointerByReference exptCoefficients = OpenMM_DoubleArray_create(4);
4294 OpenMM_DoubleArray_set(exptCoefficients, 0, -0.154);
4295 OpenMM_DoubleArray_set(exptCoefficients, 1, 0.017);
4296 OpenMM_DoubleArray_set(exptCoefficients, 2, 0.657);
4297 OpenMM_DoubleArray_set(exptCoefficients, 3, 0.475);
4298
4299
4300 OpenMM_AmoebaMultipoleForce_setExtrapolationCoefficients(amoebaMultipoleForce,
4301 exptCoefficients);
4302 OpenMM_DoubleArray_destroy(exptCoefficients);
4303 } else {
4304 OpenMM_AmoebaMultipoleForce_setPolarizationType(amoebaMultipoleForce,
4305 OpenMM_AmoebaMultipoleForce_Mutual);
4306 }
4307 }
4308
4309 PointerByReference dipoles = OpenMM_DoubleArray_create(3);
4310 PointerByReference quadrupoles = OpenMM_DoubleArray_create(9);
4311
4312 for (int i = 0; i < nAtoms; i++) {
4313 Atom atom = atoms[i];
4314 MultipoleType multipoleType = pme.getMultipoleType(i);
4315 PolarizeType polarType = pme.getPolarizeType(i);
4316
4317
4318 int axisType = switch (multipoleType.frameDefinition) {
4319 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
4320 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
4321 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
4322 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
4323 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
4324 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
4325 };
4326
4327 double useFactor = 1.0;
4328 if (!atoms[i].getUse() || !atoms[i].getElectrostatics()) {
4329 useFactor = 0.0;
4330 }
4331
4332 double lambdaScale = lambdaElec;
4333 if (!atom.applyLambda()) {
4334 lambdaScale = 1.0;
4335 }
4336
4337 useFactor *= lambdaScale;
4338
4339
4340 for (int j = 0; j < 3; j++) {
4341 OpenMM_DoubleArray_set(dipoles, j,
4342 multipoleType.dipole[j] * OpenMM_NmPerAngstrom * useFactor);
4343 }
4344 int l = 0;
4345 for (int j = 0; j < 3; j++) {
4346 for (int k = 0; k < 3; k++) {
4347 OpenMM_DoubleArray_set(quadrupoles, l++,
4348 multipoleType.quadrupole[j][k] * quadrupoleConversion * useFactor / 3.0);
4349 }
4350 }
4351
4352 int zaxis = -1;
4353 int xaxis = -1;
4354 int yaxis = -1;
4355 int[] refAtoms = axisAtom[i];
4356 if (refAtoms != null) {
4357 zaxis = refAtoms[0];
4358 if (refAtoms.length > 1) {
4359 xaxis = refAtoms[1];
4360 if (refAtoms.length > 2) {
4361 yaxis = refAtoms[2];
4362 }
4363 }
4364 } else {
4365 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
4366 }
4367
4368 double charge = multipoleType.charge * useFactor;
4369
4370
4371 OpenMM_AmoebaMultipoleForce_addMultipole(amoebaMultipoleForce, charge, dipoles, quadrupoles,
4372 axisType, zaxis, xaxis, yaxis, polarType.thole,
4373 polarType.pdamp * dampingFactorConversion,
4374 polarType.polarizability * polarityConversion * polarScale);
4375 }
4376 OpenMM_DoubleArray_destroy(dipoles);
4377 OpenMM_DoubleArray_destroy(quadrupoles);
4378
4379 Crystal crystal = getCrystal();
4380 if (!crystal.aperiodic()) {
4381 OpenMM_AmoebaMultipoleForce_setNonbondedMethod(amoebaMultipoleForce,
4382 OpenMM_AmoebaMultipoleForce_PME);
4383 OpenMM_AmoebaMultipoleForce_setCutoffDistance(amoebaMultipoleForce,
4384 pme.getEwaldCutoff() * OpenMM_NmPerAngstrom);
4385 OpenMM_AmoebaMultipoleForce_setAEwald(amoebaMultipoleForce,
4386 pme.getEwaldCoefficient() / OpenMM_NmPerAngstrom);
4387
4388 double ewaldTolerance = 1.0e-04;
4389 OpenMM_AmoebaMultipoleForce_setEwaldErrorTolerance(amoebaMultipoleForce, ewaldTolerance);
4390
4391 PointerByReference gridDimensions = OpenMM_IntArray_create(3);
4392 ReciprocalSpace recip = pme.getReciprocalSpace();
4393 OpenMM_IntArray_set(gridDimensions, 0, recip.getXDim());
4394 OpenMM_IntArray_set(gridDimensions, 1, recip.getYDim());
4395 OpenMM_IntArray_set(gridDimensions, 2, recip.getZDim());
4396 OpenMM_AmoebaMultipoleForce_setPmeGridDimensions(amoebaMultipoleForce, gridDimensions);
4397 OpenMM_IntArray_destroy(gridDimensions);
4398 } else {
4399 OpenMM_AmoebaMultipoleForce_setNonbondedMethod(amoebaMultipoleForce,
4400 OpenMM_AmoebaMultipoleForce_NoCutoff);
4401 }
4402
4403 OpenMM_AmoebaMultipoleForce_setMutualInducedMaxIterations(amoebaMultipoleForce, 500);
4404 OpenMM_AmoebaMultipoleForce_setMutualInducedTargetEpsilon(amoebaMultipoleForce,
4405 pme.getPolarEps());
4406
4407 int[][] ip11 = pme.getPolarization11();
4408
4409 PointerByReference covalentMap = OpenMM_IntArray_create(0);
4410 for (int i = 0; i < nAtoms; i++) {
4411 Atom ai = atoms[i];
4412
4413
4414 OpenMM_IntArray_resize(covalentMap, 0);
4415 for (Atom ak : ai.get12List()) {
4416 OpenMM_IntArray_append(covalentMap, ak.getIndex() - 1);
4417 }
4418 OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
4419 OpenMM_AmoebaMultipoleForce_Covalent12, covalentMap);
4420
4421
4422 OpenMM_IntArray_resize(covalentMap, 0);
4423 for (Atom ak : ai.get13List()) {
4424 OpenMM_IntArray_append(covalentMap, ak.getIndex() - 1);
4425 }
4426 OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
4427 OpenMM_AmoebaMultipoleForce_Covalent13, covalentMap);
4428
4429
4430 OpenMM_IntArray_resize(covalentMap, 0);
4431 for (Atom ak : ai.get14List()) {
4432 OpenMM_IntArray_append(covalentMap, ak.getIndex() - 1);
4433 }
4434 OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
4435 OpenMM_AmoebaMultipoleForce_Covalent14, covalentMap);
4436
4437
4438 OpenMM_IntArray_resize(covalentMap, 0);
4439 for (Atom ak : ai.get15List()) {
4440 OpenMM_IntArray_append(covalentMap, ak.getIndex() - 1);
4441 }
4442 OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
4443 OpenMM_AmoebaMultipoleForce_Covalent15, covalentMap);
4444
4445
4446 OpenMM_IntArray_resize(covalentMap, 0);
4447 for (int j = 0; j < ip11[i].length; j++) {
4448 OpenMM_IntArray_append(covalentMap, ip11[i][j]);
4449 }
4450 OpenMM_AmoebaMultipoleForce_setCovalentMap(amoebaMultipoleForce, i,
4451 OpenMM_AmoebaMultipoleForce_PolarizationCovalent11, covalentMap);
4452
4453
4454 }
4455
4456 OpenMM_IntArray_destroy(covalentMap);
4457
4458 int forceGroup = forceField.getInteger("PME_FORCE_GROUP", 1);
4459 OpenMM_Force_setForceGroup(amoebaMultipoleForce, forceGroup);
4460 OpenMM_System_addForce(system, amoebaMultipoleForce);
4461
4462 logger.log(Level.INFO, format(" AMOEBA polarizable multipole force \t%d", forceGroup));
4463
4464 GeneralizedKirkwood gk = getGK();
4465 if (gk != null) {
4466 addGeneralizedKirkwoodForce();
4467 }
4468
4469 if (scfAlgorithm == SCFAlgorithm.EPT) {
4470 logger.info(" Using extrapolated perturbation theory for polarization energy.");
4471 }
4472 }
4473
4474
4475
4476
4477
4478
4479 private void updateAmoebaMultipoleForce(Atom[] atoms) {
4480 ParticleMeshEwald pme = getPmeNode();
4481 double quadrupoleConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
4482 double polarityConversion = OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom;
4483 double dampingFactorConversion = sqrt(OpenMM_NmPerAngstrom);
4484
4485 double polarScale = 1.0;
4486 if (pme.getPolarizationType() == Polarization.NONE) {
4487 polarScale = 0.0;
4488 }
4489
4490 PointerByReference dipoles = OpenMM_DoubleArray_create(3);
4491 PointerByReference quadrupoles = OpenMM_DoubleArray_create(9);
4492
4493 for (Atom atom : atoms) {
4494 int index = atom.getXyzIndex() - 1;
4495 MultipoleType multipoleType = pme.getMultipoleType(index);
4496 PolarizeType polarizeType = pme.getPolarizeType(index);
4497 int[] axisAtoms = atom.getAxisAtomIndices();
4498
4499 double useFactor = 1.0;
4500 if (!atom.getUse() || !atom.getElectrostatics()) {
4501 useFactor = 0.0;
4502 }
4503
4504 double lambdaScale = lambdaElec;
4505 if (!atom.applyLambda()) {
4506 lambdaScale = 1.0;
4507 }
4508 useFactor *= lambdaScale;
4509
4510
4511 int axisType = switch (multipoleType.frameDefinition) {
4512 case NONE -> OpenMM_AmoebaMultipoleForce_NoAxisType;
4513 case ZONLY -> OpenMM_AmoebaMultipoleForce_ZOnly;
4514 case ZTHENX -> OpenMM_AmoebaMultipoleForce_ZThenX;
4515 case BISECTOR -> OpenMM_AmoebaMultipoleForce_Bisector;
4516 case ZTHENBISECTOR -> OpenMM_AmoebaMultipoleForce_ZBisect;
4517 case THREEFOLD -> OpenMM_AmoebaMultipoleForce_ThreeFold;
4518 };
4519
4520
4521 for (int j = 0; j < 3; j++) {
4522 OpenMM_DoubleArray_set(dipoles, j,
4523 multipoleType.dipole[j] * OpenMM_NmPerAngstrom * useFactor);
4524 }
4525 int l = 0;
4526 for (int j = 0; j < 3; j++) {
4527 for (int k = 0; k < 3; k++) {
4528 OpenMM_DoubleArray_set(quadrupoles, l++,
4529 multipoleType.quadrupole[j][k] * quadrupoleConversion / 3.0 * useFactor);
4530 }
4531 }
4532
4533 int zaxis = 1;
4534 int xaxis = 1;
4535 int yaxis = 1;
4536
4537 if (axisAtoms != null) {
4538 zaxis = axisAtoms[0];
4539 if (axisAtoms.length > 1) {
4540 xaxis = axisAtoms[1];
4541 if (axisAtoms.length > 2) {
4542 yaxis = axisAtoms[2];
4543 }
4544 }
4545 } else {
4546 axisType = OpenMM_AmoebaMultipoleForce_NoAxisType;
4547 }
4548
4549
4550 OpenMM_AmoebaMultipoleForce_setMultipoleParameters(amoebaMultipoleForce, index,
4551 multipoleType.charge * useFactor, dipoles, quadrupoles, axisType, zaxis, xaxis, yaxis,
4552 polarizeType.thole, polarizeType.pdamp * dampingFactorConversion,
4553 polarizeType.polarizability * polarityConversion * polarScale * useFactor);
4554 }
4555
4556 OpenMM_DoubleArray_destroy(dipoles);
4557 OpenMM_DoubleArray_destroy(quadrupoles);
4558
4559 if (context.contextPointer != null) {
4560 OpenMM_AmoebaMultipoleForce_updateParametersInContext(amoebaMultipoleForce,
4561 context.contextPointer);
4562 }
4563 }
4564
4565
4566 private void addGeneralizedKirkwoodForce() {
4567 GeneralizedKirkwood gk = getGK();
4568
4569 amoebaGeneralizedKirkwoodForce = OpenMM_AmoebaGeneralizedKirkwoodForce_create();
4570 OpenMM_AmoebaGeneralizedKirkwoodForce_setSolventDielectric(amoebaGeneralizedKirkwoodForce,
4571 gk.getSolventPermittivity());
4572 OpenMM_AmoebaGeneralizedKirkwoodForce_setSoluteDielectric(amoebaGeneralizedKirkwoodForce, 1.0);
4573 OpenMM_AmoebaGeneralizedKirkwoodForce_setDielectricOffset(amoebaGeneralizedKirkwoodForce,
4574 gk.getDescreenOffset() * OpenMM_NmPerAngstrom);
4575
4576 boolean usePerfectRadii = gk.getUsePerfectRadii();
4577 double perfectRadiiScale = 1.0;
4578 if (usePerfectRadii) {
4579
4580 perfectRadiiScale = 0.0;
4581 }
4582
4583
4584 int tanhRescale = 0;
4585 if (gk.getTanhCorrection() && !usePerfectRadii) {
4586 tanhRescale = 1;
4587 }
4588 double[] betas = gk.getTanhBetas();
4589 OpenMM_AmoebaGeneralizedKirkwoodForce_setTanhRescaling(amoebaGeneralizedKirkwoodForce,
4590 tanhRescale);
4591 OpenMM_AmoebaGeneralizedKirkwoodForce_setTanhParameters(amoebaGeneralizedKirkwoodForce,
4592 betas[0], betas[1], betas[2]);
4593
4594 double[] baseRadius = gk.getBaseRadii();
4595 if (usePerfectRadii) {
4596 baseRadius = gk.getPerfectRadii();
4597 }
4598
4599 double[] overlapScale = gk.getOverlapScale();
4600 double[] descreenRadius = gk.getDescreenRadii();
4601 double[] neckFactor = gk.getNeckScale();
4602
4603 if (!usePerfectRadii && logger.isLoggable(Level.FINE)) {
4604 logger.fine(" GK Base Radii Descreen Radius Overlap Scale Overlap");
4605 }
4606
4607 for (int i = 0; i < nAtoms; i++) {
4608 MultipoleType multipoleType = atoms[i].getMultipoleType();
4609 double base = baseRadius[i] * OpenMM_NmPerAngstrom;
4610 double descreen = descreenRadius[i] * OpenMM_NmPerAngstrom * perfectRadiiScale;
4611 double overlap = overlapScale[i] * perfectRadiiScale;
4612 double neck = neckFactor[i] * perfectRadiiScale;
4613 OpenMM_AmoebaGeneralizedKirkwoodForce_addParticle_1(amoebaGeneralizedKirkwoodForce,
4614 multipoleType.charge, base, overlap, descreen, neck);
4615
4616 if (!usePerfectRadii && logger.isLoggable(Level.FINE)) {
4617 logger.fine(format(" %s %8.6f %8.6f %5.3f", atoms[i].toString(), baseRadius[i],
4618 descreenRadius[i], overlapScale[i]));
4619 }
4620 }
4621
4622 OpenMM_AmoebaGeneralizedKirkwoodForce_setProbeRadius(amoebaGeneralizedKirkwoodForce,
4623 gk.getProbeRadius() * OpenMM_NmPerAngstrom);
4624
4625 NonPolarModel nonpolar = gk.getNonPolarModel();
4626 switch (nonpolar) {
4627 default -> {
4628
4629 double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal * OpenMM_AngstromsPerNm
4630 * OpenMM_AngstromsPerNm;
4631 OpenMM_AmoebaGeneralizedKirkwoodForce_setIncludeCavityTerm(amoebaGeneralizedKirkwoodForce,
4632 OpenMM_True);
4633 OpenMM_AmoebaGeneralizedKirkwoodForce_setSurfaceAreaFactor(amoebaGeneralizedKirkwoodForce,
4634 -surfaceTension);
4635 }
4636 case CAV, CAV_DISP, GAUSS_DISP, SEV_DISP, HYDROPHOBIC_PMF, NONE ->
4637
4638 OpenMM_AmoebaGeneralizedKirkwoodForce_setIncludeCavityTerm(
4639 amoebaGeneralizedKirkwoodForce, OpenMM_False);
4640 }
4641
4642 int forceGroup = forceField.getInteger("GK_FORCE_GROUP", 1);
4643 OpenMM_Force_setForceGroup(amoebaGeneralizedKirkwoodForce, forceGroup);
4644 OpenMM_System_addForce(system, amoebaGeneralizedKirkwoodForce);
4645
4646 logger.log(Level.INFO, format(" Generalized Kirkwood force \t\t%d", forceGroup));
4647
4648
4649 switch (nonpolar) {
4650 case CAV_DISP, GAUSS_DISP, SEV_DISP, BORN_CAV_DISP -> addWCAForce();
4651 default -> {
4652
4653 }
4654 }
4655
4656
4657 if (nonpolar == GAUSS_DISP) {
4658 addCavitationForce();
4659 }
4660 }
4661
4662
4663
4664
4665
4666
4667 private void updateGeneralizedKirkwoodForce(Atom[] atoms) {
4668 GeneralizedKirkwood gk = getGK();
4669
4670 for (int i = 0; i < nAtoms; i++) {
4671 gk.udpateSoluteParameters(i);
4672 }
4673
4674 boolean usePerfectRadii = gk.getUsePerfectRadii();
4675 double perfectRadiiScale = 1.0;
4676 if (usePerfectRadii) {
4677
4678 perfectRadiiScale = 0.0;
4679 }
4680
4681 double[] baseRadii = gk.getBaseRadii();
4682 if (usePerfectRadii) {
4683 baseRadii = gk.getPerfectRadii();
4684 }
4685 double[] overlapScale = gk.getOverlapScale();
4686 double[] descreenRadius = gk.getDescreenRadii();
4687 double[] neckFactors = gk.getNeckScale();
4688
4689 boolean nea = gk.getNativeEnvironmentApproximation();
4690
4691 for (Atom atom : atoms) {
4692 int index = atom.getXyzIndex() - 1;
4693 double chargeUseFactor = 1.0;
4694 if (!atom.getUse() || !atom.getElectrostatics()) {
4695 chargeUseFactor = 0.0;
4696 }
4697
4698 double lambdaScale = lambdaElec;
4699 if (!atom.applyLambda()) {
4700 lambdaScale = 1.0;
4701 }
4702
4703 double baseSize = baseRadii[index] * OpenMM_NmPerAngstrom;
4704 double descreenSize = descreenRadius[index] * OpenMM_NmPerAngstrom * perfectRadiiScale;
4705
4706 chargeUseFactor *= lambdaScale;
4707 double overlapScaleUseFactor = nea ? 1.0 : chargeUseFactor;
4708 overlapScaleUseFactor = overlapScaleUseFactor * perfectRadiiScale;
4709 double overlap = overlapScale[index] * overlapScaleUseFactor;
4710 double neckFactor = neckFactors[index] * overlapScaleUseFactor;
4711
4712 MultipoleType multipoleType = atom.getMultipoleType();
4713 OpenMM_AmoebaGeneralizedKirkwoodForce_setParticleParameters_1(amoebaGeneralizedKirkwoodForce,
4714 index, multipoleType.charge * chargeUseFactor, baseSize, overlap, descreenSize,
4715 neckFactor);
4716 }
4717
4718
4719
4720
4721
4722
4723
4724
4725
4726
4727
4728
4729
4730
4731
4732
4733
4734
4735
4736
4737
4738
4739
4740 if (context.contextPointer != null) {
4741 OpenMM_AmoebaGeneralizedKirkwoodForce_updateParametersInContext(
4742 amoebaGeneralizedKirkwoodForce, context.contextPointer);
4743 }
4744 }
4745
4746
4747 private void addWCAForce() {
4748
4749 GeneralizedKirkwood gk = getGK();
4750 DispersionRegion dispersionRegion = gk.getDispersionRegion();
4751
4752 double epso = 0.1100;
4753 double epsh = 0.0135;
4754 double rmino = 1.7025;
4755 double rminh = 1.3275;
4756 double awater = 0.033428;
4757 double slevy = 1.0;
4758 double dispoff = dispersionRegion.getDispersionOffset();
4759 double shctd = dispersionRegion.getDispersionOverlapFactor();
4760
4761 VanDerWaals vdW = getVdwNode();
4762 VanDerWaalsForm vdwForm = vdW.getVDWForm();
4763 double radScale = 1.0;
4764 if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
4765 radScale = 0.5;
4766 }
4767
4768 amoebaWcaDispersionForce = OpenMM_AmoebaWcaDispersionForce_create();
4769
4770 for (Atom atom : atoms) {
4771 VDWType vdwType = atom.getVDWType();
4772 double radius = vdwType.radius;
4773 double eps = vdwType.wellDepth;
4774 OpenMM_AmoebaWcaDispersionForce_addParticle(amoebaWcaDispersionForce,
4775 OpenMM_NmPerAngstrom * radius * radScale, OpenMM_KJPerKcal * eps);
4776 }
4777
4778 OpenMM_AmoebaWcaDispersionForce_setEpso(amoebaWcaDispersionForce, epso * OpenMM_KJPerKcal);
4779 OpenMM_AmoebaWcaDispersionForce_setEpsh(amoebaWcaDispersionForce, epsh * OpenMM_KJPerKcal);
4780 OpenMM_AmoebaWcaDispersionForce_setRmino(amoebaWcaDispersionForce,
4781 rmino * OpenMM_NmPerAngstrom);
4782 OpenMM_AmoebaWcaDispersionForce_setRminh(amoebaWcaDispersionForce,
4783 rminh * OpenMM_NmPerAngstrom);
4784 OpenMM_AmoebaWcaDispersionForce_setDispoff(amoebaWcaDispersionForce,
4785 dispoff * OpenMM_NmPerAngstrom);
4786 OpenMM_AmoebaWcaDispersionForce_setAwater(amoebaWcaDispersionForce,
4787 awater / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
4788 OpenMM_AmoebaWcaDispersionForce_setSlevy(amoebaWcaDispersionForce, slevy);
4789 OpenMM_AmoebaWcaDispersionForce_setShctd(amoebaWcaDispersionForce, shctd);
4790
4791 int forceGroup = forceField.getInteger("GK_FORCE_GROUP", 1);
4792
4793 OpenMM_Force_setForceGroup(amoebaWcaDispersionForce, forceGroup);
4794 OpenMM_System_addForce(system, amoebaWcaDispersionForce);
4795
4796 logger.log(Level.INFO, format(" WCA dispersion force \t\t\t%d", forceGroup));
4797 }
4798
4799
4800
4801
4802
4803
4804 private void updateWCAForce(Atom[] atoms) {
4805 VanDerWaals vdW = getVdwNode();
4806 VanDerWaalsForm vdwForm = vdW.getVDWForm();
4807 double radScale = 1.0;
4808 if (vdwForm.radiusSize == VanDerWaalsForm.RADIUS_SIZE.DIAMETER) {
4809 radScale = 0.5;
4810 }
4811
4812 for (Atom atom : atoms) {
4813 int index = atom.getXyzIndex() - 1;
4814 double useFactor = 1.0;
4815 if (!atom.getUse()) {
4816 useFactor = 0.0;
4817 }
4818
4819
4820
4821 double lambdaScale = lambdaElec * lambdaElec;
4822 if (!atom.applyLambda()) {
4823 lambdaScale = 1.0;
4824 }
4825 useFactor *= lambdaScale;
4826
4827 VDWType vdwType = atom.getVDWType();
4828 double radius = vdwType.radius;
4829 double eps = vdwType.wellDepth;
4830 OpenMM_AmoebaWcaDispersionForce_setParticleParameters(amoebaWcaDispersionForce, index,
4831 OpenMM_NmPerAngstrom * radius * radScale, OpenMM_KJPerKcal * eps * useFactor);
4832 }
4833
4834 if (context.contextPointer != null) {
4835 OpenMM_AmoebaWcaDispersionForce_updateParametersInContext(amoebaWcaDispersionForce,
4836 context.contextPointer);
4837 }
4838 }
4839
4840
4841 private void addCavitationForce() {
4842
4843 GeneralizedKirkwood generalizedKirkwood = getGK();
4844 ChandlerCavitation chandlerCavitation = generalizedKirkwood.getChandlerCavitation();
4845 GaussVol gaussVol = chandlerCavitation.getGaussVol();
4846 if (gaussVol == null) {
4847 return;
4848 }
4849
4850 amoebaCavitationForce = OpenMM_AmoebaGKCavitationForce_create();
4851 double surfaceTension =
4852 chandlerCavitation.getSurfaceTension() * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom
4853 / OpenMM_NmPerAngstrom;
4854 double[] rad = gaussVol.getRadii();
4855
4856 int index = 0;
4857 for (Atom atom : atoms) {
4858 int isHydrogen = OpenMM_False;
4859 double radius = rad[index++];
4860 if (atom.isHydrogen()) {
4861 isHydrogen = OpenMM_True;
4862 radius = 0.0;
4863 }
4864 OpenMM_AmoebaGKCavitationForce_addParticle(amoebaCavitationForce,
4865 radius * OpenMM_NmPerAngstrom, surfaceTension, isHydrogen);
4866 }
4867
4868 OpenMM_AmoebaGKCavitationForce_setNonbondedMethod(amoebaCavitationForce,
4869 OpenMM_AmoebaGKCavitationForce_NoCutoff);
4870
4871 int forceGroup = forceField.getInteger("GK_FORCE_GROUP", 2);
4872 OpenMM_Force_setForceGroup(amoebaCavitationForce, forceGroup);
4873 OpenMM_System_addForce(system, amoebaCavitationForce);
4874
4875 logger.log(Level.INFO, format(" GaussVol cavitation force \t\t%d", forceGroup));
4876 }
4877
4878
4879
4880
4881
4882
4883 private void updateCavitationForce(Atom[] atoms) {
4884 GeneralizedKirkwood generalizedKirkwood = getGK();
4885 ChandlerCavitation chandlerCavitation = generalizedKirkwood.getChandlerCavitation();
4886 GaussVol gaussVol = chandlerCavitation.getGaussVol();
4887 if (gaussVol == null) {
4888 return;
4889 }
4890
4891 double surfaceTension =
4892 chandlerCavitation.getSurfaceTension() * OpenMM_KJPerKcal / OpenMM_NmPerAngstrom
4893 / OpenMM_NmPerAngstrom;
4894
4895
4896
4897
4898
4899 double[] rad = gaussVol.getRadii();
4900
4901 for (Atom atom : atoms) {
4902 int index = atom.getXyzIndex() - 1;
4903 double useFactor = 1.0;
4904 if (!atom.getUse()) {
4905 useFactor = 0.0;
4906 }
4907
4908
4909 double lambdaScale = lambdaElec * lambdaElec;
4910 if (!atom.applyLambda()) {
4911 lambdaScale = 1.0;
4912 }
4913 useFactor *= lambdaScale;
4914
4915 double radius = rad[index];
4916 int isHydrogen = OpenMM_False;
4917 if (atom.isHydrogen()) {
4918 isHydrogen = OpenMM_True;
4919 radius = 0.0;
4920 }
4921
4922 OpenMM_AmoebaGKCavitationForce_setParticleParameters(amoebaCavitationForce, index,
4923 radius * OpenMM_NmPerAngstrom, surfaceTension * useFactor, isHydrogen);
4924 }
4925
4926 if (context.contextPointer != null) {
4927 OpenMM_AmoebaGKCavitationForce_updateParametersInContext(amoebaCavitationForce,
4928 context.contextPointer);
4929 }
4930 }
4931
4932
4933
4934
4935
4936
4937 private void addHarmonicRestraintForce() {
4938 int nRestraints = getCoordRestraints().size();
4939
4940 int forceGroup = forceField.getInteger("COORD_RESTRAINT_FORCE_GROUP", 0);
4941
4942 for (RestrainPosition restrainPosition : getCoordRestraints()) {
4943 double forceConstant = restrainPosition.getForceConstant();
4944 forceConstant *= OpenMM_KJPerKcal;
4945 forceConstant *= (OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm);
4946 Atom[] restAtoms = restrainPosition.getAtoms();
4947 int nRestAts = restrainPosition.getNumAtoms();
4948 double[][] oCoords = restrainPosition.getEquilibriumCoordinates();
4949 for (int i = 0; i < nRestAts; i++) {
4950 oCoords[i][0] *= OpenMM_NmPerAngstrom;
4951 oCoords[i][1] *= OpenMM_NmPerAngstrom;
4952 oCoords[i][2] *= OpenMM_NmPerAngstrom;
4953 }
4954
4955 PointerByReference theRestraint = OpenMM_CustomExternalForce_create(
4956 "k*periodicdistance(x,y,z,x0,y0,z0)^2");
4957 OpenMM_CustomExternalForce_addGlobalParameter(theRestraint, "k", forceConstant);
4958 OpenMM_CustomExternalForce_addPerParticleParameter(theRestraint, "x0");
4959 OpenMM_CustomExternalForce_addPerParticleParameter(theRestraint, "y0");
4960 OpenMM_CustomExternalForce_addPerParticleParameter(theRestraint, "z0");
4961
4962 PointerByReference xyzOrigArray = OpenMM_DoubleArray_create(3);
4963 for (int i = 0; i < nRestAts; i++) {
4964 int ommIndex = restAtoms[i].getXyzIndex() - 1;
4965 for (int j = 0; j < 3; j++) {
4966 OpenMM_DoubleArray_set(xyzOrigArray, j, oCoords[i][j]);
4967 }
4968 OpenMM_CustomExternalForce_addParticle(theRestraint, ommIndex, xyzOrigArray);
4969 }
4970 OpenMM_DoubleArray_destroy(xyzOrigArray);
4971
4972 OpenMM_Force_setForceGroup(theRestraint, forceGroup);
4973 OpenMM_System_addForce(system, theRestraint);
4974 }
4975
4976 if (nRestraints > 0) {
4977 logger.log(Level.INFO,
4978 format(" Harmonic restraint force \t%6d\t%d", nRestraints, forceGroup));
4979 }
4980 }
4981
4982
4983 private void addRestraintBondForce() {
4984 List<RestraintBond> restraintBonds = getRestraintBonds();
4985 if (restraintBonds == null || restraintBonds.isEmpty()) {
4986 return;
4987 }
4988
4989 int forceGroup = forceField.getInteger("BOND_RESTRAINT_FORCE_GROUP", 0);
4990
4991
4992 double kParameterConversion =
4993 2.0 * OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
4994
4995
4996 Map<BondType.BondFunction, PointerByReference> restraintForces = new HashMap<>();
4997
4998 for (RestraintBond restraintBond : getRestraintBonds()) {
4999 PointerByReference theForce;
5000 BondType bondType = restraintBond.bondType;
5001 BondType.BondFunction bondFunction = bondType.bondFunction;
5002 if (restraintForces.containsKey(bondFunction)) {
5003 theForce = restraintForces.get(bondFunction);
5004 } else {
5005 theForce = OpenMM_CustomBondForce_create(bondFunction.toMathematicalForm());
5006 OpenMM_CustomBondForce_addPerBondParameter(theForce, "k");
5007 OpenMM_CustomBondForce_addPerBondParameter(theForce, "r0");
5008 if (bondFunction.hasFlatBottom()) {
5009 OpenMM_CustomBondForce_addPerBondParameter(theForce, "fb");
5010 }
5011
5012 switch (bondFunction) {
5013 case QUARTIC, FLAT_BOTTOM_QUARTIC -> {
5014 OpenMM_CustomBondForce_addGlobalParameter(theForce, "cubic",
5015 bondType.cubic / OpenMM_NmPerAngstrom);
5016 OpenMM_CustomBondForce_addGlobalParameter(theForce, "quartic",
5017 bondType.quartic / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom));
5018 }
5019 default -> {
5020
5021 }
5022 }
5023
5024 OpenMM_Force_setForceGroup(theForce, forceGroup);
5025 OpenMM_System_addForce(system, theForce);
5026 restraintForces.put(bondFunction, theForce);
5027 }
5028
5029 double forceConst = bondType.forceConstant * bondType.bondUnit * kParameterConversion;
5030 double equilDist = bondType.distance * OpenMM_NmPerAngstrom;
5031 Atom[] ats = restraintBond.getAtomArray();
5032 int at1 = ats[0].getXyzIndex() - 1;
5033 int at2 = ats[1].getXyzIndex() - 1;
5034
5035 PointerByReference bondParams = OpenMM_DoubleArray_create(0);
5036 OpenMM_DoubleArray_append(bondParams, forceConst);
5037 OpenMM_DoubleArray_append(bondParams, equilDist);
5038 if (bondFunction.hasFlatBottom()) {
5039 OpenMM_DoubleArray_append(bondParams, bondType.flatBottomRadius * OpenMM_NmPerAngstrom);
5040 }
5041 OpenMM_CustomBondForce_addBond(theForce, at1, at2, bondParams);
5042 OpenMM_DoubleArray_destroy(bondParams);
5043 }
5044
5045 logger.log(Level.INFO,
5046 format(" Restraint bonds force \t%6d\t%d", restraintBonds.size(), forceGroup));
5047 }
5048
5049 private void addRestrainGroupsForce() {
5050 RestrainGroups restrainGroups = getRestrainGroups();
5051 if (restrainGroups == null) {
5052 return;
5053 }
5054
5055
5056 PointerByReference force = OpenMM_CustomCentroidBondForce_create(2,
5057 "step(distance(g1,g2)-u)*k*(distance(g1,g2)-u)^2+step(l-distance(g1,g2))*k*(distance(g1,g2)-l)^2");
5058 OpenMM_CustomCentroidBondForce_addPerBondParameter(force, "k");
5059 OpenMM_CustomCentroidBondForce_addPerBondParameter(force, "l");
5060 OpenMM_CustomCentroidBondForce_addPerBondParameter(force, "u");
5061
5062
5063 int nGroups = restrainGroups.getNumberOfGroups();
5064 for (int j = 0; j < nGroups; j++) {
5065 PointerByReference group = OpenMM_IntArray_create(0);
5066 PointerByReference weight = OpenMM_DoubleArray_create(0);
5067 int[] groupMembers = restrainGroups.getGroupMembers(j);
5068 for (int i : groupMembers) {
5069 OpenMM_IntArray_append(group, i);
5070 OpenMM_DoubleArray_append(weight, atoms[i].getMass());
5071 }
5072 OpenMM_CustomCentroidBondForce_addGroup(force, group, weight);
5073 OpenMM_IntArray_destroy(group);
5074 OpenMM_DoubleArray_destroy(weight);
5075 }
5076
5077
5078 double convert = OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
5079 int nRestraints = restrainGroups.getNumberOfRestraints();
5080 int[] group1 = restrainGroups.getGroup1();
5081 int[] group2 = restrainGroups.getGroup2();
5082 double[] forceConstants = restrainGroups.getForceConstants();
5083 double[] smallerDistance = restrainGroups.getSmallerDistance();
5084 double[] largerDistance = restrainGroups.getLargerDistance();
5085 for (int i = 0; i < nRestraints; i++) {
5086 PointerByReference group = OpenMM_IntArray_create(0);
5087 OpenMM_IntArray_append(group, group1[i]);
5088 OpenMM_IntArray_append(group, group2[i]);
5089 PointerByReference params = OpenMM_DoubleArray_create(0);
5090 OpenMM_DoubleArray_append(params, forceConstants[i] * convert);
5091 OpenMM_DoubleArray_append(params, smallerDistance[i] * OpenMM_NmPerAngstrom);
5092 OpenMM_DoubleArray_append(params, largerDistance[i] * OpenMM_NmPerAngstrom);
5093 OpenMM_CustomCentroidBondForce_addBond(force, group, params);
5094 OpenMM_IntArray_destroy(group);
5095 OpenMM_DoubleArray_destroy(params);
5096 }
5097
5098
5099 int forceGroup = forceField.getInteger("BOND_RESTRAINT_FORCE_GROUP", 0);
5100 OpenMM_Force_setForceGroup(force, forceGroup);
5101
5102 if (getCrystal().aperiodic()) {
5103 OpenMM_CustomCentroidBondForce_setUsesPeriodicBoundaryConditions(force, OpenMM_False);
5104 } else {
5105 OpenMM_CustomCentroidBondForce_setUsesPeriodicBoundaryConditions(force, OpenMM_True);
5106 }
5107 OpenMM_System_addForce(system, force);
5108 logger.log(Level.INFO, format(" Restrain Groups \t%6d\t\t%1d", nRestraints, forceGroup));
5109 }
5110
5111
5112 private void addUpBondConstraints() {
5113 Bond[] bonds = getBonds();
5114 if (bonds == null || bonds.length < 1) {
5115 return;
5116 }
5117
5118 logger.info(" Adding constraints for all bonds.");
5119 for (Bond bond : bonds) {
5120 Atom atom1 = bond.getAtom(0);
5121 Atom atom2 = bond.getAtom(1);
5122 int iAtom1 = atom1.getXyzIndex() - 1;
5123 int iAtom2 = atom2.getXyzIndex() - 1;
5124 OpenMM_System_addConstraint(system, iAtom1, iAtom2,
5125 bond.bondType.distance * OpenMM_NmPerAngstrom);
5126 }
5127 }
5128
5129
5130 private void addHydrogenConstraints() {
5131 Bond[] bonds = getBonds();
5132 if (bonds == null || bonds.length < 1) {
5133 return;
5134 }
5135
5136 logger.info(" Adding constraints for hydrogen bonds.");
5137 for (Bond bond : bonds) {
5138 Atom atom1 = bond.getAtom(0);
5139 Atom atom2 = bond.getAtom(1);
5140 if (atom1.isHydrogen() || atom2.isHydrogen()) {
5141 BondType bondType = bond.bondType;
5142 int iAtom1 = atom1.getXyzIndex() - 1;
5143 int iAtom2 = atom2.getXyzIndex() - 1;
5144 OpenMM_System_addConstraint(system, iAtom1, iAtom2,
5145 bondType.distance * OpenMM_NmPerAngstrom);
5146 }
5147 }
5148 }
5149
5150
5151 private void setUpHydrogenAngleConstraints() {
5152 Angle[] angles = getAngles();
5153 if (angles == null || angles.length < 1) {
5154 return;
5155 }
5156
5157 logger.info(" Adding hydrogen angle constraints.");
5158
5159 for (Angle angle : angles) {
5160 if (isHydrogenAngle(angle)) {
5161 Atom atom1 = angle.getAtom(0);
5162 Atom atom3 = angle.getAtom(2);
5163
5164
5165
5166 Bond bond1 = angle.getBond(0);
5167 double distance1 = bond1.bondType.distance;
5168
5169 Bond bond2 = angle.getBond(1);
5170 double distance2 = bond2.bondType.distance;
5171
5172
5173 double angleVal = angle.angleType.angle[angle.nh];
5174
5175
5176 double falseBondLength = sqrt(
5177 distance1 * distance1 + distance2 * distance2 - 2.0 * distance1 * distance2 * cos(
5178 toRadians(angleVal)));
5179
5180 int iAtom1 = atom1.getXyzIndex() - 1;
5181 int iAtom3 = atom3.getXyzIndex() - 1;
5182 OpenMM_System_addConstraint(system, iAtom1, iAtom3,
5183 falseBondLength * OpenMM_NmPerAngstrom);
5184 }
5185 }
5186 }
5187
5188
5189
5190
5191
5192
5193
5194
5195 private boolean isHydrogenAngle(Angle angle) {
5196 if (angle.containsHydrogen()) {
5197
5198 double angleVal = angle.angleType.angle[angle.nh];
5199
5200
5201
5202 if (angleVal < 160.0) {
5203 Atom atom1 = angle.getAtom(0);
5204 Atom atom2 = angle.getAtom(1);
5205 Atom atom3 = angle.getAtom(2);
5206
5207
5208 return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
5209 }
5210 }
5211 return false;
5212 }
5213
5214 }
5215
5216
5217 public class State {
5218
5219
5220 public final double potentialEnergy;
5221
5222 public final double kineticEnergy;
5223
5224 public final double totalEnergy;
5225
5226 public final double temperature;
5227
5228 private final PointerByReference state;
5229
5230 private final boolean positions;
5231
5232 private final boolean velocities;
5233
5234 private final boolean forces;
5235
5236
5237
5238
5239
5240
5241
5242
5243
5244 public State(boolean positions, boolean energies, boolean forces, boolean velocities) {
5245 this.positions = positions;
5246 this.forces = forces;
5247 this.velocities = velocities;
5248
5249
5250 int mask = 0;
5251 if (positions) {
5252 mask = OpenMM_State_Positions;
5253 }
5254 if (energies) {
5255 mask += OpenMM_State_Energy;
5256 }
5257 if (velocities) {
5258 mask += OpenMM_State_Velocities;
5259 }
5260 if (forces) {
5261 mask += OpenMM_State_Forces;
5262 }
5263
5264
5265 state = OpenMM_Context_getState(context.getContextPointer(), mask, enforcePBC);
5266
5267 if (energies) {
5268 potentialEnergy = OpenMM_State_getPotentialEnergy(state) * OpenMM_KcalPerKJ;
5269 kineticEnergy = OpenMM_State_getKineticEnergy(state) * OpenMM_KcalPerKJ;
5270 totalEnergy = potentialEnergy + kineticEnergy;
5271 double dof = system.calculateDegreesOfFreedom();
5272 temperature = 2.0 * kineticEnergy * KCAL_TO_GRAM_ANG2_PER_PS2 / (kB * dof);
5273 } else {
5274 potentialEnergy = 0.0;
5275 kineticEnergy = 0.0;
5276 totalEnergy = 0.0;
5277 temperature = 0.0;
5278 }
5279 }
5280
5281 public void free() {
5282 if (state != null && freeOpenMM) {
5283 logger.fine(" Free OpenMM State.");
5284 OpenMM_State_destroy(state);
5285 logger.fine(" Free OpenMM State completed.");
5286 }
5287 }
5288
5289
5290
5291
5292
5293
5294
5295
5296 public double[] getAccelerations(double[] a) {
5297 if (!forces) {
5298 return a;
5299 }
5300
5301 int n = getNumberOfVariables();
5302 if (a == null || a.length != n) {
5303 a = new double[n];
5304 }
5305
5306 PointerByReference forcePointer = OpenMM_State_getForces(state);
5307
5308 for (int i = 0, index = 0; i < nAtoms; i++) {
5309 Atom atom = atoms[i];
5310 if (atom.isActive()) {
5311 double mass = atom.getMass();
5312 OpenMM_Vec3 acc = OpenMM_Vec3Array_get(forcePointer, i);
5313 double xx = acc.x * OpenMM_AngstromsPerNm / mass;
5314 double yy = acc.y * OpenMM_AngstromsPerNm / mass;
5315 double zz = acc.z * OpenMM_AngstromsPerNm / mass;
5316 a[index++] = xx;
5317 a[index++] = yy;
5318 a[index++] = zz;
5319 atom.setAcceleration(xx, yy, zz);
5320 }
5321 }
5322 return a;
5323 }
5324
5325
5326
5327
5328
5329
5330
5331
5332 public double[] getGradient(double[] g) {
5333 if (!forces) {
5334 return g;
5335 }
5336
5337 int n = getNumberOfVariables();
5338 if (g == null || g.length != n) {
5339 g = new double[n];
5340 }
5341
5342 PointerByReference forcePointer = OpenMM_State_getForces(state);
5343 for (int i = 0, index = 0; i < nAtoms; i++) {
5344 Atom atom = atoms[i];
5345 if (atom.isActive()) {
5346 OpenMM_Vec3 force = OpenMM_Vec3Array_get(forcePointer, i);
5347 double xx = -force.x * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
5348 double yy = -force.y * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
5349 double zz = -force.z * OpenMM_NmPerAngstrom * OpenMM_KcalPerKJ;
5350 if (isNaN(xx) || isInfinite(xx) || isNaN(yy) || isInfinite(yy) || isNaN(zz) || isInfinite(
5351 zz)) {
5352 StringBuilder sb = new StringBuilder(
5353 format(" The gradient of atom %s is (%8.3f,%8.3f,%8.3f).", atom, xx, yy, zz));
5354 double[] vals = new double[3];
5355 atom.getVelocity(vals);
5356 sb.append(format("\n Velocities: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
5357 atom.getAcceleration(vals);
5358 sb.append(format("\n Accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
5359 atom.getPreviousAcceleration(vals);
5360 sb.append(
5361 format("\n Previous accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
5362 throw new EnergyException(sb.toString());
5363 }
5364 g[index++] = xx;
5365 g[index++] = yy;
5366 g[index++] = zz;
5367 atom.setXYZGradient(xx, yy, zz);
5368 }
5369 }
5370 return g;
5371 }
5372
5373
5374
5375
5376
5377
5378 public void getPeriodicBoxVectors() {
5379 if (!positions) {
5380 return;
5381 }
5382
5383 Crystal crystal = getCrystal();
5384 if (!crystal.aperiodic()) {
5385 OpenMM_Vec3 a = new OpenMM_Vec3();
5386 OpenMM_Vec3 b = new OpenMM_Vec3();
5387 OpenMM_Vec3 c = new OpenMM_Vec3();
5388 OpenMM_State_getPeriodicBoxVectors(state, a, b, c);
5389 double[][] latticeVectors = new double[3][3];
5390 latticeVectors[0][0] = a.x * OpenMM_AngstromsPerNm;
5391 latticeVectors[0][1] = a.y * OpenMM_AngstromsPerNm;
5392 latticeVectors[0][2] = a.z * OpenMM_AngstromsPerNm;
5393 latticeVectors[1][0] = b.x * OpenMM_AngstromsPerNm;
5394 latticeVectors[1][1] = b.y * OpenMM_AngstromsPerNm;
5395 latticeVectors[1][2] = b.z * OpenMM_AngstromsPerNm;
5396 latticeVectors[2][0] = c.x * OpenMM_AngstromsPerNm;
5397 latticeVectors[2][1] = c.y * OpenMM_AngstromsPerNm;
5398 latticeVectors[2][2] = c.z * OpenMM_AngstromsPerNm;
5399 crystal.setCellVectors(latticeVectors);
5400 setCrystal(crystal);
5401 }
5402 }
5403
5404
5405
5406
5407
5408
5409
5410
5411 public double[] getPositions(double[] x) {
5412 if (!positions) {
5413 return x;
5414 }
5415
5416 int n = getNumberOfVariables();
5417 if (x == null || x.length != n) {
5418 x = new double[n];
5419 }
5420
5421 PointerByReference positionsPointer = OpenMM_State_getPositions(state);
5422 for (int i = 0, index = 0; i < nAtoms; i++) {
5423 Atom atom = atoms[i];
5424 if (atom.isActive()) {
5425 OpenMM_Vec3 pos = OpenMM_Vec3Array_get(positionsPointer, i);
5426 double xx = pos.x * OpenMM_AngstromsPerNm;
5427 double yy = pos.y * OpenMM_AngstromsPerNm;
5428 double zz = pos.z * OpenMM_AngstromsPerNm;
5429 x[index++] = xx;
5430 x[index++] = yy;
5431 x[index++] = zz;
5432 atom.moveTo(xx, yy, zz);
5433 }
5434 }
5435 return x;
5436 }
5437
5438
5439
5440
5441
5442
5443
5444
5445 public double[] getVelocities(double[] v) {
5446 if (!velocities) {
5447 return v;
5448 }
5449
5450 int n = getNumberOfVariables();
5451 if (v == null || v.length != n) {
5452 v = new double[n];
5453 }
5454
5455 PointerByReference velocitiesPointer = OpenMM_State_getVelocities(state);
5456 for (int i = 0, index = 0; i < nAtoms; i++) {
5457 Atom atom = atoms[i];
5458 if (atom.isActive()) {
5459 OpenMM_Vec3 vel = OpenMM_Vec3Array_get(velocitiesPointer, i);
5460 double xx = vel.x * OpenMM_AngstromsPerNm;
5461 double yy = vel.y * OpenMM_AngstromsPerNm;
5462 double zz = vel.z * OpenMM_AngstromsPerNm;
5463 v[index++] = xx;
5464 v[index++] = yy;
5465 v[index++] = zz;
5466 atom.setVelocity(xx, yy, zz);
5467 }
5468 }
5469 return v;
5470 }
5471 }
5472 }