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.abs;
42
43 import edu.rit.pj.IntegerSchedule;
44 import edu.rit.pj.WorkerIntegerForLoop;
45 import edu.rit.pj.WorkerRegion;
46 import ffx.algorithms.optimize.RotamerOptimization;
47 import ffx.potential.bonded.Residue;
48 import java.util.List;
49 import java.util.Map;
50 import java.util.Set;
51 import java.util.logging.Logger;
52 import java.util.stream.DoubleStream;
53
54
55 public class FourBodyEnergyRegion extends WorkerRegion {
56
57 private static final Logger logger = Logger.getLogger(FourBodyEnergyRegion.class.getName());
58 private final Residue[] residues;
59 private final RotamerOptimization rO;
60 private final DistanceMatrix dM;
61 private final EnergyExpansion eE;
62 private final EliminatedRotamers eR;
63
64
65
66
67 private final List<Residue> allResiduesList;
68
69 private final Map<Integer, Integer[]> fourBodyEnergyMap;
70
71
72
73
74 private final double superpositionThreshold;
75
76 private Set<Integer> keySet;
77
78 public FourBodyEnergyRegion(RotamerOptimization rotamerOptimization, DistanceMatrix dM,
79 EnergyExpansion eE, EliminatedRotamers eR, Residue[] residues, List<Residue> allResiduesList,
80 double superpositionThreshold) {
81 this.rO = rotamerOptimization;
82 this.dM = dM;
83 this.eE = eE;
84 this.eR = eR;
85 this.residues = residues;
86 this.allResiduesList = allResiduesList;
87 this.superpositionThreshold = superpositionThreshold;
88
89 this.fourBodyEnergyMap = eE.getFourBodyEnergyMap();
90 logger.info(format(" Running quads: %d jobs.", fourBodyEnergyMap.size()));
91 }
92
93 @Override
94 public void run() throws Exception {
95 if (!keySet.isEmpty()) {
96 execute(0, keySet.size() - 1, new QuadsEnergyLoop());
97 }
98 }
99
100 @Override
101 public void start() {
102 keySet = fourBodyEnergyMap.keySet();
103 }
104
105 private class QuadsEnergyLoop extends WorkerIntegerForLoop {
106
107 @Override
108 public void run(int lb, int ub) {
109 for (int key = lb; key <= ub; key++) {
110 if (!fourBodyEnergyMap.containsKey(key)) {
111 continue;
112 }
113
114 Integer[] job = fourBodyEnergyMap.get(key);
115 int i = job[0];
116 int ri = job[1];
117 int j = job[2];
118 int rj = job[3];
119 int k = job[4];
120 int rk = job[5];
121 int l = job[6];
122 int rl = job[7];
123
124 if (eR.check(i, ri) || eR.check(j, rj) || eR.check(k, rk) || eR.check(l, rl)
125 || eR.check(i, ri, j, rj) || eR.check(i, ri, k, rk) || eR.check(i, ri, l, rl)
126 || eR.check(j, rj, k, rk) || eR.check(j, rj, l, rl) || eR.check(k, rk, l, rl)) {
127
128 continue;
129 }
130
131 Residue resI = residues[i];
132 Residue resJ = residues[j];
133 Residue resK = residues[k];
134 Residue resL = residues[l];
135
136 int indexI = allResiduesList.indexOf(residues[i]);
137 int indexJ = allResiduesList.indexOf(residues[j]);
138 int indexK = allResiduesList.indexOf(residues[k]);
139 int indexL = allResiduesList.indexOf(residues[l]);
140
141 double rawDist = dM.getRawNBodyDistance(indexI, ri, indexJ, rj, indexK, rk, indexL, rl);
142 double dIJ = dM.checkDistMatrix(indexI, ri, indexJ, rj);
143 double dIK = dM.checkDistMatrix(indexI, ri, indexK, rk);
144 double dIL = dM.checkDistMatrix(indexI, ri, indexL, rl);
145 double dJK = dM.checkDistMatrix(indexJ, rj, indexK, rk);
146 double dJL = dM.checkDistMatrix(indexJ, rj, indexL, rl);
147 double dKL = dM.checkDistMatrix(indexK, rk, indexL, rl);
148
149 double minDist = DoubleStream.of(dIJ, dIK, dIL, dJK, dJL, dKL).min().getAsDouble();
150
151 String distString = " large";
152 if (rawDist < Double.MAX_VALUE) {
153 distString = format("%10.3f", rawDist);
154 }
155
156 double resDist = dM.get4BodyResidueDistance(indexI, ri, indexJ, rj, indexK, rk, indexL, rl);
157 String resDistString = " large";
158 if (resDist < Double.MAX_VALUE) {
159 resDistString = format("%5.3f", resDist);
160 }
161
162 double fourBodyEnergy;
163 if (minDist < superpositionThreshold) {
164 logger.info(format(
165 " Quad %8s %-2d, %8s %-2d, %8s %-2d, %8s %-2d: set to NaN at %13.6f Ang (%s Ang by residue) < %5.3f Ang.",
166 residues[i], ri, residues[j].toFormattedString(false, true), rj,
167 residues[k].toFormattedString(false, true), rk,
168 residues[l].toFormattedString(false, true), rl, minDist, resDistString,
169 superpositionThreshold));
170 } else if (dM.checkQuadDistThreshold(indexI, ri, indexJ, rj, indexK, rk, indexL, rl)) {
171
172 fourBodyEnergy = 0.0;
173 logger.info(format(
174 " Quad %8s %-2d, %8s %-2d, %8s %-2d, %8s %-2d: %s at %s Ang (%s Ang by residue).",
175 resI.toFormattedString(false, true), ri, resJ.toFormattedString(false, true), rj,
176 resK.toFormattedString(false, true), rk, resL.toFormattedString(false, true), rl,
177 rO.formatEnergy(fourBodyEnergy), distString, resDistString));
178 } else {
179 try {
180 fourBodyEnergy = eE.compute4BodyEnergy(residues, i, ri, j, rj, k, rk, l, rl);
181 logger.info(format(
182 " Quad %8s %-2d, %8s %-2d, %8s %-2d, %8s %-2d: %s at %s Ang (%s Ang by residue).",
183 resI.toFormattedString(false, true), ri, resJ.toFormattedString(false, true), rj,
184 resK.toFormattedString(false, true), rk, resL.toFormattedString(false, true), rl,
185 rO.formatEnergy(fourBodyEnergy), distString, resDistString));
186 if (abs(fourBodyEnergy) > 1.0) {
187 StringBuilder sb = new StringBuilder();
188 sb.append(format(
189 " Quad %8s %-2d, %8s %-2d, %8s %-2d, %8s %-2d: %s at %s Ang (%s Ang by residue).\n",
190 resI.toFormattedString(false, true), ri, resJ.toFormattedString(false, true), rj,
191 resK.toFormattedString(false, true), rk, resL.toFormattedString(false, true), rl,
192 rO.formatEnergy(fourBodyEnergy), distString, resDistString));
193 sb.append(format(" Explain: (ref %d) \n", key));
194 sb.append(
195 format(" Self %3d %3d: %.3f\n", i, ri, eE.getSelf(i, ri)));
196 sb.append(
197 format(" Self %3d %3d: %.3f\n", j, rj, eE.getSelf(j, rj)));
198 sb.append(
199 format(" Self %3d %3d: %.3f\n", k, rk, eE.getSelf(k, rk)));
200 sb.append(
201 format(" Self %3d %3d: %.3f\n", l, rl, eE.getSelf(l, rl)));
202 sb.append(format(" Pair %3d %3d %3d %3d: %.3f\n", i, ri, j, rj,
203 eE.get2Body(i, ri, j, rj)));
204 sb.append(format(" Pair %3d %3d %3d %3d: %.3f\n", i, ri, k, rk,
205 eE.get2Body(i, ri, k, rk)));
206 sb.append(format(" Pair %3d %3d %3d %3d: %.3f\n", i, ri, l, rl,
207 eE.get2Body(i, ri, l, rl)));
208 sb.append(format(" Pair %3d %3d %3d %3d: %.3f\n", j, rj, k, rk,
209 eE.get2Body(j, rj, k, rk)));
210 sb.append(format(" Pair %3d %3d %3d %3d: %.3f\n", j, rj, l, rl,
211 eE.get2Body(j, rj, l, rl)));
212 sb.append(format(" Pair %3d %3d %3d %3d: %.3f\n", k, rk, l, rl,
213 eE.get2Body(k, rk, l, rl)));
214 sb.append(format(" Tri %3d %3d %3d %3d %3d %3d: %.3f\n", i, ri, j, rj, k, rk,
215 eE.get3Body(residues, i, ri, j, rj, k, rk)));
216 sb.append(format(" Tri %3d %3d %3d %3d %3d %3d: %.3f\n", i, ri, j, rj, l, rl,
217 eE.get3Body(residues, i, ri, j, rj, l, rl)));
218 sb.append(format(" Tri %3d %3d %3d %3d %3d %3d: %.3f\n", i, ri, k, rk, l, rl,
219 eE.get3Body(residues, i, ri, k, rk, l, rl)));
220 sb.append(format(" Tri %3d %3d %3d %3d %3d %3d: %.3f\n", j, rj, k, rk, l, rl,
221 eE.get3Body(residues, j, rj, k, rk, l, rl)));
222 sb.append(
223 format(" backbone: %.3f\n", rO.getBackboneEnergy()));
224 sb.append(format(" quadEnergy: %.3f\n", fourBodyEnergy));
225 sb.append(" --s--\n");
226 sb.append(" Active residues:\n");
227 for (Residue residue : residues) {
228 if (residue.getSideChainAtoms().get(0).getUse()) {
229 sb.append(format(" %s\n", residue));
230 }
231 }
232 sb.append(" --f--\n");
233 logger.info(sb.toString());
234 }
235 } catch (ArithmeticException ex) {
236 logger.info(format(
237 " Quad %8s %-2d, %8s %-2d, %8s %-2d, %8s %-2d: NaN at %s Ang (%s Ang by residue).",
238 resI.toFormattedString(false, true), ri, resJ.toFormattedString(false, true), rj,
239 resK.toFormattedString(false, true), rk, resL.toFormattedString(false, true), rl,
240 distString, resDistString));
241 }
242 }
243 }
244 }
245
246 @Override
247 public IntegerSchedule schedule() {
248 return IntegerSchedule.fixed();
249 }
250 }
251 }