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 InducedDipoleFieldReduceRegion extends ParallelRegion {
61
62 private static final Logger logger = Logger.getLogger(DirectRegion.class.getName());
63 private final InducedDipoleFieldReduceLoop[] inducedDipoleFieldReduceLoop;
64
65 public double[][][] inducedDipole;
66
67 public double[][][] inducedDipoleCR;
68
69 private Atom[] atoms;
70
71 private double[][] cartesianDipolePhi;
72 private double[][] cartesianDipolePhiCR;
73
74 private AtomicDoubleArray3D field;
75
76 private AtomicDoubleArray3D fieldCR;
77
78 private boolean generalizedKirkwoodTerm;
79
80 private GeneralizedKirkwood generalizedKirkwood;
81 private double aewald;
82 private double aewald3;
83
84 private double dielectric;
85
86 public InducedDipoleFieldReduceRegion(int nt) {
87 inducedDipoleFieldReduceLoop = new InducedDipoleFieldReduceLoop[nt];
88 }
89
90
91
92
93
94
95 public void executeWith(ParallelTeam parallelTeam) {
96 try {
97 parallelTeam.execute(this);
98 } catch (Exception e) {
99 String message = " Exception computing induced dipole field.\n";
100 logger.log(Level.WARNING, message, e);
101 }
102 }
103
104 public void init(
105 Atom[] atoms,
106 double[][][] inducedDipole,
107 double[][][] inducedDipoleCR,
108 boolean generalizedKirkwoodTerm,
109 GeneralizedKirkwood generalizedKirkwood,
110 EwaldParameters ewaldParameters,
111 double dielectric,
112 double[][] cartesianDipolePhi,
113 double[][] cartesianDipolePhiCR,
114 AtomicDoubleArray3D field,
115 AtomicDoubleArray3D fieldCR) {
116
117 this.atoms = atoms;
118 this.inducedDipole = inducedDipole;
119 this.inducedDipoleCR = inducedDipoleCR;
120 this.generalizedKirkwoodTerm = generalizedKirkwoodTerm;
121 this.generalizedKirkwood = generalizedKirkwood;
122 this.aewald = ewaldParameters.aewald;
123 this.aewald3 = ewaldParameters.aewald3;
124 this.dielectric = dielectric;
125 this.cartesianDipolePhi = cartesianDipolePhi;
126 this.cartesianDipolePhiCR = cartesianDipolePhiCR;
127
128 this.field = field;
129 this.fieldCR = fieldCR;
130 }
131
132 @Override
133 public void run() throws Exception {
134 try {
135 int threadID = getThreadIndex();
136 if (inducedDipoleFieldReduceLoop[threadID] == null) {
137 inducedDipoleFieldReduceLoop[threadID] = new InducedDipoleFieldReduceLoop();
138 }
139 int nAtoms = atoms.length;
140 execute(0, nAtoms - 1, inducedDipoleFieldReduceLoop[threadID]);
141 } catch (Exception e) {
142 String message =
143 " Fatal exception computing the mutual induced dipoles in thread "
144 + getThreadIndex()
145 + "\n";
146 logger.log(Level.SEVERE, message, e);
147 }
148 }
149
150 private class InducedDipoleFieldReduceLoop extends IntegerForLoop {
151
152 @Override
153 public void run(int lb, int ub) throws Exception {
154 int threadID = getThreadIndex();
155
156 final double[][] induced0 = inducedDipole[0];
157 final double[][] inducedCR0 = inducedDipoleCR[0];
158
159 if (aewald > 0.0) {
160 for (int i = lb; i <= ub; i++) {
161 double[] dipolei = induced0[i];
162 double[] dipoleCRi = inducedCR0[i];
163 final double[] phii = cartesianDipolePhi[i];
164 final double[] phiCRi = cartesianDipolePhiCR[i];
165 double fx = aewald3 * dipolei[0] - phii[t100];
166 double fy = aewald3 * dipolei[1] - phii[t010];
167 double fz = aewald3 * dipolei[2] - phii[t001];
168 double fxCR = aewald3 * dipoleCRi[0] - phiCRi[t100];
169 double fyCR = aewald3 * dipoleCRi[1] - phiCRi[t010];
170 double fzCR = aewald3 * dipoleCRi[2] - phiCRi[t001];
171 field.add(threadID, i, fx, fy, fz);
172 fieldCR.add(threadID, i, fxCR, fyCR, fzCR);
173 }
174 }
175
176
177 field.reduce(lb, ub);
178 fieldCR.reduce(lb, ub);
179
180
181 if (dielectric > 1.0) {
182 double inverseDielectric = 1.0 / dielectric;
183 for (int i = lb; i <= ub; i++) {
184 field.scale(0, i, inverseDielectric);
185 fieldCR.scale(0, i, inverseDielectric);
186 }
187 }
188
189
190 if (generalizedKirkwoodTerm) {
191 AtomicDoubleArray3D fieldGK = generalizedKirkwood.getFieldGK();
192 AtomicDoubleArray3D fieldGKCR = generalizedKirkwood.getFieldGKCR();
193 for (int i = lb; i <= ub; i++) {
194 field.add(0, i, fieldGK.getX(i), fieldGK.getY(i), fieldGK.getZ(i));
195 fieldCR.add(0, i, fieldGKCR.getX(i), fieldGKCR.getY(i), fieldGKCR.getZ(i));
196 }
197 }
198 }
199
200 @Override
201 public IntegerSchedule schedule() {
202 return IntegerSchedule.fixed();
203 }
204 }
205 }