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