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.potential.nonbonded.pme;
39
40 import edu.rit.pj.IntegerForLoop;
41 import edu.rit.pj.IntegerSchedule;
42 import edu.rit.pj.ParallelRegion;
43 import edu.rit.pj.ParallelTeam;
44 import ffx.numerics.atomic.AtomicDoubleArray3D;
45 import ffx.potential.bonded.Atom;
46 import ffx.potential.parameters.ForceField;
47 import ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition;
48 import ffx.potential.utils.EnergyException;
49
50 import java.util.Arrays;
51 import java.util.logging.Level;
52 import java.util.logging.Logger;
53
54 import static java.lang.Double.isInfinite;
55 import static java.lang.Double.isNaN;
56 import static java.lang.String.format;
57
58
59
60
61
62
63
64 public class ReduceRegion extends ParallelRegion {
65
66 private static final Logger logger = Logger.getLogger(ReduceRegion.class.getName());
67
68
69
70
71
72
73 private final boolean rotateMultipoles;
74
75 private final TorqueLoop[] torqueLoop;
76 private final ReduceLoop[] reduceLoop;
77
78
79
80
81 private boolean lambdaTerm;
82
83
84
85 private boolean gradient;
86
87
88
89 private Atom[] atoms;
90
91
92
93 private double[][][] coordinates;
94
95
96
97 private MultipoleFrameDefinition[] frame;
98
99
100
101 private int[][] axisAtom;
102
103
104
105 private AtomicDoubleArray3D grad;
106
107
108
109 private AtomicDoubleArray3D torque;
110
111
112
113 private AtomicDoubleArray3D lambdaGrad;
114
115
116
117 private AtomicDoubleArray3D lambdaTorque;
118
119 public ReduceRegion(int threadCount, ForceField forceField) {
120 torqueLoop = new TorqueLoop[threadCount];
121 reduceLoop = new ReduceLoop[threadCount];
122 rotateMultipoles = forceField.getBoolean("ROTATE_MULTIPOLES", true);
123 }
124
125
126
127
128
129
130 public void executeWith(ParallelTeam parallelTeam) {
131 try {
132 parallelTeam.execute(this);
133 } catch (Exception e) {
134 String message = "Exception calculating torques.";
135 logger.log(Level.SEVERE, message, e);
136 }
137 }
138
139 public void init(
140 boolean lambdaTerm,
141 boolean gradient,
142 Atom[] atoms,
143 double[][][] coordinates,
144 MultipoleFrameDefinition[] frame,
145 int[][] axisAtom,
146 AtomicDoubleArray3D grad,
147 AtomicDoubleArray3D torque,
148 AtomicDoubleArray3D lambdaGrad,
149 AtomicDoubleArray3D lambdaTorque) {
150 this.lambdaTerm = lambdaTerm;
151 this.gradient = gradient;
152 this.atoms = atoms;
153 this.coordinates = coordinates;
154 this.frame = frame;
155 this.axisAtom = axisAtom;
156 this.grad = grad;
157 this.torque = torque;
158 this.lambdaGrad = lambdaGrad;
159 this.lambdaTorque = lambdaTorque;
160 }
161
162 @Override
163 public void run() throws EnergyException {
164 int nAtoms = atoms.length;
165 try {
166 int threadIndex = getThreadIndex();
167 if (torqueLoop[threadIndex] == null) {
168 torqueLoop[threadIndex] = new TorqueLoop();
169 reduceLoop[threadIndex] = new ReduceLoop();
170 }
171 if (rotateMultipoles) {
172 execute(0, nAtoms - 1, torqueLoop[threadIndex]);
173 }
174 execute(0, nAtoms - 1, reduceLoop[threadIndex]);
175 } catch (Exception e) {
176 if (e instanceof EnergyException) {
177 throw (EnergyException) e;
178 }
179 String message = "Fatal exception computing torque in thread " + getThreadIndex() + "\n";
180 logger.log(Level.SEVERE, message, e);
181 }
182 }
183
184 private class TorqueLoop extends IntegerForLoop {
185
186 private Torque torques;
187 private int threadID;
188 double[] trq = new double[3];
189 double[][] g = new double[4][3];
190 int[] frameIndex = new int[4];
191
192 TorqueLoop() {
193 torques = new Torque();
194 }
195
196 @Override
197 public void run(int lb, int ub) throws EnergyException {
198 if (gradient) {
199 torque.reduce(lb, ub);
200 for (int i = lb; i <= ub; i++) {
201
202 Arrays.fill(frameIndex, -1);
203 trq[0] = torque.getX(i);
204 trq[1] = torque.getY(i);
205 trq[2] = torque.getZ(i);
206
207
208 if (isNaN(trq[0]) || isInfinite(trq[0])
209 || isNaN(trq[1]) || isInfinite(trq[1])
210 || isNaN(trq[2]) || isInfinite(trq[2])) {
211 Atom a = atoms[i];
212 throw new EnergyException(
213 format(" Undefined torque (%8.3f,%8.3f,%8.3f) for atom %s.", trq[0], trq[1], trq[2], a));
214 }
215
216 torques.torque(i, 0, trq, frameIndex, g);
217 for (int j = 0; j < 4; j++) {
218 int index = frameIndex[j];
219 if (index >= 0) {
220 double[] gj = g[j];
221
222
223 if (isNaN(gj[0]) || isInfinite(gj[0])
224 || isNaN(gj[1]) || isInfinite(gj[1])
225 || isNaN(gj[2]) || isInfinite(gj[2])) {
226 Atom ai = atoms[i];
227 Atom aj = atoms[index];
228 throw new EnergyException(
229 format(" Undefined gradient (%8.3f,%8.3f,%8.3f)\n For atom: %s\n From torque of atom %s",
230 gj[0], gj[1], gj[2], aj, ai));
231 }
232
233
234 grad.add(threadID, index, gj[0], gj[1], gj[2]);
235 }
236 }
237 }
238 }
239 if (lambdaTerm) {
240 lambdaTorque.reduce(lb, ub);
241 for (int i = lb; i <= ub; i++) {
242
243 Arrays.fill(frameIndex, -1);
244 trq[0] = lambdaTorque.getX(i);
245 trq[1] = lambdaTorque.getY(i);
246 trq[2] = lambdaTorque.getZ(i);
247 torques.torque(i, 0, trq, frameIndex, g);
248 for (int j = 0; j < 4; j++) {
249 int index = frameIndex[j];
250 if (index >= 0) {
251 double[] gj = g[j];
252 lambdaGrad.add(threadID, index, gj[0], gj[1], gj[2]);
253 }
254 }
255 }
256 }
257 }
258
259 @Override
260 public IntegerSchedule schedule() {
261 return IntegerSchedule.fixed();
262 }
263
264 @Override
265 public void start() {
266 threadID = getThreadIndex();
267 torques.init(axisAtom, frame, coordinates);
268 }
269 }
270
271 private class ReduceLoop extends IntegerForLoop {
272
273 @Override
274 public void run(int lb, int ub) throws EnergyException {
275 if (gradient) {
276 grad.reduce(lb, ub);
277 for (int i = lb; i <= ub; i++) {
278 Atom ai = atoms[i];
279 ai.addToXYZGradient(grad.getX(i), grad.getY(i), grad.getZ(i));
280 }
281 }
282 if (lambdaTerm) {
283 lambdaGrad.reduce(lb, ub);
284 }
285 }
286
287 @Override
288 public IntegerSchedule schedule() {
289 return IntegerSchedule.fixed();
290 }
291 }
292 }