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-2024.
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.cli;
39  
40  import static java.lang.String.format;
41  
42  import ffx.algorithms.AlgorithmFunctions;
43  import ffx.algorithms.AlgorithmListener;
44  import ffx.algorithms.optimize.RotamerOptimization;
45  import ffx.crystal.CrystalPotential;
46  import ffx.numerics.Potential;
47  import ffx.potential.DualTopologyEnergy;
48  import ffx.potential.MolecularAssembly;
49  import ffx.potential.QuadTopologyEnergy;
50  import ffx.potential.bonded.LambdaInterface;
51  import ffx.potential.bonded.Polymer;
52  import ffx.potential.bonded.Residue;
53  import ffx.potential.bonded.RotamerLibrary;
54  import ffx.potential.cli.AlchemicalOptions;
55  import ffx.potential.cli.TopologyOptions;
56  import java.io.File;
57  import java.util.ArrayList;
58  import java.util.List;
59  import java.util.logging.Logger;
60  import java.util.regex.Matcher;
61  import java.util.regex.Pattern;
62  import org.apache.commons.configuration2.CompositeConfiguration;
63  import picocli.CommandLine.ArgGroup;
64  import picocli.CommandLine.Option;
65  
66  /**
67   * Represents command line options for scripts that can create multiple walkers, such as multi-walker
68   * OST. Should be kept agnostic to whether it is an MD-based algorithm, or some other flavor of Monte
69   * Carlo.
70   *
71   * @author Michael J. Schnieders
72   * @author Jacob M. Litman
73   * @since 1.0
74   */
75  public class MultiDynamicsOptions {
76  
77    private static final Logger logger = Logger.getLogger(MultiDynamicsOptions.class.getName());
78  
79    /**
80     * The ArgGroup keeps the MultiDynamicsOptionGroup together when printing help.
81     */
82    @ArgGroup(heading = "%n Multiple Walker Options for MPI Simulations%n", validate = false)
83    private final MultiDynamicsOptionGroup group = new MultiDynamicsOptionGroup();
84  
85    /**
86     * If residues selected for distributing initial configurations, performs many-body optimization
87     * for this distribution.
88     *
89     * @param molecularAssemblies an array of {@link ffx.potential.MolecularAssembly} objects.
90     * @param potentials ForceFieldEnergy for each topology.
91     * @param crystalPotential Overall CrystalPotential in use.
92     * @param algorithmFunctions a {@link ffx.algorithms.AlgorithmFunctions} object.
93     * @param rank The MPI rank of this process.
94     * @param worldSize The number of MPI processes.
95     */
96    public void distribute(MolecularAssembly[] molecularAssemblies, Potential[] potentials,
97        CrystalPotential crystalPotential, AlgorithmFunctions algorithmFunctions, int rank,
98        int worldSize) {
99      if (!group.distributeWalkersString.equalsIgnoreCase("AUTO")
100         && !group.distributeWalkersString.equalsIgnoreCase("OFF")) {
101       logger.info(" Distributing walker conformations.");
102       int nSys = molecularAssemblies.length;
103       assert nSys == potentials.length;
104       switch (nSys) {
105         case 1 -> optStructure(molecularAssemblies[0], crystalPotential, algorithmFunctions, rank,
106             worldSize);
107         case 2 -> {
108           DualTopologyEnergy dte = (DualTopologyEnergy) crystalPotential;
109           if (dte.getNumSharedVariables() == dte.getNumberOfVariables()) {
110             logger.info(" Generating starting structures based on dual-topology:");
111             optStructure(molecularAssemblies[0], dte, algorithmFunctions, rank, worldSize);
112           } else {
113             logger.info(
114                 " Generating separate starting structures for each topology of the dual topology:");
115             optStructure(molecularAssemblies[0], potentials[0], algorithmFunctions, rank, worldSize);
116             optStructure(molecularAssemblies[1], potentials[1], algorithmFunctions, rank, worldSize);
117           }
118         }
119         case 4 -> {
120           QuadTopologyEnergy qte = (QuadTopologyEnergy) crystalPotential;
121           optStructure(molecularAssemblies[0], qte.getDualTopA(), algorithmFunctions, rank,
122               worldSize);
123           optStructure(molecularAssemblies[3], qte.getDualTopB(), algorithmFunctions, rank,
124               worldSize);
125         }
126         // Oct-topology is deprecated on account of not working as intended.
127         default -> logger.severe(" First: must have 1, 2, or 4 topologies.");
128       }
129     } else {
130       logger.finer(" Skipping RO-based distribution of initial configurations.");
131     }
132   }
133 
134   /**
135    * If residues selected for distributing initial configurations, performs many-body optimization
136    * for this distribution.
137    *
138    * @param molecularAssemblies an array of {@link ffx.potential.MolecularAssembly} objects.
139    * @param crystalPotential Overall CrystalPotential in use.
140    * @param algorithmFunctions a {@link ffx.algorithms.AlgorithmFunctions} object.
141    * @param rank The MPI rank of this process.
142    * @param worldSize The number of MPI processes.
143    */
144   public void distribute(MolecularAssembly[] molecularAssemblies, CrystalPotential crystalPotential,
145       AlgorithmFunctions algorithmFunctions, int rank, int worldSize) {
146     int ntops = molecularAssemblies.length;
147     Potential[] energies = new Potential[ntops];
148     for (int i = 0; i < ntops; i++) {
149       energies[i] = molecularAssemblies[i].getPotentialEnergy();
150     }
151     distribute(molecularAssemblies, energies, crystalPotential, algorithmFunctions, rank, worldSize);
152   }
153 
154   /**
155    * Synchronous walker communication.
156    *
157    * <p>isSynchronous.
158    *
159    * @return a boolean.
160    */
161   public boolean isSynchronous() {
162     return group.synchronous;
163   }
164 
165   public void setSynchronous(boolean synchronous) {
166     group.synchronous = synchronous;
167   }
168 
169   /**
170    * Opens a file and processes it. Extends the behavior of AlchemicalOptions.openFile by permitting
171    * use of a rank-dependent File.
172    *
173    * @param algorithmFunctions AlgorithmFunctions object.
174    * @param topologyOptions Topology Options.
175    * @param threadsPer Threads to use per system.
176    * @param toOpen Filename to open.
177    * @param topNum Number of the topology to open.
178    * @param alchemicalOptions Alchemical Options.
179    * @param rank Rank in the world communicator.
180    * @param structureFile a {@link java.io.File} object.
181    * @return a {@link ffx.potential.MolecularAssembly} object.
182    */
183   public MolecularAssembly openFile(AlgorithmFunctions algorithmFunctions,
184       TopologyOptions topologyOptions, int threadsPer, String toOpen, int topNum,
185       AlchemicalOptions alchemicalOptions, File structureFile, int rank) {
186     boolean autoDist = group.distributeWalkersString.equalsIgnoreCase("AUTO");
187 
188     if (autoDist) {
189       String openName = format("%s_%d", toOpen, rank + 1);
190       File testFile = new File(openName);
191       if (testFile.exists()) {
192         toOpen = openName;
193       } else {
194         logger.warning(format(" File %s does not exist; using default %s", openName, toOpen));
195       }
196     }
197     MolecularAssembly assembly = alchemicalOptions.openFile(algorithmFunctions, topologyOptions,
198         threadsPer, toOpen, topNum);
199     assembly.setFile(structureFile);
200     return assembly;
201   }
202 
203   /**
204    * Parses --dw into optimization tokens if it's not "OFF", "AUTO", or null.
205    *
206    * @return An array of Strings from splitting the distributed flag.
207    */
208   private String[] parseDistributed() {
209     String distributeWalkersString = group.distributeWalkersString;
210     if (distributeWalkersString.equalsIgnoreCase("OFF")
211         || distributeWalkersString.equalsIgnoreCase("AUTO") || distributeWalkersString.isEmpty()) {
212       return null;
213     }
214     return distributeWalkersString.split("\\.");
215   }
216 
217   /**
218    * Allows walkers to start from multiple conformations; AUTO picks up per-walker conformations as
219    * filename.pdb_(walker number), and specifying a residue starts a rotamer optimization to generate
220    * side-chain configurations to start from.
221    *
222    * @return Returns the Distribute Walkers string.
223    */
224   public String getDistributeWalkersString() {
225     return group.distributeWalkersString;
226   }
227 
228   /**
229    * Distribute side-chain conformations of molecularAssembly.
230    *
231    * @param molecularAssembly To distribute
232    * @param potential Potential to use
233    */
234   private void optStructure(MolecularAssembly molecularAssembly, Potential potential,
235       AlgorithmFunctions algorithmFunctions, int rank, int worldSize) {
236     RotamerLibrary rLib = new RotamerLibrary(false);
237     String[] distribRes = parseDistributed();
238 
239     if (distribRes == null || distribRes.length == 0) {
240       throw new IllegalArgumentException(
241           " Programming error: Must have list of residues to split on!");
242     }
243 
244     LambdaInterface lambdaInterface =
245         (potential instanceof LambdaInterface) ? (LambdaInterface) potential : null;
246     double initLam = -1.0;
247     if (lambdaInterface != null) {
248       initLam = lambdaInterface.getLambda();
249       lambdaInterface.setLambda(0.5);
250     }
251 
252     Pattern chainMatcher = Pattern.compile("^([a-zA-Z])?([0-9]+)$");
253 
254     List<Residue> residueList = new ArrayList<>(distribRes.length);
255 
256     for (String ts : distribRes) {
257       Matcher m = chainMatcher.matcher(ts);
258       Character chainID;
259       int resNum;
260       if (m.find()) {
261         if (m.groupCount() == 2) {
262           chainID = m.group(1).charAt(0);
263           resNum = Integer.parseInt(m.group(2));
264         } else {
265           chainID = ' ';
266           resNum = Integer.parseInt(m.group(1));
267         }
268       } else {
269         logger.warning(format(" Could not parse %s as a valid residue!", ts));
270         continue;
271       }
272       logger.info(format(" Looking for chain %c residue %d", chainID, resNum));
273 
274       for (Polymer p : molecularAssembly.getChains()) {
275         if (p.getChainID() == chainID) {
276           for (Residue r : p.getResidues()) {
277             if (r.getResidueNumber() == resNum && r.setRotamers(rLib) != null) {
278               residueList.add(r);
279             }
280           }
281         }
282       }
283     }
284 
285     if (residueList.isEmpty()) {
286       throw new IllegalArgumentException(" No valid entries for distWalkers!");
287     }
288 
289     AlgorithmListener alist = algorithmFunctions.getDefaultListener();
290     RotamerOptimization ropt = new RotamerOptimization(molecularAssembly, potential, alist);
291     ropt.setRotamerLibrary(rLib);
292 
293     ropt.setThreeBodyEnergy(false);
294 
295     CompositeConfiguration properties = molecularAssembly.getProperties();
296     if (!properties.containsKey("ro-ensembleNumber") && !properties.containsKey(
297         "ro-ensembleEnergy")) {
298       logger.info(format(" Setting ensemble to default of number of walkers %d", worldSize));
299       ropt.setEnsemble(worldSize);
300     }
301 
302     ropt.setPrintFiles(false);
303     ropt.setResiduesIgnoreNull(residueList);
304 
305     RotamerLibrary.measureRotamers(residueList, false);
306 
307     boolean lazyMat = properties.getBoolean("ro-lazyMatrix", false);
308     if (!lazyMat) {
309       properties.clearProperty("ro-lazyMatrix");
310       properties.addProperty("ro-lazyMatrix", true);
311     }
312 
313     ropt.optimize(RotamerOptimization.Algorithm.ALL);
314     ropt.setCoordinatesToEnsemble(rank);
315 
316     // One final energy call to ensure the coordinates are properly set at the
317     // end of rotamer optimization.
318     double[] xyz = new double[potential.getNumberOfVariables()];
319     potential.getCoordinates(xyz);
320     logger.info(" Final Optimized Energy:");
321     potential.energy(xyz, true);
322 
323     if (lambdaInterface != null) {
324       lambdaInterface.setLambda(initLam);
325     }
326 
327     if (!lazyMat) {
328       properties.clearProperty("ro-lazyMatrix");
329     }
330   }
331 
332   public void setDistributeWalkersString(String distributeWalkersString) {
333     group.distributeWalkersString = distributeWalkersString;
334   }
335 
336   public int getFirstDir() {
337     return group.firstDir;
338   }
339 
340   public void setFirstDir(int firstDir) {
341     group.firstDir = firstDir;
342   }
343 
344   /**
345    * Collection of Multi-Walker Dynamics Options.
346    */
347   private static class MultiDynamicsOptionGroup {
348 
349     /** --firstDir The first directory to use for multiple walker jobs. */
350     @Option(names = {"--firstDir"}, defaultValue = "0", paramLabel = "0",
351         description = "The first directory to use for multiple walker jobs.")
352     private int firstDir = 0;
353 
354     /** -y or --synchronous sets synchronous walker communication (not recommended) */
355     @Option(names = {"-y", "--synchronous"}, defaultValue = "false",
356         description = "Walker communication is synchronous")
357     private boolean synchronous = false;
358 
359     /**
360      * -dw or --distributeWalkers allows walkers to start from multiple conformations; AUTO picks up
361      * per-walker conformations as filename.pdb_(walker number), and specifying a residue starts a
362      * rotamer optimization to generate side-chain configurations to start from.
363      */
364     @Option(names = {"--dw", "--distributeWalkers"}, paramLabel = "OFF", defaultValue = "OFF",
365         description = "AUTO: Pick up per-walker configurations as [filename.pdb]_[num], or specify a residue to distribute on.")
366     private String distributeWalkersString = "OFF";
367   }
368 }