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.nonbonded.pme;
39
40 import edu.rit.pj.IntegerSchedule;
41 import edu.rit.util.Range;
42 import ffx.crystal.Crystal;
43 import ffx.potential.parameters.ForceField;
44
45 import java.util.logging.Logger;
46
47 import static java.lang.String.format;
48 import static org.apache.commons.math3.util.FastMath.pow;
49
50 public class AlchemicalParameters {
51
52 private static final Logger logger = Logger.getLogger(AlchemicalParameters.class.getName());
53
54
55
56
57 private final Polarization polarization;
58
59
60
61 public final AlchemicalMode mode;
62
63
64
65 public double permLambdaAlpha = 1.0;
66
67
68
69 public double permLambdaExponent = 3.0;
70
71
72
73 public double permLambdaStart = 0.4;
74
75
76
77 public double permLambdaEnd = 1.0;
78
79
80
81
82 public double polLambdaStart = 0.75;
83 public double polLambdaEnd = 1.0;
84
85
86
87 public double polLambdaExponent = 3.0;
88
89
90
91 public boolean doLigandVaporElec = true;
92
93
94
95 public boolean doLigandGKElec = false;
96
97
98
99
100 public boolean doNoLigandCondensedSCF = true;
101
102
103
104 public double lAlpha = 0.0;
105 public double dlAlpha = 0.0;
106 public double d2lAlpha = 0.0;
107 public double dEdLSign = 1.0;
108
109
110
111 public double lPowPerm = 1.0;
112 public double dlPowPerm = 0.0;
113 public double d2lPowPerm = 0.0;
114 public boolean doPermanentRealSpace = true;
115 public double permanentScale = 1.0;
116
117
118
119 public double lPowPol = 1.0;
120 public double dlPowPol = 0.0;
121 public double d2lPowPol = 0.0;
122 public boolean doPolarization = true;
123
124
125
126
127
128
129
130
131
132
133
134
135
136 public double polarizationScale = 1.0;
137
138
139
140
141 public double polLambda = 1.0;
142
143
144
145
146 public double permLambda = 1.0;
147
148
149
150 public Crystal vaporCrystal = null;
151 public int[][][] vaporLists = null;
152 public Range[] vacuumRanges = null;
153 public IntegerSchedule vaporPermanentSchedule = null;
154 public IntegerSchedule vaporEwaldSchedule = null;
155
156 public AlchemicalParameters(ForceField forceField, boolean lambdaTerm,
157 boolean nnTerm, Polarization polarization) {
158 this.polarization = polarization;
159
160 AlchemicalMode tempMode;
161 try {
162 tempMode = AlchemicalMode.valueOf(forceField.getString("ALCHEMICAL_MODE", "OST").toUpperCase());
163 } catch (IllegalArgumentException e) {
164 logger.info(" Invalid value for alchemical-mode; reverting to OST");
165 tempMode = AlchemicalMode.OST;
166 }
167 mode = tempMode;
168
169 if (lambdaTerm) {
170
171 permLambdaAlpha = forceField.getDouble("PERMANENT_LAMBDA_ALPHA", 2.0);
172 if (permLambdaAlpha < 0.0 || permLambdaAlpha > 3.0) {
173 logger.warning("Invalid value for permanent-lambda-alpha (<0.0 || >3.0); reverting to 2.0");
174 permLambdaAlpha = 2.0;
175 }
176
177
178
179
180
181
182
183
184
185 permLambdaExponent = forceField.getDouble("PERMANENT_LAMBDA_EXPONENT", 3.0);
186 if (permLambdaExponent < 0.0) {
187 logger.warning("Invalid value for permanent-lambda-exponent (<0.0); reverting to 3.0");
188 permLambdaExponent = 3.0;
189 }
190
191
192
193
194
195
196
197
198
199 polLambdaExponent = forceField.getDouble("POLARIZATION_LAMBDA_EXPONENT", 3.0);
200 if (polLambdaExponent < 0.0) {
201 logger.warning("Invalid value for polarization-lambda-exponent (<0.0); reverting to 3.0");
202 polLambdaExponent = 3.0;
203 }
204
205
206 permLambdaStart = forceField.getDouble("PERMANENT_LAMBDA_START", 0.4);
207 if (permLambdaStart < 0.0 || permLambdaStart > 1.0) {
208 logger.warning("Invalid value for perm-lambda-start (<0.0 || >1.0); reverting to 0.4");
209 permLambdaStart = 0.4;
210 }
211
212
213 permLambdaEnd = forceField.getDouble("PERMANENT_LAMBDA_END", 1.0);
214 if (permLambdaEnd < permLambdaStart || permLambdaEnd > 1.0) {
215 logger.warning("Invalid value for perm-lambda-end (<start || >1.0); reverting to 1.0");
216 permLambdaEnd = 1.0;
217 }
218
219
220
221
222
223
224
225
226
227
228 polLambdaStart = forceField.getDouble("POLARIZATION_LAMBDA_START", 0.75);
229 if (polLambdaStart < 0.0 || polLambdaStart > 1.0) {
230 logger.warning("Invalid value for polarization-lambda-start; reverting to 0.75");
231 polLambdaStart = 0.75;
232 }
233
234
235
236
237
238
239 polLambdaEnd = forceField.getDouble("POLARIZATION_LAMBDA_END", 1.0);
240 if (polLambdaEnd < polLambdaStart || polLambdaEnd > 1.0) {
241 logger.warning(
242 "Invalid value for polarization-lambda-end (<start || >1.0); reverting to 1.0");
243 polLambdaEnd = 1.0;
244 }
245
246
247
248 if (mode == AlchemicalMode.OST) {
249 doLigandVaporElec = forceField.getBoolean("LIGAND_VAPOR_ELEC", true);
250 doLigandGKElec = forceField.getBoolean("LIGAND_GK_ELEC", false);
251 doNoLigandCondensedSCF = forceField.getBoolean("NO_LIGAND_CONDENSED_SCF", true);
252 } else {
253 doLigandVaporElec = false;
254 doLigandGKElec = false;
255 doNoLigandCondensedSCF = false;
256 }
257 } else if (nnTerm) {
258 permLambdaAlpha = 0.0;
259 permLambdaExponent = 1.0;
260 polLambdaExponent = 1.0;
261 permLambdaStart = 0.0;
262 permLambdaEnd = 1.0;
263 polLambdaStart = 0.0;
264 polLambdaEnd = 1.0;
265
266
267 if (mode == AlchemicalMode.OST) {
268 doLigandVaporElec = forceField.getBoolean("LIGAND_VAPOR_ELEC", true);
269 } else {
270 doLigandVaporElec = false;
271 }
272 doLigandGKElec = false;
273 doNoLigandCondensedSCF = false;
274 }
275 }
276
277 public String toString() {
278 StringBuilder sb = new StringBuilder(" Alchemical Parameters\n");
279 sb.append(
280 format(
281 " Permanent Multipole Range: %5.3f-%5.3f\n", permLambdaStart, permLambdaEnd));
282 sb.append(format(" Permanent Multipole Softcore Alpha: %5.3f\n", permLambdaAlpha));
283 sb.append(format(" Permanent Multipole Lambda Exponent: %5.3f\n", permLambdaExponent));
284 if (polarization != Polarization.NONE) {
285 sb.append(format(" Polarization Lambda Exponent: %5.3f\n", polLambdaExponent));
286 sb.append(
287 format(
288 " Polarization Range: %5.3f-%5.3f\n", polLambdaStart, polLambdaEnd));
289 sb.append(format(" Condensed SCF Without Ligand: %B\n", doNoLigandCondensedSCF));
290 }
291 if (!doLigandGKElec) {
292 sb.append(format(" Vapor Electrostatics: %B\n", doLigandVaporElec));
293 } else {
294 sb.append(format(" GK Electrostatics at L=0: %B\n", doLigandGKElec));
295 }
296 return sb.toString();
297 }
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316 public void update(double lambda) {
317 lPowPerm = 1.0;
318 permLambda = 1.0;
319 dlPowPerm = 0.0;
320 d2lPowPerm = 0.0;
321 lAlpha = 0.0;
322 dlAlpha = 0.0;
323 d2lAlpha = 0.0;
324 if (lambda < permLambdaStart) {
325 lPowPerm = 0.0;
326 permLambda = 0.0;
327 } else if (lambda <= permLambdaEnd) {
328 double permWindow = permLambdaEnd - permLambdaStart;
329 double permLambdaScale = 1.0 / permWindow;
330 permLambda = permLambdaScale * (lambda - permLambdaStart);
331 if (mode == AlchemicalMode.OST) {
332 lAlpha = permLambdaAlpha * (1.0 - permLambda) * (1.0 - permLambda);
333 dlAlpha = permLambdaAlpha * (1.0 - permLambda);
334 d2lAlpha = -permLambdaAlpha;
335
336 lPowPerm = pow(permLambda, permLambdaExponent);
337 dlPowPerm = permLambdaExponent * pow(permLambda, permLambdaExponent - 1.0);
338 d2lPowPerm = 0.0;
339 if (permLambdaExponent >= 2.0) {
340 d2lPowPerm =
341 permLambdaExponent
342 * (permLambdaExponent - 1.0)
343 * pow(permLambda, permLambdaExponent - 2.0);
344 }
345
346 dlAlpha *= permLambdaScale;
347 d2lAlpha *= (permLambdaScale * permLambdaScale);
348 dlPowPerm *= permLambdaScale;
349 d2lPowPerm *= (permLambdaScale * permLambdaScale);
350 }
351 }
352
353
354 lPowPol = 1.0;
355 polLambda = 1.0;
356 dlPowPol = 0.0;
357 d2lPowPol = 0.0;
358 if (lambda < polLambdaStart) {
359 lPowPol = 0.0;
360 polLambda = 0.0;
361 } else if (lambda <= polLambdaEnd) {
362 double polWindow = polLambdaEnd - polLambdaStart;
363 double polLambdaScale = 1.0 / polWindow;
364 polLambda = polLambdaScale * (lambda - polLambdaStart);
365 if (mode == AlchemicalMode.OST) {
366 if (polLambdaExponent > 0.0) {
367 lPowPol = pow(polLambda, polLambdaExponent);
368 if (polLambdaExponent >= 1.0) {
369 dlPowPol = polLambdaExponent * pow(polLambda, polLambdaExponent - 1.0);
370 if (polLambdaExponent >= 2.0) {
371 d2lPowPol =
372 polLambdaExponent
373 * (polLambdaExponent - 1.0)
374 * pow(polLambda, polLambdaExponent - 2.0);
375 }
376 }
377 }
378
379 dlPowPol *= polLambdaScale;
380 d2lPowPol *= (polLambdaScale * polLambdaScale);
381 }
382 }
383 }
384
385
386
387
388
389
390
391
392 public enum AlchemicalMode {
393 OST, SCALE
394 }
395
396 }