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.nonbonded.GeneralizedKirkwood;
47
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50
51 import static ffx.potential.parameters.MultipoleType.t001;
52 import static ffx.potential.parameters.MultipoleType.t010;
53 import static ffx.potential.parameters.MultipoleType.t100;
54
55
56
57
58
59
60
61 public class DirectRegion extends ParallelRegion {
62
63 private static final Logger logger = Logger.getLogger(DirectRegion.class.getName());
64 private final DirectLoop[] directLoop;
65
66 public double[][][] inducedDipole;
67
68 public double[][][] inducedDipoleCR;
69
70 public double[][] directDipole;
71
72 public double[][] directDipoleCR;
73 public double[][] directField;
74 public double[][] directFieldCR;
75
76 private Atom[] atoms;
77
78 private double[] polarizability;
79
80 private double[][][] globalMultipole;
81
82 private double[][] cartMultipolePhi;
83
84 private AtomicDoubleArray3D field;
85
86 private AtomicDoubleArray3D fieldCR;
87
88 private boolean generalizedKirkwoodTerm;
89
90 private GeneralizedKirkwood generalizedKirkwood;
91 private double aewald;
92 private double aewald3;
93
94 private double soluteDielectric;
95
96 public DirectRegion(int nt) {
97 directLoop = new DirectLoop[nt];
98 }
99
100
101
102
103
104
105 public void executeWith(ParallelTeam parallelTeam) {
106 try {
107 parallelTeam.execute(this);
108 } catch (Exception e) {
109 String message = " Exception computing direct induced dipoles.\n";
110 logger.log(Level.WARNING, message, e);
111 }
112 }
113
114 public void init(
115 Atom[] atoms,
116 double[] polarizability,
117 double[][][] globalMultipole,
118 double[][] cartMultipolePhi,
119 AtomicDoubleArray3D field,
120 AtomicDoubleArray3D fieldCR,
121 boolean generalizedKirkwoodTerm,
122 GeneralizedKirkwood generalizedKirkwood,
123 EwaldParameters ewaldParameters,
124 double soluteDielectric,
125 double[][][] inducedDipole,
126 double[][][] inducedDipoleCR,
127 double[][] directDipole,
128 double[][] directDipoleCR,
129 double[][] directField,
130 double[][] directFieldCR) {
131
132 this.atoms = atoms;
133 this.polarizability = polarizability;
134 this.globalMultipole = globalMultipole;
135 this.cartMultipolePhi = cartMultipolePhi;
136 this.field = field;
137 this.fieldCR = fieldCR;
138 this.generalizedKirkwoodTerm = generalizedKirkwoodTerm;
139 this.generalizedKirkwood = generalizedKirkwood;
140 this.aewald = ewaldParameters.aewald;
141 this.aewald3 = ewaldParameters.aewald3;
142 this.soluteDielectric = soluteDielectric;
143
144 this.inducedDipole = inducedDipole;
145 this.inducedDipoleCR = inducedDipoleCR;
146 this.directDipole = directDipole;
147 this.directDipoleCR = directDipoleCR;
148 this.directField = directField;
149 this.directFieldCR = directFieldCR;
150 }
151
152 @Override
153 public void run() throws Exception {
154 int ti = getThreadIndex();
155 if (directLoop[ti] == null) {
156 directLoop[ti] = new DirectLoop();
157 }
158 try {
159 int nAtoms = atoms.length;
160 execute(0, nAtoms - 1, directLoop[ti]);
161 } catch (Exception e) {
162 String message =
163 "Fatal exception computing the direct induced dipoles in thread "
164 + getThreadIndex()
165 + "\n";
166 logger.log(Level.SEVERE, message, e);
167 }
168 }
169
170 private class DirectLoop extends IntegerForLoop {
171
172 @Override
173 public void run(int lb, int ub) throws Exception {
174 int threadID = getThreadIndex();
175
176 if (aewald > 0.0) {
177
178 for (int i = lb; i <= ub; i++) {
179 double[] mpolei = globalMultipole[0][i];
180 double[] phii = cartMultipolePhi[i];
181 double fx = aewald3 * mpolei[t100] - phii[t100];
182 double fy = aewald3 * mpolei[t010] - phii[t010];
183 double fz = aewald3 * mpolei[t001] - phii[t001];
184 field.add(threadID, i, fx, fy, fz);
185 fieldCR.add(threadID, i, fx, fy, fz);
186 }
187 }
188
189
190 field.reduce(lb, ub);
191 fieldCR.reduce(lb, ub);
192
193
194 if (soluteDielectric > 1.0) {
195 double inverseDielectric = 1.0 / soluteDielectric;
196 for (int i = lb; i <= ub; i++) {
197 field.scale(0, i, inverseDielectric);
198 fieldCR.scale(0, i, inverseDielectric);
199 }
200 }
201
202 if (generalizedKirkwoodTerm) {
203
204 AtomicDoubleArray3D fieldGK = generalizedKirkwood.getFieldGK();
205 for (int i = lb; i <= ub; i++) {
206 double fx = fieldGK.getX(i);
207 double fy = fieldGK.getY(i);
208 double fz = fieldGK.getZ(i);
209 field.add(0, i, fx, fy, fz);
210 fieldCR.add(0, i, fx, fy, fz);
211 }
212 }
213
214
215
216 final double[][] induced0 = inducedDipole[0];
217 final double[][] inducedCR0 = inducedDipoleCR[0];
218 for (int i = lb; i <= ub; i++) {
219 final double polar = polarizability[i];
220 final double[] ind = induced0[i];
221 final double[] directi = directDipole[i];
222 directField[i][0] = field.getX(i);
223 directField[i][1] = field.getY(i);
224 directField[i][2] = field.getZ(i);
225 ind[0] = polar * directField[i][0];
226 ind[1] = polar * directField[i][1];
227 ind[2] = polar * directField[i][2];
228 directi[0] = ind[0];
229 directi[1] = ind[1];
230 directi[2] = ind[2];
231 final double[] indCR = inducedCR0[i];
232 final double[] directCRi = directDipoleCR[i];
233 directFieldCR[i][0] = fieldCR.getX(i);
234 directFieldCR[i][1] = fieldCR.getY(i);
235 directFieldCR[i][2] = fieldCR.getZ(i);
236 indCR[0] = polar * directFieldCR[i][0];
237 indCR[1] = polar * directFieldCR[i][1];
238 indCR[2] = polar * directFieldCR[i][2];
239 directCRi[0] = indCR[0];
240 directCRi[1] = indCR[1];
241 directCRi[2] = indCR[2];
242 }
243 }
244
245 @Override
246 public IntegerSchedule schedule() {
247 return IntegerSchedule.fixed();
248 }
249 }
250 }