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 setSoluteDielectric(1.0);
70 setDielectricOffset(gk.getDescreenOffset() * OpenMM_NmPerAngstrom);
71
72 boolean usePerfectRadii = gk.getUsePerfectRadii();
73 double perfectRadiiScale = 1.0;
74 if (usePerfectRadii) {
75
76 perfectRadiiScale = 0.0;
77 }
78
79
80 int tanhRescale = 0;
81 if (gk.getTanhCorrection() && !usePerfectRadii) {
82 tanhRescale = 1;
83 }
84 double[] betas = gk.getTanhBetas();
85 setTanhRescaling(tanhRescale);
86 setTanhParameters(betas[0], betas[1], betas[2]);
87
88 double[] baseRadius = gk.getBaseRadii();
89 if (usePerfectRadii) {
90 baseRadius = gk.getPerfectRadii();
91 }
92
93 double[] overlapScale = gk.getOverlapScale();
94 double[] descreenRadius = gk.getDescreenRadii();
95 double[] neckFactor = gk.getNeckScale();
96
97 if (!usePerfectRadii && logger.isLoggable(Level.FINE)) {
98 logger.fine(" GK Base Radii Descreen Radius Overlap Scale Overlap");
99 }
100
101 Atom[] atoms = openMMEnergy.getMolecularAssembly().getAtomArray();
102 int nAtoms = atoms.length;
103 for (int i = 0; i < nAtoms; i++) {
104 MultipoleType multipoleType = atoms[i].getMultipoleType();
105 double base = baseRadius[i] * OpenMM_NmPerAngstrom;
106 double descreen = descreenRadius[i] * OpenMM_NmPerAngstrom * perfectRadiiScale;
107 double overlap = overlapScale[i] * perfectRadiiScale;
108 double neck = neckFactor[i] * perfectRadiiScale;
109 addParticle_1(multipoleType.charge, base, overlap, descreen, neck);
110 if (!usePerfectRadii && logger.isLoggable(Level.FINE)) {
111 logger.fine(format(" %s %8.6f %8.6f %5.3f", atoms[i].toString(), baseRadius[i], descreenRadius[i], overlapScale[i]));
112 }
113 }
114
115 setProbeRadius(gk.getProbeRadius() * OpenMM_NmPerAngstrom);
116
117 GeneralizedKirkwood.NonPolarModel nonpolar = gk.getNonPolarModel();
118 switch (nonpolar) {
119 case BORN_CAV_DISP, BORN_SOLV -> {
120
121 double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal * OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm;
122 setIncludeCavityTerm(OpenMM_True);
123 setSurfaceAreaFactor(-surfaceTension);
124 }
125
126 default -> setIncludeCavityTerm(OpenMM_False);
127 }
128
129 ForceField forceField = openMMEnergy.getMolecularAssembly().getForceField();
130 int forceGroup = forceField.getInteger("GK_FORCE_GROUP", 1);
131 setForceGroup(forceGroup);
132 logger.log(Level.INFO, format(" Generalized Kirkwood force \t\t%d", forceGroup));
133 }
134
135
136
137
138
139
140
141 public static Force constructForce(OpenMMEnergy openMMEnergy) {
142 GeneralizedKirkwood gk = openMMEnergy.getGK();
143 if (gk == null) {
144 return null;
145 }
146 return new AmoebaGeneralizedKirkwoodForce(openMMEnergy);
147 }
148
149
150
151
152
153
154
155
156 public void updateForce(Atom[] atoms, OpenMMEnergy openMMEnergy) {
157 GeneralizedKirkwood gk = openMMEnergy.getGK();
158 if (gk == null || pointer == null) {
159 return;
160 }
161
162
163 int nAtoms = openMMEnergy.getMolecularAssembly().getAtomArray().length;
164 for (int i = 0; i < nAtoms; i++) {
165 gk.udpateSoluteParameters(i);
166 }
167
168 boolean usePerfectRadii = gk.getUsePerfectRadii();
169 double perfectRadiiScale = 1.0;
170 if (usePerfectRadii) {
171
172 perfectRadiiScale = 0.0;
173 }
174
175 double[] baseRadii = gk.getBaseRadii();
176 if (usePerfectRadii) {
177 baseRadii = gk.getPerfectRadii();
178 }
179 double[] overlapScale = gk.getOverlapScale();
180 double[] descreenRadius = gk.getDescreenRadii();
181 double[] neckFactors = gk.getNeckScale();
182
183 boolean nea = gk.getNativeEnvironmentApproximation();
184 double lambdaElec = openMMEnergy.getSystem().getLambdaElec();
185
186 for (Atom atom : atoms) {
187 int index = atom.getXyzIndex() - 1;
188 double chargeUseFactor = 1.0;
189 if (!atom.getUse() || !atom.getElectrostatics()) {
190 chargeUseFactor = 0.0;
191 }
192
193 double lambdaScale = lambdaElec;
194 if (!atom.applyLambda()) {
195 lambdaScale = 1.0;
196 }
197
198 double baseSize = baseRadii[index] * OpenMM_NmPerAngstrom;
199 double descreenSize = descreenRadius[index] * OpenMM_NmPerAngstrom * perfectRadiiScale;
200
201 chargeUseFactor *= lambdaScale;
202 double overlapScaleUseFactor = nea ? 1.0 : chargeUseFactor;
203 overlapScaleUseFactor = overlapScaleUseFactor * perfectRadiiScale;
204 double overlap = overlapScale[index] * overlapScaleUseFactor;
205 double neckFactor = neckFactors[index] * overlapScaleUseFactor;
206
207 MultipoleType multipoleType = atom.getMultipoleType();
208 setParticleParameters_1(index, multipoleType.charge * chargeUseFactor, baseSize, overlap, descreenSize, neckFactor);
209 }
210
211
212 GeneralizedKirkwood.NonPolarModel nonpolar = gk.getNonPolarModel();
213 switch (nonpolar) {
214 case BORN_CAV_DISP, BORN_SOLV -> {
215
216 double surfaceTension = gk.getSurfaceTension() * OpenMM_KJPerKcal * OpenMM_AngstromsPerNm * OpenMM_AngstromsPerNm;
217 setSurfaceAreaFactor(-surfaceTension);
218 }
219 }
220
221 updateParametersInContext(openMMEnergy.getContext());
222 }
223 }