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