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 GKEnergyQI {
51
52 private final GKSource gkSource;
53 private final GKTensorQI gkMonopole;
54 private final GKTensorQI gkDipole;
55 private final GKTensorQI gkQuadrupole;
56
57
58
59
60
61
62
63
64
65 public GKEnergyQI(double soluteDielectric, double solventDielectric, double gkc, boolean gradient) {
66 int monopoleOrder = 2;
67 int dipoleOrder = 3;
68 int quadrupoleOrder = 4;
69 if (gradient) {
70 monopoleOrder = 3;
71 dipoleOrder = 4;
72 quadrupoleOrder = 5;
73 }
74 gkSource = new GKSource(quadrupoleOrder, gkc);
75 gkMonopole = new GKTensorQI(MONOPOLE, monopoleOrder, gkSource, soluteDielectric, solventDielectric);
76 gkDipole = new GKTensorQI(DIPOLE, dipoleOrder, gkSource, soluteDielectric, solventDielectric);
77 gkQuadrupole = new GKTensorQI(QUADRUPOLE, quadrupoleOrder, gkSource, soluteDielectric, solventDielectric);
78 }
79
80
81
82
83
84
85
86
87
88 public void initPotential(double[] r, double r2, double rbi, double rbk) {
89 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, rbi, rbk);
90 gkMonopole.setR(r);
91 gkDipole.setR(r);
92 gkQuadrupole.setR(r);
93 gkMonopole.generateTensor();
94 gkDipole.generateTensor();
95 gkQuadrupole.generateTensor();
96 }
97
98
99
100
101
102
103
104
105
106 public void initBorn(double[] r, double r2, double rbi, double rbk) {
107 gkSource.generateSource(BORN, QUADRUPOLE, r2, rbi, rbk);
108 gkMonopole.setR(r);
109 gkDipole.setR(r);
110 gkQuadrupole.setR(r);
111 gkMonopole.generateTensor();
112 gkDipole.generateTensor();
113 gkQuadrupole.generateTensor();
114 }
115
116
117
118
119
120
121
122
123 public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
124 double em = gkMonopole.multipoleEnergy(mI, mK);
125 double ed = gkDipole.multipoleEnergy(mI, mK);
126 double eq = gkQuadrupole.multipoleEnergy(mI, mK);
127 return em + ed + eq;
128 }
129
130
131
132
133
134
135
136
137 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
138 double emp = gkMonopole.polarizationEnergy(mI, mK);
139 double edp = gkDipole.polarizationEnergy(mI, mK);
140 double eqp = gkQuadrupole.polarizationEnergy(mI, mK);
141 return emp + edp + eqp;
142 }
143
144
145
146
147
148
149
150
151
152
153
154 public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
155 double[] gradI, double[] torqueI, double[] torqueK) {
156 double[] gI = new double[3];
157 double[] gK = new double[3];
158 double[] tI = new double[3];
159 double[] tK = new double[3];
160 double em = gkMonopole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
161 for (int j = 0; j < 3; j++) {
162 gradI[j] = gI[j];
163 torqueI[j] = tI[j];
164 torqueK[j] = tK[j];
165 }
166 fill(gI, 0.0);
167 fill(gK, 0.0);
168 fill(tI, 0.0);
169 fill(tK, 0.0);
170 double ed = gkDipole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
171 for (int j = 0; j < 3; j++) {
172 gradI[j] += gI[j];
173 torqueI[j] += tI[j];
174 torqueK[j] += tK[j];
175 }
176 fill(gI, 0.0);
177 fill(gK, 0.0);
178 fill(tI, 0.0);
179 fill(tK, 0.0);
180 double eq = gkQuadrupole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
181 for (int j = 0; j < 3; j++) {
182 gradI[j] += gI[j];
183 torqueI[j] += tI[j];
184 torqueK[j] += tK[j];
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,
201 double mutualMask, 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 }