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.optimize.manybody;
39  
40  import edu.rit.mp.DoubleBuf;
41  import edu.rit.pj.Comm;
42  import edu.rit.pj.IntegerSchedule;
43  import edu.rit.pj.WorkerIntegerForLoop;
44  import edu.rit.pj.WorkerRegion;
45  import ffx.algorithms.optimize.RotamerOptimization;
46  import ffx.potential.bonded.Residue;
47  import ffx.potential.bonded.Rotamer;
48  
49  import java.io.BufferedWriter;
50  import java.io.IOException;
51  import java.util.List;
52  import java.util.Map;
53  import java.util.Set;
54  import java.util.logging.Level;
55  import java.util.logging.Logger;
56  
57  import static java.lang.String.format;
58  import static org.apache.commons.math3.util.FastMath.min;
59  
60  /** Compute 3-Body energy values in parallel across nodes. */
61  public class ThreeBodyEnergyRegion extends WorkerRegion {
62  
63    private static final Logger logger = Logger.getLogger(ThreeBodyEnergyRegion.class.getName());
64    private final Residue[] residues;
65    private final RotamerOptimization rO;
66    private final DistanceMatrix dM;
67    private final EnergyExpansion eE;
68    private final EliminatedRotamers eR;
69    /**
70     * A list of all residues being optimized. Note that Box and Window optimizations operate on
71     * subsets of this list.
72     */
73    private final List<Residue> allResiduesList;
74    /** Map of 3-body energy values to compute. */
75    private final Map<Integer, Integer[]> threeBodyEnergyMap;
76    /** Writes energies to restart file. */
77    private final BufferedWriter energyWriter;
78    /** World Parallel Java communicator. */
79    private final Comm world;
80    /** Number of Parallel Java processes. */
81    private final int numProc;
82    /**
83     * If a pair of residues have two atoms closer together than the superposition threshold, the
84     * energy is set to NaN.
85     */
86    private final double superpositionThreshold;
87    /** Flag to indicate if this is the master process. */
88    private final boolean master;
89    /** Rank of this process. */
90    private final int rank;
91    /** Flag to indicate verbose logging. */
92    private final boolean verbose;
93    /** If true, write out an energy restart file. */
94    private final boolean writeEnergyRestart;
95    /**
96     * Sets whether files should be printed; true for standalone applications, false for some
97     * applications which use rotamer optimization as part of a larger process.
98     */
99    private final boolean printFiles;
100 
101   private Set<Integer> keySet;
102 
103   public ThreeBodyEnergyRegion(RotamerOptimization rotamerOptimization, DistanceMatrix dM,
104       EnergyExpansion eE, EliminatedRotamers eR, Residue[] residues, List<Residue> allResiduesList,
105       BufferedWriter energyWriter, Comm world, int numProc, double superpositionThreshold,
106       boolean master, int rank, boolean verbose, boolean writeEnergyRestart, boolean printFiles) {
107     this.rO = rotamerOptimization;
108     this.dM = dM;
109     this.eE = eE;
110     this.eR = eR;
111     this.residues = residues;
112     this.allResiduesList = allResiduesList;
113     this.energyWriter = energyWriter;
114     this.world = world;
115     this.numProc = numProc;
116     this.superpositionThreshold = superpositionThreshold;
117     this.master = master;
118     this.rank = rank;
119     this.verbose = verbose;
120     this.writeEnergyRestart = writeEnergyRestart;
121     this.printFiles = printFiles;
122 
123     this.threeBodyEnergyMap = eE.getThreeBodyEnergyMap();
124     logger.info(format(" Number of 3-Body energies to calculate: %d", threeBodyEnergyMap.size()));
125   }
126 
127   @Override
128   public void finish() {
129     // Print what we've got so far.
130     if (master && verbose) {
131       for (int i = 0; i < residues.length; i++) {
132         Residue resI = residues[i];
133         Rotamer[] rotI = resI.getRotamers();
134         for (int ri = 0; ri < rotI.length; ri++) {
135           if (eR.check(i, ri)) {
136             continue;
137           }
138           for (int j = i + 1; j < residues.length; j++) {
139             Residue resJ = residues[j];
140             Rotamer[] rotJ = resJ.getRotamers();
141             for (int rj = 0; rj < rotJ.length; rj++) {
142               if (eR.check(j, rj) || eR.check(i, ri, j, rj)) {
143                 continue;
144               }
145               for (int k = j + 1; k < residues.length; k++) {
146                 Residue resK = residues[k];
147                 Rotamer[] rotK = resK.getRotamers();
148                 for (int rk = 0; rk < rotK.length; rk++) {
149                   if (eR.check(k, rk) || eR.check(i, ri, k, rk) || eR.check(j, rj, k, rk)) {
150                     continue;
151                   }
152                   logger.info(format(" 3-Body energy %8s %-2d, %8s %-2d, %8s %-2d: %s",
153                       resI.toString(rotI[ri]), ri, resJ.toString(rotJ[rj]), rj,
154                       resK.toString(rotK[rk]), rk,
155                       rO.formatEnergy(eE.get3Body(residues, i, ri, j, rj, k, rk))));
156                 }
157               }
158             }
159           }
160         }
161       }
162     }
163   }
164 
165   @Override
166   public void run() throws Exception {
167     if (!keySet.isEmpty()) {
168       execute(0, keySet.size() - 1, new ThreeBodyEnergyLoop());
169     }
170   }
171 
172   @Override
173   public void start() {
174     int numTriple = threeBodyEnergyMap.size();
175     int remainder = numTriple % numProc;
176     // Set padded residue and rotamer to less than zero.
177     Integer[] padding = {-1, -1, -1, -1, -1, -1};
178 
179     int padKey = numTriple;
180     while (remainder != 0) {
181       threeBodyEnergyMap.put(padKey++, padding);
182       remainder = threeBodyEnergyMap.size() % numProc;
183     }
184 
185     numTriple = threeBodyEnergyMap.size();
186     if (numTriple % numProc != 0) {
187       logger.severe(" Logic error padding triple energies.");
188     }
189 
190     // Load the keySet of triple energies.
191     keySet = threeBodyEnergyMap.keySet();
192   }
193 
194   private class ThreeBodyEnergyLoop extends WorkerIntegerForLoop {
195 
196     final DoubleBuf[] resultBuffer;
197     final DoubleBuf myBuffer;
198 
199     ThreeBodyEnergyLoop() {
200       resultBuffer = new DoubleBuf[numProc];
201       for (int i = 0; i < numProc; i++) {
202         resultBuffer[i] = DoubleBuf.buffer(new double[7]);
203       }
204       myBuffer = resultBuffer[rank];
205     }
206 
207     @Override
208     public void run(int lb, int ub) {
209       for (int key = lb; key <= ub; key++) {
210         long time = -System.nanoTime();
211         Integer[] job = threeBodyEnergyMap.get(key);
212         int i = job[0];
213         int ri = job[1];
214         int j = job[2];
215         int rj = job[3];
216         int k = job[4];
217         int rk = job[5];
218 
219         myBuffer.put(0, i);
220         myBuffer.put(1, ri);
221         myBuffer.put(2, j);
222         myBuffer.put(3, rj);
223         myBuffer.put(4, k);
224         myBuffer.put(5, rk);
225         myBuffer.put(6, 0.0);
226 
227         // Initialize result.
228         if (i >= 0 && ri >= 0 && j >= 0 && rj >= 0 && k >= 0 && rk >= 0) {
229           if ((!eR.check(i, ri) || !eR.check(j, rj) || !eR.check(k, rk) || !eR.check(i, ri, j, rj)
230               || !eR.check(i, ri, k, rk) || !eR.check(j, rj, k, rk))) {
231 
232             Residue residueI = residues[i];
233             Residue residueJ = residues[j];
234             Residue residueK = residues[k];
235 
236             Rotamer[] rotI = residueI.getRotamers();
237             Rotamer[] rotJ = residueJ.getRotamers();
238             Rotamer[] rotK = residueK.getRotamers();
239 
240             int indexI = allResiduesList.indexOf(residueI);
241             int indexJ = allResiduesList.indexOf(residueJ);
242             int indexK = allResiduesList.indexOf(residueK);
243 
244             double rawDist = dM.getRawNBodyDistance(indexI, ri, indexJ, rj, indexK, rk);
245             double dIJ = dM.checkDistMatrix(indexI, ri, indexJ, rj);
246             double dIK = dM.checkDistMatrix(indexI, ri, indexK, rk);
247             double dJK = dM.checkDistMatrix(indexJ, rj, indexK, rk);
248             double minDist = min(min(dIJ, dIK), dJK);
249 
250             double resDist = dM.get3BodyResidueDistance(indexI, ri, indexJ, rj, indexK, rk);
251             String resDistString = "     large";
252             if (resDist < Double.MAX_VALUE) {
253               resDistString = format("%5.3f", resDist);
254             }
255 
256             String distString = "     large";
257             if (rawDist < Double.MAX_VALUE) {
258               distString = format("%10.3f", rawDist);
259             }
260 
261             double threeBodyEnergy;
262             if (minDist < superpositionThreshold) {
263               threeBodyEnergy = Double.NaN;
264               logger.info(format(
265                   " 3-Body %8s %-2d, %8s %-2d, %8s %-2d:\t    NaN      at %13.6f Ang (%s Ang by residue) < %5.3f Ang.",
266                   residueI.toString(rotI[ri]), ri, residueJ.toString(rotJ[rj]), rj,
267                   residueK.toString(rotK[rk]), rk, minDist, resDistString, superpositionThreshold));
268             } else if (dM.checkTriDistThreshold(indexI, ri, indexJ, rj, indexK, rk)) {
269               // Set the two-body energy to 0.0 for separation distances larger than the two-body
270               // cutoff.
271               threeBodyEnergy = 0.0;
272               time += System.nanoTime();
273               logger.fine(format(
274                   " 3-Body %8s %-2d, %8s %-2d, %8s %-2d: %s at %s Ang (%s Ang by residue) in %6.4f (sec).",
275                   residueI.toString(rotI[ri]), ri, residueJ.toString(rotJ[rj]), rj,
276                   residueK.toString(rotK[rk]), rk, rO.formatEnergy(threeBodyEnergy), distString,
277                   resDistString, time * 1.0e-9));
278             } else {
279               try {
280                 threeBodyEnergy = eE.compute3BodyEnergy(residues, i, ri, j, rj, k, rk);
281                 time += System.nanoTime();
282                 logger.info(format(
283                     " 3-Body %8s %-2d, %8s %-2d, %8s %-2d: %s at %s Ang (%s Ang by residue) in %6.4f (sec).",
284                     residueI.toString(rotI[ri]), ri, residueJ.toString(rotJ[rj]), rj,
285                     residueK.toString(rotK[rk]), rk, rO.formatEnergy(threeBodyEnergy), distString,
286                     resDistString, time * 1.0e-9));
287               } catch (ArithmeticException ex) {
288                 threeBodyEnergy = Double.NaN;
289                 time += System.nanoTime();
290                 logger.info(format(
291                     " 3-Body %8s %-2d, %8s %-2d, %8s %-2d:\t    NaN      at %s Ang (%s Ang by residue) in %6.4f (sec).",
292                     residueI.toString(rotI[ri]), ri, residueJ.toString(rotJ[rj]), rj,
293                     residueK.toString(rotK[rk]), rk, distString, resDistString, time * 1.0e-9));
294               }
295             }
296             myBuffer.put(6, threeBodyEnergy);
297           }
298         }
299 
300         // All to All communication
301         if (numProc > 1) {
302           try {
303             world.allGather(myBuffer, resultBuffer);
304           } catch (Exception e) {
305             logger.log(Level.SEVERE, " Exception communicating pair energies.", e);
306           }
307         }
308 
309         // Process the three-body energy received from each process.
310         for (DoubleBuf doubleBuf : resultBuffer) {
311           int resI = (int) doubleBuf.get(0);
312           int rotI = (int) doubleBuf.get(1);
313           int resJ = (int) doubleBuf.get(2);
314           int rotJ = (int) doubleBuf.get(3);
315           int resK = (int) doubleBuf.get(4);
316           int rotK = (int) doubleBuf.get(5);
317           double energy = doubleBuf.get(6);
318           // Skip for padded result.
319           if (resI >= 0 && rotI >= 0 && resJ >= 0 && rotJ >= 0 && resK >= 0 && rotK >= 0) {
320             if (!Double.isFinite(energy)) {
321               logger.info(
322                   " Rotamer pair eliminated: " + resI + ", " + rotI + ", " + resJ + ", " + rotJ);
323               eR.eliminateRotamerPair(residues, resI, rotI, resJ, rotJ, false);
324             }
325             eE.set3Body(residues, resI, rotI, resJ, rotJ, resK, rotK, energy);
326             if (rank == 0 && writeEnergyRestart && printFiles) {
327               try {
328                 energyWriter.append(
329                     format("Triple %d %d, %d %d, %d %d: %16.8f", resI, rotI, resJ, rotJ, resK, rotK,
330                         energy));
331                 energyWriter.newLine();
332                 energyWriter.flush();
333               } catch (IOException ex) {
334                 logger.log(Level.SEVERE, " Exception writing energy restart file.", ex);
335               }
336             }
337           }
338         }
339       }
340     }
341 
342     @Override
343     public IntegerSchedule schedule() {
344       // The schedule must be fixed.
345       return IntegerSchedule.fixed();
346     }
347   }
348 }