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 java.util.Arrays;
49 import java.util.logging.Level;
50 import java.util.logging.Logger;
51
52
53
54
55
56
57
58 public class ReduceRegion extends ParallelRegion {
59
60 private static final Logger logger = Logger.getLogger(ReduceRegion.class.getName());
61
62
63
64
65
66
67 private final boolean rotateMultipoles;
68
69 private final TorqueLoop[] torqueLoop;
70 private final ReduceLoop[] reduceLoop;
71
72
73
74
75 private boolean lambdaTerm;
76
77 private boolean gradient;
78
79 private Atom[] atoms;
80
81 private double[][][] coordinates;
82
83 private MultipoleFrameDefinition[] frame;
84
85 private int[][] axisAtom;
86
87 private AtomicDoubleArray3D grad;
88
89 private AtomicDoubleArray3D torque;
90
91 private AtomicDoubleArray3D lambdaGrad;
92
93 private AtomicDoubleArray3D lambdaTorque;
94
95 public ReduceRegion(int threadCount, ForceField forceField) {
96 torqueLoop = new TorqueLoop[threadCount];
97 reduceLoop = new ReduceLoop[threadCount];
98 rotateMultipoles = forceField.getBoolean("ROTATE_MULTIPOLES", true);
99 }
100
101
102
103
104
105
106 public void excuteWith(ParallelTeam parallelTeam) {
107 try {
108 parallelTeam.execute(this);
109 } catch (Exception e) {
110 String message = "Exception calculating torques.";
111 logger.log(Level.SEVERE, message, e);
112 }
113 }
114
115 public void init(
116 boolean lambdaTerm,
117 boolean gradient,
118 Atom[] atoms,
119 double[][][] coordinates,
120 MultipoleFrameDefinition[] frame,
121 int[][] axisAtom,
122 AtomicDoubleArray3D grad,
123 AtomicDoubleArray3D torque,
124 AtomicDoubleArray3D lambdaGrad,
125 AtomicDoubleArray3D lambdaTorque) {
126 this.lambdaTerm = lambdaTerm;
127 this.gradient = gradient;
128 this.atoms = atoms;
129 this.coordinates = coordinates;
130 this.frame = frame;
131 this.axisAtom = axisAtom;
132 this.grad = grad;
133 this.torque = torque;
134 this.lambdaGrad = lambdaGrad;
135 this.lambdaTorque = lambdaTorque;
136 }
137
138 @Override
139 public void run() {
140 int nAtoms = atoms.length;
141 try {
142 int threadIndex = getThreadIndex();
143 if (torqueLoop[threadIndex] == null) {
144 torqueLoop[threadIndex] = new TorqueLoop();
145 reduceLoop[threadIndex] = new ReduceLoop();
146 }
147 if (rotateMultipoles) {
148 execute(0, nAtoms - 1, torqueLoop[threadIndex]);
149 }
150 execute(0, nAtoms - 1, reduceLoop[threadIndex]);
151 } catch (Exception e) {
152 String message = "Fatal exception computing torque in thread " + getThreadIndex() + "\n";
153 logger.log(Level.SEVERE, message, e);
154 }
155 }
156
157 private class TorqueLoop extends IntegerForLoop {
158
159 private Torque torques;
160 private int threadID;
161 double[] trq = new double[3];
162 double[][] g = new double[4][3];
163 int[] frameIndex = new int[4];
164
165 TorqueLoop() {
166 torques = new Torque();
167 }
168
169 @Override
170 public void run(int lb, int ub) {
171 if (gradient) {
172 torque.reduce(lb, ub);
173 for (int i = lb; i <= ub; i++) {
174
175 Arrays.fill(frameIndex, -1);
176 trq[0] = torque.getX(i);
177 trq[1] = torque.getY(i);
178 trq[2] = torque.getZ(i);
179 torques.torque(i, 0, trq, frameIndex, g);
180 for (int j = 0; j < 4; j++) {
181 int index = frameIndex[j];
182 if (index >= 0) {
183 double[] gj = g[j];
184 grad.add(threadID, index, gj[0], gj[1], gj[2]);
185 }
186 }
187 }
188 }
189 if (lambdaTerm) {
190 lambdaTorque.reduce(lb, ub);
191 for (int i = lb; i <= ub; i++) {
192
193 Arrays.fill(frameIndex, -1);
194 trq[0] = lambdaTorque.getX(i);
195 trq[1] = lambdaTorque.getY(i);
196 trq[2] = lambdaTorque.getZ(i);
197 torques.torque(i, 0, trq, frameIndex, g);
198 for (int j = 0; j < 4; j++) {
199 int index = frameIndex[j];
200 if (index >= 0) {
201 double[] gj = g[j];
202 lambdaGrad.add(threadID, index, gj[0], gj[1], gj[2]);
203 }
204 }
205 }
206 }
207 }
208
209 @Override
210 public IntegerSchedule schedule() {
211 return IntegerSchedule.fixed();
212 }
213
214 @Override
215 public void start() {
216 threadID = getThreadIndex();
217 torques.init(axisAtom, frame, coordinates);
218 }
219 }
220
221 private class ReduceLoop extends IntegerForLoop {
222
223 @Override
224 public void run(int lb, int ub) throws Exception {
225 if (gradient) {
226 grad.reduce(lb, ub);
227 for (int i = lb; i <= ub; i++) {
228 Atom ai = atoms[i];
229 ai.addToXYZGradient(grad.getX(i), grad.getY(i), grad.getZ(i));
230 }
231 }
232 if (lambdaTerm) {
233 lambdaGrad.reduce(lb, ub);
234 }
235 }
236
237 @Override
238 public IntegerSchedule schedule() {
239 return IntegerSchedule.fixed();
240 }
241 }
242 }