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