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 static ffx.potential.parameters.MultipoleType.t001;
41 import static ffx.potential.parameters.MultipoleType.t010;
42 import static ffx.potential.parameters.MultipoleType.t100;
43
44 import edu.rit.pj.IntegerForLoop;
45 import edu.rit.pj.IntegerSchedule;
46 import edu.rit.pj.ParallelRegion;
47 import ffx.numerics.atomic.AtomicDoubleArray3D;
48 import ffx.potential.bonded.Atom;
49 import ffx.potential.nonbonded.GeneralizedKirkwood;
50 import java.util.logging.Level;
51 import java.util.logging.Logger;
52
53
54
55
56
57
58
59 public class OPTRegion extends ParallelRegion {
60
61 private static final Logger logger = Logger.getLogger(OPTRegion.class.getName());
62 public final double[] optCoefficients;
63
64 public final int optOrder = 2;
65
66 private final OPTLoop[] optLoop;
67 private final double[] optCoefficientsSum;
68
69 public double[][][] inducedDipole;
70
71 public double[][][] inducedDipoleCR;
72 public double[][][] optDipole;
73 public double[][][] optDipoleCR;
74
75 private Atom[] atoms;
76
77 private double[] polarizability;
78 private double[][] cartesianDipolePhi;
79 private double[][] cartesianDipolePhiCR;
80
81 private AtomicDoubleArray3D field;
82
83 private AtomicDoubleArray3D fieldCR;
84
85 private boolean generalizedKirkwoodTerm;
86
87 private GeneralizedKirkwood generalizedKirkwood;
88 private double aewald;
89 private double aewald3;
90 private double dielectric;
91 private int currentOptOrder;
92
93 public OPTRegion(int nt) {
94 optLoop = new OPTLoop[nt];
95 optCoefficients = new double[optOrder + 1];
96 optCoefficientsSum = new double[optOrder + 1];
97 switch (optOrder) {
98 case 1:
99 optCoefficients[0] = 0.530;
100 optCoefficients[1] = 0.604;
101 break;
102 case 2:
103 optCoefficients[0] = 0.042;
104 optCoefficients[1] = 0.635;
105 optCoefficients[2] = 0.414;
106 break;
107 case 3:
108 optCoefficients[0] = -0.132;
109 optCoefficients[1] = 0.218;
110 optCoefficients[2] = 0.637;
111 optCoefficients[3] = 0.293;
112 break;
113 case 4:
114 optCoefficients[0] = -0.071;
115 optCoefficients[1] = -0.096;
116 optCoefficients[2] = 0.358;
117 optCoefficients[3] = 0.587;
118 optCoefficients[4] = 0.216;
119 break;
120 case 5:
121 optCoefficients[0] = -0.005;
122 optCoefficients[1] = -0.129;
123 optCoefficients[2] = -0.026;
124 optCoefficients[3] = 0.465;
125 optCoefficients[4] = 0.528;
126 optCoefficients[5] = 0.161;
127 break;
128 case 6:
129 optCoefficients[0] = 0.014;
130 optCoefficients[1] = -0.041;
131 optCoefficients[2] = -0.172;
132 optCoefficients[3] = 0.073;
133 optCoefficients[4] = 0.535;
134 optCoefficients[5] = 0.467;
135 optCoefficients[6] = 0.122;
136 break;
137 default:
138 logger.severe(" Unsupported OPT order.");
139 }
140
141 for (int i = 0; i <= optOrder; i++) {
142 for (int j = optOrder; j >= i; j--) {
143 optCoefficientsSum[i] += optCoefficients[j];
144 }
145 }
146 }
147
148 public void init(
149 int currentOptOrder,
150 Atom[] atoms,
151 double[] polarizability,
152 double[][][] inducedDipole,
153 double[][][] inducedDipoleCR,
154 double[][] cartesianDipolePhi,
155 double[][] cartesianDipolePhiCR,
156 AtomicDoubleArray3D field,
157 AtomicDoubleArray3D fieldCR,
158 boolean generalizedKirkwoodTerm,
159 GeneralizedKirkwood generalizedKirkwood,
160 EwaldParameters ewaldParameters,
161 double dielectric) {
162 this.currentOptOrder = currentOptOrder;
163 this.atoms = atoms;
164 this.polarizability = polarizability;
165 this.inducedDipole = inducedDipole;
166 this.inducedDipoleCR = inducedDipoleCR;
167 this.cartesianDipolePhi = cartesianDipolePhi;
168 this.cartesianDipolePhiCR = cartesianDipolePhiCR;
169 this.field = field;
170 this.fieldCR = fieldCR;
171 this.generalizedKirkwoodTerm = generalizedKirkwoodTerm;
172 this.generalizedKirkwood = generalizedKirkwood;
173 this.aewald = ewaldParameters.aewald;
174 this.aewald3 = ewaldParameters.aewald3;
175 this.dielectric = dielectric;
176 }
177
178 @Override
179 public void run() throws Exception {
180 try {
181 int ti = getThreadIndex();
182 if (optLoop[ti] == null) {
183 optLoop[ti] = new OPTRegion.OPTLoop();
184 }
185 int nAtoms = atoms.length;
186 execute(0, nAtoms - 1, optLoop[ti]);
187 } catch (RuntimeException ex) {
188 logger.warning(
189 "Fatal exception computing the opt induced dipoles in thread " + getThreadIndex());
190 throw ex;
191 } catch (Exception e) {
192 String message =
193 "Fatal exception computing the opt induced dipoles in thread " + getThreadIndex() + "\n";
194 logger.log(Level.SEVERE, message, e);
195 }
196 }
197
198 private class OPTLoop extends IntegerForLoop {
199
200 @Override
201 public void run(int lb, int ub) throws Exception {
202 int threadID = getThreadIndex();
203
204 final double[][] induced0 = inducedDipole[0];
205 final double[][] inducedCR0 = inducedDipoleCR[0];
206 if (aewald > 0.0) {
207
208 for (int i = lb; i <= ub; i++) {
209 double[] dipolei = induced0[i];
210 double[] dipoleCRi = inducedCR0[i];
211 final double[] phii = cartesianDipolePhi[i];
212 final double[] phiCRi = cartesianDipolePhiCR[i];
213 double fx = aewald3 * dipolei[0] - phii[t100];
214 double fy = aewald3 * dipolei[1] - phii[t010];
215 double fz = aewald3 * dipolei[2] - phii[t001];
216 double fxCR = aewald3 * dipoleCRi[0] - phiCRi[t100];
217 double fyCR = aewald3 * dipoleCRi[1] - phiCRi[t010];
218 double fzCR = aewald3 * dipoleCRi[2] - phiCRi[t001];
219 field.add(threadID, i, fx, fy, fz);
220 fieldCR.add(threadID, i, fxCR, fyCR, fzCR);
221 }
222 }
223
224
225 field.reduce(lb, ub);
226 fieldCR.reduce(lb, ub);
227
228
229 if (dielectric > 1.0) {
230 double inverseDielectric = 1.0 / dielectric;
231 for (int i = lb; i <= ub; i++) {
232 field.scale(0, i, inverseDielectric);
233 fieldCR.scale(0, i, inverseDielectric);
234 }
235 }
236
237 if (generalizedKirkwoodTerm) {
238 AtomicDoubleArray3D fieldGK = generalizedKirkwood.getFieldGK();
239 AtomicDoubleArray3D fieldGKCR = generalizedKirkwood.getFieldGKCR();
240
241 for (int i = lb; i <= ub; i++) {
242 field.add(0, i, fieldGK.getX(i), fieldGK.getY(i), fieldGK.getZ(i));
243 fieldCR.add(0, i, fieldGKCR.getX(i), fieldGKCR.getY(i), fieldGKCR.getZ(i));
244 }
245 }
246
247
248 for (int i = lb; i <= ub; i++) {
249 final double[] ind = induced0[i];
250 final double[] indCR = inducedCR0[i];
251 final double polar = polarizability[i];
252 for (int j = 0; j < 3; j++) {
253 optDipole[currentOptOrder][i][j] = polar * field.get(j, i);
254 optDipoleCR[currentOptOrder][i][j] = polar * fieldCR.get(j, i);
255 ind[j] = optDipole[currentOptOrder][i][j];
256 indCR[j] = optDipoleCR[currentOptOrder][i][j];
257 }
258 }
259 }
260
261 @Override
262 public IntegerSchedule schedule() {
263 return IntegerSchedule.fixed();
264 }
265 }
266 }