View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2023.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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  * Compute the potential energy and derivatives using OpenMM.
383  *
384  * @author Michael J. Schnieders
385  * @since 1.0
386  */
387 @SuppressWarnings("deprecation")
388 public class ForceFieldEnergyOpenMM extends ForceFieldEnergy {
389 
390   private static final Logger logger = Logger.getLogger(ForceFieldEnergyOpenMM.class.getName());
391   /** Whether to enforce periodic boundary conditions when obtaining new States. */
392   public final int enforcePBC;
393   /** The atoms this ForceFieldEnergyOpenMM operates on. */
394   private final Atom[] atoms;
395   /** Number of particles. */
396   private final int nAtoms;
397   /** Lambda step size for finite difference dU/dL. */
398   private final double finiteDifferenceStepSize;
399   /** OpenMM Context. */
400   private Context context;
401   /** OpenMM System. */
402   private System system;
403   /**
404    * Truncate the normal OpenMM Lambda Path from 0 ... 1 to Lambda_Start ... 1. This is useful for
405    * conformational optimization if full removal of vdW interactions is not desired (i.e. lambdaStart
406    * = ~0.2).
407    */
408   private double lambdaStart = 0.0;
409   /** Use two-sided finite difference dU/dL. */
410   private boolean twoSidedFiniteDifference = true;
411 
412   private final boolean freeOpenMM;
413 
414   /**
415    * ForceFieldEnergyOpenMM constructor; offloads heavy-duty computation to an OpenMM Platform while
416    * keeping track of information locally.
417    *
418    * @param molecularAssembly Assembly to construct energy for.
419    * @param requestedPlatform requested OpenMM platform to be used.
420    * @param nThreads Number of threads to use in the super class ForceFieldEnergy instance.
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    * Gets the default coprocessor device, ignoring any CUDA_DEVICE over-ride. This is either
453    * determined by process rank and the availableDevices/CUDA_DEVICES property, or just 0 if neither
454    * property is sets.
455    *
456    * @param props Properties in use.
457    * @return Pre-override device index.
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         // Format the host as a CharacterBuf of length 100.
484         int messageLen = 100;
485         String host = world.host();
486         // Truncate to max 100 characters.
487         host = host.substring(0, Math.min(messageLen, host.length()));
488         // Pad to 100 characters.
489         host = String.format("%-100s", host);
490         char[] messageOut = host.toCharArray();
491         CharacterBuf out = CharacterBuf.buffer(messageOut);
492 
493         // Now create CharacterBuf array for all incoming messages.
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       // Behavior is just to keep index = 0.
531     }
532     return devs[index];
533   }
534 
535   /**
536    * Create a JNA Pointer to a String.
537    *
538    * @param string WARNING: assumes ascii-only string
539    * @return pointer.
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    * Returns the platform array as a String
549    *
550    * @param stringArray The OpenMM String array.
551    * @param i The index of the String to return.
552    * @return String The requested String.
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    * Create an OpenMM Context.
564    *
565    * <p>Context.free() must be called to free OpenMM memory.
566    *
567    * @param integratorString Integrator to use.
568    * @param timeStep Time step.
569    * @param temperature Temperature (K).
570    * @param forceCreation Force a new Context to be created, even if the existing one matches the
571    *     request.
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    * Create an immutable OpenMM State.
580    *
581    * <p>State.free() must be called to free OpenMM memory.
582    *
583    * @param positions Retrieve positions.
584    * @param energies Retrieve energies.
585    * @param forces Retrieve forces.
586    * @param velocities Retrieve velocities.
587    * @return Returns the State.
588    */
589   public State createState(boolean positions, boolean energies, boolean forces, boolean velocities) {
590     return new State(positions, energies, forces, velocities);
591   }
592 
593   /** {@inheritDoc} */
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   /** {@inheritDoc} */
605   @Override
606   public double energy(double[] x) {
607     return energy(x, false);
608   }
609 
610   /** {@inheritDoc} */
611   @Override
612   public double energy(double[] x, boolean verbose) {
613 
614     if (lambdaBondedTerms) {
615       return 0.0;
616     }
617 
618     // Make sure a context has been created.
619     context.getContextPointer();
620 
621     updateParameters(atoms);
622 
623     // Unscale the coordinates.
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     // Rescale the coordinates.
646     scaleCoordinates(x);
647 
648     return e;
649   }
650 
651   /** {@inheritDoc} */
652   @Override
653   public double energyAndGradient(double[] x, double[] g) {
654     return energyAndGradient(x, g, false);
655   }
656 
657   /** {@inheritDoc} */
658   @Override
659   public double energyAndGradient(double[] x, double[] g, boolean verbose) {
660     if (lambdaBondedTerms) {
661       return 0.0;
662     }
663 
664     // ZE BUG: updateParameters only gets called for energy(), not energyAndGradient().
665 
666     // Un-scale the coordinates.
667     unscaleCoordinates(x);
668 
669     // Make sure a context has been created.
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     // if (vdwLambdaTerm) {
689     //    PointerByReference parameterArray = OpenMM_State_getEnergyParameterDerivatives(state);
690     //    int numDerives = OpenMM_ParameterArray_getSize(parameterArray);
691     //    if (numDerives > 0) {
692     //        double vdwdUdL = OpenMM_ParameterArray_get(parameterArray,
693     // pointerForString("vdw_lambda")) / OpenMM_KJPerKcal;
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     // Scale the coordinates and gradients.
722     scaleCoordinatesAndGradient(x, g);
723 
724     return e;
725   }
726 
727   /**
728    * Compute the energy and gradient using the pure Java code path.
729    *
730    * @param x Input atomic coordinates
731    * @param g Storage for the gradient vector.
732    * @return The energy (kcal/mol)
733    */
734   public double energyAndGradientFFX(double[] x, double[] g) {
735     return super.energyAndGradient(x, g, false);
736   }
737 
738   /**
739    * Compute the energy and gradient using the pure Java code path.
740    *
741    * @param x Input atomic coordinates
742    * @param g Storage for the gradient vector.
743    * @param verbose Use verbose logging.
744    * @return The energy (kcal/mol)
745    */
746   public double energyAndGradientFFX(double[] x, double[] g, boolean verbose) {
747     return super.energyAndGradient(x, g, verbose);
748   }
749 
750   /**
751    * Compute the energy using the pure Java code path.
752    *
753    * @param x Atomic coordinates.
754    * @return The energy (kcal/mol)
755    */
756   public double energyFFX(double[] x) {
757     return super.energy(x, false);
758   }
759 
760   /**
761    * Compute the energy using the pure Java code path.
762    *
763    * @param x Input atomic coordinates
764    * @param verbose Use verbose logging.
765    * @return The energy (kcal/mol)
766    */
767   public double energyFFX(double[] x, boolean verbose) {
768     return super.energy(x, verbose);
769   }
770 
771   /**
772    * Returns the Context instance.
773    *
774    * @return context
775    */
776   public Context getContext() {
777     return context;
778   }
779 
780   /**
781    * Re-compute the gradient using OpenMM and return it.
782    *
783    * @param g Gradient array.
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   /** {@inheritDoc} */
794   @Override
795   public Platform getPlatform() {
796     return context.platform;
797   }
798 
799   /**
800    * Get a reference to the System instance.
801    *
802    * @return Java wrapper to an OpenMM system.
803    */
804   public System getSystem() {
805     return system;
806   }
807 
808   /** {@inheritDoc} */
809   @Override
810   public double getd2EdL2() {
811     return 0.0;
812   }
813 
814   /** {@inheritDoc} */
815   @Override
816   public double getdEdL() {
817     // No lambda dependence.
818     if (!lambdaTerm) {
819       return 0.0;
820     }
821 
822     // Small optimization to only create the x array once.
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         // Two sided finite difference estimate of dE/dL.
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       // One-sided finite difference estimates of dE/dL
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     // Compute the finite difference derivative.
867     double dEdL = (ePlus - eMinus) / width;
868 
869     // logger.info(format(" getdEdL currentLambda: CL=%8.6f L=%8.6f dEdL=%12.6f", currentLambda,
870     // lambda, dEdL));
871     return dEdL;
872   }
873 
874   /** {@inheritDoc} */
875   @Override
876   public void getdEdXdL(double[] gradients) {
877     // Note for ForceFieldEnergyOpenMM this method is not implemented.
878   }
879 
880   /** Update active atoms. */
881   public void setActiveAtoms() {
882     system.updateAtomMass();
883     // Tests show reinitialization of the OpenMM Context is not necessary to pick up mass changes.
884     // context.reinitContext();
885   }
886 
887   /**
888    * Set FFX and OpenMM coordinates for active atoms.
889    *
890    * @param x Atomic coordinates.
891    */
892   @Override
893   public void setCoordinates(double[] x) {
894     // Set both OpenMM and FFX coordinates to x.
895     context.setOpenMMPositions(x);
896   }
897 
898   /** {@inheritDoc} */
899   @Override
900   public void setCrystal(Crystal crystal) {
901     super.setCrystal(crystal);
902     context.setPeriodicBoxVectors();
903   }
904 
905   /** {@inheritDoc} */
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     // Check for lambda outside the range [0 .. 1].
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     // Remove the beginning of the normal Lambda path.
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         // Update force field parameters based on defined lambda values.
940         updateParameters(atomList.toArray(new Atom[0]));
941       } else {
942         updateParameters(null);
943       }
944     }
945   }
946 
947   /**
948    * Update parameters if the Use flags and/or Lambda value has changed.
949    *
950    * @param atoms Atoms in this list are considered.
951    */
952   public void updateParameters(Atom[] atoms) {
953     if (atoms == null) {
954       atoms = this.atoms;
955     }
956     system.updateParameters(atoms);
957   }
958 
959   /** Free OpenMM memory for the Context, Integrator and System. */
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    * Creates and manage an OpenMM Context.
973    *
974    * <p>A Context stores the complete state of a simulation. More specifically, it includes: The
975    * current time The position of each particle The velocity of each particle The values of
976    * configurable parameters defined by Force objects in the System
977    *
978    * <p>You can retrieve a snapshot of the current state at any time by calling getState(). This
979    * allows you to record the state of the simulation at various points, either for analysis or for
980    * checkpointing. getState() can also be used to retrieve the current forces on each particle and
981    * the current energy of the System.
982    */
983   public class Context {
984 
985     /** Requested Platform (i.e. Java or an OpenMM platform). */
986     private final Platform platform;
987     /** Instance of the OpenMM Integrator class. */
988     private final Integrator integrator;
989     /** Constraint tolerance as a fraction of the constrained bond length. */
990     private final double constraintTolerance = ForceFieldEnergy.DEFAULT_CONSTRAINT_TOLERANCE;
991     /** OpenMM Context pointer. */
992     private PointerByReference contextPointer = null;
993     /** Integrator string (default = VERLET). */
994     private String integratorString = "VERLET";
995     /** Time step (default = 0.001 psec). */
996     private double timeStep = 0.001;
997     /** OpenMM Platform pointer. */
998     private PointerByReference platformPointer = null;
999     /** Temperature (default = 298.15). */
1000     private double temperature = 298.15;
1001 
1002     /**
1003      * Create an OpenMM Context.
1004      *
1005      * @param forceField ForceField to used to create an integrator.
1006      * @param requestedPlatform Platform requested.
1007      */
1008     Context(ForceField forceField, Platform requestedPlatform) {
1009       loadPlatform(requestedPlatform);
1010       platform = requestedPlatform;
1011       integrator = new Integrator(forceField, constraintTolerance);
1012     }
1013 
1014     /**
1015      * Get a Pointer to the OpenMM Context. A Context is created if none has been instantiated yet.
1016      *
1017      * @return Context pointer.
1018      */
1019     public PointerByReference getContextPointer() {
1020       if (contextPointer == null) {
1021         create(integratorString, timeStep, temperature, true);
1022       }
1023       return contextPointer;
1024     }
1025 
1026     /**
1027      * Get a Pointer to the OpenMM Integrator.
1028      *
1029      * @return Integrator pointer.
1030      */
1031     public PointerByReference getIntegrator() {
1032       return integrator.getIntegratorPointer();
1033     }
1034 
1035     /**
1036      * Use the Context / Integrator combination to take the requested number of steps.
1037      *
1038      * @param numSteps Number of steps to take.
1039      */
1040     public void integrate(int numSteps) {
1041       OpenMM_Integrator_step(integrator.getIntegratorPointer(), numSteps);
1042     }
1043 
1044     /**
1045      * Use the Context to optimize the system to the requested tolerance.
1046      *
1047      * @param eps Convergence criteria (kcal/mole/A).
1048      * @param maxIterations Maximum number of iterations.
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      * The array x contains atomic coordinates only for active atoms.
1057      *
1058      * @param x Atomic coordinate array for only active atoms.
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           // OpenMM requires coordinates for even "inactive" atoms with mass of zero.
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      * The array v contains velocity values for active atomic coordinates.
1091      *
1092      * @param v Velocity array for active atoms.
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           // OpenMM requires velocities for even "inactive" atoms with mass of zero.
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     /** Set the periodic box vectors for a context based on the crystal instance. */
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     /** {@inheritDoc} */
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     /** Free OpenMM memory for the current Context and Integrator. */
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      * Construct a new Context in which to run a simulation.
1169      *
1170      * @param integratorString Requested integrator.
1171      * @param timeStep Time step (psec).
1172      * @param temperature Temperature (K).
1173      * @param forceCreation Force creation of a new context, even if the current one matches.
1174      * @return Pointer to the created OpenMM context.
1175      */
1176     Context create(String integratorString, double timeStep, double temperature,
1177         boolean forceCreation) {
1178       // Check if the current context is consistent with the requested context.
1179       if (contextPointer != null && !forceCreation) {
1180         if (this.temperature == temperature && this.timeStep == timeStep
1181             && this.integratorString.equalsIgnoreCase(integratorString)) {
1182           // All requested features agree.
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       // Set lambda to 1.0 when creating a context to avoid OpenMM compiling out any terms.
1204       double currentLambda = getLambda();
1205 
1206       if (lambdaTerm) {
1207         ForceFieldEnergyOpenMM.this.setLambda(1.0);
1208       }
1209 
1210       // Create a context.
1211       contextPointer = OpenMM_Context_create_2(system.getSystem(), integratorPointer,
1212           platformPointer);
1213 
1214       // Revert to the current lambda value.
1215       if (lambdaTerm) {
1216         ForceFieldEnergyOpenMM.this.setLambda(currentLambda);
1217       }
1218 
1219       // Get initial positions and velocities for active atoms.
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           // X-axis
1229           x[index] = a.getX();
1230           v[index++] = vel3[0];
1231           // Y-axis
1232           x[index] = a.getY();
1233           v[index++] = vel3[1];
1234           // Z-axis
1235           x[index] = a.getZ();
1236           v[index++] = vel3[2];
1237         }
1238       }
1239 
1240       // Load the current periodic box vectors.
1241       setPeriodicBoxVectors();
1242 
1243       // Load current atomic positions.
1244       setOpenMMPositions(x);
1245 
1246       // Load current velocities.
1247       setOpenMMVelocities(v);
1248 
1249       // Apply constraints starting from current atomic positions.
1250       OpenMM_Context_applyConstraints(contextPointer, constraintTolerance);
1251 
1252       // Application of constraints can change coordinates and velocities.
1253       // Retrieve them for consistency.
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      * Reinitialize the context.
1264      *
1265      * <p>When a Context is created, it may cache information about the System being simulated and
1266      * the Force objects contained in it. This means that, if the System or Forces are then modified,
1267      * the Context might not see all changes. Call reinitialize() to force the Context to rebuild its
1268      * internal representation of the System and pick up any changes that have been made.
1269      *
1270      * <p>This is an expensive operation, so you should try to avoid calling it too frequently.
1271      */
1272     void reinitContext() {
1273       if (contextPointer != null) {
1274         int preserveState = 1;
1275         OpenMM_Context_reinitialize(contextPointer, preserveState);
1276       }
1277     }
1278 
1279     /**
1280      * Set the value of an adjustable parameter defined by a Force object in the System.
1281      *
1282      * @param name the name of the parameter to set.
1283      * @param value the value of the parameter.
1284      */
1285     void setParameter(String name, double value) {
1286       if (contextPointer != null) {
1287         OpenMM_Context_setParameter(contextPointer, name, value);
1288       }
1289     }
1290 
1291     /** Load an OpenMM Platform */
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       // Print out the OpenMM Version.
1299       Pointer version = OpenMM_Platform_getOpenMMVersion();
1300       logger.log(Level.INFO, " Version: {0}", version.getString(0));
1301 
1302       // Print out the OpenMM lib directory.
1303       logger.log(Level.FINE, " Lib Directory:       {0}", OpenMMUtils.getLibDirectory());
1304       // Load platforms and print out their names.
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       // Print out the OpenMM plugin directory.
1316       logger.log(Level.INFO, "\n Plugin Directory:  {0}", OpenMMUtils.getPluginDirectory());
1317       // Load plugins and print out their names.
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       // Extra logging to print out plugins that failed to load.
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    * Create and manage an OpenMM Integrator.
1397    *
1398    * <p>An Integrator defines a method for simulating a System by integrating the equations of
1399    * motion.
1400    *
1401    * <p>Each Integrator object is bound to a particular Context which it integrates. This connection
1402    * is specified by passing the Integrator as an argument to the constructor of the Context.
1403    */
1404   private class Integrator {
1405 
1406     /** Constraint tolerance as a fraction of the constrained bond length. */
1407     private final double constraintTolerance;
1408     /** Langevin friction coefficient. */
1409     private final double frictionCoeff;
1410     /** OpenMM Integrator pointer. */
1411     private PointerByReference integratorPointer = null;
1412 
1413     /**
1414      * Create an Integrator instance.
1415      *
1416      * @param forceField the ForceField instance containing integrator parameters.
1417      * @param constraintTolerance The integrator constraint tolerance.
1418      */
1419     Integrator(ForceField forceField, double constraintTolerance) {
1420       this.constraintTolerance = constraintTolerance;
1421       frictionCoeff = forceField.getDouble("FRICTION_COEFF", 91.0);
1422     }
1423 
1424     /**
1425      * Return a reference to the integrator.
1426      *
1427      * @return Integrator reference.
1428      */
1429     public PointerByReference getIntegratorPointer() {
1430       return integratorPointer;
1431     }
1432 
1433     /**
1434      * Create a integrator.
1435      *
1436      * @param integratorString Name of the integrator to use.
1437      * @param timeStep Time step (psec).
1438      * @param temperature Target temperature (kelvin).
1439      * @return Integrator reference.
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      * Create a Langevin integrator.
1456      *
1457      * @param dt Time step (psec).
1458      * @param temperature Temperature (K).
1459      * @param frictionCoeff Frictional coefficient.
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      * Create a Custom MTS Integrator.
1480      *
1481      * @param dt The outer time step (psec).
1482      */
1483     private void createCustomMTSIntegrator(double dt) {
1484       createCustomIntegrator(dt);
1485 
1486       int n = 4;
1487       // Force group 1 contains slowly varying forces.
1488       // Force group 0 contains the fast varying forces.
1489       int[] forceGroups = {1, 0};
1490       // There will be 1 force evaluation per outer step, and 4 per inner step.
1491       int[] subSteps = {1, 4};
1492       if (system.amoebaCavitationForce != null) {
1493         n = 8;
1494         // Force group 2 contains the cavitation force.
1495         // Force group 1 contains slowly varying forces.
1496         // Force group 0 contains the fast varying forces.
1497         forceGroups = new int[] {2, 1, 0};
1498         // There will be 1 force evaluation per outer step.
1499         // There will be 2 force evaluations per middle step.
1500         // There will be 8 force evaluations per inner step.
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      * Create substeps for the MTS CustomIntegrator.
1517      *
1518      * @param parentSubsteps The number of substeps for the previous force group.
1519      * @param forceGroups The force groups to be evaluated.
1520      * @param subSteps The number of substeps for each force group.
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      * Create a Custom MTS Langevin integrator.
1554      *
1555      * @param dt The outer time step (psec).
1556      * @param temperature The target temperature (K).
1557      * @param frictionCoeff The friction coefficient (1/psec).
1558      */
1559     private void createCustomMTSLangevinIntegrator(double dt, double temperature,
1560         double frictionCoeff) {
1561       createCustomIntegrator(dt);
1562 
1563       int n = 4;
1564       // Force group 1 contains slowly varying forces.
1565       // Force group 0 contains the fast varying forces.
1566       int[] forceGroups = {1, 0};
1567       // There will be 1 force evaluation per outer step, and 4 per inner step.
1568       int[] subSteps = {1, 4};
1569       if (system.amoebaCavitationForce != null) {
1570         n = 8;
1571         // Force group 2 contains the cavitation force.
1572         // Force group 1 contains slowly varying forces.
1573         // Force group 0 contains the fast varying forces.
1574         forceGroups = new int[] {2, 1, 0};
1575         // There will be 1 force evaluation per outer step.
1576         // There will be 2 force evaluations per middle step.
1577         // There will be 8 force evaluations per inner step.
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       // Log the substeps.
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      * Create substeps for the MTS Langevin CustomIntegrator.
1604      *
1605      * @param parentSubsteps The number of substeps for the previous force group.
1606      * @param forceGroups The force groups to be evaluated.
1607      * @param subSteps The number of substeps for each force group.
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         // String step;
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      * Create a Verlet integrator.
1664      *
1665      * @param dt Time step (psec).
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      * Create a Custom integrator.
1678      *
1679      * @param dt Time step (psec).
1680      */
1681     private void createCustomIntegrator(double dt) {
1682       free();
1683       integratorPointer = OpenMM_CustomIntegrator_create(dt);
1684       OpenMM_Integrator_setConstraintTolerance(integratorPointer, constraintTolerance);
1685     }
1686 
1687     /** Destroy the integrator instance. */
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    * Create and manage an OpenMM System.
1700    *
1701    * <p>The definition of a System involves four elements:
1702    *
1703    * <p>The particles and constraints are defined directly by the System object, while forces are
1704    * defined by objects that extend the Force class. After creating a System, call addParticle() once
1705    * for each particle, addConstraint() for each constraint, and addForce() for each Force.
1706    *
1707    * <p>In addition, particles may be designated as "virtual sites". These are particles whose
1708    * positions are computed automatically based on the positions of other particles. To define a
1709    * virtual site, call setVirtualSite(), passing in a VirtualSite object that defines the rules for
1710    * computing its position.
1711    */
1712   public class System {
1713 
1714     private static final double DEFAULT_MELD_SCALE_FACTOR = -1.0;
1715     private final double meldScaleFactor;
1716     /** Andersen thermostat collision frequency. */
1717     private final double collisionFreq;
1718     /**
1719      * When using MELD, our goal will be to scale down the potential by this factor. A negative value
1720      * indicates we're not using MELD.
1721      */
1722     private final boolean useMeld;
1723     /** The Force Field in use. */
1724     ForceField forceField;
1725     /** Array of atoms in the sytem. */
1726     Atom[] atoms;
1727     /** OpenMM System. */
1728     private PointerByReference system;
1729     /** Barostat to be added if NPT (isothermal-isobaric) dynamics is requested. */
1730     private PointerByReference ommBarostat = null;
1731     /** OpenMM center-of-mass motion remover. */
1732     private PointerByReference commRemover = null;
1733     /**
1734      * This flag indicates bonded force constants and equilibria are updated (e.g. during ManyBody
1735      * titration).
1736      */
1737     private boolean updateBondedTerms = false;
1738     /**
1739      * If true, all torsions are treated as 6-fold, and all angles are treated as possibly changing
1740      * between normal and in-plane types.
1741      */
1742     private final boolean manyBodyTitration;
1743     /** OpenMM Custom Bond Force */
1744     private PointerByReference bondForce = null;
1745     /** OpenMM Custom Angle Force */
1746     private PointerByReference angleForce = null;
1747     /** OpenMM Custom Stretch-Bend Force */
1748     private PointerByReference stretchBendForce = null;
1749     /** OpenMM Custom In-Plane Angle Force */
1750     private PointerByReference inPlaneAngleForce = null;
1751     /** OpenMM Custom Urey-Bradley Force */
1752     private PointerByReference ureyBradleyForce = null;
1753     /** OpenMM Custom Out-of-Plane Bend Force */
1754     private PointerByReference outOfPlaneBendForce = null;
1755     /** OpenMM Custom Pi-Torsion Force */
1756     private PointerByReference piTorsionForce = null;
1757     /** OpenMM AMOEBA Torsion Force. */
1758     private PointerByReference torsionForce = null;
1759     private PointerByReference[] restraintTorsions = null;
1760     /** OpenMM Improper Torsion Force. */
1761     private PointerByReference improperTorsionForce = null;
1762     /** OpenMM AMOEBA van der Waals Force. */
1763     private PointerByReference amoebaVDWForce = null;
1764     /** OpenMM AMOEBA Multipole Force. */
1765     private PointerByReference amoebaMultipoleForce = null;
1766     /** OpenMM Generalized Kirkwood Force. */
1767     private PointerByReference amoebaGeneralizedKirkwoodForce = null;
1768     /** OpenMM AMOEBA WCA Dispersion Force. */
1769     private PointerByReference amoebaWcaDispersionForce = null;
1770     /** OpenMM AMOEBA WCA Cavitation Force. */
1771     private PointerByReference amoebaCavitationForce = null;
1772     /** OpenMM Custom GB Force. */
1773     private PointerByReference customGBForce = null;
1774     /** OpenMM Fixed Charge Non-Bonded Force. */
1775     private PointerByReference fixedChargeNonBondedForce = null;
1776     /** Fixed charge softcore vdW force boolean. */
1777     private boolean softcoreCreated = false;
1778     /** Boolean array, holds charge exclusion list. */
1779     private boolean[] chargeExclusion;
1780     /** Boolean array, holds van Der Waals exclusion list. */
1781     private boolean[] vdWExclusion;
1782     /** Double array, holds charge quantity value for exceptions. */
1783     private double[] exceptionChargeProd;
1784     /** Double array, holds epsilon quantity value for exceptions. */
1785     private double[] exceptionEps;
1786     /**
1787      * A map from vdW class values to OpenMM vdW types.
1788      */
1789     private Map<Integer, Integer> vdwClassToOpenMMType;
1790     /**
1791      * A class for a special vdW type that specifies zero energy (eps = 0.0; sigma = 1.0) for use
1792      * with the FFX "use" flag (e.g. use = false should give zero vdW energy for a many-body
1793      * expansion).
1794      */
1795     private int vdWClassForNoInteraction;
1796     /**
1797      * Lambda flag to indicate control of electrostatic scaling. If both elec and vdW are being
1798      * scaled, then vdW is scaled first, followed by elec.
1799      */
1800     private boolean elecLambdaTerm;
1801     /**
1802      * Lambda flag to indicate control of vdW scaling. If both elec and vdW are being scaled, then
1803      * vdW is scaled first, followed by elec.
1804      */
1805     private boolean vdwLambdaTerm;
1806     /**
1807      * Lambda flag to indicate control of torsional force constants (L=0 corresponds to torsions
1808      * being off, and L=1 to torsions at full strength).
1809      */
1810     private boolean torsionLambdaTerm;
1811     /** Value of the van der Waals lambda state variable. */
1812     private double lambdaVDW = 1.0;
1813     /** Value of the electrostatics lambda state variable. */
1814     private double lambdaElec = 1.0;
1815     /** Value of the electrostatics lambda state variable. */
1816     private double lambdaTorsion = 1.0;
1817     /**
1818      * The lambda value that defines when the electrostatics will start to turn on for full path
1819      * non-bonded term scaling.
1820      *
1821      * <p>A value of 0.6 works well for Chloride ion solvation, which is a difficult case due to the
1822      * ion having a formal negative charge and a large polarizability.
1823      */
1824     private double electrostaticStart = 0.6;
1825     /** Electrostatics lambda is raised to this power. */
1826     private double electrostaticLambdaPower;
1827     /** van der Waals softcore alpha. */
1828     private double vdWSoftcoreAlpha = 0.25;
1829     /** OpenMM thermostat. Currently, an Andersen thermostat is supported. */
1830     private PointerByReference ommThermostat = null;
1831     /** van der Waals softcore beta. */
1832     private double vdwSoftcorePower = 3.0;
1833     /** Torsional lambda power. */
1834     private double torsionalLambdaPower = 2.0;
1835 
1836     /**
1837      * OpenMMSystem constructor.
1838      *
1839      * @param molecularAssembly MolecularAssembly used to construct the OpenMM system.
1840      */
1841     System(MolecularAssembly molecularAssembly) {
1842       // Create the OpenMM System
1843       system = OpenMM_System_create();
1844       logger.info("\n System created");
1845 
1846       forceField = molecularAssembly.getForceField();
1847       atoms = molecularAssembly.getAtomArray();
1848 
1849       // Load atoms.
1850       try {
1851         addAtoms();
1852       } catch (Exception e) {
1853         logger.severe(" Atom without mass encountered.");
1854       }
1855 
1856       // Check for MELD use. If we're using MELD, set all lambda terms to true.
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       // Read alchemical information -- this needs to be done before creating forces.
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       // Expand the path [lambda-start .. 1.0] to the interval [0.0 .. 1.0].
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         // lambda path starts at 0.0
1907         lambdaStart = 0.0;
1908         // electrostaticStart is ignored for MELD.
1909         electrostaticStart = 0.0;
1910         // electrostaticLambdaPower is ignored for MELD.
1911         electrostaticLambdaPower = 1.0;
1912         // vdW is linearly scaled for MELD.
1913         vdwSoftcorePower = 1;
1914         // No softcore offset for MELD.
1915         vdWSoftcoreAlpha = 0.0;
1916         // Torsions are linearly scaled for MELD.
1917         torsionalLambdaPower = 1.0;
1918         // Only need single-sided dU/dL
1919         twoSidedFiniteDifference = false;
1920       }
1921 
1922       collisionFreq = forceField.getDouble("COLLISION_FREQ", 0.1);
1923 
1924       // Set up rigid constraints. These flags need to be set before bonds and angles are created
1925       // below.
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       // Add Bond Force.
1942       if (rigidBonds) {
1943         logger.info(" Not creating AmoebaBondForce because bonds are constrained.");
1944       } else {
1945         addBondForce();
1946       }
1947 
1948       // Add Angle Force.
1949       addAngleForce();
1950       addInPlaneAngleForce();
1951 
1952       // Add Stretch-Bend Force.
1953       addStretchBendForce();
1954 
1955       // Add Urey-Bradley Force.
1956       addUreyBradleyForce();
1957 
1958       // Out-of Plane Bend Force.
1959       addOutOfPlaneBendForce();
1960 
1961       // Add Torsion Force.
1962       addTorsionForce();
1963 
1964       // Add Improper Torsion Force.
1965       addImproperTorsionForce();
1966 
1967       // Add Pi-Torsion Force.
1968       addPiTorsionForce();
1969 
1970       // Add Torsion-Torsion Force.
1971       addTorsionTorsionForce();
1972 
1973       // Add coordinate restraints.
1974       addHarmonicRestraintForce();
1975 
1976       // Add bond restraints.
1977       addRestraintBondForce();
1978 
1979       // Add Restrain Groups.
1980       addRestrainGroupsForce();
1981 
1982       // Add stretch-torsion coupling terms.
1983       addStretchTorsionForce();
1984 
1985       // Add angle-torsion coupling terms.
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           // Add vdW Force.
1999           addAmoebaVDWForce();
2000 
2001           // Add Multipole Force.
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      * Add an Andersen thermostat to the system.
2032      *
2033      * @param targetTemp Target temperature in Kelvins.
2034      */
2035     public void addAndersenThermostatForce(double targetTemp) {
2036       addAndersenThermostatForce(targetTemp, collisionFreq);
2037     }
2038 
2039     /**
2040      * Add an Andersen thermostat to the system.
2041      *
2042      * @param targetTemp Target temperature in Kelvins.
2043      * @param collisionFreq Collision frequency in 1/psec.
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     /** Adds a force that removes center-of-mass motion. */
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      * Add a Monte Carlo Barostat to the system.
2074      *
2075      * @param targetPressure The target pressure (in atm).
2076      * @param targetTemp The target temperature.
2077      * @param frequency The frequency to apply the barostat.
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      * Calculate the number of degrees of freedom.
2104      *
2105      * @return Number of degrees of freedom.
2106      */
2107     public int calculateDegreesOfFreedom() {
2108       // Begin from the 3 times the number of active atoms.
2109       int dof = getNumberOfVariables();
2110       // Remove OpenMM constraints.
2111       dof = dof - OpenMM_System_getNumConstraints(system);
2112       // Remove center of mass motion.
2113       if (commRemover != null) {
2114         dof -= 3;
2115       }
2116       return dof;
2117     }
2118 
2119     /** Destroy the system. */
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     /** Print current lambda values. */
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      * Set the overall lambda value for the system.
2137      *
2138      * @param lambda Current lambda value.
2139      */
2140     public void setLambda(double lambda) {
2141 
2142       // Initially set all lambda values to 1.0.
2143       lambdaTorsion = 1.0;
2144 
2145       // Applied to softcore vdW forces.
2146       lambdaVDW = 1.0;
2147 
2148       // Applied to normal electrostatic parameters for alchemical atoms.
2149       lambdaElec = 1.0;
2150 
2151       if (torsionLambdaTerm) {
2152         // Multiply torsional potentials by L^2 (dU/dL = 0 at L=0).
2153         lambdaTorsion = pow(lambda, torsionalLambdaPower);
2154         if (useMeld) {
2155           lambdaTorsion = meldScaleFactor + lambda * (1.0 - meldScaleFactor);
2156         }
2157       }
2158 
2159       if (elecLambdaTerm && vdwLambdaTerm) {
2160         // Lambda effects both vdW and electrostatics.
2161         if (lambda < electrostaticStart) {
2162           // Begin turning vdW on with electrostatics off.
2163           lambdaElec = 0.0;
2164         } else {
2165           // Turn electrostatics on during the latter part of the path.
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         // Lambda effects vdW, with electrostatics turned off.
2177         lambdaElec = 0.0;
2178         lambdaVDW = lambda;
2179         if (useMeld) {
2180           lambdaVDW = meldScaleFactor + lambda * (1.0 - meldScaleFactor);
2181         }
2182 
2183       } else if (elecLambdaTerm) {
2184         // Lambda effects electrostatics, but not vdW.
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      * Return a reference to the System.
2198      *
2199      * @return System referenece.
2200      */
2201     PointerByReference getSystem() {
2202       return system;
2203     }
2204 
2205     /**
2206      * Set the default values of the vectors defining the axes of the periodic box (measured in nm).
2207      *
2208      * <p>Any newly created Context will have its box vectors set to these. They will affect any
2209      * Force added to the System that uses periodic boundary conditions.
2210      *
2211      * <p>Triclinic boxes are supported, but the vectors must satisfy certain requirements. In
2212      * particular, a must point in the x direction, b must point "mostly" in the y direction, and c
2213      * must point "mostly" in the z direction. See the documentation for details.
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      * Update parameters if the Use flags and/or Lambda value has changed.
2237      *
2238      * @param atoms Atoms in this list are considered.
2239      */
2240     void updateParameters(Atom[] atoms) {
2241 
2242       if (vdwLambdaTerm) {
2243         if (fixedChargeNonBondedForce != null) {
2244           if (!softcoreCreated) {
2245             addCustomNonbondedSoftcoreForce();
2246             // Re-initialize the context.
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             // Avoid any updateParametersInContext calls if vdwLambdaTerm is true, but not other
2255             // alchemical terms.
2256             if (!torsionLambdaTerm && !elecLambdaTerm) {
2257               return;
2258             }
2259           } else {
2260             softcoreCreated = true;
2261           }
2262         }
2263       }
2264 
2265       // Note Stretch-Torsion and Angle-Torsion terms (for nucleic acids)
2266       // and Torsion-Torsion terms (for protein backbones) are not udpated yet.
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       // Update fixed charge non-bonded parameters.
2310       if (fixedChargeNonBondedForce != null) {
2311         updateFixedChargeNonBondedForce(atoms);
2312       }
2313 
2314       // Update fixed charge GB parameters.
2315       if (customGBForce != null) {
2316         updateCustomGBForce(atoms);
2317       }
2318 
2319       // Update AMOEBA vdW parameters.
2320       if (amoebaVDWForce != null) {
2321         updateAmoebaVDWForce(atoms);
2322       }
2323 
2324       // Update AMOEBA polarizable multipole parameters.
2325       if (amoebaMultipoleForce != null) {
2326         updateAmoebaMultipoleForce(atoms);
2327       }
2328 
2329       // Update GK force.
2330       if (amoebaGeneralizedKirkwoodForce != null) {
2331         updateGeneralizedKirkwoodForce(atoms);
2332       }
2333 
2334       // Update WCA Force.
2335       if (amoebaWcaDispersionForce != null) {
2336         updateWCAForce(atoms);
2337       }
2338 
2339       // Update WCA Force.
2340       if (amoebaCavitationForce != null) {
2341         updateCavitationForce(atoms);
2342       }
2343     }
2344 
2345     /**
2346      * Adds atoms from the molecular assembly to the OpenMM System and reports to the user the number
2347      * of particles added.
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     /** This method sets the mass of inactive atoms to zero. */
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     /** Add a bond force to the OpenMM System. */
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     /** Update an existing bond force for the OpenMM System. */
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     /** Add an angle force to the OpenMM System. */
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           // Skip In-Plane angles unless this is ManyBody Titration.
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             // This is a place-holder Angle, in case the In-Plane Angle is swtiched to a
2488             // Normal Angle during in the udpateAngleForce.
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     /** Update the angle force. */
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           // Skip In-Plane angles unless this is ManyBody Titration.
2519         }
2520         // Update angles that do not involve rigid hydrogen atoms.
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             // Zero the force constant for In-Plane Angles.
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     /** Add an in-plane angle force to the OpenMM System. */
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           // Skip Normal angles unless this is ManyBody Titration.
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             // This is a place-holder Angle, in case the Normal Angle is switched to an
2595             // In-Plane Angle during in the udpateInPlaneAngleForce.
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     /** Update the in-plane angle force. */
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           // Skip Normal angles unless this is ManyBody Titration.
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           // There is no 4th atom for normal angles, so set the index to first atom.
2648           int i4 = 0;
2649           if (angleMode == AngleMode.NORMAL) {
2650             // Zero the force constant for Normal Angles.
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     /** Add the Urey-Bradley force to the OpenMM System. */
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         // The implementation of UreyBradley in FFX & Tinker: k x^2
2700         // The implementation of Harmonic Bond Force in OpenMM:  k x^2 / 2
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     /** Update the Urey-Bradley force. */
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         // The implementation of UreyBradley in FFX & Tinker: k x^2
2729         // The implementation of Harmonic Bond Force in OpenMM:  k x^2 / 2
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     /** Add an out-of-plane bend force to the OpenMM System. */
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     /** Update the Out-of-Plane bend force. */
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     /** Add a stretch-bend force to the OpenMM System. */
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     /** Update the Stretch-Bend force. */
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     /** Add a torsion force to the OpenMM System. */
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         // Enforce 6-fold torsions since TorsionType instances can have different lengths
2941         // when side-chain protonation changes.
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     /** Update the Torsion force. */
2957     private void updateTorsionForce() {
2958       // Check if this system has torsions.
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         // Enforce 6-fold torsions since TorsionType instances can have different lengths
2979         // when side-chain protonation changes.
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     /** Add an improper-torsion force to the OpenMM System. */
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     /** Update the Improper Torsion force. */
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     /** Add a Pi-Torsion force to the OpenMM System. */
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     /** Update the Pi-Torsion force. */
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     /** Add a Torsion-Torsion force to the OpenMM System. */
3145     private void addTorsionTorsionForce() {
3146       TorsionTorsion[] torsionTorsions = getTorsionTorsions();
3147       if (torsionTorsions == null || torsionTorsions.length < 1) {
3148         return;
3149       }
3150 
3151       // Load the torsion-torsions.
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         // Check if the TorTor parameters have already been added to the Hash.
3166         int gridIndex = 0;
3167         if (torTorTypes.containsKey(key)) {
3168 
3169           // If the TorTor has been added, get its (ordered) index in the Hash.
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           // Add the new TorTor.
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       // Load the Torsion-Torsion parameters.
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         // Create the 3D grid.
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     /** Adds stretch-torsion couplings (as defined for the 2017 AMOEBA nucleic acid model). */
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     /** Adds stretch-torsion couplings (as defined for the 2017 AMOEBA nucleic acid model). */
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         // One thing that concerns me is whether it's correct to get angle[0] instead of angle[num
3342         // hydrogens].
3343         // This is the way it is in FFX, but that may be a bug.
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       // update restraint torsions ONLY here.
3398       // Only update parameters if torsions are being scaled by lambda.
3399 
3400       // Check if this system has torsions.
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     /** Uses arithmetic mean to define sigma and geometric mean for epsilon. */
3429     private void addFixedChargeNonBondedForce() {
3430       VanDerWaals vdW = getVdwNode();
3431       if (vdW == null) {
3432         return;
3433       }
3434 
3435       /*
3436        Only 6-12 LJ with arithmetic mean to define sigma and geometric mean
3437        for epsilon is supported.
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       // OpenMM vdW force requires a diameter (i.e. not radius).
3452       double radScale = 1.0;
3453       if (vdwForm.radiusSize == RADIUS) {
3454         radScale = 2.0;
3455       }
3456 
3457       // OpenMM vdw force requires atomic sigma values (i.e. not r-min).
3458       if (vdwForm.radiusType == R_MIN) {
3459         radScale /= 1.122462048309372981;
3460       }
3461 
3462       // Add particles.
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       // Define 1-4 scale factors.
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           // Units of the Ewald coefficient are A^-1; Multiply by AngstromsPerNM to convert to
3536           // (Nm^-1).
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      * Updates the fixed-charge non-bonded force for changes in Use flags or Lambda.
3589      *
3590      * @param atoms Array of atoms to update.
3591      */
3592     private void updateFixedChargeNonBondedForce(Atom[] atoms) {
3593       VanDerWaals vdW = getVdwNode();
3594       // Only 6-12 LJ with arithmetic mean to define sigma and geometric mean for epsilon is
3595       // supported.
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       // OpenMM vdW force requires a diameter (i.e. not radius).
3604       double radScale = 1.0;
3605       if (vdwForm.radiusSize == RADIUS) {
3606         radScale = 2.0;
3607       }
3608 
3609       // OpenMM vdw force requires atomic sigma values (i.e. not r-min).
3610       if (vdwForm.radiusType == R_MIN) {
3611         radScale /= 1.122462048309372981;
3612       }
3613 
3614       // Update parameters.
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           // If we're using vdwLambdaTerm, this atom's vdW interactions are handled by the Custom
3631           // Non-Bonded force.
3632           if (vdwLambdaTerm) {
3633             eps = 0.0;
3634           }
3635           // Always scale the charge by lambdaElec
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       // Update Exceptions.
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         // Only update exceptions.
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         Note that the minimum epsilon value cannot be zero, or OpenMM may
3677         report an error that the number of Exceptions has changed.
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      * 1. Handle interactions between non-alchemical atoms with our default OpenMM NonBondedForce.
3714      * Note that alchemical atoms must have eps=0 to turn them off in this force.
3715      * <p>
3716      * 2. Handle interactions between alchemical atoms and mixed non-alchemical <-> alchemical
3717      * interactions with an OpenMM CustomNonBondedForce.
3718      */
3719     private void addCustomNonbondedSoftcoreForce() {
3720       VanDerWaals vdW = getVdwNode();
3721       if (vdW == null) {
3722         return;
3723       }
3724 
3725       /*
3726        Only 6-12 LJ with arithmetic mean to define sigma and geometric mean
3727        for epsilon is supported.
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       // Sterics mixing rules.
3740       String stericsMixingRules = " epsilon = sqrt(epsilon1*epsilon2);";
3741       stericsMixingRules += " rmin = 0.5 * (sigma1 + sigma2) * 1.122462048309372981;";
3742 
3743       // Softcore Lennard-Jones, with a form equivalent to that used in FFX VanDerWaals class.
3744       String stericsEnergyExpression = "(vdw_lambda^beta)*epsilon*x*(x-2.0);";
3745       // Effective softcore distance for sterics.
3746       stericsEnergyExpression += " x = 1.0 / (alpha*(1.0-vdw_lambda)^2.0 + (r/rmin)^6.0);";
3747       // Define energy expression for sterics.
3748       String energyExpression = stericsEnergyExpression + stericsMixingRules;
3749 
3750       PointerByReference fixedChargeSoftcore = OpenMM_CustomNonbondedForce_create(energyExpression);
3751 
3752       // Get the Alpha and Beta constants from the VanDerWaals instance.
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       // Add particles.
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         // Handle cases where sigma is 0.0; for example Amber99 tyrosine hydrogen atoms.
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       // Add energy parameter derivative
3835       // OpenMM_CustomNonbondedForce_addEnergyParameterDerivative(fixedChargeSoftcore,
3836       // "vdw_lambda");
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       // Alchemical with Alchemical could be either softcore or normal interactions (softcore here).
3844       PointerByReference alchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(
3845           stericsEnergyExpression);
3846 
3847       // Non-Alchemical with Alchemical is essentially always softcore.
3848       PointerByReference nonAlchemicalAlchemicalStericsForce = OpenMM_CustomBondForce_create(
3849           stericsEnergyExpression);
3850 
3851       // Currently both are treated the same (so we could condense the code below).
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         // Omit both Exclusions (1-2, 1-3) and Exceptions (scaled 1-4) from the
3876         // CustomNonbondedForce.
3877         OpenMM_CustomNonbondedForce_addExclusion(fixedChargeSoftcore, atomi.getValue(),
3878             atomj.getValue());
3879 
3880         // Deal with scaled 1-4 torsions using the CustomBondForce
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       // OpenMM_CustomBondForce_addEnergyParameterDerivative(alchemicalAlchemicalStericsForce,
3928       // "vdw_lambda");
3929       // OpenMM_CustomBondForce_addEnergyParameterDerivative(nonAlchemicalAlchemicalStericsForce,
3930       // "vdw_lambda");
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     /** Add a custom GB force to the OpenMM System. */
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; // 100 square Angstroms per square nanometer.
3952         // logger.info(String.format(" FFX surface tension: %9.5g kcal/mol/Ang^2", sTens));
3953         // logger.info(String.format(" OpenMM surface tension: %9.5g kJ/mol/nm^2", sTens));
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); // Factor of 0.1 for Ang to nm.
3967       OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "probeRadius",
3968           gk.getProbeRadius() * OpenMM_NmPerAngstrom);
3969 
3970       OpenMM_CustomGBForce_addComputedValue(customGBForce, "I",
3971           // "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r-sr2*sr2/r)+0.5*log(L/U)/r+C);"
3972           // "step(r+sr2-or1)*0.5*((1/L^3-1/U^3)/3+(1/U^4-1/L^4)/8*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);"
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               // + "C=2*(1/or1-1/L)*step(sr2-r-or1);"
3976               + "C=2/3*(1/or1^3-1/L^3)*step(sr2-r-or1);"
3977               // + "L=step(or1-D)*or1 + (1-step(or1-D))*D;"
3978               // + "D=step(r-sr2)*(r-sr2) + (1-step(r-sr2))*(sr2-r);"
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           // "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
3986           // "psi=I*or; or=radius-0.009"
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       // Particle pair term is the generalized Born cross term.
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      * Updates the custom GB force for changes in Use flags or Lambda.
4027      *
4028      * @param atoms Array of atoms to update.
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; // 100 square Angstroms per square nanometer.
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         // Reset the double array for the next atom.
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     /** Add an AMOEBA van der Waals force to the OpenMM System. */
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         The vdW class used to specify no vdW interactions for an atom will be Zero
4101         if all atom classes are greater than zero. Otherwise:
4102         vdWClassForNoInteraction = min(atomClass) - 1
4103        */
4104       vdWClassForNoInteraction = 0;
4105       // Add vdW parameters to the force and record their type.
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       // Add a special vdW type for zero vdW energy and forces (e.g. to support the FFX "use" flag).
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       // Create exclusion lists.
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      * Updates the AMOEBA van der Waals force for changes in Use flags or Lambda.
4210      *
4211      * @param atoms Array of atoms to update.
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         // Get the OpenMM index for this vdW type.
4230         int type = vdwClassToOpenMMType.get(vdwType.atomClass);
4231         if (!atom.getUse()) {
4232           // Get the OpenMM index for a special vdW type that has no interactions.
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     /** Add an AMOEBA polarizable multipole force to the OpenMM System. */
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            * Citation:
4289            * Simmonett, A. C.;  Pickard, F. C. t.;  Shao, Y.;  Cheatham, T. E., 3rd; Brooks, B. R., Efficient treatment of induced dipoles. The Journal of chemical physics 2015, 143 (7), 074115-074115.
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           // PointerByReference exptCoefficients = OpenMM_DoubleArray_create(1);
4299           // OpenMM_DoubleArray_set(exptCoefficients, 0, 1.044);
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         // Define the frame definition.
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; // Should be 1.0 at this point.
4333         if (!atom.applyLambda()) {
4334           lambdaScale = 1.0;
4335         }
4336 
4337         useFactor *= lambdaScale;
4338 
4339         // Load local multipole coefficients.
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         // Add the multipole.
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         // 1-2 Mask
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         // 1-3 Mask
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         // 1-4 Mask
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         // 1-5 Mask
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         // 1-1 Polarization Groups.
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         // AMOEBA does not scale between 1-2, 1-3, etc. polarization groups.
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      * Updates the Amoeba electrostatic multipolar force for changes in Use flags or Lambda.
4476      *
4477      * @param atoms Array of atoms to update.
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         // Define the frame definition.
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         // Load local multipole coefficients.
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         // Set the multipole parameters.
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     /** Add a Generalized Kirkwood force to the OpenMM System. */
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         // No descreening when using perfect radii (OpenMM will just load the base radii).
4580         perfectRadiiScale = 0.0;
4581       }
4582 
4583       // Turn on tanh rescaling only when not using perfect radii.
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           // Configure a Born Radii based surface area term.
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           // This NonPolar model does not use a Born Radii based surface area term.
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       // Add dispersion
4649       switch (nonpolar) {
4650         case CAV_DISP, GAUSS_DISP, SEV_DISP, BORN_CAV_DISP -> addWCAForce();
4651         default -> {
4652           // WCA force is not being used.
4653         }
4654       }
4655 
4656       // Add cavitation
4657       if (nonpolar == GAUSS_DISP) {
4658         addCavitationForce();
4659       }
4660     }
4661 
4662     /**
4663      * Updates the AMOEBA Generalized Kirkwood force for changes in Use flags or Lambda.
4664      *
4665      * @param atoms Array of atoms to update.
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         // No descreening when using perfect radii (OpenMM will just load the base radii).
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       //        OpenMM Bug: Surface Area is not Updated by "updateParametersInContext"
4719       //
4720       //        NonPolar nonpolar = gk.getNonPolarModel();
4721       //        switch (nonpolar) {
4722       //            case BORN_SOLV:
4723       //            case BORN_CAV_DISP:
4724       //            default:
4725       //                // Configure a Born Radii based surface area term.
4726       //                double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal
4727       //                        * OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm * lambdaElec;
4728       //
4729       // OpenMM_AmoebaGeneralizedKirkwoodForce_setSurfaceAreaFactor(amoebaGeneralizedKirkwoodForce,
4730       // -surfaceTension);
4731       //                break;
4732       //            case CAV:
4733       //            case CAV_DISP:
4734       //            case HYDROPHOBIC_PMF:
4735       //            case NONE:
4736       //                // This NonPolar model does not use a Born Radii based surface area term.
4737       //                break;
4738       //        }
4739 
4740       if (context.contextPointer != null) {
4741         OpenMM_AmoebaGeneralizedKirkwoodForce_updateParametersInContext(
4742             amoebaGeneralizedKirkwoodForce, context.contextPointer);
4743       }
4744     }
4745 
4746     /** Add a nonpolar Weeks-Chandler-Andersen dispersion force to the OpenMM System. */
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      * Updates the WCA force for changes in Use flags or Lambda.
4801      *
4802      * @param atoms Array of atoms to update.
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         // Scale all implicit solvent terms with the square of electrostatics lambda
4820         // (so dUdisp / dL is 0 at lambdaElec = 0).
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     /** Add a GaussVol cavitation force to the OpenMM System. */
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      * Updates the Cavitation force for changes in Use flags or Lambda.
4880      *
4881      * @param atoms Array of atoms to update.
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       // Changing cavitation radii is not supported.
4896       // for (int i=0; i<nAtoms; i++) {
4897       //  gaussVol.updateAtom(i);
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         // Scale all implicit solvent terms with the square of electrostatics lambda
4908         // (so dUcav / dL is 0 at lambdaElec = 0).
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      * Adds harmonic restraints (CoordRestraint objects) to OpenMM as a custom external force.
4934      *
4935      * <p>TODO: Make robust to flat-bottom restraints.
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     /** Adds restraint bonds, if any. */
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       // OpenMM's HarmonicBondForce class uses k, not 1/2*k as does FFX.
4992       double kParameterConversion =
4993           2.0 * OpenMM_KJPerKcal / (OpenMM_NmPerAngstrom * OpenMM_NmPerAngstrom);
4994 
4995       // Map from bond functional forms to the restraint-bonds using that functional form.
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               // Do nothing.
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       // In the expression below, u and l are the upper and lower threshold
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       // Create the Restrain Groups.
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       // Add the restraints between groups.
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       // Add the constraint force.
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     /** Add a constraint to every bond. */
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     /** Add a constraint to every bond that includes a hydrogen atom. */
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     /** Add a constraint to every angle that includes two hydrogen atoms. */
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           // Calculate a "false bond" length between atoms 1 and 3 to constrain the angle using the
5165           // law of cosines.
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           // Equilibrium angle value in degrees.
5173           double angleVal = angle.angleType.angle[angle.nh];
5174 
5175           // Law of cosines.
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      * Check to see if an angle is a hydrogen angle. This method only returns true for hydrogen
5190      * angles that are less than 160 degrees.
5191      *
5192      * @param angle Angle to check.
5193      * @return boolean indicating whether an angle is a hydrogen angle that is less than 160 degrees.
5194      */
5195     private boolean isHydrogenAngle(Angle angle) {
5196       if (angle.containsHydrogen()) {
5197         // Equilibrium angle value in degrees.
5198         double angleVal = angle.angleType.angle[angle.nh];
5199         // Make sure angle is less than 160 degrees because greater than 160 degrees will not be
5200         // constrained
5201         // well using the law of cosines.
5202         if (angleVal < 160.0) {
5203           Atom atom1 = angle.getAtom(0);
5204           Atom atom2 = angle.getAtom(1);
5205           Atom atom3 = angle.getAtom(2);
5206           // Setting constraints only on angles where atom1 or atom3 is a hydrogen while atom2 is
5207           // not a hydrogen.
5208           return atom1.isHydrogen() && atom3.isHydrogen() && !atom2.isHydrogen();
5209         }
5210       }
5211       return false;
5212     }
5213 
5214   }
5215 
5216   /** Retrieve state information from an OpenMM Simulation. */
5217   public class State {
5218 
5219     /** Potential energy (kcal/mol). */
5220     public final double potentialEnergy;
5221     /** Kinetic energy (kcal/mol). */
5222     public final double kineticEnergy;
5223     /** Total energy (kcal/mol). */
5224     public final double totalEnergy;
5225     /** Temperature (K). */
5226     public final double temperature;
5227     /** Pointer to an OpenMM state. */
5228     private final PointerByReference state;
5229     /** If true, positions are available. */
5230     private final boolean positions;
5231     /** If true, velocities are available. */
5232     private final boolean velocities;
5233     /** If true. forces are available. */
5234     private final boolean forces;
5235 
5236     /**
5237      * Construct an OpenMM State with the requested information.
5238      *
5239      * @param positions Retrieve positions.
5240      * @param energies Retrieve energies.
5241      * @param forces Retrieve forces.
5242      * @param velocities Retrieve velocities.
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       // Construct the mask.
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       // Retrieve the state.
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      * The force array contains the OpenMM force information for all atoms. The returned array a
5291      * contains accelerations for active atoms only.
5292      *
5293      * @param a Acceleration components for only active atomic coordinates.
5294      * @return The acceleration for each active atomic coordinate.
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      * The force array contains the OpenMM force information for all atoms. The returned array g
5327      * contains components for active atoms only.
5328      *
5329      * @param g Gradient components for only active atomic coordinates.
5330      * @return g The gradient includes only active atoms
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      * Read the periodic lattice vectors from a state.
5375      *
5376      * <p>The crystal instance will be updated, and passed to the ForceFieldEnergy instance.
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      * The positions array contains the OpenMM atomic position information for all atoms. The
5406      * returned array x contains coordinates only for active atoms.
5407      *
5408      * @param x Atomic coordinates only for active atoms.
5409      * @return x The atomic coordinates for only active atoms.
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      * The positions array contains the OpenMM atomic position information for all atoms. The
5440      * returned array x contains coordinates for active atoms only.
5441      *
5442      * @param v Velocity only for active atomic coordinates.
5443      * @return v The velocity for each active atomic coordinate.
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 }