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.potential.cli;
39  
40  import static java.lang.Double.parseDouble;
41  import static java.lang.String.format;
42  
43  import ffx.numerics.Potential;
44  import ffx.numerics.switching.MultiplicativeSwitch;
45  import ffx.numerics.switching.PowerSwitch;
46  import ffx.numerics.switching.SquaredTrigSwitch;
47  import ffx.numerics.switching.UnivariateSwitchingFunction;
48  import ffx.potential.DualTopologyEnergy;
49  import ffx.potential.MolecularAssembly;
50  import ffx.potential.QuadTopologyEnergy;
51  import ffx.potential.bonded.Atom;
52  import java.util.Arrays;
53  import java.util.Collections;
54  import java.util.HashSet;
55  import java.util.List;
56  import java.util.Set;
57  import java.util.logging.Level;
58  import java.util.logging.Logger;
59  import java.util.regex.Matcher;
60  import java.util.stream.Collectors;
61  import picocli.CommandLine.ArgGroup;
62  import picocli.CommandLine.Option;
63  
64  /**
65   * Represents command line options for scripts that utilize multiple physical topologies.
66   *
67   * @author Michael J. Schnieders
68   * @author Jacob M. Litman
69   * @since 1.0
70   */
71  public class TopologyOptions {
72  
73    /** The logger for this class. */
74    public static final Logger logger = Logger.getLogger(TopologyOptions.class.getName());
75  
76    /**
77     * The ArgGroup keeps the TopologyOptions together when printing help.
78     */
79    @ArgGroup(heading = "%n Alchemical Options for Dual and Quad Topologies%n", validate = false)
80    private final TopologyOptionGroup group = new TopologyOptionGroup();
81  
82    /**
83     * --ac2 or --alchemicalAtoms2 Specify alchemical atoms [ALL, NONE, Range(s): 1-3,6-N].
84     *
85     * @return Returns alchemical atoms for the 2nd topology.
86     */
87    public String getAlchemicalAtoms2() {
88      return group.alchemicalAtoms2;
89    }
90  
91    /**
92     * --uc2 or --unchargedAtoms2 Specify atoms without electrostatics [ALL, NONE, Range(s): 1-3,6-N].
93     *
94     * @return Returns atoms without electrostatics for the 2nd topology.
95     */
96    public String getUnchargedAtoms2() {
97      return group.unchargedAtoms2;
98    }
99  
100   /**
101    * -np or --nParallel sets the number of topologies to evaluate in parallel; currently 1, 2, or 4.
102    *
103    * @return Returns Number of topologies to evaluate in parallel.
104    */
105   public int getNPar() {
106     return group.nPar;
107   }
108 
109   /**
110    * --uaA or --unsharedA sets atoms unique to the A dual-topology, as period-separated hyphenated
111    * ranges or singletons.
112    *
113    * @return Return atoms unique to the A dual-topology.
114    */
115   public String getUnsharedA() {
116     return group.unsharedA;
117   }
118 
119   /**
120    * --uaB or --unsharedB sets atoms unique to the A dual-topology, as period-separated hyphenated
121    * ranges or singletons.
122    *
123    * @return Return atoms unique to the B dual-topology.
124    */
125   public String getUnsharedB() {
126     return group.unsharedB;
127   }
128 
129   /**
130    * -sf or --switchingFunction
131    *
132    * <p>Sets the switching function to be used by dual topologies.
133    *
134    * <ul>
135    *   <li>TRIG produces the function sin^2(pi/2*lambda)*E1(lambda) + cos^2(pi/2*lambda)*E2(1-lambda)</li>
136    *   <li>MULT uses a 5th-order polynomial switching function with zero first and second derivatives at
137    *   the end (same function as used for van der Waals switch)</li>
138    *   <li>A number uses the original function, of l^beta*E1(lambda) + (1-lambda)^beta*E2(1-lambda).</li>
139    * </ul>
140    *
141    * <p>All of these are generalizations of <code>Udt = f(l)*E1(l) + f(1-l)*E2(1-lambda)</code>,
142    * where f(l) is a continuous switching function such that f(0) = 0, f(1) = 1, and 0 &lt;= f(l) &lt;= 1
143    * for lambda 0-1. The trigonometric switch can be restated thusly, since cos^2(pi/2*lambda) is
144    * identical to sin^2(pi/2*(1-lambda)), f(1-l).
145    *
146    * @return Returns the lambda function.
147    */
148   public String getLambdaFunction() {
149     return group.lambdaFunction;
150   }
151 
152   /**
153    * The number of topologies to run in parallel.
154    *
155    * @param nThreads The number of threads available.
156    * @param nTopologies The number of topologies to run in parallel.
157    * @return Number of topologies to run in parallel.
158    */
159   public int getNumParallel(int nThreads, int nTopologies) {
160     int numParallel = getNPar();
161     if (nThreads % numParallel != 0) {
162       logger.warning(format(
163           " Number of threads available %d not evenly divisible by np %d; reverting to sequential.",
164           nThreads, numParallel));
165       numParallel = 1;
166     } else if (nTopologies % numParallel != 0) {
167       logger.warning(
168           format(" Number of topologies %d not evenly divisible by np %d; reverting to sequential.",
169               nTopologies, numParallel));
170       numParallel = 1;
171     }
172     return numParallel;
173   }
174 
175   /**
176    * Performs the bulk of the work of setting up a multi-topology system.
177    *
178    * <p>The sb StringBuilder is often something like "Timing energy and gradients for". The method
179    * will append the exact type of Potential being assembled.
180    *
181    * @param assemblies Opened MolecularAssembly(s).
182    * @param threadsAvail Number of available threads.
183    * @param sb A StringBuilder describing what is to be done.
184    * @return a {@link ffx.numerics.Potential} object.
185    */
186   public Potential assemblePotential(MolecularAssembly[] assemblies, int threadsAvail,
187       StringBuilder sb) {
188     int nargs = assemblies.length;
189     int numPar = getNumParallel(threadsAvail, nargs);
190     UnivariateSwitchingFunction sf = nargs > 1 ? getSwitchingFunction() : null;
191     List<Integer> uniqueA;
192     List<Integer> uniqueB;
193     if (assemblies.length >= 4) {
194       uniqueA = getUniqueAtomsA(assemblies[0]);
195       uniqueB = getUniqueAtomsB(assemblies[2]);
196     } else {
197       uniqueA = Collections.emptyList();
198       uniqueB = Collections.emptyList();
199     }
200     return getTopology(assemblies, sf, uniqueA, uniqueB, numPar, sb);
201   }
202 
203   /**
204    * Return the switching function between topology energies.
205    *
206    * @return the switching function.
207    */
208   public UnivariateSwitchingFunction getSwitchingFunction() {
209     UnivariateSwitchingFunction sf;
210     String lambdaFunction = getLambdaFunction();
211     if (!lambdaFunction.equalsIgnoreCase("1.0")) {
212       String lf = lambdaFunction.toUpperCase();
213       switch (lf) {
214         case "TRIG" -> sf = new SquaredTrigSwitch(false);
215         case "MULT" -> sf = new MultiplicativeSwitch(1.0, 0.0);
216         default -> {
217           try {
218             double beta = parseDouble(lf);
219             sf = new PowerSwitch(1.0, beta);
220           } catch (NumberFormatException ex) {
221             logger.warning(format(
222                 "Argument to option -sf %s could not be properly parsed; using default linear switch",
223                 lambdaFunction));
224             sf = new PowerSwitch(1.0, 1.0);
225           }
226         }
227       }
228     } else {
229       sf = new PowerSwitch(1.0, 1.0);
230     }
231     return sf;
232   }
233 
234   /**
235    * Collect unique atoms for the A dual-topology.
236    *
237    * @param topology A MolecularAssembly from dual-topology A.
238    * @return A List of Integers.
239    */
240   public List<Integer> getUniqueAtomsA(MolecularAssembly topology) {
241     return getUniqueAtoms(topology, "A", getUnsharedA());
242   }
243 
244   /**
245    * Configure a Dual-, Quad- or Oct- Topology.
246    *
247    * @param topologies The topologies.
248    * @param sf The Potential switching function.
249    * @param uniqueA The unique atoms of topology A.
250    * @param uniqueB The unique atoms of topology B.
251    * @param numParallel The number of energies to evaluate in parallel.
252    * @param sb A StringBuilder for logging.
253    * @return The Potential for the Topology.
254    */
255   public Potential getTopology(MolecularAssembly[] topologies, UnivariateSwitchingFunction sf,
256       List<Integer> uniqueA, List<Integer> uniqueB, int numParallel, StringBuilder sb) {
257     Potential potential = null;
258     switch (topologies.length) {
259       case 1 -> {
260         sb.append("Single Topology ");
261         potential = topologies[0].getPotentialEnergy();
262       }
263       case 2 -> {
264         sb.append("Dual Topology ");
265         DualTopologyEnergy dte = new DualTopologyEnergy(topologies[0], topologies[1], sf);
266         if (numParallel == 2) {
267           dte.setParallel(true);
268         }
269         potential = dte;
270       }
271       case 4 -> {
272         sb.append("Quad Topology ");
273         DualTopologyEnergy dta = new DualTopologyEnergy(topologies[0], topologies[1], sf);
274         DualTopologyEnergy dtb = new DualTopologyEnergy(topologies[3], topologies[2], sf);
275         QuadTopologyEnergy qte = new QuadTopologyEnergy(dta, dtb, uniqueA, uniqueB);
276         if (numParallel >= 2) {
277           qte.setParallel(true);
278           if (numParallel == 4) {
279             dta.setParallel(true);
280             dtb.setParallel(true);
281           }
282         }
283         potential = qte;
284       }
285       default -> logger.severe(" Must have 2, 4, or 8 topologies!");
286     }
287     sb.append(Arrays.stream(topologies).map(MolecularAssembly::toString)
288         .collect(Collectors.joining(", ", " [", "] ")));
289     return potential;
290   }
291 
292   /**
293    * Collect unique atoms for a dual-topology. List MUST be sorted at the end.
294    *
295    * @param assembly A MolecularAssembly from the dual topology.
296    * @param label Either 'A' or 'B'.
297    * @param unshared Atoms this dual topology isn't sharing.
298    * @return A sorted List of Integers.
299    */
300   public static List<Integer> getUniqueAtoms(MolecularAssembly assembly, String label,
301       String unshared) {
302     if (!unshared.isEmpty()) {
303       logger.info(" Finding unique atoms for dual topology " + label);
304       Set<Integer> indices = new HashSet<>();
305       String[] toks = unshared.split("\\.");
306       Atom[] atoms1 = assembly.getAtomArray();
307       for (String range : toks) {
308         Matcher m = AlchemicalOptions.rangeRegEx.matcher(range);
309         if (m.find()) {
310           int rangeStart = Integer.parseInt(m.group(1));
311           int rangeEnd = (m.group(2) != null) ? Integer.parseInt(m.group(2)) : rangeStart;
312           if (rangeStart > rangeEnd) {
313             logger.severe(format(" Range %s was invalid start was greater than end", range));
314           }
315           logger.info(
316               format(" Range %s for %s, start %d end %d", range, label, rangeStart, rangeEnd));
317           if(logger.isLoggable(Level.FINE)) {
318             logger.fine(format(" First atom in range: %s", atoms1[rangeStart - 1]));
319             if (rangeEnd > rangeStart) {
320               logger.fine(format(" Last atom in range: %s", atoms1[rangeEnd - 1]));
321             }
322           }
323           for (int i = rangeStart; i <= rangeEnd; i++) {
324             indices.add(i - 1);
325           }
326         }
327       }
328       int counter = 0;
329       Set<Integer> adjustedIndices = new HashSet<>(); // Indexed by common variables in dtA.
330       int atomLength = atoms1.length;
331       for (int i = 0; i < atomLength; i++) {
332         Atom ai = atoms1[i];
333         if (indices.contains(i)) {
334           if (ai.applyLambda()) {
335             logger.warning(format(
336                 " Ranges defined in %s should not overlap with ligand atoms they are assumed to not be shared.",
337                 label));
338           } else {
339             logger.fine(format(" Unshared %s: %d variables %d-%d", label, i, counter, counter + 2));
340             for (int j = 0; j < 3; j++) {
341               adjustedIndices.add(counter + j);
342             }
343           }
344         }
345         if (!ai.applyLambda()) {
346           counter += 3;
347         }
348       }
349       return adjustedIndices.stream().sorted().collect(Collectors.toList());
350     } else {
351       return Collections.emptyList();
352     }
353   }
354 
355   /**
356    * Collect unique atoms for the B dual-topology.
357    *
358    * @param topology A MolecularAssembly from dual-topology B.
359    * @return A List of Integers.
360    */
361   public List<Integer> getUniqueAtomsB(MolecularAssembly topology) {
362     return getUniqueAtoms(topology, "B", getUnsharedB());
363   }
364 
365   /**
366    * If any softcore Atoms have been detected.
367    *
368    * @return Presence of softcore Atoms.
369    */
370   public boolean hasSoftcore() {
371     String alchemicalAtoms2 = getAlchemicalAtoms2();
372     return (alchemicalAtoms2 != null && !alchemicalAtoms2.equalsIgnoreCase("NONE")
373         && !alchemicalAtoms2.equalsIgnoreCase(""));
374   }
375 
376   /**
377    * Set the alchemical atoms for this topology.
378    *
379    * @param topology a {@link ffx.potential.MolecularAssembly} object.
380    */
381   public void setSecondSystemAlchemistry(MolecularAssembly topology) {
382     String alchemicalAtoms2 = getAlchemicalAtoms2();
383     AlchemicalOptions.setAlchemicalAtoms(topology, alchemicalAtoms2);
384   }
385 
386   /**
387    * Set uncharged atoms for this topology.
388    *
389    * @param topology a {@link ffx.potential.MolecularAssembly} object.
390    */
391   public void setSecondSystemUnchargedAtoms(MolecularAssembly topology) {
392     String unchargedAtoms2 = getUnchargedAtoms2();
393     AlchemicalOptions.setUnchargedAtoms(topology, unchargedAtoms2);
394   }
395 
396   /**
397    * Collection of Topology Options.
398    */
399   private static class TopologyOptionGroup {
400 
401     /** --ac2 or --alchemicalAtoms2 Specify alchemical atoms [ALL, NONE, Range(s): 1-3,6-N]. */
402     @Option(names = {"--ac2", "--alchemicalAtoms2"}, paramLabel = "<selection>", defaultValue = "",
403         description = "Specify alchemical atoms for the 2nd topology [ALL, NONE, Range(s): 1-3,6-N].")
404     String alchemicalAtoms2 = "";
405 
406     /**
407      * --uc2 or --unchargedAtoms2 Specify atoms without electrostatics [ALL, NONE, Range(s):
408      * 1-3,6-N]."
409      */
410     @Option(names = {"--uc2", "--unchargedAtoms2"}, paramLabel = "<selection>", defaultValue = "",
411         description = "Specify atoms without electrostatics for the 2nd topology [ALL, NONE, Range(s): 1-3,6-N].")
412     String unchargedAtoms2 = "";
413 
414     /**
415      * -np or --nParallel sets the number of topologies to evaluate in parallel; currently 1, 2, or
416      * 4.
417      */
418     @Option(names = {"--np", "--nParallel"}, paramLabel = "1",
419         description = "Number of topologies to evaluate in parallel")
420     int nPar = 1;
421 
422     /**
423      * --uaA or --unsharedA sets atoms unique to the A dual-topology, as period-separated hyphenated
424      * ranges or singletons.
425      */
426     @Option(names = {"--uaA", "--unsharedA"}, paramLabel = "-1",
427         description = "Unshared atoms in the A dual topology (e.g. 1-24.32-65).")
428     String unsharedA = null;
429 
430     /**
431      * --uaB or --unsharedB sets atoms unique to the B dual-topology, as period-separated hyphenated
432      * ranges or singletons.
433      */
434     @Option(names = {"--uaB", "--unsharedB"}, paramLabel = "-1",
435         description = "Unshared atoms in the B dual topology (e.g. 1-24.32-65).")
436     String unsharedB = null;
437 
438     /**
439      * -sf or --switchingFunction
440      *
441      * <p>Sets the switching function to be used by dual topologies.
442      *
443      * <p>
444      *
445      * <ul>
446      *   TRIG produces the function sin^2(pi/2*lambda)*E1(lambda) + cos^2(pi/2*lambda)*E2(1-lambda)
447      * </ul>
448      *
449      * <ul>
450      *   MULT uses a 5th-order polynomial switching function with zero first and second derivatives at
451      *   the end (same function as used for van der Waals switch)
452      * </ul>
453      *
454      * <ul>
455      *   A number uses the original function, of l^beta*E1(lambda) + (1-lambda)^beta*E2(1-lambda).
456      * </ul>
457      *
458      * <p>All of these are generalizations of <code>Udt = f(l)*E1(l) + f(1-l)*E2(1-lambda)</code>,
459      * where f(l) is a continuous switching function such that f(0) = 0, f(1) = 1, and 0 <= f(l) <= 1
460      * for lambda 0-1. The trigonometric switch can be restated thusly, since cos^2(pi/2*lambda) is
461      * identical to sin^2(pi/2*(1-lambda)), f(1-l).
462      */
463     @Option(names = {"--sf", "--switchingFunction"}, paramLabel = "1.0",
464         description = "Switching function to use for dual topology: options are TRIG, MULT, or a number (original behavior with specified lambda exponent)")
465     String lambdaFunction = "1.0";
466   }
467 }