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