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