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