1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.openmm;
39
40 import ffx.openmm.Force;
41 import ffx.openmm.amoeba.GeneralizedKirkwoodForce;
42 import ffx.potential.bonded.Atom;
43 import ffx.potential.nonbonded.GeneralizedKirkwood;
44 import ffx.potential.parameters.ForceField;
45 import ffx.potential.parameters.MultipoleType;
46
47 import java.util.logging.Level;
48 import java.util.logging.Logger;
49
50 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_AngstromsPerNm;
51 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_KJPerKcal;
52 import static edu.uiowa.jopenmm.OpenMMAmoebaLibrary.OpenMM_NmPerAngstrom;
53 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_False;
54 import static edu.uiowa.jopenmm.OpenMMLibrary.OpenMM_Boolean.OpenMM_True;
55 import static java.lang.String.format;
56
57 public class AmoebaGeneralizedKirkwoodForce extends GeneralizedKirkwoodForce {
58
59 private static final Logger logger = Logger.getLogger(AmoebaGeneralizedKirkwoodForce.class.getName());
60
61 public AmoebaGeneralizedKirkwoodForce(OpenMMEnergy openMMEnergy) {
62 GeneralizedKirkwood gk = openMMEnergy.getGK();
63 if (gk == null) {
64 destroy();
65 return;
66 }
67
68 setSolventDielectric(gk.getSolventPermittivity());
69 double soluteDielectric = gk.getSolutePermittivity();
70 if (soluteDielectric != 1.0) {
71 logger.severe(" Solute dielectric is not 1.0, which is not supported by OpenMM.");
72 }
73 setSoluteDielectric(1.0);
74 setDielectricOffset(gk.getDescreenOffset() * OpenMM_NmPerAngstrom);
75
76
77 boolean usePerfectRadii = gk.getUsePerfectRadii();
78 double perfectRadiiScale = 1.0;
79 if (usePerfectRadii) {
80
81 perfectRadiiScale = 0.0;
82 }
83
84
85 int tanhRescale = 0;
86 if (gk.getTanhCorrection() && !usePerfectRadii) {
87 tanhRescale = 1;
88 }
89 double[] betas = gk.getTanhBetas();
90 setTanhRescaling(tanhRescale);
91 setTanhParameters(betas[0], betas[1], betas[2]);
92
93 double[] baseRadius = gk.getBaseRadii();
94 if (usePerfectRadii) {
95 baseRadius = gk.getPerfectRadii();
96 }
97
98 double[] overlapScale = gk.getOverlapScale();
99 double[] descreenRadius = gk.getDescreenRadii();
100 double[] neckFactor = gk.getNeckScale();
101
102 if (!usePerfectRadii && logger.isLoggable(Level.FINE)) {
103 logger.fine(" GK Base Radii Descreen Radius Overlap Scale Overlap");
104 }
105
106 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
107 int nAtoms = atoms.length;
108 for (int i = 0; i < nAtoms; i++) {
109 MultipoleType multipoleType = atoms[i].getMultipoleType();
110 double base = baseRadius[i] * OpenMM_NmPerAngstrom;
111 double descreen = descreenRadius[i] * OpenMM_NmPerAngstrom * perfectRadiiScale;
112 double overlap = overlapScale[i] * perfectRadiiScale;
113 double neck = neckFactor[i] * perfectRadiiScale;
114 addParticle_1(multipoleType.charge, base, overlap, descreen, neck);
115 if (!usePerfectRadii && logger.isLoggable(Level.FINE)) {
116 logger.fine(format(" %s %8.6f %8.6f %5.3f", atoms[i].toString(), baseRadius[i], descreenRadius[i], overlapScale[i]));
117 }
118 }
119
120 setProbeRadius(gk.getProbeRadius() * OpenMM_NmPerAngstrom);
121
122 GeneralizedKirkwood.NonPolarModel nonpolar = gk.getNonPolarModel();
123 switch (nonpolar) {
124 case BORN_CAV_DISP, BORN_SOLV -> {
125
126 double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal * OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm;
127 setIncludeCavityTerm(OpenMM_True);
128 setSurfaceAreaFactor(-surfaceTension);
129 }
130
131 default -> setIncludeCavityTerm(OpenMM_False);
132 }
133
134 ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
135 int forceGroup = forceField.getInteger("GK_FORCE_GROUP", 1);
136 setForceGroup(forceGroup);
137 logger.log(Level.INFO, format(" Generalized Kirkwood force \t\t%d", forceGroup));
138 }
139
140
141
142
143
144
145
146 public static Force constructForce(OpenMMEnergy openMMEnergy) {
147 GeneralizedKirkwood gk = openMMEnergy.getGK();
148 if (gk == null) {
149 return null;
150 }
151 return new AmoebaGeneralizedKirkwoodForce(openMMEnergy);
152 }
153
154
155
156
157
158
159
160
161 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
162 GeneralizedKirkwood gk = openMMEnergy.getGK();
163 if (gk == null || pointer == null) {
164 return;
165 }
166
167
168 int nAtoms = openMMEnergy.getMolecularAssembly().getAtomArray().length;
169 for (int i = 0; i < nAtoms; i++) {
170 gk.udpateSoluteParameters(i);
171 }
172
173 boolean usePerfectRadii = gk.getUsePerfectRadii();
174 double perfectRadiiScale = 1.0;
175 if (usePerfectRadii) {
176
177 perfectRadiiScale = 0.0;
178 }
179
180 double[] baseRadii = gk.getBaseRadii();
181 if (usePerfectRadii) {
182 baseRadii = gk.getPerfectRadii();
183 }
184 double[] overlapScale = gk.getOverlapScale();
185 double[] descreenRadius = gk.getDescreenRadii();
186 double[] neckFactors = gk.getNeckScale();
187
188 boolean nea = gk.getNativeEnvironmentApproximation();
189 double lambdaElec = openMMEnergy.getPmeNode().getAlchemicalParameters().permLambda;
190
191 for (Atom atom : atoms) {
192 int index = atom.getXyzIndex() - 1;
193 double chargeUseFactor = 1.0;
194 if (!atom.getUse() || !atom.getElectrostatics()) {
195 chargeUseFactor = 0.0;
196 }
197
198 double lambdaScale = lambdaElec;
199 if (!atom.applyLambda()) {
200 lambdaScale = 1.0;
201 }
202
203 double baseSize = baseRadii[index] * OpenMM_NmPerAngstrom;
204 double descreenSize = descreenRadius[index] * OpenMM_NmPerAngstrom * perfectRadiiScale;
205
206 chargeUseFactor *= lambdaScale;
207 double overlapScaleUseFactor = nea ? 1.0 : chargeUseFactor;
208 overlapScaleUseFactor = overlapScaleUseFactor * perfectRadiiScale;
209 double overlap = overlapScale[index] * overlapScaleUseFactor;
210 double neckFactor = neckFactors[index] * overlapScaleUseFactor;
211
212 MultipoleType multipoleType = atom.getMultipoleType();
213 setParticleParameters_1(index, multipoleType.charge * chargeUseFactor, baseSize, overlap, descreenSize, neckFactor);
214 }
215
216
217 GeneralizedKirkwood.NonPolarModel nonpolar = gk.getNonPolarModel();
218 switch (nonpolar) {
219 case BORN_CAV_DISP, BORN_SOLV -> {
220
221 double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal * OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm;
222 setSurfaceAreaFactor(-surfaceTension);
223 }
224 }
225
226 updateParametersInContext(openMMEnergy.getContext());
227 }
228 }