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 final 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 public double getEps() {
111 double eps = sharedEps.get();
112 double epsCR = sharedEpsCR.get();
113 return max(eps, epsCR);
114 }
115
116 public double getSOR() {
117 return polsor;
118 }
119
120 public void init(
121 Atom[] atoms,
122 double[] polarizability,
123 double[][][] inducedDipole,
124 double[][][] inducedDipoleCR,
125 double[][] directDipole,
126 double[][] directDipoleCR,
127 double[][] cartesianDipolePhi,
128 double[][] cartesianDipolePhiCR,
129 AtomicDoubleArray3D field,
130 AtomicDoubleArray3D fieldCR,
131 boolean generalizedKirkwoodTerm,
132 GeneralizedKirkwood generalizedKirkwood,
133 EwaldParameters ewaldParameters) {
134 this.atoms = atoms;
135 this.polarizability = polarizability;
136 this.inducedDipole = inducedDipole;
137 this.inducedDipoleCR = inducedDipoleCR;
138 this.directDipole = directDipole;
139 this.directDipoleCR = directDipoleCR;
140 this.cartesianDipolePhi = cartesianDipolePhi;
141 this.cartesianDipolePhiCR = cartesianDipolePhiCR;
142 this.field = field;
143 this.fieldCR = fieldCR;
144 this.generalizedKirkwoodTerm = generalizedKirkwoodTerm;
145 this.generalizedKirkwood = generalizedKirkwood;
146 this.aewald = ewaldParameters.aewald;
147 this.aewald3 = ewaldParameters.aewald3;
148 }
149
150 @Override
151 public void run() throws Exception {
152 try {
153 int ti = getThreadIndex();
154 if (sorLoop[ti] == null) {
155 sorLoop[ti] = new SORLoop();
156 }
157 int nAtoms = atoms.length;
158 execute(0, nAtoms - 1, sorLoop[ti]);
159 } catch (RuntimeException ex) {
160 logger.warning(
161 "Fatal exception computing the mutual induced dipoles in thread " + getThreadIndex());
162 throw ex;
163 } catch (Exception e) {
164 String message =
165 "Fatal exception computing the mutual induced dipoles in thread "
166 + getThreadIndex()
167 + "\n";
168 logger.log(Level.SEVERE, message, e);
169 }
170 }
171
172 @Override
173 public void start() {
174 sharedEps.set(0.0);
175 sharedEpsCR.set(0.0);
176 }
177
178 private class SORLoop extends IntegerForLoop {
179
180 private double eps, epsCR;
181
182 @Override
183 public void finish() {
184 sharedEps.addAndGet(eps);
185 sharedEpsCR.addAndGet(epsCR);
186 }
187
188 @Override
189 public void run(int lb, int ub) throws Exception {
190 int threadID = getThreadIndex();
191 final double[][] induced0 = inducedDipole[0];
192 final double[][] inducedCR0 = inducedDipoleCR[0];
193 if (aewald > 0.0) {
194
195 for (int i = lb; i <= ub; i++) {
196 double[] dipolei = induced0[i];
197 double[] dipoleCRi = inducedCR0[i];
198 final double[] phii = cartesianDipolePhi[i];
199 final double[] phiCRi = cartesianDipolePhiCR[i];
200 double fx = aewald3 * dipolei[0] - phii[t100];
201 double fy = aewald3 * dipolei[1] - phii[t010];
202 double fz = aewald3 * dipolei[2] - phii[t001];
203 double fxCR = aewald3 * dipoleCRi[0] - phiCRi[t100];
204 double fyCR = aewald3 * dipoleCRi[1] - phiCRi[t010];
205 double fzCR = aewald3 * dipoleCRi[2] - phiCRi[t001];
206 field.add(threadID, i, fx, fy, fz);
207 fieldCR.add(threadID, i, fxCR, fyCR, fzCR);
208 }
209 }
210
211 if (generalizedKirkwoodTerm) {
212 AtomicDoubleArray3D fieldGK = generalizedKirkwood.getFieldGK();
213 AtomicDoubleArray3D fieldGKCR = generalizedKirkwood.getFieldGKCR();
214
215 for (int i = lb; i <= ub; i++) {
216 field.add(threadID, i, fieldGK.getX(i), fieldGK.getY(i), fieldGK.getZ(i));
217 fieldCR.add(threadID, i, fieldGKCR.getX(i), fieldGKCR.getY(i), fieldGKCR.getZ(i));
218 }
219 }
220
221
222 field.reduce(lb, ub);
223 fieldCR.reduce(lb, ub);
224
225
226 for (int i = lb; i <= ub; i++) {
227 final double[] ind = induced0[i];
228 final double[] indCR = inducedCR0[i];
229 final double[] direct = directDipole[i];
230 final double[] directCR = directDipoleCR[i];
231 final double polar = polarizability[i];
232 for (int j = 0; j < 3; j++) {
233 double previous = ind[j];
234 double mutual = polar * field.get(j, i);
235 ind[j] = direct[j] + mutual;
236 double delta = polsor * (ind[j] - previous);
237 ind[j] = previous + delta;
238 eps += delta * delta;
239 previous = indCR[j];
240 mutual = polar * fieldCR.get(j, i);
241 indCR[j] = directCR[j] + mutual;
242 delta = polsor * (indCR[j] - previous);
243 indCR[j] = previous + delta;
244 epsCR += delta * delta;
245 }
246 }
247 }
248
249 @Override
250 public IntegerSchedule schedule() {
251 return IntegerSchedule.fixed();
252 }
253
254 @Override
255 public void start() {
256 eps = 0.0;
257 epsCR = 0.0;
258 }
259 }
260 }