View Javadoc
1   //******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  //******************************************************************************
38  package ffx.algorithms.commands;
39  
40  import ffx.algorithms.cli.AlgorithmsCommand;
41  import ffx.algorithms.cli.MinimizeOptions;
42  import ffx.crystal.Crystal;
43  import ffx.numerics.Potential;
44  import ffx.potential.ForceFieldEnergy;
45  import ffx.potential.MolecularAssembly;
46  import ffx.potential.bonded.LambdaInterface;
47  import ffx.potential.cli.AlchemicalOptions;
48  import ffx.potential.cli.AtomSelectionOptions;
49  import ffx.potential.cli.TopologyOptions;
50  import ffx.potential.parsers.PDBFilter;
51  import ffx.potential.parsers.SystemFilter;
52  import ffx.potential.parsers.XYZFilter;
53  import ffx.utilities.FFXBinding;
54  import ffx.utilities.FileUtils;
55  import org.apache.commons.io.FilenameUtils;
56  import picocli.CommandLine.Command;
57  import picocli.CommandLine.Mixin;
58  import picocli.CommandLine.Parameters;
59  
60  import java.io.File;
61  import java.util.ArrayList;
62  import java.util.Collections;
63  import java.util.List;
64  
65  import static java.lang.String.format;
66  
67  /**
68   * The Minimize script uses a limited-memory BFGS algorithm to minimize the energy of a molecular system.
69   * <br>
70   * Usage:
71   * <br>
72   * ffxc Minimize [options] &lt;filename&gt; [file2...]
73   */
74  @Command(description = " Run L-BFGS minimization on a system.", name = "Minimize")
75  public class Minimize extends AlgorithmsCommand {
76  
77    @Mixin
78    private MinimizeOptions minimizeOptions;
79  
80    @Mixin
81    private AtomSelectionOptions atomSelectionOptions;
82  
83    @Mixin
84    private AlchemicalOptions alchemicalOptions;
85  
86    @Mixin
87    private TopologyOptions topologyOptions;
88  
89    /**
90     * The final argument(s) should be one or more filenames.
91     */
92    @Parameters(arity = "1..*", paramLabel = "files", description = "Atomic coordinate files in PDB or XYZ format.")
93    private List<String> filenames = null;
94  
95    private MolecularAssembly[] topologies;
96    private Potential potential;
97  
98    /**
99     * The minimization algorithm.
100    */
101   private ffx.algorithms.optimize.Minimize minimize;
102 
103   /**
104    * Minimize Constructor.
105    */
106   public Minimize() {
107     super();
108   }
109 
110   /**
111    * Minimize Constructor.
112    *
113    * @param binding The Binding to use.
114    */
115   public Minimize(FFXBinding binding) {
116     super(binding);
117   }
118 
119   /**
120    * Minimize constructor that sets the command line arguments.
121    *
122    * @param args Command line arguments.
123    */
124   public Minimize(String[] args) {
125     super(args);
126   }
127 
128   /**
129    * {@inheritDoc}
130    */
131   @Override
132   public Minimize run() {
133 
134     if (!init()) {
135       return this;
136     }
137 
138     // Determine the number of topologies to be read and allocate the array.
139     int numTopologies = topologyOptions.getNumberOfTopologies(filenames);
140     int threadsPerTopology = topologyOptions.getThreadsPerTopology(numTopologies);
141     topologies = new MolecularAssembly[numTopologies];
142 
143     // Turn on computation of lambda derivatives if softcore atoms exist.
144     alchemicalOptions.setAlchemicalProperties();
145     topologyOptions.setAlchemicalProperties(numTopologies);
146 
147     // Read in files.
148     if (filenames == null || filenames.isEmpty()) {
149       activeAssembly = getActiveAssembly(null);
150       if (activeAssembly == null) {
151         logger.info(helpString());
152         return this;
153       }
154       filenames = new ArrayList<>();
155       filenames.add(activeAssembly.getFile().getName());
156       topologies[0] = alchemicalOptions.processFile(topologyOptions, activeAssembly, 0);
157     } else {
158       logger.info(format(" Initializing %d topologies...", numTopologies));
159       for (int i = 0; i < numTopologies; i++) {
160         topologies[i] = alchemicalOptions.openFile(algorithmFunctions,
161             topologyOptions, threadsPerTopology, filenames.get(i), i);
162       }
163       activeAssembly = topologies[0];
164     }
165 
166     if (topologies.length == 1) {
167       atomSelectionOptions.setActiveAtoms(topologies[0]);
168     }
169 
170     // Configure the potential to test.
171     StringBuilder sb = new StringBuilder("\n Minimizing energy of ");
172     potential = topologyOptions.assemblePotential(topologies, sb);
173     logger.info(sb.toString());
174 
175     LambdaInterface linter = (potential instanceof LambdaInterface) ? (LambdaInterface) potential : null;
176     if (linter != null) {
177       linter.setLambda(alchemicalOptions.getInitialLambda());
178     }
179 
180     SystemFilter systemFilter = algorithmFunctions.getFilter();
181 
182     double[] x = new double[potential.getNumberOfVariables()];
183     potential.getCoordinates(x);
184     potential.energy(x, true);
185 
186     minimize = new ffx.algorithms.optimize.Minimize(topologies[0], potential, algorithmListener);
187     minimize.minimize(minimizeOptions.getNBFGS(), minimizeOptions.getEps(), minimizeOptions.getIterations());
188 
189     potential.getCoordinates(x);
190     activeAssembly = systemFilter.getActiveMolecularSystem();
191     updateTitle(potential.energy(x, true));
192     if (topologies.length > 1) {
193       // Handle Multiple Topology Cases.
194       for (MolecularAssembly molecularAssembly : topologies) {
195         String modelFilename = molecularAssembly.getFile().getAbsolutePath();
196 
197         if (baseDir == null || !baseDir.exists() || !baseDir.isDirectory() || !baseDir.canWrite()) {
198           baseDir = new File(FilenameUtils.getFullPath(modelFilename));
199         }
200 
201         String dirName = baseDir.toString() + File.separator;
202         String fileName = FilenameUtils.getName(modelFilename);
203         String ext = FilenameUtils.getExtension(fileName);
204         fileName = FilenameUtils.removeExtension(fileName);
205 
206         if (ext.toUpperCase().contains("XYZ")) {
207           algorithmFunctions.saveAsXYZ(molecularAssembly, new File(dirName + fileName + ".xyz"));
208         } else {
209           algorithmFunctions.saveAsPDB(molecularAssembly, new File(dirName + fileName + ".pdb"));
210         }
211       }
212     } else {
213       // Handle Single Topology Cases.
214       setActiveAssembly(topologies[0]);
215       String modelFilename = activeAssembly.getFile().getAbsolutePath();
216       if (baseDir == null || !baseDir.exists() || !baseDir.isDirectory() || !baseDir.canWrite()) {
217         baseDir = new File(FilenameUtils.getFullPath(modelFilename));
218       }
219 
220       String dirName = baseDir.toString() + File.separator;
221       String fileName = FilenameUtils.getName(modelFilename);
222       String ext = FilenameUtils.getExtension(fileName);
223       fileName = FilenameUtils.removeExtension(fileName);
224       File saveFile;
225       SystemFilter writeFilter;
226       if (ext.toUpperCase().contains("XYZ")) {
227         saveFile = new File(dirName + fileName + ".xyz");
228         writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
229             activeAssembly.getProperties());
230         algorithmFunctions.saveAsXYZ(activeAssembly, saveFile);
231       } else if (ext.toUpperCase().contains("ARC")) {
232         saveFile = new File(dirName + fileName + ".arc");
233         saveFile = algorithmFunctions.versionFile(saveFile);
234         writeFilter = new XYZFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
235             activeAssembly.getProperties());
236         algorithmFunctions.saveAsXYZ(activeAssembly, saveFile);
237       } else {
238         saveFile = new File(dirName + fileName + ".pdb");
239         saveFile = algorithmFunctions.versionFile(saveFile);
240         PDBFilter pdbFilter = new PDBFilter(saveFile, activeAssembly, activeAssembly.getForceField(),
241             activeAssembly.getProperties());
242         writeFilter = pdbFilter;
243         int numModels = systemFilter.countNumModels();
244         if (numModels > 1) {
245           pdbFilter.setModelNumbering(0);
246         }
247         pdbFilter.writeFile(saveFile, true, false, false);
248       }
249 
250       if (systemFilter instanceof XYZFilter || systemFilter instanceof PDBFilter) {
251         while (systemFilter.readNext()) {
252           Crystal crystal = activeAssembly.getCrystal();
253           ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy();
254           forceFieldEnergy.setCrystal(crystal);
255           if (systemFilter instanceof PDBFilter) {
256             FileUtils.append(saveFile, "ENDMDL\n");
257             minimize.minimize(minimizeOptions.getEps(), minimizeOptions.getIterations());
258             PDBFilter pdbFilter = (PDBFilter) systemFilter;
259             pdbFilter.writeFile(saveFile, true, false, false);
260           } else if (systemFilter instanceof XYZFilter) {
261             minimize.minimize(minimizeOptions.getEps(), minimizeOptions.getIterations());
262             writeFilter.writeFile(saveFile, true);
263           }
264         }
265         if (systemFilter instanceof PDBFilter) {
266           FileUtils.append(saveFile, "END\n");
267         }
268       }
269     }
270 
271     return this;
272   }
273 
274   /**
275    * {@inheritDoc}
276    */
277   @Override
278   public String toString() {
279     StringBuilder sb = new StringBuilder();
280     if (minimize != null) {
281       sb.append(minimize.toString());
282     } else {
283       sb.append("No minimization has been performed.");
284     }
285     return sb.toString();
286   }
287 
288   /**
289    * Get the list of energies recorded during minimization.
290    *
291    * @return A List of Double energies, or null if minimization was not run.
292    */
293   public List<Double> getEnergyList() {
294     if (minimize != null) {
295       return minimize.getEnergyList();
296     }
297     return null;
298   }
299 
300   /**
301    * Get the final RMS gradient after minimization.
302    *
303    * @return The final RMS gradient or NaN if minimization was not run.
304    */
305   public double getRMSGradient() {
306     if (minimize != null) {
307       return minimize.getRMSGradient();
308     }
309     return Double.NaN;
310   }
311 
312   /**
313    * Get the final energy after minimization.
314    *
315    * @return The final energy or NaN if minimization was not run.
316    */
317   public double getEnergy() {
318     if (minimize != null) {
319       return minimize.getEnergy();
320     }
321     return Double.NaN;
322   }
323 
324   /**
325    * {@inheritDoc}
326    */
327   @Override
328   public List<Potential> getPotentials() {
329     List<Potential> potentials;
330     if (potential == null) {
331       potentials = Collections.emptyList();
332     } else {
333       potentials = Collections.singletonList(potential);
334     }
335     return potentials;
336   }
337 
338 }