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 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
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
71
72
73 private final List<Residue> allResiduesList;
74
75 private final Map<Integer, Integer[]> threeBodyEnergyMap;
76
77 private final BufferedWriter energyWriter;
78
79 private final Comm world;
80
81 private final int numProc;
82
83
84
85
86 private final double superpositionThreshold;
87
88 private final boolean master;
89
90 private final int rank;
91
92 private final boolean verbose;
93
94 private final boolean writeEnergyRestart;
95
96
97
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
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
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
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
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
270
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
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
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
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
345 return IntegerSchedule.fixed();
346 }
347 }
348 }