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