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.LambdaInterface;
47 import ffx.potential.bonded.Torsion;
48 import ffx.potential.cli.AlchemicalOptions;
49 import ffx.potential.cli.AtomSelectionOptions;
50 import ffx.potential.cli.PotentialCommand;
51 import ffx.potential.cli.TopologyOptions;
52 import ffx.potential.parsers.PDBFilter;
53 import ffx.potential.parsers.SystemFilter;
54 import ffx.potential.parsers.XYZFilter;
55 import ffx.utilities.FFXBinding;
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.File;
62 import java.nio.file.Files;
63 import java.nio.file.StandardOpenOption;
64 import java.util.ArrayList;
65 import java.util.Collections;
66 import java.util.List;
67 import java.util.PriorityQueue;
68 import java.util.logging.Level;
69
70 import static ffx.potential.utils.StructureMetrics.momentsOfInertia;
71 import static ffx.potential.utils.StructureMetrics.radiusOfGyration;
72 import static ffx.utilities.StringUtils.parseAtomRanges;
73 import static java.lang.String.format;
74 import static org.apache.commons.io.FilenameUtils.getExtension;
75 import static org.apache.commons.io.FilenameUtils.getName;
76 import static org.apache.commons.io.FilenameUtils.removeExtension;
77
78
79
80
81
82
83
84 @Command(name = "Energy", description = " Compute the force field potential energy.")
85 public class Energy extends PotentialCommand {
86
87 @Mixin
88 private AtomSelectionOptions atomSelectionOptions = new AtomSelectionOptions();
89
90 @Mixin
91 private AlchemicalOptions alchemicalOptions = new AlchemicalOptions();
92
93 @Mixin
94 private TopologyOptions topologyOptions = new TopologyOptions();
95
96 @Option(names = {"-m", "--moments"}, paramLabel = "false", defaultValue = "false",
97 description = "Print out electrostatic moments.")
98 private boolean moments = false;
99
100 @Option(names = {"--rg", "--gyrate"}, paramLabel = "false", defaultValue = "false",
101 description = "Print out the radius of gyration.")
102 private boolean gyrate = false;
103
104 @Option(names = {"--in", "--inertia"}, paramLabel = "false", defaultValue = "false",
105 description = "Print out the moments of inertia.")
106 private boolean inertia = false;
107
108 @Option(names = {"-g", "--gradient"}, paramLabel = "false", defaultValue = "false",
109 description = "Compute the atomic gradient as well as energy.")
110 private boolean gradient = false;
111
112 @Option(names = {"--fl", "--findLowest"}, paramLabel = "0", defaultValue = "0",
113 description = "Return the n lowest energies from an ARC/PDB file.")
114 private int fl = 0;
115
116 @Option(names = {"-v", "--verbose"}, paramLabel = "false", defaultValue = "false",
117 description = "Print out all energy components for each snapshot.")
118 private boolean verbose = false;
119
120 @Option(names = {"--dc", "--densityCutoff"}, paramLabel = "0.0", defaultValue = "0.0",
121 description = "Create ARC file of structures above a specified density.")
122 private double dCutoff = 0.0;
123
124 @Option(names = {"--ec", "--energyCutoff"}, paramLabel = "0.0", defaultValue = "0.0",
125 description = "Create ARC file of structures within a specified energy of the lowest energy structure.")
126 private double eCutoff = 0.0;
127
128 @Option(names = {"--pd", "--printDihedral"}, paramLabel = "", defaultValue = "",
129 description = "Atom indices to print dihedral angle values.")
130 private String dihedralAtoms = "";
131
132 @Option(names = {"--pb", "--printBondedTerms"}, paramLabel = "", defaultValue = "false",
133 description = "Print all bonded energy terms.")
134 private boolean printBondedTerms = false;
135
136 @Parameters(arity = "1..4", paramLabel = "file",
137 description = "The atomic coordinate file in PDB or XYZ format.")
138 private List<String> filenames = null;
139
140 public double energy = 0.0;
141 public ForceFieldEnergy forceFieldEnergy = null;
142 private AssemblyState assemblyState = null;
143
144 private Potential potential;
145 MolecularAssembly[] topologies;
146
147 public Energy() {
148 super();
149 }
150
151 public Energy(FFXBinding binding) {
152 super(binding);
153 }
154
155 public Energy(String[] args) {
156 super(args);
157 }
158
159 @Override
160 public Energy run() {
161 if (!init()) {
162 return this;
163 }
164
165 int numTopologies = topologyOptions.getNumberOfTopologies(filenames);
166 int threadsPerTopology = topologyOptions.getThreadsPerTopology(numTopologies);
167 topologies = new MolecularAssembly[numTopologies];
168
169 alchemicalOptions.setAlchemicalProperties();
170 topologyOptions.setAlchemicalProperties(numTopologies);
171
172 if (filenames == null || filenames.isEmpty()) {
173 activeAssembly = getActiveAssembly(filenames);
174 if (activeAssembly == null) {
175 logger.info(helpString());
176 return this;
177 }
178 filenames = new ArrayList<>();
179 filenames.add(activeAssembly.getFile().getName());
180 topologies[0] = alchemicalOptions.processFile(topologyOptions, activeAssembly, 0);
181 } else {
182 logger.info(format(" Initializing %d topologies...", numTopologies));
183 for (int i = 0; i < numTopologies; i++) {
184 topologies[i] = alchemicalOptions.openFile(potentialFunctions, topologyOptions, threadsPerTopology, filenames.get(i), i);
185 }
186 activeAssembly = topologies[0];
187 }
188
189 if (topologies.length == 1) {
190 atomSelectionOptions.setActiveAtoms(topologies[0]);
191 }
192
193 StringBuilder sb = new StringBuilder("\n Calculating energy of ");
194 potential = topologyOptions.assemblePotential(topologies, sb);
195 logger.info(sb.toString());
196
197 LambdaInterface linter = (potential instanceof LambdaInterface) ? (LambdaInterface) potential : null;
198 if (linter != null) {
199 linter.setLambda(alchemicalOptions.getInitialLambda());
200 }
201
202 String filename = activeAssembly.getFile().getAbsolutePath();
203 logger.info("\n Running Energy on " + filename);
204
205 forceFieldEnergy = activeAssembly.getPotentialEnergy();
206 int nVars = potential.getNumberOfVariables();
207 double[] x = new double[nVars];
208 potential.getCoordinates(x);
209
210 if (gradient) {
211 double[] g = new double[nVars];
212 int nAts = nVars / 3;
213 energy = potential.energyAndGradient(x, g, true);
214 logger.info(format(" Atom X, Y and Z Gradient Components (kcal/mol/A)"));
215 for (int i = 0; i < nAts; i++) {
216 int i3 = 3 * i;
217 logger.info(format(" %7d %16.8f %16.8f %16.8f", i + 1, g[i3], g[i3 + 1], g[i3 + 2]));
218 }
219 } else {
220 energy = potential.energy(x, true);
221 }
222
223 if (moments) {
224 forceFieldEnergy.getPmeNode().computeMoments(activeAssembly.getActiveAtomArray(), false);
225 }
226
227 if (gyrate) {
228 double rg = radiusOfGyration(activeAssembly.getActiveAtomArray());
229 logger.info(format(" Radius of gyration: %10.5f A", rg));
230 }
231
232 if (inertia) {
233 momentsOfInertia(activeAssembly.getActiveAtomArray(), false, true, true);
234 }
235
236 ArrayList<Integer> unique = null;
237 if (dihedralAtoms != null && !dihedralAtoms.isEmpty()) {
238 unique = new ArrayList<>(parseAtomRanges("Dihedral Atoms", dihedralAtoms, activeAssembly.getAtomList().size()));
239 printDihedral(activeAssembly, unique);
240 }
241
242 if (printBondedTerms) {
243 forceFieldEnergy.logBondedTermsAndRestraints();
244 }
245
246 SystemFilter systemFilter = potentialFunctions.getFilter();
247 if (systemFilter instanceof XYZFilter || systemFilter instanceof PDBFilter) {
248 int numSnaps = fl;
249 double lowestEnergy = energy;
250 assemblyState = new AssemblyState(activeAssembly);
251 int index = 1;
252 PriorityQueue<StateContainer> lowestEnergyQueue = null;
253 if (fl > 0) {
254 lowestEnergyQueue = new PriorityQueue<>(numSnaps, Collections.reverseOrder());
255 lowestEnergyQueue.add(new StateContainer(assemblyState, lowestEnergy));
256 }
257
258 int numModels = systemFilter.countNumModels();
259 double[] densities = new double[numModels];
260 double[] energies = new double[numModels];
261 energies[0] = energy;
262
263 while (systemFilter.readNext()) {
264 index++;
265 Crystal crystal = activeAssembly.getCrystal();
266 densities[index - 1] = crystal.getDensity(activeAssembly.getMass());
267 forceFieldEnergy.setCrystal(crystal);
268 forceFieldEnergy.getCoordinates(x);
269 if (verbose) {
270 logger.info(format(" Snapshot %4d", index));
271 if (!crystal.aperiodic()) {
272 logger.info(format("\n Density: %6.3f (g/cc)", crystal.getDensity(activeAssembly.getMass())));
273 if (logger.isLoggable(Level.FINE)) {
274 logger.fine(crystal.toString());
275 }
276 }
277 energy = forceFieldEnergy.energy(x, true);
278 } else {
279 energy = forceFieldEnergy.energy(x, false);
280 logger.info(format(" Snapshot %4d: %16.8f (kcal/mol)", index, energy));
281 }
282
283 energies[index - 1] = energy;
284 if (energy < lowestEnergy) {
285 lowestEnergy = energy;
286 }
287
288 if (fl > 0) {
289 StateContainer sc = new StateContainer(new AssemblyState(activeAssembly), energy);
290 if (lowestEnergyQueue.size() < numSnaps) {
291 lowestEnergyQueue.add(sc);
292 } else {
293 StateContainer worst = lowestEnergyQueue.peek();
294 if (worst != null && energy < worst.getEnergy()) {
295 lowestEnergyQueue.poll();
296 lowestEnergyQueue.add(sc);
297 }
298 }
299 }
300
301 if (moments) {
302 forceFieldEnergy.getPmeNode().computeMoments(activeAssembly.getActiveAtomArray(), false);
303 }
304 if (gyrate) {
305 double rg = radiusOfGyration(activeAssembly.getActiveAtomArray());
306 logger.info(format(" Radius of gyration: %10.5f A", rg));
307 }
308 if (inertia) {
309 momentsOfInertia(activeAssembly.getActiveAtomArray(), false, true, true);
310 }
311 if (dihedralAtoms != null && !dihedralAtoms.isEmpty()) {
312 printDihedral(activeAssembly, unique);
313 }
314 if (printBondedTerms) {
315 forceFieldEnergy.logBondedTermsAndRestraints();
316 }
317 }
318
319 String dirString = getBaseDirString(filename);
320
321 if ((eCutoff > 0.0 || dCutoff > 0.0) && numModels > 1) {
322 systemFilter.readNext(true);
323 activeAssembly = systemFilter.getActiveMolecularSystem();
324
325 String name = getName(filename);
326 String ext = getExtension(name);
327 name = removeExtension(name);
328
329 File saveFile;
330 SystemFilter writeFilter;
331 if (ext.toUpperCase().contains("XYZ")) {
332 saveFile = new File(dirString + name + ".xyz");
333 writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(), activeAssembly.getProperties());
334 potentialFunctions.saveAsXYZ(activeAssembly, saveFile);
335 } else if (ext.toUpperCase().contains("ARC")) {
336 saveFile = new File(dirString + name + ".arc");
337 saveFile = potentialFunctions.versionFile(saveFile);
338 writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(), activeAssembly.getProperties());
339 logger.info(" Saving to " + saveFile.getAbsolutePath());
340 try { saveFile.createNewFile(); } catch (Exception e) { }
341 } else {
342 saveFile = new File(dirString + name + ".pdb");
343 saveFile = potentialFunctions.versionFile(saveFile);
344 writeFilter = new PDBFilter(saveFile, activeAssembly, activeAssembly.getForceField(), activeAssembly.getProperties());
345 if (numModels > 1 && writeFilter instanceof PDBFilter) {
346 ((PDBFilter) writeFilter).setModelNumbering(0);
347 }
348 try { saveFile.createNewFile(); } catch (Exception e) { }
349 }
350
351 int structNum = 0;
352 do {
353 if (dCutoff > 0.0 && eCutoff > 0.0) {
354 if (energies[structNum] < lowestEnergy + eCutoff && densities[structNum] > dCutoff) {
355 if (systemFilter instanceof PDBFilter pdbFilter) {
356 pdbFilter.writeFile(saveFile, true, false, false);
357 try { Files.writeString(saveFile.toPath(), "ENDMDL\n", StandardOpenOption.APPEND); } catch (Exception e) { }
358 } else if (systemFilter instanceof XYZFilter) {
359 writeFilter.writeFile(saveFile, true);
360 }
361 }
362 } else if (dCutoff > 0.0) {
363 if (densities[structNum] > dCutoff) {
364 if (systemFilter instanceof PDBFilter pdbFilter) {
365 pdbFilter.writeFile(saveFile, true, false, false);
366 try { Files.writeString(saveFile.toPath(), "ENDMDL\n", StandardOpenOption.APPEND); } catch (Exception e) { }
367 } else if (systemFilter instanceof XYZFilter) {
368 writeFilter.writeFile(saveFile, true);
369 }
370 }
371 } else if (eCutoff > 0.0) {
372 if (energies[structNum] < lowestEnergy + eCutoff) {
373 if (systemFilter instanceof PDBFilter pdbFilter) {
374 pdbFilter.writeFile(saveFile, true, false, false);
375 try { Files.writeString(saveFile.toPath(), "ENDMDL\n", StandardOpenOption.APPEND); } catch (Exception e) { }
376 } else if (systemFilter instanceof XYZFilter) {
377 writeFilter.writeFile(saveFile, true);
378 }
379 }
380 }
381 structNum++;
382 } while (systemFilter.readNext());
383 if (systemFilter instanceof PDBFilter) {
384 try { Files.writeString(saveFile.toPath(), "END\n", StandardOpenOption.APPEND); } catch (Exception e) { }
385 }
386 }
387
388 if (fl > 0) {
389 if (numSnaps > index) {
390 logger.warning(format(" Requested %d snapshots, but file %s has only %d snapshots. All %d energies will be reported", numSnaps, filenames, index, index));
391 numSnaps = index;
392 }
393
394 String name = getName(filenames.get(0));
395 logger.info(" Saving the " + numSnaps + " lowest energy snapshots to " + dirString + name);
396 SystemFilter saveFilter = potentialFunctions.getFilter();
397 File saveFile = potentialFunctions.versionFile(new File(dirString + name));
398 int toSave = Math.min(numSnaps, lowestEnergyQueue.size());
399 for (int i = 0; i < toSave; i++) {
400 StateContainer savedState = lowestEnergyQueue.poll();
401 AssemblyState finalAssembly = savedState.getState();
402 finalAssembly.revertState();
403 String remark = format("Potential Energy: %16.8f (kcal/mol)", savedState.getEnergy());
404 logger.info(format(" Snapshot %d %s", i + 1, remark));
405 saveFilter.writeFile(saveFile, true, new String[]{remark});
406 }
407 }
408 }
409
410 return this;
411 }
412
413 private static void printDihedral(MolecularAssembly molecularAssembly, ArrayList<Integer> indices) {
414 if (indices != null && indices.size() == 4) {
415 Atom[] atoms = molecularAssembly.getAtomArray();
416 Atom atom0 = atoms[indices.get(0)];
417 Atom atom1 = atoms[indices.get(1)];
418 Atom atom2 = atoms[indices.get(2)];
419 Atom atom3 = atoms[indices.get(3)];
420 try {
421 for (Torsion torsion : atom0.getTorsions()) {
422 if (torsion.compare(atom0, atom1, atom2, atom3)) {
423 logger.info("\n Torsion: " + torsion);
424 return;
425 }
426 }
427 } catch (Exception e) {
428 logger.info(" Exception during dihedral print " + e);
429 }
430 logger.info("\n No torsion between atoms:\n " + atom0 + "\n " + atom1 + "\n " + atom2 + "\n " + atom3);
431 }
432 }
433
434 @Override
435 public List<Potential> getPotentials() {
436 if (forceFieldEnergy == null) {
437 return Collections.emptyList();
438 } else {
439 return Collections.singletonList(forceFieldEnergy);
440 }
441 }
442
443 public ForceFieldEnergy getForceFieldEnergy() {
444 return forceFieldEnergy;
445 }
446
447 public double getEnergy() {
448 return energy;
449 }
450
451 private static class StateContainer implements Comparable<StateContainer> {
452 private final AssemblyState state;
453 private final double e;
454 StateContainer(AssemblyState state, double e) { this.state = state; this.e = e; }
455 AssemblyState getState() { return state; }
456 double getEnergy() { return e; }
457 @Override
458 public int compareTo(StateContainer o) { return Double.compare(e, o.getEnergy()); }
459 }
460 }