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.misc;
39  
40  import ffx.algorithms.AlgorithmListener;
41  import ffx.numerics.Potential;
42  import ffx.potential.MolecularAssembly;
43  import ffx.potential.bonded.AminoAcidUtils;
44  import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
45  import ffx.potential.bonded.Atom;
46  import ffx.potential.bonded.NucleicAcidUtils.NucleicAcid3;
47  import ffx.potential.bonded.Residue;
48  import ffx.potential.bonded.Rotamer;
49  import ffx.potential.bonded.RotamerLibrary;
50  import ffx.potential.parsers.PDBFilter;
51  
52  import java.io.BufferedWriter;
53  import java.io.File;
54  import java.io.FileWriter;
55  import java.io.IOException;
56  import java.util.Arrays;
57  import java.util.logging.Level;
58  import java.util.logging.Logger;
59  import java.util.regex.Matcher;
60  import java.util.regex.Pattern;
61  
62  import static java.lang.String.format;
63  
64  /**
65   * The GenerateRotamers class helps generate a rotamer library (particularly for nonstandard amino
66   * acids) for a Residue.
67   *
68   * @author Jacob M. Litman
69   * @author Michael J. Schnieders
70   */
71  public class GenerateRotamers {
72  
73    private static final Logger logger = Logger.getLogger(GenerateRotamers.class.getName());
74    private static final Pattern atRangePatt = Pattern.compile("(\\d+)-(\\d+)");
75    private final Residue residue;
76    private final int nChi;
77    private final double[] currentChi;
78    private final File outFile;
79    private final MolecularAssembly molecularAssembly;
80    private final AlgorithmListener listener;
81    private final RotamerLibrary library;
82    private final Potential potential;
83    private double[] x;
84    private AminoAcid3 baselineAAres;
85    private Rotamer[] baselineRotamers;
86    private double incr = 10.0; // Degrees to rotate each torsion by.
87    private boolean aroundLibrary = false;
88    private double width = 180.0; // Left 180 + right 180 = 360 degree coverage.
89    private int startDepth = 0;
90    private int endDepth;
91    private boolean print = false; // If true, log energy at each torsion.
92    private int nEval = 0;
93    private boolean writeVideo = false;
94    private File videoFile;
95    private PDBFilter videoFilter;
96  
97    /**
98     * Intended to create rotamer sets for nonstandard amino acids.
99     *
100    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
101    * @param potential a {@link ffx.numerics.Potential} object.
102    * @param residue a {@link ffx.potential.bonded.Residue} object.
103    * @param file Output file
104    * @param nChi Number of rotameric torsions in the residue
105    * @param listener a {@link ffx.algorithms.AlgorithmListener} object.
106    */
107   public GenerateRotamers(MolecularAssembly molecularAssembly, Potential potential, Residue residue,
108       File file, int nChi, AlgorithmListener listener) {
109     this(molecularAssembly, potential, residue, file, nChi, listener,
110         RotamerLibrary.getDefaultLibrary());
111   }
112 
113   /**
114    * Intended to create rotamer sets for nonstandard amino acids.
115    *
116    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
117    * @param potential a {@link ffx.numerics.Potential} object.
118    * @param residue a {@link ffx.potential.bonded.Residue} object.
119    * @param file Output file
120    * @param nChi Number of rotameric torsions in the residue
121    * @param listener a {@link ffx.algorithms.AlgorithmListener} object.
122    * @param library Rotamer library to use
123    */
124   public GenerateRotamers(MolecularAssembly molecularAssembly, Potential potential, Residue residue,
125       File file, int nChi, AlgorithmListener listener, RotamerLibrary library) {
126     this.residue = residue;
127     this.nChi = nChi;
128     endDepth = nChi - 1;
129     this.potential = potential;
130     this.molecularAssembly = molecularAssembly;
131     this.listener = listener;
132     this.currentChi = new double[nChi];
133     this.library = library;
134     Arrays.fill(currentChi, 0.0);
135     File outputFile = file;
136     if (outputFile.exists()) {
137       String outName = outputFile.getName();
138       for (int i = 1; i < 1000; i++) {
139         outputFile = new File(format("%s_%d", outName, i));
140         if (!outputFile.exists()) {
141           break;
142         }
143       }
144       if (outputFile.exists()) {
145         logger.severe(format(" Could not version file %s", outName));
146       }
147     }
148     outFile = outputFile;
149     this.baselineAAres = residue.getAminoAcid3();
150     baselineRotamers = library.getRotamers(baselineAAres);
151   }
152 
153   /**
154    * Accessory method for more simplistic saving of specific torsion states.
155    *
156    * @param torSets an array of {@link java.lang.String} objects.
157    */
158   public void applyAndSaveTorsions(String[] torSets) {
159     for (String torSet : torSets) {
160       String[] torsions = torSet.split(",");
161       double[] values = new double[nChi * 2];
162       Arrays.fill(values, 0.0);
163       for (int i = 0; i < (Math.min(torsions.length, nChi)); i++) {
164         double chival = Double.parseDouble(torsions[i]);
165         currentChi[i] = chival;
166         values[2 * i] = chival;
167       }
168       Rotamer newRot = generateRotamer(values);
169       RotamerLibrary.applyRotamer(residue, newRot);
170       writeSnapshot();
171     }
172   }
173 
174   /**
175    * Sets a standard amino acid to be the baseline for rotamer generation. For example, use TYR as a
176    * baseline for phospho-tyrosine.
177    *
178    * @param aa3 a {@link AminoAcid3} object.
179    */
180   public void setBaselineAARes(AminoAcid3 aa3) {
181     this.baselineAAres = aa3;
182     if (aa3 == AminoAcidUtils.AminoAcid3.UNK) {
183       boolean orig = library.getUsingOrigCoordsRotamer();
184       library.setUseOrigCoordsRotamer(false);
185       baselineRotamers = residue.setRotamers(library);
186       library.setUseOrigCoordsRotamer(orig);
187     } else {
188       baselineRotamers = library.getRotamers(aa3);
189     }
190     aroundLibrary = true;
191   }
192 
193   /**
194    * Set which torsions to work on. Negative end values set the final depth to be the total number of
195    * torsions.
196    *
197    * @param start The first torsion to work on.
198    * @param end The last torsion to work on.
199    */
200   public void setDepth(int start, int end) {
201     startDepth = start;
202     if (end < 1 || end > nChi) {
203       endDepth = nChi - 1;
204     } else {
205       endDepth = end;
206     }
207   }
208 
209   /**
210    * Inactivates electrostatics for atom sets defined by 'start-end,start-end,...'.
211    *
212    * @param electrostatics Input string
213    */
214   public void setElectrostatics(String electrostatics) {
215     if (electrostatics != null) {
216       String[] tokens = electrostatics.split(",");
217       Atom[] atoms = molecularAssembly.getAtomArray();
218       for (String tok : tokens) {
219         Matcher m = atRangePatt.matcher(tok);
220         if (m.matches()) {
221           int begin = Integer.parseInt(m.group(1));
222           int end = Integer.parseInt(m.group(2));
223           logger.info(format(" Inactivating electrostatics for atoms %d-%d", begin, end));
224           for (int i = begin; i <= end; i++) {
225             Atom ai = atoms[i - 1];
226             ai.setElectrostatics(false);
227             ai.print(Level.FINE);
228           }
229         } else {
230           logger.info(format(" Discarding electrostatics input %s", tok));
231         }
232       }
233     }
234   }
235 
236   /**
237    * Inactivates atom sets defined by 'start-end,start-end,...'.
238    *
239    * @param iatoms Input string
240    */
241   public void setInactiveAtoms(String iatoms) {
242     if (iatoms != null) {
243       String[] tokens = iatoms.split(",");
244       Atom[] atoms = molecularAssembly.getAtomArray();
245       for (String tok : tokens) {
246         Matcher m = atRangePatt.matcher(tok);
247         if (m.matches()) {
248           int begin = Integer.parseInt(m.group(1));
249           int end = Integer.parseInt(m.group(2));
250           logger.info(format(" Inactivating atoms %d-%d", begin, end));
251           for (int i = begin; i <= end; i++) {
252             Atom ai = atoms[i - 1];
253             ai.setUse(false);
254             ai.print(Level.FINE);
255           }
256         } else {
257           logger.info(format(" Discarding inactive atoms input %s", tok));
258         }
259       }
260     }
261   }
262 
263   /**
264    * Sets the angle to change torsions by.
265    *
266    * @param incr a double.
267    */
268   public void setIncrement(double incr) {
269     this.incr = incr;
270   }
271 
272   /**
273    * Sets algorithm to log all torsions/energies (not just to file).
274    *
275    * @param print a boolean.
276    */
277   public void setPrint(boolean print) {
278     this.print = print;
279   }
280 
281   /**
282    * Sets the width around each torsion to search (+/-, so 10 degree width will search a 20 degree
283    * arc).
284    *
285    * @param width a double.
286    */
287   public void setSearchWidth(double width) {
288     this.width = width;
289   }
290 
291   /**
292    * Null file indicates to not write a video.
293    *
294    * @param videoFile Filename for video or null
295    */
296   public void setVideo(String videoFile) {
297     if (videoFile != null) {
298       File vidFile = new File(videoFile);
299       if (vidFile.exists()) {
300         for (int i = 0; i < 1000; i++) {
301           vidFile = new File(format("%s_%d", videoFile, i));
302           if (!vidFile.exists()) {
303             this.videoFile = vidFile;
304             writeVideo = true;
305             videoFilter = new PDBFilter(this.videoFile, molecularAssembly,
306                 molecularAssembly.getForceField(), null);
307             videoFilter.setLogWrites(false);
308             break;
309           }
310         }
311         if (vidFile.exists()) {
312           logger.warning(format(" Could not version video file %s", videoFile));
313         }
314       } else {
315         this.videoFile = vidFile;
316         writeVideo = true;
317         videoFilter = new PDBFilter(this.videoFile, molecularAssembly,
318             molecularAssembly.getForceField(), null);
319         videoFilter.setLogWrites(false);
320       }
321     } else {
322       writeVideo = false;
323     }
324   }
325 
326   /** Main driver method; spins torsions, evaluates energy, and prints to file. */
327   public void tryRotamers() {
328     try (BufferedWriter bw = new BufferedWriter(new FileWriter(outFile))) {
329       if (aroundLibrary) {
330         for (Rotamer rotamer : baselineRotamers) {
331           applyLibRotamer(rotamer);
332           turnChi(startDepth, bw);
333         }
334       } else {
335         RotamerLibrary.measureRotamer(residue, currentChi, false);
336         turnChi(startDepth, bw);
337       }
338     } catch (IOException ex) {
339       logger.warning(format(" IO exception in rotamer generation: %s", ex));
340     }
341   }
342 
343   /**
344    * Applies a library rotamer, filling in the end with 0s as necessary. For example, PTR requires
345    * more torsions than TYR, so one cannot simply apply a TYR rotamer to a PTR residue.
346    *
347    * @param rotamer The rotamer to apply.
348    */
349   private void applyLibRotamer(Rotamer rotamer) {
350     double[] rotValues = new double[nChi * 2];
351     Arrays.fill(rotValues, 0.0);
352     Arrays.fill(currentChi, 0.0);
353     int fillTo = Math.min(startDepth, (rotamer.length - 1));
354     for (int i = 0; i < fillTo; i++) {
355       int ii = 2 * i;
356       rotValues[ii] = rotamer.angles[i];
357       currentChi[i] = rotamer.angles[i];
358       rotValues[ii + 1] = rotamer.sigmas[i];
359     }
360     Rotamer newRot = generateRotamer(rotValues);
361     RotamerLibrary.applyRotamer(residue, newRot);
362   }
363 
364   /**
365    * Generates an aa/na/unk Rotamer for the selected residue given values.
366    *
367    * @param values The values to use.
368    * @return The Rotamer.
369    */
370   private Rotamer generateRotamer(double[] values) {
371     switch (residue.getResidueType()) {
372       case AA -> {
373         AminoAcid3 aa3 = residue.getAminoAcid3();
374         return new Rotamer(aa3, values);
375       }
376       case NA -> {
377         NucleicAcid3 na3 = residue.getNucleicAcid3();
378         return new Rotamer(na3, values);
379       }
380       default -> {
381         return new Rotamer(values);
382       }
383     }
384   }
385 
386   /**
387    * Recursive method for turning the torsions.
388    *
389    * @param depth Current depth of the recursion
390    * @param bw The BufferedWriter to write to.
391    * @throws IOException If there is an IO exception.
392    */
393   private void turnChi(int depth, BufferedWriter bw) throws IOException {
394     double chi = currentChi[depth];
395     for (double i = chi - width; i <= chi + width; i += incr) {
396       currentChi[depth] = i;
397       if (depth == endDepth) {
398         evaluateChi(bw);
399         if (listener != null) {
400           listener.algorithmUpdate(molecularAssembly);
401         }
402         if (writeVideo) {
403           writeSnapshot();
404         }
405       } else {
406         turnChi(depth + 1, bw);
407       }
408     }
409     currentChi[depth] = chi; // Reset the chi value to where it was.
410   }
411 
412   private void writeSnapshot() {
413     try (BufferedWriter bw = new BufferedWriter(new FileWriter(videoFile, true))) {
414       bw.write(format("MODEL %d", nEval));
415       bw.newLine();
416       bw.write(format("REMARK 301 TORSIONS %s", formatChi()));
417       bw.newLine();
418       bw.flush();
419       videoFilter.writeFile(videoFile, true);
420       bw.write("ENDMDL");
421       bw.newLine();
422       bw.flush();
423     } catch (IOException ex) {
424       logger.warning(
425           format(" Exception writing to video file %s: " + "%s", videoFile.getName(), ex));
426     }
427   }
428 
429   /**
430    * Called at a leaf of the recursion.
431    *
432    * @param bw The BufferedWriter to write to.
433    * @throws IOException If there is an IO exception.
434    */
435   private void evaluateChi(BufferedWriter bw) throws IOException {
436     try {
437       applyChi();
438       double e = currentEnergy();
439       String result = format("%s,%12f", formatChi(), e);
440       ++nEval;
441       if (print) {
442         logger.info(format(" Evaluation %10d %s, energy %10.5f kcal/mol", nEval, formatChi(), e));
443       }
444       bw.write(result);
445       bw.newLine();
446       if (nEval % 1000 == 0) {
447         logger.info(format(" %12.7e states evaluated", (double) nEval));
448       }
449     } catch (ArithmeticException ex) {
450       logger.info(format(" Force field exception at chi %s", formatChi()));
451     }
452   }
453 
454   /** Converts the chi array to a Rotamer and applies it. */
455   private void applyChi() {
456     double[] rotValues = new double[nChi * 2];
457     for (int i = 0; i < nChi; i++) {
458       int ii = 2 * i;
459       rotValues[ii] = currentChi[i];
460       rotValues[ii + 1] = 0.0;
461     }
462     RotamerLibrary.applyRotamer(residue, generateRotamer(rotValues));
463   }
464 
465   /**
466    * Returns a formatted String with chi values.
467    *
468    * @return A formatted String with chi values.
469    */
470   private String formatChi() {
471     StringBuilder sb = new StringBuilder(format("%8f", currentChi[0]));
472     for (int i = 1; i < nChi; i++) {
473       // sb.append(",").append(currentChi[i]);
474       sb.append(format(",%8f", currentChi[i]));
475     }
476     return sb.toString();
477   }
478 
479   /**
480    * Calculates the energy at the current state.
481    *
482    * @return Energy of the current state.
483    */
484   private double currentEnergy() throws ArithmeticException {
485     if (x == null) {
486       int nVar = potential.getNumberOfVariables();
487       x = new double[nVar * 3];
488     }
489     potential.getCoordinates(x);
490     return potential.energy(x);
491   }
492 }