1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
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
70
71
72 private final List<Residue> allResiduesList;
73
74 private final Map<Integer, Integer[]> threeBodyEnergyMap;
75
76 private final BufferedWriter energyWriter;
77
78 private final Comm world;
79
80 private final int numProc;
81
82
83
84
85 private final double superpositionThreshold;
86
87 private final boolean master;
88
89 private final int rank;
90
91 private final boolean verbose;
92
93 private final boolean writeEnergyRestart;
94
95
96
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
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
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
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
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
269
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
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
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
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
344 return IntegerSchedule.fixed();
345 }
346 }
347 }