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.reduction.SharedDouble;
44 import ffx.numerics.atomic.AtomicDoubleArray3D;
45 import ffx.potential.bonded.Atom;
46 import ffx.potential.nonbonded.GeneralizedKirkwood;
47 import ffx.potential.parameters.ForceField;
48 import ffx.utilities.FFXProperty;
49
50 import java.util.logging.Level;
51 import java.util.logging.Logger;
52
53 import static ffx.potential.parameters.MultipoleType.t001;
54 import static ffx.potential.parameters.MultipoleType.t010;
55 import static ffx.potential.parameters.MultipoleType.t100;
56 import static ffx.utilities.PropertyGroup.ElectrostaticsFunctionalForm;
57 import static org.apache.commons.math3.util.FastMath.max;
58
59
60
61
62
63
64
65 public class SORRegion extends ParallelRegion {
66
67 private static final Logger logger = Logger.getLogger(SORRegion.class.getName());
68 private final int maxThreads;
69
70 private static final double DEFAULT_POLAR_SOR = 0.7;
71 @FFXProperty(name = "polar-sor", propertyGroup = ElectrostaticsFunctionalForm, defaultValue = "0.7",
72 description = "The induced dipole successive over-relaxation convergence acceleration factor.")
73 private double polsor;
74 private final SORLoop[] sorLoop;
75 private final SharedDouble sharedEps;
76 private final SharedDouble sharedEpsCR;
77
78 public double[][][] inducedDipole;
79
80 public double[][][] inducedDipoleCR;
81
82 public double[][] directDipole;
83
84 public double[][] directDipoleCR;
85
86 private Atom[] atoms;
87
88 private double[] polarizability;
89 private double[][] cartesianDipolePhi;
90 private double[][] cartesianDipolePhiCR;
91
92 private AtomicDoubleArray3D field;
93
94 private AtomicDoubleArray3D fieldCR;
95
96 private boolean generalizedKirkwoodTerm;
97
98 private GeneralizedKirkwood generalizedKirkwood;
99 private double aewald;
100 private double aewald3;
101
102 public SORRegion(int nt, ForceField forceField) {
103 maxThreads = nt;
104 sorLoop = new SORLoop[nt];
105 sharedEps = new SharedDouble();
106 sharedEpsCR = new SharedDouble();
107 polsor = forceField.getDouble("POLAR_SOR", DEFAULT_POLAR_SOR);
108 }
109
110
111
112
113
114 public double getPolarSOR() {
115 return polsor;
116 }
117
118
119
120
121
122 public void setPolarSOR(double polarSOR) {
123 this.polsor = polarSOR;
124 }
125
126 public double getEps() {
127 double eps = sharedEps.get();
128 double epsCR = sharedEpsCR.get();
129 return max(eps, epsCR);
130 }
131
132 public double getSOR() {
133 return polsor;
134 }
135
136 public void init(
137 Atom[] atoms,
138 double[] polarizability,
139 double[][][] inducedDipole,
140 double[][][] inducedDipoleCR,
141 double[][] directDipole,
142 double[][] directDipoleCR,
143 double[][] cartesianDipolePhi,
144 double[][] cartesianDipolePhiCR,
145 AtomicDoubleArray3D field,
146 AtomicDoubleArray3D fieldCR,
147 boolean generalizedKirkwoodTerm,
148 GeneralizedKirkwood generalizedKirkwood,
149 EwaldParameters ewaldParameters) {
150 this.atoms = atoms;
151 this.polarizability = polarizability;
152 this.inducedDipole = inducedDipole;
153 this.inducedDipoleCR = inducedDipoleCR;
154 this.directDipole = directDipole;
155 this.directDipoleCR = directDipoleCR;
156 this.cartesianDipolePhi = cartesianDipolePhi;
157 this.cartesianDipolePhiCR = cartesianDipolePhiCR;
158 this.field = field;
159 this.fieldCR = fieldCR;
160 this.generalizedKirkwoodTerm = generalizedKirkwoodTerm;
161 this.generalizedKirkwood = generalizedKirkwood;
162 this.aewald = ewaldParameters.aewald;
163 this.aewald3 = ewaldParameters.aewald3;
164 }
165
166 @Override
167 public void run() throws Exception {
168 try {
169 int ti = getThreadIndex();
170 if (sorLoop[ti] == null) {
171 sorLoop[ti] = new SORLoop();
172 }
173 int nAtoms = atoms.length;
174 execute(0, nAtoms - 1, sorLoop[ti]);
175 } catch (RuntimeException ex) {
176 logger.warning("Fatal exception computing the mutual induced dipoles in thread " + getThreadIndex());
177 throw ex;
178 } catch (Exception e) {
179 String message = "Fatal exception computing the mutual induced dipoles in thread "
180 + getThreadIndex() + "\n";
181 logger.log(Level.SEVERE, message, e);
182 }
183 }
184
185 @Override
186 public void start() {
187 sharedEps.set(0.0);
188 sharedEpsCR.set(0.0);
189 }
190
191 private class SORLoop extends IntegerForLoop {
192
193 private double eps, epsCR;
194
195 @Override
196 public void finish() {
197 sharedEps.addAndGet(eps);
198 sharedEpsCR.addAndGet(epsCR);
199 }
200
201 @Override
202 public void run(int lb, int ub) throws Exception {
203 int threadID = getThreadIndex();
204 final double[][] induced0 = inducedDipole[0];
205 final double[][] inducedCR0 = inducedDipoleCR[0];
206 if (aewald > 0.0) {
207
208 for (int i = lb; i <= ub; i++) {
209 double[] dipolei = induced0[i];
210 double[] dipoleCRi = inducedCR0[i];
211 final double[] phii = cartesianDipolePhi[i];
212 final double[] phiCRi = cartesianDipolePhiCR[i];
213 double fx = aewald3 * dipolei[0] - phii[t100];
214 double fy = aewald3 * dipolei[1] - phii[t010];
215 double fz = aewald3 * dipolei[2] - phii[t001];
216 double fxCR = aewald3 * dipoleCRi[0] - phiCRi[t100];
217 double fyCR = aewald3 * dipoleCRi[1] - phiCRi[t010];
218 double fzCR = aewald3 * dipoleCRi[2] - phiCRi[t001];
219 field.add(threadID, i, fx, fy, fz);
220 fieldCR.add(threadID, i, fxCR, fyCR, fzCR);
221 }
222 }
223
224 if (generalizedKirkwoodTerm) {
225 AtomicDoubleArray3D fieldGK = generalizedKirkwood.getFieldGK();
226 AtomicDoubleArray3D fieldGKCR = generalizedKirkwood.getFieldGKCR();
227
228 for (int i = lb; i <= ub; i++) {
229 field.add(threadID, i, fieldGK.getX(i), fieldGK.getY(i), fieldGK.getZ(i));
230 fieldCR.add(threadID, i, fieldGKCR.getX(i), fieldGKCR.getY(i), fieldGKCR.getZ(i));
231 }
232 }
233
234
235 field.reduce(lb, ub);
236 fieldCR.reduce(lb, ub);
237
238
239 for (int i = lb; i <= ub; i++) {
240 final double[] ind = induced0[i];
241 final double[] indCR = inducedCR0[i];
242 final double[] direct = directDipole[i];
243 final double[] directCR = directDipoleCR[i];
244 final double polar = polarizability[i];
245 for (int j = 0; j < 3; j++) {
246 double previous = ind[j];
247 double mutual = polar * field.get(j, i);
248 ind[j] = direct[j] + mutual;
249 double delta = polsor * (ind[j] - previous);
250 ind[j] = previous + delta;
251 eps += delta * delta;
252 previous = indCR[j];
253 mutual = polar * fieldCR.get(j, i);
254 indCR[j] = directCR[j] + mutual;
255 delta = polsor * (indCR[j] - previous);
256 indCR[j] = previous + delta;
257 epsCR += delta * delta;
258 }
259 }
260 }
261
262 @Override
263 public IntegerSchedule schedule() {
264 return IntegerSchedule.fixed();
265 }
266
267 @Override
268 public void start() {
269 eps = 0.0;
270 epsCR = 0.0;
271 }
272 }
273 }