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.numerics.multipole;
39
40 import static ffx.numerics.multipole.GKMultipoleOrder.DIPOLE;
41 import static ffx.numerics.multipole.GKMultipoleOrder.MONOPOLE;
42 import static ffx.numerics.multipole.GKMultipoleOrder.QUADRUPOLE;
43 import static ffx.numerics.multipole.GKTensorMode.BORN;
44 import static ffx.numerics.multipole.GKTensorMode.POTENTIAL;
45 import static java.util.Arrays.fill;
46
47
48
49
50 public class GKEnergyGlobal {
51
52 private final GKSource gkSource;
53 private final GKTensorGlobal gkMonopole;
54 private final GKTensorGlobal gkDipole;
55 private final GKTensorGlobal gkQuadrupole;
56
57
58
59
60
61
62
63
64 public GKEnergyGlobal(double gkc, double epsilon, boolean gradient) {
65 int monopoleOrder = 2;
66 int dipoleOrder = 3;
67 int quadrupoleOrder = 4;
68 if (gradient) {
69 monopoleOrder = 3;
70 dipoleOrder = 4;
71 quadrupoleOrder = 5;
72 }
73 gkSource = new GKSource(quadrupoleOrder, gkc);
74 gkMonopole = new GKTensorGlobal(MONOPOLE, monopoleOrder, gkSource, 1.0, epsilon);
75 gkDipole = new GKTensorGlobal(DIPOLE, dipoleOrder, gkSource, 1.0, epsilon);
76 gkQuadrupole = new GKTensorGlobal(QUADRUPOLE, quadrupoleOrder, gkSource, 1.0, epsilon);
77 }
78
79
80
81
82
83
84
85
86
87 public void initPotential(double[] r, double r2, double rbi, double rbk) {
88 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, rbi, rbk);
89 gkMonopole.setR(r);
90 gkDipole.setR(r);
91 gkQuadrupole.setR(r);
92 gkMonopole.generateTensor();
93 gkDipole.generateTensor();
94 gkQuadrupole.generateTensor();
95 }
96
97
98
99
100
101
102
103
104
105 public void initBorn(double[] r, double r2, double rbi, double rbk) {
106 gkSource.generateSource(BORN, QUADRUPOLE, r2, rbi, rbk);
107 gkMonopole.setR(r);
108 gkDipole.setR(r);
109 gkQuadrupole.setR(r);
110 gkMonopole.generateTensor();
111 gkDipole.generateTensor();
112 gkQuadrupole.generateTensor();
113 }
114
115
116
117
118
119
120
121
122 public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
123 double em = gkMonopole.multipoleEnergy(mI, mK);
124 double ed = gkDipole.multipoleEnergy(mI, mK);
125 double eq = gkQuadrupole.multipoleEnergy(mI, mK);
126 return em + ed + eq;
127 }
128
129
130
131
132
133
134
135
136 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
137 double emp = gkMonopole.polarizationEnergy(mI, mK);
138 double edp = gkDipole.polarizationEnergy(mI, mK);
139 double eqp = gkQuadrupole.polarizationEnergy(mI, mK);
140 return emp + edp + eqp;
141 }
142
143
144
145
146
147
148
149
150
151
152
153 public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
154 double[] gradI, double[] torqueI, double[] torqueK) {
155 double[] gI = new double[3];
156 double[] gK = new double[3];
157 double[] tI = new double[3];
158 double[] tK = new double[3];
159 double em = gkMonopole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
160 for (int j = 0; j < 3; j++) {
161 gradI[j] = gI[j];
162 torqueI[j] = tI[j];
163 torqueK[j] = tK[j];
164 }
165 fill(gI, 0.0);
166 fill(gK, 0.0);
167 fill(tI, 0.0);
168 fill(tK, 0.0);
169 double ed = gkDipole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
170 for (int j = 0; j < 3; j++) {
171 gradI[j] += gI[j];
172 torqueI[j] += tI[j];
173 torqueK[j] += tK[j];
174 }
175 fill(gI, 0.0);
176 fill(gK, 0.0);
177 fill(tI, 0.0);
178 fill(tK, 0.0);
179 double eq = gkQuadrupole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
180 for (int j = 0; j < 3; j++) {
181 gradI[j] += gI[j];
182 torqueI[j] += tI[j];
183 torqueK[j] += tK[j];
184 }
185
186 return em + ed + eq;
187 }
188
189
190
191
192
193
194
195
196
197
198
199
200 public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK, double mutualMask,
201 double[] gradI, double[] torqueI, double[] torqueK) {
202 double[] gI = new double[3];
203 double[] tI = new double[3];
204 double[] tK = new double[3];
205 double emp = gkMonopole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
206 for (int j = 0; j < 3; j++) {
207 gradI[j] = gI[j];
208 torqueI[j] = tI[j];
209 torqueK[j] = tK[j];
210 }
211 fill(gI, 0.0);
212 fill(tI, 0.0);
213 fill(tK, 0.0);
214 double edp = gkDipole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
215 for (int j = 0; j < 3; j++) {
216 gradI[j] += gI[j];
217 torqueI[j] += tI[j];
218 torqueK[j] += tK[j];
219 }
220 fill(gI, 0.0);
221 fill(tI, 0.0);
222 fill(tK, 0.0);
223 double eqp = gkQuadrupole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
224 for (int j = 0; j < 3; j++) {
225 gradI[j] += gI[j];
226 torqueI[j] += tI[j];
227 torqueK[j] += tK[j];
228 }
229
230
231 return emp + edp + eqp;
232 }
233
234
235
236
237
238
239
240
241 public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
242 double db = gkMonopole.multipoleEnergyBornGrad(mI, mK);
243 db += gkDipole.multipoleEnergyBornGrad(mI, mK);
244 db += gkQuadrupole.multipoleEnergyBornGrad(mI, mK);
245 return db;
246 }
247
248
249
250
251
252
253
254
255
256 public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK,
257 boolean mutual) {
258
259 double db = gkMonopole.polarizationEnergyBornGrad(mI, mK);
260 db += gkDipole.polarizationEnergyBornGrad(mI, mK);
261 db += gkQuadrupole.polarizationEnergyBornGrad(mI, mK);
262
263 if (mutual) {
264 db += gkDipole.mutualPolarizationEnergyBornGrad(mI, mK);
265 }
266 return db;
267 }
268
269 }