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.GKSource.GK_MULTIPOLE_ORDER.DIPOLE;
41 import static ffx.numerics.multipole.GKSource.GK_MULTIPOLE_ORDER.MONOPOLE;
42 import static ffx.numerics.multipole.GKSource.GK_MULTIPOLE_ORDER.QUADRUPOLE;
43 import static ffx.numerics.multipole.GKSource.GK_TENSOR_MODE.BORN;
44 import static ffx.numerics.multipole.GKSource.GK_TENSOR_MODE.POTENTIAL;
45 import static java.util.Arrays.fill;
46
47 public class GKEnergyQI {
48
49 private final GKSource gkSource;
50 private final GKTensorQI gkMonopole;
51 private final GKTensorQI gkDipole;
52 private final GKTensorQI gkQuadrupole;
53
54
55
56
57
58
59
60
61
62 public GKEnergyQI(double soluteDielectric, double solventDielectric, double gkc, boolean gradient) {
63 int monopoleOrder = 2;
64 int dipoleOrder = 3;
65 int quadrupoleOrder = 4;
66 if (gradient) {
67 monopoleOrder = 3;
68 dipoleOrder = 4;
69 quadrupoleOrder = 5;
70 }
71 gkSource = new GKSource(quadrupoleOrder, gkc);
72 gkMonopole = new GKTensorQI(MONOPOLE, monopoleOrder, gkSource, soluteDielectric, solventDielectric);
73 gkDipole = new GKTensorQI(DIPOLE, dipoleOrder, gkSource, soluteDielectric, solventDielectric);
74 gkQuadrupole = new GKTensorQI(QUADRUPOLE, quadrupoleOrder, gkSource, soluteDielectric, solventDielectric);
75 }
76
77 public void initPotential(double[] r, double r2, double rbi, double rbk) {
78 gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, rbi, rbk);
79 gkMonopole.setR(r);
80 gkDipole.setR(r);
81 gkQuadrupole.setR(r);
82 gkMonopole.generateTensor();
83 gkDipole.generateTensor();
84 gkQuadrupole.generateTensor();
85 }
86
87 public void initBorn(double[] r, double r2, double rbi, double rbk) {
88 gkSource.generateSource(BORN, 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 public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
98 double em = gkMonopole.multipoleEnergy(mI, mK);
99 double ed = gkDipole.multipoleEnergy(mI, mK);
100 double eq = gkQuadrupole.multipoleEnergy(mI, mK);
101 return em + ed + eq;
102 }
103
104 public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
105 double emp = gkMonopole.polarizationEnergy(mI, mK);
106 double edp = gkDipole.polarizationEnergy(mI, mK);
107 double eqp = gkQuadrupole.polarizationEnergy(mI, mK);
108 return emp + edp + eqp;
109 }
110
111 public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
112 double[] gradI, double[] torqueI, double[] torqueK) {
113 double[] gI = new double[3];
114 double[] gK = new double[3];
115 double[] tI = new double[3];
116 double[] tK = new double[3];
117 double em = gkMonopole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
118 for (int j = 0; j < 3; j++) {
119 gradI[j] = gI[j];
120 torqueI[j] = tI[j];
121 torqueK[j] = tK[j];
122 }
123 fill(gI, 0.0);
124 fill(gK, 0.0);
125 fill(tI, 0.0);
126 fill(tK, 0.0);
127 double ed = gkDipole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
128 for (int j = 0; j < 3; j++) {
129 gradI[j] += gI[j];
130 torqueI[j] += tI[j];
131 torqueK[j] += tK[j];
132 }
133 fill(gI, 0.0);
134 fill(gK, 0.0);
135 fill(tI, 0.0);
136 fill(tK, 0.0);
137 double eq = gkQuadrupole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
138 for (int j = 0; j < 3; j++) {
139 gradI[j] += gI[j];
140 torqueI[j] += tI[j];
141 torqueK[j] += tK[j];
142 }
143 return em + ed + eq;
144 }
145
146 public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
147 double mutualMask, double[] gradI, double[] torqueI, double[] torqueK) {
148 double[] gI = new double[3];
149 double[] tI = new double[3];
150 double[] tK = new double[3];
151 double emp = gkMonopole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
152 for (int j = 0; j < 3; j++) {
153 gradI[j] = gI[j];
154 torqueI[j] = tI[j];
155 torqueK[j] = tK[j];
156 }
157 fill(gI, 0.0);
158 fill(tI, 0.0);
159 fill(tK, 0.0);
160 double edp = gkDipole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, 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(tI, 0.0);
168 fill(tK, 0.0);
169 double eqp = gkQuadrupole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, 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
176
177 return emp + edp + eqp;
178 }
179
180 public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
181 double db = gkMonopole.multipoleEnergyBornGrad(mI, mK);
182 db += gkDipole.multipoleEnergyBornGrad(mI, mK);
183 db += gkQuadrupole.multipoleEnergyBornGrad(mI, mK);
184 return db;
185 }
186
187 public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK,
188 boolean mutual) {
189
190 double db = gkMonopole.polarizationEnergyBornGrad(mI, mK);
191 db += gkDipole.polarizationEnergyBornGrad(mI, mK);
192 db += gkQuadrupole.polarizationEnergyBornGrad(mI, mK);
193
194 if (mutual) {
195 db += gkDipole.mutualPolarizationEnergyBornGrad(mI, mK);
196 }
197 return db;
198 }
199
200 }