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.commands;
39
40 import ffx.crystal.Crystal;
41 import ffx.numerics.Potential;
42 import ffx.potential.ForceFieldEnergy;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.bonded.Residue;
45 import ffx.potential.cli.AtomSelectionOptions;
46 import ffx.potential.cli.PotentialCommand;
47 import ffx.potential.extended.ExtendedSystem;
48 import ffx.potential.parsers.PDBFilter;
49 import ffx.potential.parsers.SystemFilter;
50 import ffx.potential.parsers.XPHFilter;
51 import ffx.potential.parsers.XYZFilter;
52 import ffx.utilities.FFXBinding;
53 import picocli.CommandLine.Command;
54 import picocli.CommandLine.Mixin;
55 import picocli.CommandLine.Option;
56 import picocli.CommandLine.Parameters;
57
58 import java.io.BufferedWriter;
59 import java.io.File;
60 import java.io.FileWriter;
61 import java.io.IOException;
62 import java.math.BigDecimal;
63 import java.nio.file.Files;
64 import java.util.ArrayList;
65 import java.util.Arrays;
66 import java.util.Collections;
67 import java.util.List;
68 import java.util.regex.Matcher;
69 import java.util.regex.Pattern;
70
71 import static ffx.potential.utils.StructureMetrics.momentsOfInertia;
72 import static ffx.potential.utils.StructureMetrics.radiusOfGyration;
73 import static java.lang.String.format;
74 import static org.apache.commons.io.FilenameUtils.getBaseName;
75 import static org.apache.commons.io.FilenameUtils.removeExtension;
76 import static org.apache.commons.math3.util.FastMath.pow;
77
78
79
80
81
82
83
84
85 @Command(description = " Compute the force field potential energy for a CpHMD system.", name = "PhEnergy")
86 public class PhEnergy extends PotentialCommand {
87
88 @Mixin
89 private AtomSelectionOptions atomSelectionOptions = new AtomSelectionOptions();
90
91
92
93
94 @Option(names = {"-m", "--moments"}, paramLabel = "false", defaultValue = "false",
95 description = "Print out electrostatic moments.")
96 private boolean moments = false;
97
98
99
100
101 @Option(names = {"--rg", "--gyrate"}, paramLabel = "false", defaultValue = "false",
102 description = "Print out the radius of gyration.")
103 private boolean gyrate = false;
104
105
106
107
108 @Option(names = {"--in", "--inertia"}, paramLabel = "false", defaultValue = "false",
109 description = "Print out the moments of inertia.")
110 private boolean inertia = false;
111
112
113
114
115 @Option(names = {"-g", "--gradient"}, paramLabel = "false", defaultValue = "false",
116 description = "Compute the atomic gradient as well as energy.")
117 private boolean gradient = false;
118
119
120
121
122
123 @Option(names = {"-v", "--verbose"}, paramLabel = "false", defaultValue = "false",
124 description = "Print out all energy components for each snapshot.")
125 private boolean verbose = false;
126
127
128
129
130 @Option(names = {"--pH", "--constantPH"}, paramLabel = "7.4",
131 description = "pH value for the energy evaluation. (Only applies when esvTerm is true)")
132 double pH = 7.4;
133
134 @Option(names = {"--aFi", "--arcFile"}, paramLabel = "traj",
135 description = "A file containing snapshots to evaluate on when using a PDB as a reference to build from. There is currently no default.")
136 private String arcFileName = null;
137
138 @Option(names = {"--bar", "--mbar"}, paramLabel = "false",
139 description = "Run (restartable) energy evaluations for MBAR. Requires an ARC file to be passed in. Set the tautomer flag to true for tautomer parameterization.")
140 boolean mbar = false;
141
142 @Option(names = {"--numLambda", "--nL", "--nw"}, paramLabel = "-1",
143 description = "Required for lambda energy evaluations. Ensure numLambda is consistent with the trajectory lambdas, i.e. gaps between traj can be filled easily. nL >> nTraj is recommended.")
144 int numLambda = -1;
145
146 @Option(names = {"--outputDir", "--oD"}, paramLabel = "",
147 description = "Where to place MBAR files. Default is ../mbarFiles/energy_(window#).mbar. Will write out a file called energy_0.mbar.")
148 String outputDirectory = "";
149
150 @Option(names = {"--lambdaDerivative", "--lD"}, paramLabel = "false",
151 description = "Perform dU/dL evaluations and save to mbarFiles.")
152 boolean derivatives = false;
153
154 @Option(names = {"--perturbTautomer"}, paramLabel = "false",
155 description = "Change tautomer instead of lambda state for MBAR energy evaluations.")
156 boolean tautomer = false;
157
158 @Option(names = {"--testEndStateEnergies"}, paramLabel = "false",
159 description = "Test both ESV energy end states as if the polarization damping factor is initialized from the respective protonated or deprotonated state")
160 boolean testEndstateEnergies = false;
161
162 @Option(names = {"--recomputeAverage"}, paramLabel = "false",
163 description = "Recompute average position and spit out said structure from trajectory")
164 boolean recomputeAverage = false;
165
166
167
168
169 @Parameters(arity = "1", paramLabel = "file",
170 description = "The atomic coordinate file in PDB or XPH format.")
171 private String filename = null;
172
173 public double energy = 0.0;
174 public ForceFieldEnergy forceFieldEnergy = null;
175
176
177
178
179 public PhEnergy() {
180 super();
181 }
182
183
184
185
186
187
188 public PhEnergy(FFXBinding binding) {
189 super(binding);
190 }
191
192
193
194
195
196
197 public PhEnergy(String[] args) {
198 super(args);
199 }
200
201
202
203
204 public PhEnergy run() {
205
206 if (!init()) {
207 return this;
208 }
209 if (mbar) {
210 System.setProperty("lock.esv.states", "false");
211 }
212
213
214 activeAssembly = getActiveAssembly(filename);
215 if (activeAssembly == null) {
216 logger.info(helpString());
217 return this;
218 }
219
220
221 filename = activeAssembly.getFile().getAbsolutePath();
222
223 logger.info("\n Running Energy on " + filename);
224 forceFieldEnergy = activeAssembly.getPotentialEnergy();
225
226
227 File esv = new File(removeExtension(filename) + ".esv");
228 if (!esv.exists()) {
229 esv = null;
230 }
231
232 ExtendedSystem esvSystem = new ExtendedSystem(activeAssembly, pH, esv);
233 if (testEndstateEnergies && BigDecimal.valueOf(esvSystem.getExtendedLambdas()[0]).compareTo(BigDecimal.ZERO) == 0) {
234 for (Atom atom : esvSystem.getExtendedAtoms()) {
235 int atomIndex = atom.getArrayIndex();
236 if (esvSystem.isTitratingHeavy(atomIndex)) {
237 double endstatePolar = esvSystem.getTitrationUtils().getPolarizability(atom, 0.0, 0.0, atom.getPolarizeType().polarizability);
238 double sixth = 1.0 / 6.0;
239 atom.getPolarizeType().pdamp = pow(endstatePolar, sixth);
240 }
241 }
242 } else if (testEndstateEnergies && BigDecimal.valueOf(esvSystem.getExtendedLambdas()[0]).compareTo(BigDecimal.ONE) == 0) {
243 for (Atom atom : esvSystem.getExtendedAtoms()) {
244 int atomIndex = atom.getArrayIndex();
245 if (esvSystem.isTitratingHeavy(atomIndex)) {
246 double endstatePolar = esvSystem.getTitrationUtils().getPolarizability(atom, 1.0, 0.0, atom.getPolarizeType().polarizability);
247 double sixth = 1.0 / 6.0;
248 atom.getPolarizeType().pdamp = pow(endstatePolar, sixth);
249 }
250 }
251 }
252 esvSystem.setConstantPh(pH);
253 int numESVs = esvSystem.getExtendedResidueList().size();
254 forceFieldEnergy.attachExtendedSystem(esvSystem);
255 logger.info(format(" Attached extended system with %d residues.", numESVs));
256
257
258 atomSelectionOptions.setActiveAtoms(activeAssembly);
259
260 int nVars = forceFieldEnergy.getNumberOfVariables();
261 double[] x = new double[nVars];
262 forceFieldEnergy.getCoordinates(x);
263 double[] averageCoordinates = Arrays.copyOf(x, x.length);
264 if (gradient) {
265 double[] g = new double[nVars];
266 int nAts = nVars / 3;
267 energy = forceFieldEnergy.energyAndGradient(x, g, true);
268 logger.info(" Atom X, Y and Z Gradient Components (kcal/mol/A)");
269 for (int i = 0; i < nAts; i++) {
270 int i3 = 3 * i;
271 logger.info(format(" %7d %16.8f %16.8f %16.8f", i + 1, g[i3], g[i3 + 1], g[i3 + 2]));
272 }
273 } else {
274 energy = forceFieldEnergy.energy(x, true);
275 }
276
277 if (moments) {
278 forceFieldEnergy.getPmeNode().computeMoments(activeAssembly.getActiveAtomArray(), false);
279 }
280
281 if (gyrate) {
282 double rg = radiusOfGyration(activeAssembly.getActiveAtomArray());
283 logger.info(format(" Radius of gyration: %10.5f A", rg));
284 }
285
286 if (inertia) {
287 momentsOfInertia(activeAssembly.getActiveAtomArray(), false, true, true);
288 }
289
290 SystemFilter systemFilter;
291 if (arcFileName != null) {
292 File arcFile = new File(arcFileName);
293 systemFilter = new XPHFilter(arcFile, activeAssembly, activeAssembly.getForceField(), activeAssembly.getProperties(), esvSystem);
294 } else {
295 systemFilter = potentialFunctions.getFilter();
296 if (systemFilter instanceof XYZFilter) {
297 systemFilter = new XPHFilter(activeAssembly.getFile(), activeAssembly, activeAssembly.getForceField(), activeAssembly.getProperties(), esvSystem);
298 systemFilter.readFile();
299 logger.info("Reading ESV lambdas from XPH file");
300 forceFieldEnergy.getCoordinates(x);
301 forceFieldEnergy.energy(x, true);
302 }
303 }
304
305 if (mbar) {
306
307 computeESVEnergiesAndWriteFile(systemFilter, esvSystem);
308 return this;
309 }
310
311 if (systemFilter instanceof XPHFilter || systemFilter instanceof PDBFilter) {
312 int index = 1;
313 while (systemFilter.readNext()) {
314 index++;
315 Crystal crystal = activeAssembly.getCrystal();
316 forceFieldEnergy.setCrystal(crystal);
317 forceFieldEnergy.getCoordinates(x);
318 if (recomputeAverage) {
319 for (int i = 0; i < x.length; i++) {
320 averageCoordinates[i] += x[i];
321 }
322 }
323 if (verbose) {
324 logger.info(format(" Snapshot %4d", index));
325 if (!crystal.aperiodic()) {
326 logger.info(format("\n Density: %6.3f (g/cc)",
327 crystal.getDensity(activeAssembly.getMass())));
328 }
329 energy = forceFieldEnergy.energy(x, true);
330 } else {
331 energy = forceFieldEnergy.energy(x, false);
332 logger.info(format(" Snapshot %4d: %16.8f (kcal/mol)", index, energy));
333 }
334 }
335 if (recomputeAverage) {
336 for (int i = 0; i < x.length; i++) {
337 x[i] = averageCoordinates[i] / index;
338 }
339 forceFieldEnergy.setCoordinates(x);
340 }
341 }
342 if (recomputeAverage) {
343 saveByOriginalExtension(activeAssembly, filename);
344 }
345
346 return this;
347 }
348
349 @Override
350 public List<Potential> getPotentials() {
351 List<Potential> potentials;
352 if (forceFieldEnergy == null) {
353 potentials = Collections.emptyList();
354 } else {
355 potentials = Collections.singletonList(forceFieldEnergy);
356 }
357 return potentials;
358 }
359
360 void computeESVEnergiesAndWriteFile(SystemFilter systemFilter, ExtendedSystem esvSystem) {
361
362 File mbarFile;
363 File mbarGradFile;
364 double[] lambdas;
365 if (outputDirectory.isEmpty()) {
366 File dir = new File(filename).getParentFile();
367 File parentDir = dir.getParentFile();
368 int thisRung = -1;
369 Pattern pattern = Pattern.compile("(\\d+)");
370 Matcher matcher = pattern.matcher(dir.getName());
371 if (matcher.find()) {
372 thisRung = Integer.parseInt(matcher.group(1));
373 }
374 assert thisRung != -1 : "Could not determine the rung number from the directory name.";
375 mbarFile = new File(parentDir.getAbsolutePath() + File.separator + "mbarFiles" + File.separator + "energy_"
376 + thisRung + ".mbar");
377 mbarGradFile = new File(parentDir.getAbsolutePath() + File.separator + "mbarFiles" + File.separator + "derivative_"
378 + thisRung + ".mbar");
379 mbarFile.getParentFile().mkdir();
380 File[] lsFiles = parentDir.listFiles();
381 List<File> rungFiles = new ArrayList<>();
382 for (File file : lsFiles) {
383 if (file.isDirectory() && file.getName().matches("\\d+")) {
384 rungFiles.add(file);
385 }
386 }
387 if (numLambda == -1) {
388 numLambda = rungFiles.size();
389 }
390 lambdas = new double[numLambda];
391 for (int i = 0; i < numLambda; i++) {
392 double dL = 1.0 / (numLambda - 1);
393 lambdas[i] = i * dL;
394 }
395
396
397 logger.info(" Computing energies for each lambda state for generation of mbar file.");
398 logger.info(" MBAR File: " + mbarFile);
399 logger.info(" Lambda States: " + Arrays.toString(lambdas));
400 } else {
401 mbarFile = new File(outputDirectory + File.separator + "energy_0.mbar");
402 mbarGradFile = new File(outputDirectory + File.separator + "derivative_0.mbar");
403 lambdas = new double[numLambda];
404 if (numLambda == -1) {
405 logger.severe("numLambda must be set when outputDirectory is set.");
406 }
407 for (int i = 0; i < numLambda; i++) {
408 double dL = 1.0 / (numLambda - 1);
409 lambdas[i] = i * dL;
410 }
411 }
412
413 int progress = 1;
414 if (mbarFile.exists()) {
415
416 try {
417 progress = (int) Files.lines(mbarFile.toPath()).count() - 1;
418 } catch (IOException e) {
419 logger.severe("Error reading MBAR file for restart.");
420 progress = 1;
421 }
422 for (int i = 0; i < progress; i++) {
423 systemFilter.readNext();
424 }
425 progress += 1;
426 logger.info("\n Restarting MBAR file at snapshot " + progress);
427 }
428
429 if (systemFilter instanceof XPHFilter || systemFilter instanceof PDBFilter) {
430 int index = progress;
431 double[] x = new double[forceFieldEnergy.getNumberOfVariables()];
432 try (FileWriter fw = new FileWriter(mbarFile, mbarFile.exists());
433 BufferedWriter writer = new BufferedWriter(fw);
434 FileWriter fwGrad = new FileWriter(mbarGradFile, mbarGradFile.exists());
435 BufferedWriter writerGrad = new BufferedWriter(fwGrad)
436 ) {
437
438 StringBuilder sb = new StringBuilder(systemFilter.countNumModels() + "\t" + "298.0" + "\t" + getBaseName(filename));
439 StringBuilder sbGrad = new StringBuilder(systemFilter.countNumModels() + "\t" + "298.0" + "\t" + getBaseName(filename));
440 logger.info(" MBAR file temp is hardcoded to 298.0 K. Please change if necessary.");
441 sb.append("\n");
442 sbGrad.append("\n");
443 if (progress == 1) {
444 writer.write(sb.toString());
445 writer.flush();
446 logger.info(" Header: " + sb);
447 if (derivatives) {
448 writerGrad.write(sbGrad.toString());
449 writerGrad.flush();
450 logger.info(" Header: " + sbGrad);
451 }
452 }
453 while (systemFilter.readNext()) {
454
455 sb = new StringBuilder("\t" + index + "\t");
456 sbGrad = new StringBuilder("\t" + index + "\t");
457 index++;
458 Crystal crystal = activeAssembly.getCrystal();
459 forceFieldEnergy.setCrystal(crystal);
460 for (double lambda : lambdas) {
461 if (tautomer) {
462 setESVTautomer(lambda, esvSystem);
463 } else {
464 setESVLambda(lambda, esvSystem);
465 }
466 forceFieldEnergy.getCoordinates(x);
467 if (derivatives) {
468 energy = forceFieldEnergy.energyAndGradient(x, new double[x.length * 3]);
469 double grad = esvSystem.getDerivatives()[0];
470 sbGrad.append(grad).append(" ");
471 } else {
472 energy = forceFieldEnergy.energy(x, false);
473 }
474 sb.append(energy).append(" ");
475 }
476 sb.append("\n");
477 writer.write(sb.toString());
478 writer.flush();
479 logger.info(sb.toString());
480 if (derivatives) {
481 sbGrad.append("\n");
482 writerGrad.write(sbGrad.toString());
483 writerGrad.flush();
484 logger.info(sbGrad.toString());
485 }
486 }
487 } catch (IOException e) {
488 logger.severe("Error writing to MBAR file.");
489 }
490 }
491 }
492
493
494
495
496
497
498
499
500 public static void setESVLambda(double lambda, ExtendedSystem extendedSystem) {
501 List<Residue> residueList = extendedSystem.getExtendedResidueList();
502 if (residueList.size() == 1 || (residueList.size() == 2 && extendedSystem.isTautomer(residueList.getFirst()))) {
503 extendedSystem.setTitrationLambda(residueList.getFirst(), lambda, false);
504 } else {
505 if (residueList.isEmpty()) {
506 logger.severe(" No residues found in the extended system.");
507 } else {
508 logger.severe(" Only one lambda path is allowed for MBAR energy evaluations.");
509 }
510 }
511 }
512
513
514
515
516
517
518
519
520 public static void setESVTautomer(double tautomer, ExtendedSystem extendedSystem) {
521 List<Residue> residueList = extendedSystem.getExtendedResidueList();
522 if (residueList.size() == 1 || (residueList.size() == 2 && extendedSystem.isTautomer(residueList.getFirst()))) {
523 extendedSystem.setTautomerLambda(residueList.getFirst(), tautomer, false);
524 } else {
525 if (residueList.isEmpty()) {
526 logger.severe(" No residues found in the extended system.");
527 } else {
528 logger.severe(" Only one lambda path is allowed for MBAR energy evaluations.");
529 }
530 }
531 }
532 }