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