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