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.cli;
39  
40  import static java.lang.Integer.parseInt;
41  
42  import ffx.algorithms.optimize.RotamerOptimization;
43  import ffx.algorithms.optimize.RotamerOptimization.Algorithm;
44  import ffx.numerics.math.DoubleMath;
45  import ffx.potential.MolecularAssembly;
46  import ffx.potential.bonded.*;
47  
48  import java.io.File;
49  import java.util.*;
50  import java.util.logging.Logger;
51  
52  import picocli.CommandLine.ArgGroup;
53  import picocli.CommandLine.Option;
54  
55  /**
56   * Represents command line options for scripts that use a many-body expansion for global
57   * optimization.
58   *
59   * @author Michael J. Schnieders
60   * @author Mallory R. Tollefson
61   * @since 1.0
62   */
63  public class ManyBodyOptions {
64  
65    private static final Logger logger = Logger.getLogger(ManyBodyOptions.class.getName());
66  
67    /**
68     * The ArgGroup keeps the ManyBodyOptionGroup together when printing help.
69     */
70    @ArgGroup(heading = "%n Many-Body Optimization Options%n", validate = false)
71    private final ManyBodyOptionGroup group = new ManyBodyOptionGroup();
72  
73    /**
74     * The ArgGroup keeps the BoxOptionGroup together when printing help.
75     */
76    @ArgGroup(heading = "%n Many-Body Box Optimization Options%n", validate = false)
77    private final BoxOptionGroup boxGroup = new BoxOptionGroup();
78  
79    /**
80     * The ArgGroup keeps the WindowOptionGroup together when printing help.
81     */
82    @ArgGroup(heading = "%n Many-Body Window Optimization Options%n", validate = false)
83    private final WindowOptionGroup windowGroup = new WindowOptionGroup();
84  
85    /**
86     * The ArgGroup keeps the WindowOptionGroup together when printing help.
87     */
88    @ArgGroup(heading = "%n Many-Body Energy Expansion and Cut-off Options%n", validate = false)
89    private final EnergyOptionGroup energyGroup = new EnergyOptionGroup();
90  
91    /**
92     * The ArgGroup keeps the ResidueOptionGroup together when printing help.
93     */
94    @ArgGroup(heading = "%n Many-Body Residue Selection Options%n", validate = false)
95    private final ResidueOptionGroup residueGroup = new ResidueOptionGroup();
96  
97    private RotamerOptimization rotamerOptimization;
98    private RotamerLibrary rotamerLibrary;
99  
100   /**
101    * initRotamerOptimization.
102    *
103    * @param rotamerOptimization a {@link RotamerOptimization} object.
104    * @param activeAssembly a {@link ffx.potential.MolecularAssembly} object.
105    */
106   public void initRotamerOptimization(RotamerOptimization rotamerOptimization,
107       MolecularAssembly activeAssembly) {
108     this.rotamerOptimization = rotamerOptimization;
109 
110     // Collect the residues to optimize.
111     List<Residue> residues = collectResidues(activeAssembly);
112     rotamerOptimization.setResidues(residues);
113     rotamerOptimization.setRotamerLibrary(rotamerLibrary);
114 
115     // If the user has not selected an algorithm, it will be chosen based on the number of residues.
116     Algorithm algorithm = getAlgorithm(residues.size());
117 
118     // Configure general options.
119     rotamerOptimization.setDecomposeOriginal(group.decompose);
120     rotamerOptimization.setUseGoldstein(!group.dee);
121     rotamerOptimization.setRevert(group.revert);
122     boolean monteCarloBool = group.monteCarlo > 1;
123     rotamerOptimization.setMonteCarlo(monteCarloBool, group.monteCarlo);
124     File energyRestartFile;
125     if (!group.energyRestart.equalsIgnoreCase("none")) {
126       energyRestartFile = new File(group.energyRestart);
127       rotamerOptimization.setEnergyRestartFile(energyRestartFile);
128     }
129 
130     // Configure Energy Expansion and Pruning Options.
131     rotamerOptimization.setTwoBodyCutoff(energyGroup.twoBodyCutoff);
132     rotamerOptimization.setThreeBodyEnergy(energyGroup.threeBody);
133     rotamerOptimization.setThreeBodyCutoff(energyGroup.threeBodyCutoff);
134     rotamerOptimization.setDistanceCutoff(energyGroup.cutoff);
135     rotamerOptimization.setPruning(energyGroup.prune);
136     rotamerOptimization.setSingletonClashThreshold(energyGroup.clashThreshold);
137     rotamerOptimization.setPairClashThreshold(energyGroup.pairClashThreshold);
138 
139     // Window
140     if (algorithm == Algorithm.WINDOW) {
141       rotamerOptimization.setWindowSize(windowGroup.window);
142       rotamerOptimization.setIncrement(windowGroup.increment);
143     } else if (algorithm == Algorithm.BOX) {
144       // Box
145       parseBoxSelection();
146       if (boxGroup.approxBoxLength < 0) {
147         logger.info(" Negative box length value changed to -1 * input.");
148         boxGroup.approxBoxLength *= -1;
149       }
150       rotamerOptimization.setBoxBorderSize(boxGroup.boxBorderSize);
151       rotamerOptimization.setApproxBoxLength(boxGroup.approxBoxLength);
152       rotamerOptimization.setNumXYZBoxes(boxGroup.numXYZBoxes);
153       rotamerOptimization.setBoxInclusionCriterion(boxGroup.boxInclusionCriterion);
154       rotamerOptimization.setBoxStart(boxGroup.initialBox);
155       rotamerOptimization.setBoxEnd(boxGroup.finalBox);
156       rotamerOptimization.setTitrationBoxes(boxGroup.boxTitration);
157       rotamerOptimization.setTitrationBoxSize(boxGroup.approxBoxLength);
158     }
159   }
160 
161   /**
162    * Collect residues based on residue selection flags.
163    *
164    * @param activeAssembly a {@link ffx.potential.MolecularAssembly} object.
165    */
166   public List<Residue> collectResidues(MolecularAssembly activeAssembly) {
167 
168     // Force re-initialization RotamerLibrary prior to collecting residues and rotamers.
169     initRotamerLibrary(true);
170 
171     // First, interpret the residueGroup.listResidues flag if its set.
172     if (!residueGroup.listResidues.isEmpty() && !residueGroup.listResidues.equalsIgnoreCase("none")) {
173       List<String> stringList = new ArrayList<>();
174       String[] tok = residueGroup.listResidues.split(",");
175       Collections.addAll(stringList, tok);
176       List<Residue> residueList = new ArrayList<>();
177       Polymer[] polymers = activeAssembly.getChains();
178       for (String s : stringList) {
179         Character chainID = s.charAt(0);
180         int i = parseInt(s.substring(1));
181         for (Polymer polymer : polymers) {
182           if (polymer.getChainID() == chainID) {
183             List<Residue> residues = polymer.getResidues();
184             for (Residue residue : residues) {
185               if (residue.getResidueNumber() == i) {
186                 Rotamer[] rotamers = residue.setRotamers(rotamerLibrary);
187                 if (rotamers != null && rotamers.length > 0) {
188                   residueList.add(residue);
189                 }
190               }
191             }
192           }
193         }
194       }
195       return residueList;
196     }
197 
198     // Check that the finish flag is greater than the start flag.
199     if (residueGroup.finish < residueGroup.start) {
200       residueGroup.finish = Integer.MAX_VALUE;
201     }
202     Character chainID = null;
203     if (!residueGroup.chain.equalsIgnoreCase("-1")) {
204       chainID = residueGroup.chain.charAt(0);
205     }
206 
207     // Otherwise, collect all residues with a rotamer.
208     List<Residue> residueList = new ArrayList<>();
209     Polymer[] polymers = activeAssembly.getChains();
210     for (Polymer polymer : polymers) {
211       // Enforce requested chainID.
212       if (chainID != null && chainID != polymer.getChainID()) {
213         continue;
214       }
215       List<Residue> residues = polymer.getResidues();
216       for (Residue residue : residues) {
217         int resID = residue.getResidueNumber();
218         // Enforce requested residue range.
219         if (resID >= residueGroup.start && resID <= residueGroup.finish) {
220           Rotamer[] rotamers = residue.setRotamers(rotamerLibrary);
221           if (rotamers != null) {
222             residueList.add(residue);
223           }
224         }
225       }
226     }
227 
228     return residueList;
229   }
230 
231   /**
232    * Returns the user selected algorithm or one chosen based on number of residues.
233    *
234    * @return Returns the algorithm choice.
235    */
236   public Algorithm getAlgorithm(int numResidues) {
237     if (group.algorithm == 0) {
238       if (numResidues < 100) {
239         return Algorithm.ALL;
240       } else {
241         return Algorithm.BOX;
242       }
243     }
244     return Algorithm.getAlgorithm(group.algorithm);
245   }
246 
247   public boolean getUsingOriginalCoordinates() {
248     return !group.noOriginal;
249   }
250 
251   public void setOriginalCoordinates(boolean useOrig) {
252     group.noOriginal = !useOrig;
253   }
254 
255   public double getApproximate() {
256     return rotamerOptimization.getApproximate();
257   }
258 
259   /**
260    * Gets the restart file created during rotamer optimization.
261    *
262    * @return The restart file.
263    */
264   public File getRestartFile() {
265     return rotamerOptimization.getRestartFile();
266   }
267 
268   public RotamerLibrary getRotamerLibrary(boolean reinit) {
269     initRotamerLibrary(reinit);
270     return rotamerLibrary;
271   }
272 
273   private void initRotamerLibrary(boolean reinit) {
274     if (rotamerLibrary == null || reinit) {
275       boolean useOrigCoordsRotamer = !group.noOriginal;
276       if (group.decompose) {
277         useOrigCoordsRotamer = true;
278       }
279       rotamerLibrary = new RotamerLibrary(
280           RotamerLibrary.ProteinLibrary.intToProteinLibrary(group.library), useOrigCoordsRotamer);
281     }
282   }
283 
284 
285   /** Set allStartResID, boxStart and boxEnd */
286   private void parseBoxSelection() {
287     // Parse the numBoxes flag.
288     String input = boxGroup.numBoxes;
289     Scanner boxNumInput = new java.util.Scanner(input);
290     boxNumInput.useDelimiter(",");
291     int inputLoopCounter = 0;
292     int[] numXYZBoxes = new int[3];
293     numXYZBoxes[0] = 3; // Default
294     while (inputLoopCounter < 3) {
295       if (boxNumInput.hasNextInt()) {
296         numXYZBoxes[inputLoopCounter] = boxNumInput.nextInt();
297         inputLoopCounter++;
298       } else if (boxNumInput.hasNextDouble()) {
299         numXYZBoxes[inputLoopCounter] = (int) Math.floor(boxNumInput.nextDouble());
300         inputLoopCounter++;
301         logger.info(" Double input to nB truncated to integer.");
302       } else if (boxNumInput.hasNext()) {
303         logger.info(" Non-numeric input to nB discarded");
304         boxNumInput.next();
305       } else {
306         logger.info(
307             " Insufficient input to nB. Non-input values assumed either equal to X or default to 3");
308         break;
309       }
310     }
311     boxNumInput.close();
312 
313     // Initialize dimensions not provided.
314     for (int i = inputLoopCounter; i < 3; i++) {
315       numXYZBoxes[i] = numXYZBoxes[0];
316     }
317 
318     // Correct input errors.
319     int totalCount = 1;
320     for (int i = 0; i < 3; i++) {
321       if (numXYZBoxes[i] == 0) {
322         numXYZBoxes[i] = 3;
323         logger.info(" Input of 0 to nB reset to default of 3.");
324       } else if (numXYZBoxes[i] < 0) {
325         numXYZBoxes[i] = -1 * numXYZBoxes[i];
326         logger.info(" Input of negative number to nB reset to positive number");
327       }
328       totalCount *= numXYZBoxes[i];
329     }
330 
331     boxGroup.numXYZBoxes = numXYZBoxes;
332 
333     if (boxGroup.initialBox < 0) {
334       boxGroup.initialBox = 0;
335     }
336     if (boxGroup.finalBox < 0 || boxGroup.finalBox > totalCount) {
337       boxGroup.finalBox = totalCount;
338     }
339 
340   }
341 
342   /** Sets the standard values for properties in rotamer optimization. */
343   private void setRotOptProperties(Algorithm algorithm) {
344 
345   }
346 
347   public void setAlgorithm(int algorithm) {
348     group.algorithm = algorithm;
349   }
350 
351   /**
352    * Ponder and Richards (1) or Richardson (2) rotamer library.
353    *
354    * @return Returns the Rotamer library.
355    */
356   public int getLibrary() {
357     return group.library;
358   }
359 
360   public void setLibrary(int library) {
361     group.library = library;
362   }
363 
364   /**
365    * Nucleic acid library: currently only Richardson available.
366    *
367    * @return Returns a String for the nucleic acid library.
368    */
369   public String getNaLibraryName() {
370     return group.naLibraryName;
371   }
372 
373   public void setNaLibraryName(String naLibraryName) {
374     group.naLibraryName = naLibraryName;
375   }
376 
377   /**
378    * Use dead-end elimination criteria instead of Goldstein criteria.
379    *
380    * @return Returns true if using DEE instead of Goldstein.
381    */
382   public boolean isDee() {
383     return group.dee;
384   }
385 
386   public void setDee(boolean dee) {
387     group.dee = dee;
388   }
389 
390   /**
391    * Single character chain ID of the residues to optimize.
392    *
393    * @return Returns the Chain name.
394    */
395   public String getChain() {
396     return residueGroup.chain;
397   }
398 
399   public void setChain(String chain) {
400     residueGroup.chain = chain;
401   }
402 
403   public boolean getOnlyTitration() {
404     return residueGroup.onlyTitration;
405   }
406 
407   public void setOnlyTitration(boolean onlyTitration) {
408     residueGroup.onlyTitration = onlyTitration;
409   }
410 
411   public int getInterestedResidue() {
412     return residueGroup.interestedResidue;
413   }
414 
415   public void setInterestedResidue(int interestedResidue) {
416     residueGroup.interestedResidue = interestedResidue;
417   }
418 
419   public double getInclusionCutoff() {
420     return residueGroup.inclusionCutoff;
421   }
422 
423   public void setInclusionCutoff(double inclusionCutoff) {
424     residueGroup.inclusionCutoff = inclusionCutoff;
425   }
426 
427   /**
428    * Starting residue to perform the optimization on (-1 exits). For box optimization, first box to
429    * optimize.
430    *
431    * @return Returns the starting index.
432    */
433   public int getStart() {
434     return residueGroup.start;
435   }
436 
437   public void setStart(int start) {
438     residueGroup.start = start;
439   }
440 
441   /**
442    * Final residue to perform the optimization on (-1 exits). For box optimization, final box to
443    * optimize.
444    *
445    * @return Returns the finish index.
446    */
447   public int getFinish() {
448     return residueGroup.finish;
449   }
450 
451   public void setFinish(int finish) {
452     residueGroup.finish = finish;
453   }
454 
455   /**
456    * Cutoff distance for two-body interactions.
457    *
458    * @return Returns the 2-body cutoff.
459    */
460   public double getTwoBodyCutoff() {
461     return energyGroup.twoBodyCutoff;
462   }
463 
464   public void setTwoBodyCutoff(double twoBodyCutoff) {
465     energyGroup.twoBodyCutoff = twoBodyCutoff;
466   }
467 
468   /**
469    * -T or --threeBody Include 3-Body interactions in the elimination criteria.
470    *
471    * @return Returns true if 3-body interactions are being used.
472    */
473   public boolean isThreeBody() {
474     return energyGroup.threeBody;
475   }
476 
477   public void setThreeBody(boolean threeBody) {
478     energyGroup.threeBody = threeBody;
479   }
480 
481   /**
482    * Cutoff distance for three-body interactions.
483    *
484    * @return Returns the 3-body cutoff.
485    */
486   public double getThreeBodyCutoff() {
487     return energyGroup.threeBodyCutoff;
488   }
489 
490   public void setThreeBodyCutoff(double threeBodyCutoff) {
491     energyGroup.threeBodyCutoff = threeBodyCutoff;
492   }
493 
494   /**
495    * Prune no clashes (0), only single clashes (1), or all clashes (2).
496    *
497    * @return Returns the pruning condition.
498    */
499   public int getPrune() {
500     return energyGroup.prune;
501   }
502 
503   public void setPrune(int prune) {
504     energyGroup.prune = prune;
505   }
506 
507   /**
508    * Revert unfavorable changes.
509    *
510    * @return Returns true if unfavorable changes are reverted.
511    */
512   public boolean isRevert() {
513     return group.revert;
514   }
515 
516   public void setRevert(boolean revert) {
517     group.revert = revert;
518   }
519 
520   /**
521    * Energy restart file from a previous run (requires that all parameters are the same).
522    *
523    * @return Returns the energy restart file to use.
524    */
525   public String getEnergyRestart() {
526     return group.energyRestart;
527   }
528 
529   public void setEnergyRestart(String energyRestart) {
530     group.energyRestart = energyRestart;
531   }
532 
533   /**
534    * Do not include starting coordinates as their own rotamer.
535    *
536    * @return Returns true if original side-chain coordinates should not be used as a rotamer.
537    */
538   public boolean isNoOriginal() {
539     return group.noOriginal;
540   }
541 
542   public void setNoOriginal(boolean noOriginal) {
543     group.noOriginal = noOriginal;
544   }
545 
546   /**
547    * -E or --decompose Print energy decomposition for the input structure (no optimization).
548    *
549    * @return Returns true if the input structure should undergo an energy decomposition.
550    */
551   public boolean isDecompose() {
552     return group.decompose;
553   }
554 
555   public void setDecompose(boolean decompose) {
556     group.decompose = decompose;
557   }
558 
559   /**
560    * Choose a list of individual residues to optimize (eg. A11,A24,B40).
561    *
562    * @return Returns the list of selected residues.
563    */
564   public String getListResidues() {
565     return residueGroup.listResidues;
566   }
567 
568   public void setListResidues(String listResidues) {
569     residueGroup.listResidues = listResidues;
570   }
571 
572   /**
573    * Follow elimination criteria with 'n' Monte Carlo steps, or enumerate all remaining
574    * conformations, whichever is smaller.
575    *
576    * @return Returns the number of Monte Carlo optimization steps to apply.
577    */
578   public int getMonteCarlo() {
579     return group.monteCarlo;
580   }
581 
582   public void setMonteCarlo(int monteCarlo) {
583     group.monteCarlo = monteCarlo;
584   }
585 
586   /**
587    * Save eliminated singles and eliminated pairs to a text file (global and box optimization).
588    *
589    * @return Returns true to Save eliminated rotamers to a file.
590    */
591   public boolean isSaveOutput() {
592     return group.saveOutput;
593   }
594 
595   public void setSaveOutput(boolean saveOutput) {
596     group.saveOutput = saveOutput;
597   }
598 
599   /**
600    * Size of the sliding window with respect to adjacent residues (default = 7).
601    *
602    * @return Returns the sliding window size.
603    */
604   public int getWindow() {
605     return windowGroup.window;
606   }
607 
608   public void setWindow(int window) {
609     windowGroup.window = window;
610   }
611 
612   /**
613    * Sliding window increment (default = 3).
614    *
615    * @return Returns the sliding window increment.
616    */
617   public int getIncrement() {
618     return windowGroup.increment;
619   }
620 
621   public void setIncrement(int increment) {
622     windowGroup.increment = increment;
623   }
624 
625   /**
626    * The sliding window and box cutoff radius (Angstroms).
627    *
628    * @return Returns the sliding window cutoff radius.
629    */
630   public double getCutoff() {
631     return energyGroup.cutoff;
632   }
633 
634   public void setCutoff(double cutoff) {
635     energyGroup.cutoff = cutoff;
636   }
637 
638   /**
639    * The threshold for pruning clashes. If two self-energies on the same residue have an energy
640    * difference greater than 25 kcal/mol, the high energy rotamers get pruned.
641    *
642    * @return Returns the clash threshold for self energies.
643    */
644   public double getClashThreshold() {
645     return energyGroup.clashThreshold;
646   }
647 
648   public void setClashThreshold(double clashThreshold) {
649     energyGroup.clashThreshold = clashThreshold;
650   }
651 
652   /**
653    * The threshold for pruning clashes. If two pair-energies on the same residues have an energy
654    * difference greater than 25 kcal/mol, the high energy rotamers get pruned.
655    *
656    * @return Returns the clash threshold for pair energies.
657    */
658   public double getPairClashThreshold() {
659     return energyGroup.pairClashThreshold;
660   }
661 
662   public void setPairClashThreshold(double pairClashThreshold) {
663     energyGroup.pairClashThreshold = pairClashThreshold;
664   }
665 
666   /**
667    * The number of boxes along X, Y, and Z (default: '3,3,3').
668    *
669    * @return Returns the number of boxes.
670    */
671   public String getNumBoxes() {
672     return boxGroup.numBoxes;
673   }
674 
675   public void setNumBoxes(String numBoxes) {
676     boxGroup.numBoxes = numBoxes;
677   }
678 
679   /**
680    * Extent of overlap between optimization boxes (default: 0.0 A).
681    *
682    * @return Returns the overlap between optimization boxes.
683    */
684   public double getBoxBorderSize() {
685     return boxGroup.boxBorderSize;
686   }
687 
688   public void setBoxBorderSize(double boxBorderSize) {
689     boxGroup.boxBorderSize = boxBorderSize;
690   }
691 
692   /**
693    * Approximate side lengths of boxes to be constructed (over-rides numXYZBoxes). Box sizes are
694    * rounded up to make a whole number of boxes along each axis (default of 0 disables this
695    * function).
696    *
697    * @return Returns the approximate box length.
698    */
699   public double getApproxBoxLength() {
700     return boxGroup.approxBoxLength;
701   }
702 
703   public void setApproxBoxLength(double approxBoxLength) {
704     boxGroup.approxBoxLength = approxBoxLength;
705   }
706 
707   /**
708    * Criterion to use for adding residues to boxes. (1) uses C alpha only (N1/9 for nucleic acids)
709    * (2) uses any atom. (3) uses any rotamer.
710    *
711    * @return Returns the Box inclusion criteria.
712    */
713   public int getBoxInclusionCriterion() {
714     return boxGroup.boxInclusionCriterion;
715   }
716 
717   public void setBoxInclusionCriterion(int boxInclusionCriterion) {
718     boxGroup.boxInclusionCriterion = boxInclusionCriterion;
719   }
720 
721   public void setBoxTitration(boolean boxTitration){boxGroup.boxTitration = boxTitration;}
722 
723   public boolean getBoxTitration(){return boxGroup.boxTitration;}
724 
725   public void setTitrationPH(double pH) {
726     group.titrationPH = pH;
727   }
728 
729   public double getTitrationPH() {
730     return group.titrationPH;
731   }
732 
733   public void setPHRestraint(double pHRestraint) {
734     energyGroup.pHRestraint = pHRestraint;
735   }
736 
737   public double getPHRestraint() {
738     return energyGroup.pHRestraint;
739   }
740 
741   public boolean isTitrating() {
742     return group.titrationPH == 0;
743   }
744 
745   public boolean getTitration() {
746     return group.titration;
747   }
748 
749   public String selectInclusionResidues(final List<Residue> residueList, int mutatingResidue, boolean onlyTitration,
750                                        double inclusionCutoff){
751     String listResidues = "";
752     if (mutatingResidue != -1 && inclusionCutoff != -1) {
753       List<Integer> residueNumber = new ArrayList<>();
754       for (Residue residue : residueList) {
755         residueNumber.add(residue.getResidueNumber());
756       }
757       double[] mutatingResCoor = new double[3];
758       int index = residueNumber.indexOf(mutatingResidue);
759       mutatingResCoor = residueList.get(index).getAtomByName("CA", true).getXYZ(mutatingResCoor);
760       for (Residue residue: residueList) {
761         double[] currentResCoor = new double[3];
762         currentResCoor = residue.getAtomByName("CA", true).getXYZ(currentResCoor);
763         double dist = DoubleMath.dist(mutatingResCoor, currentResCoor);
764         if (dist < inclusionCutoff) {
765           listResidues += "," + residue.getChainID() + residue.getResidueNumber();
766         }
767       }
768       listResidues = listResidues.substring(1);
769     } else if (onlyTitration){
770       String[] titratableResidues = new String[]{"HIS", "HIE", "HID", "GLU", "GLH", "ASP", "ASH", "LYS", "LYD", "CYS", "CYD"};
771       List<String> titratableResiudesList = Arrays.asList(titratableResidues);
772       for (Residue residue : residueList) {
773         if (titratableResiudesList.contains(residue.getName())) {
774           String titrateResNum = Integer.toString(residue.getResidueNumber());
775           if(!listResidues.contains(titrateResNum)){
776             listResidues += "," + residue.getChainID()+ titrateResNum;
777           }
778           if (inclusionCutoff != -1){
779             for (Residue residue2: residueList) {
780               boolean includeResidue = evaluateAllRotDist(residue, residue2, inclusionCutoff);
781               if(includeResidue){
782                 String residue2Number = Integer.toString(residue2.getResidueNumber());
783                 if(!listResidues.contains(residue2Number)){
784                   listResidues += "," + residue2.getChainID()+ residue2Number;
785                 }
786               }
787             }
788           }
789 
790         }
791 
792       }
793 
794       listResidues = listResidues.substring(1);
795     }
796     return listResidues;
797   }
798   private static boolean evaluateAllRotDist(Residue residueA, Residue residueB, double inclusionCutoff){
799     residueA.setRotamers(RotamerLibrary.getDefaultLibrary());
800     residueB.setRotamers(RotamerLibrary.getDefaultLibrary());
801     Rotamer[] rotamersA = residueA.getRotamers();
802     Rotamer[] rotamersB = residueB.getRotamers();
803     double[] aCoor = new double[3];
804     double[] bCoor = new double[3];
805     try {
806       int a = rotamersA.length;
807       int b = rotamersB.length;
808     } catch (Exception e){
809       return false;
810     }
811 
812     for(Rotamer rotamerA: rotamersA){
813       residueA.setRotamer(rotamerA);
814       for(Rotamer rotamerB: rotamersB){
815         residueB.setRotamer(rotamerB);
816         for(Atom atomA: residueA.getAtomList()){
817           for(Atom atomB: residueB.getAtomList()){
818             double dist = DoubleMath.dist(atomA.getXYZ(aCoor), atomB.getXYZ(bCoor));
819             if(dist <= inclusionCutoff){
820               return true;
821             }
822           }
823         }
824       }
825     }
826 
827     return false;
828   }
829 
830 
831   /**
832    * Collection of ManyBody Options.
833    */
834   private static class ManyBodyOptionGroup {
835 
836     /**
837      * -a or --algorithm Choices are independent residues (1), all with rotamer elimination (2), all
838      * brute force (3), sliding window (4), or box optimization (5).
839      */
840     @Option(names = {"-a",
841         "--algorithm"}, paramLabel = "0", defaultValue = "0", description = "Algorithm: default automatic settings (0), independent residues (1), all with rotamer elimination (2), all brute force (3), sliding window (4), or box optimization (5)")
842     private int algorithm;
843 
844     /** -L or --library Choose either Ponder and Richards (1) or Richardson (2) rotamer library. */
845     @Option(names = {"-L",
846         "--library"}, paramLabel = "2", defaultValue = "2", description = "Ponder and Richards (1) or Richardson (2) rotamer library.")
847     private int library;
848 
849     /** -nl or --nucleicLibrary Choose a nucleic acid library: currently only Richardson available. */
850     // @Option(
851     //    names = {"--nl", "--nucleiclibrary"},
852     //    paramLabel = "Richardson",
853     //    defaultValue = "Richardson",
854     //    description = "Nucleic acid library to select: [Richardson]")
855     private String naLibraryName = "Richardson";
856 
857     /** --dee or --deadEnd Use dead-end elimination criteria instead of Goldstein criteria. */
858     @Option(names = {"--dee",
859         "--deadEnd"}, paramLabel = "false", defaultValue = "false", description = "Use dead-end elimination criteria instead of Goldstein criteria.")
860     private boolean dee;
861 
862     /** -z or --revert Revert unfavorable changes. */
863     @Option(names = {"-z",
864         "--revert"}, defaultValue = "true", description = "Revert unfavorable changes.")
865     private boolean revert;
866 
867     /**
868      * --eR or --energyRestart Load energy restart file from a previous run (requires that all
869      * parameters are the same).
870      */
871     @Option(names = {"--eR",
872         "--energyRestart"}, paramLabel = "none", defaultValue = "none", description = "Load energy restart file from a previous run (requires that all parameters are the same).")
873     private String energyRestart;
874 
875     /** -O or --noOriginal Do not include starting coordinates as their own rotamer. */
876     @Option(names = {"-O",
877         "--noOriginal"}, defaultValue = "false", description = "Do not include starting coordinates as their own rotamer.")
878     private boolean noOriginal;
879 
880     /** -E or --decompose Print energy decomposition for the input structure (no optimization). */
881     @Option(names = {"-E",
882         "--decompose"}, defaultValue = "false", description = "Print energy decomposition for the input structure (no optimization!).")
883     private boolean decompose;
884 
885     /**
886      * --pH or --titrationPH Optimize the titration state of ASP, GLU, HIS and LYS residues.
887      */
888     @Option(names = {"--pH",
889         "--titrationPH"}, paramLabel = "0", defaultValue = "0", description = " Optimize the titration state of ASP, GLU, HIS and LYS residues at the given pH (pH = 0 turns off titration")
890     private double titrationPH;
891 
892     /**
893      * --pH or --titrationPH Optimize the titration state of ASP, GLU, HIS and LYS residues.
894      */
895     @Option(names = {"--tR",
896             "--titration"}, paramLabel = "false", defaultValue = "false", description = " Turn on titration state optimization")
897     private boolean titration;
898 
899     /**
900      * --mC or --monteCarlo Follow elimination criteria with 'n' Monte Carlo steps, or enumerate all
901      * remaining conformations, whichever is smaller.
902      */
903     @Option(names = {"--mC",
904         "--monteCarlo"}, paramLabel = "-1", defaultValue = "-1", description = "Follow elimination criteria with (n) Monte Carlo steps, or enumerate all remaining conformations, whichever is smaller.")
905     private int monteCarlo;
906 
907     /**
908      * -out or --output Save eliminated singles and eliminated pairs to a text file (global and box
909      * optimization).
910      */
911     //  @Option(
912     //      names = {"--out", "--output"},
913     //      defaultValue = "false",
914     //      description = "Save eliminated singles and eliminated pairs to a text file.")
915     private boolean saveOutput;
916 
917   }
918 
919   /**
920    * Collection of ManyBody Box Optimization Options.
921    */
922   private static class BoxOptionGroup {
923 
924     /** -nB or --numBoxes Specify number of boxes along X, Y, and Z (default: '3,3,3'). */
925     @Option(names = {"--nB",
926         "--numBoxes"}, paramLabel = "3,3,3", defaultValue = "3,3,3", description = "Specify number of boxes along X, Y, and Z (default: 3,3,3)")
927     private String numBoxes;
928 
929     /**
930      * Result of parsing numBoxes flag.
931      */
932     private int[] numXYZBoxes;
933 
934     /** -bB or --boxBorderSize Extent of overlap between optimization boxes (default: 0.0 A). */
935     @Option(names = {"--bB",
936         "--boxBorderSize"}, paramLabel = "0.0", defaultValue = "0.0", description = "Extent of overlap between optimization boxes in Angstroms.")
937     private double boxBorderSize;
938 
939     /**
940      * -bL or --approxBoxLength Approximate side lengths of boxes to be constructed (over-rides
941      * numXYZBoxes). Box sizes are rounded up to make a whole number of boxes along each axis
942      * (default of 0 disables this function).
943      */
944     @Option(names = {"--bL",
945         "--approxBoxLength"}, paramLabel = "20.0", defaultValue = "20.0", description = "Approximate side lengths of boxes to be constructed (over-rides numXYZBoxes).")
946     private double approxBoxLength;
947 
948     /**
949      * -bC or --boxInclusionCriterion Criterion to use for adding residues to boxes. (1) uses C alpha
950      * only (N1/9 for nucleic acids) (2) uses any atom. (3) uses any rotamer
951      */
952     @Option(names = {"--bC",
953         "--boxInclusionCriterion"}, paramLabel = "1", defaultValue = "1", description = "Criterion to use for adding a residue to a box: (1) uses C alpha only (N1/9 for nucleic acids), (2) uses any atom, and (3) uses any rotamer")
954     private int boxInclusionCriterion;
955 
956     /**
957      * --iB or --initialBox Initial box to optimize.
958      */
959     @Option(names = {"--iB", "--initialBox"}, defaultValue = "0", description = "Initial box to optimize.")
960     private int initialBox;
961 
962     /**
963      * --fB or --boxFinal Final box to optimize.
964      */
965     @Option(names = {"--fB", "--finalBox"}, defaultValue = "2147483647", // Integer.MAX_VALUE
966         description = "Final box to optimize.")
967     private int finalBox;
968 
969     /**
970      * --bT or --boxTitration Center boxes around titratable residues.
971      */
972     @Option(names = {"--bT", "--boxTitration"}, defaultValue = "false", // Integer.MAX_VALUE
973             description = "Center boxes around titratable residues.")
974     private boolean boxTitration;
975 
976   }
977 
978   /**
979    * Collection of ManyBody Window Optimization Options.
980    */
981   private static class WindowOptionGroup {
982 
983     /** --window Size of the sliding window with respect to adjacent residues (default = 7). */
984     @Option(names = {
985         "--window"}, paramLabel = "7", defaultValue = "7", description = "Size of the sliding window with respect to adjacent residues.")
986     private int window;
987 
988     /** --increment Sliding window increment (default = 3). */
989     @Option(names = {
990         "--increment"}, paramLabel = "3", defaultValue = "3", description = "Sliding window increment.")
991     private int increment;
992 
993   }
994 
995   /**
996    * Collection of ManyBody Energy Optimization Options.
997    */
998   private static class EnergyOptionGroup {
999 
1000     /** --radius The sliding window and box cutoff radius (Angstroms). */
1001     @Option(names = {
1002         "--radius"}, paramLabel = "2.0", defaultValue = "2.0", description = "The sliding box and window cutoff radius (Angstroms).")
1003     private double cutoff;
1004 
1005     /** --tC or --twoBodyCutoff Cutoff distance for two-body interactions. */
1006     @Option(names = {"--tC",
1007         "--twoBodyCutoff"}, paramLabel = "3.0", defaultValue = "3.0", description = "Cutoff distance for two body interactions.")
1008     private double twoBodyCutoff;
1009 
1010     /** -T or --threeBody Include 3-Body interactions in the elimination criteria. */
1011     @Option(names = {"-T",
1012         "--threeBody"}, defaultValue = "false", description = "Include 3-Body interactions in the elimination criteria.")
1013     private boolean threeBody;
1014 
1015     /** --thC or --threeBodyCutoff Cutoff distance for three-body interactions. */
1016     @Option(names = {"--thC",
1017         "--threeBodyCutoff"}, paramLabel = "3.0", defaultValue = "3.0", description = "Cutoff distance for three-body interactions.")
1018     private double threeBodyCutoff;
1019 
1020     /** --pr or --prune Prune no clashes (0), only single clashes (1), or all clashes (2). */
1021     @Option(names = {"--pr",
1022         "--prune"}, paramLabel = "1", defaultValue = "1", description = "Prune no clashes (0), only single clashes (1), or all clashes (2)")
1023     private int prune;
1024 
1025     /**
1026      * --clashThreshold The threshold for pruning clashes. If two self-energies on the same residue
1027      * have an energy difference greater than 25 kcal/mol, the high energy rotamers get pruned.
1028      */
1029     @Option(names = {
1030         "--clashThreshold"}, paramLabel = "25.0", defaultValue = "25.0", description = "The threshold for pruning clashes.")
1031     private double clashThreshold;
1032 
1033     /**
1034      * --pairClashThreshold The threshold for pruning clashes. If two pair-energies on the same
1035      * residues have an energy difference greater than 25 kcal/mol, the high energy rotamers get
1036      * pruned.
1037      */
1038     @Option(names = {
1039         "--pairClashThreshold"}, paramLabel = "25.0", defaultValue = "25.0", description = "The threshold for pruning pair clashes.")
1040     private double pairClashThreshold;
1041 
1042     /** --radius The sliding window and box cutoff radius (Angstroms). */
1043     @Option(names = {
1044             "--kPH", "--pHRestraint"}, paramLabel = "0.0", defaultValue = "0.0", description = "Only allow titration state to change from" +
1045             "standard state is self energy exceeds the restraint.")
1046     private double pHRestraint = 0;
1047   }
1048 
1049   /**
1050    * Collection of ManyBody Residue Selection Options.
1051    */
1052   private static class ResidueOptionGroup {
1053 
1054     /** --ch or --chain Single character chain ID of the residues to optimize. */
1055     @Option(names = {"--ch",
1056         "--chain"}, paramLabel = "<A>", defaultValue = "-1", description = "Include only specified chain ID (default: all chains).")
1057     private String chain;
1058 
1059     /**
1060      * --sR or --start Starting residue to perform the optimization on.
1061      */
1062     @Option(names = {"--sR", "--start"}, defaultValue = "-2147483648", // Integer.MIN_VALUE
1063         description = "Starting residue to optimize (default: all residues).")
1064     private int start;
1065 
1066     /**
1067      * --fi or --final Final residue to perform the optimization.
1068      */
1069     @Option(names = {"--fR",
1070         "--final"}, paramLabel = "<final>", defaultValue = "2147483647", // Integer.MAX_VALUE
1071         description = "Final residue to optimize (default: all residues).")
1072     private int finish;
1073 
1074     /** --lR or --listResidues Choose a list of individual residues to optimize (eg. A11,A24,B40). */
1075     @Option(names = {"--lR",
1076         "--listResidues"}, paramLabel = "<list>", defaultValue = "none", description = "Select a list of residues to optimize (eg. A11,A24,B40).")
1077     private String listResidues;
1078 
1079     /** --oT or --onlyTitrtaion Rotamer optimize only titratable residues. */
1080     @Option(names = {"--oT",
1081             "--onlyTitration"}, paramLabel = "", defaultValue = "false", description = "Rotamer optimize only titratable residues.")
1082     private boolean onlyTitration;
1083 
1084     /** --iR or --interestedResidue Optimize rotamers within some distance of a specific residue. */
1085     @Option(names = {"--iR",
1086             "--interestedResidue"}, paramLabel = "", defaultValue = "-1", description = "Optimize rotamers within some distance of a specific residue.")
1087     private int interestedResidue = -1;
1088 
1089     /** --iC or --inclusionCutoff Distance which rotamers will be included when using only protons, titratable residues, or interested residue. */
1090     @Option(names = {"--iC",
1091             "--inclusionCutoff"}, paramLabel = "", defaultValue = "-1", description = "Distance which rotamers will be included when using only protons, titratable residues, or interested residue.")
1092     private double inclusionCutoff = -1;
1093 
1094   }
1095 
1096 }