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;
39
40 import static ffx.numerics.math.DoubleMath.length2;
41 import static java.util.Arrays.fill;
42 import static org.apache.commons.math3.util.FastMath.pow;
43
44 import ffx.potential.MolecularAssembly;
45 import ffx.potential.bonded.Atom;
46 import ffx.potential.bonded.LambdaInterface;
47 import ffx.potential.bonded.MSNode;
48 import ffx.potential.bonded.Polymer;
49 import ffx.potential.parameters.ForceField;
50 import java.util.List;
51 import java.util.logging.Logger;
52
53
54
55
56
57
58
59 public class COMRestraint implements LambdaInterface {
60
61 private static final Logger logger = Logger.getLogger(COMRestraint.class.getName());
62 private final Atom[] atoms;
63 private final int nAtoms;
64 private final Polymer[] polymers;
65 private final List<MSNode> molecules;
66 private final List<MSNode> water;
67 private final List<MSNode> ions;
68 private final int nMolecules;
69
70 private final double forceConstant;
71
72 private final double[][] initialCOM;
73 private final double[][] currentCOM;
74 private final double[] dx = new double[3];
75 private final double[] dcomdx;
76 private final double lambdaExp = 1.0;
77 private final double[] lambdaGradient;
78 private double lambda = 1.0;
79 private double lambdaPow = pow(lambda, lambdaExp);
80 private double dLambdaPow = lambdaExp * pow(lambda, lambdaExp - 1.0);
81 private double d2LambdaPow = 0;
82 private double dEdL = 0.0;
83 private double d2EdL2 = 0.0;
84 private boolean lambdaTerm;
85
86
87
88
89
90
91
92 public COMRestraint(MolecularAssembly molecularAssembly) {
93 this.atoms = molecularAssembly.getAtomArray();
94 nAtoms = atoms.length;
95 this.polymers = molecularAssembly.getChains();
96 this.molecules = molecularAssembly.getMolecules();
97 this.water = molecularAssembly.getWater();
98 this.ions = molecularAssembly.getIons();
99 ForceField forceField = molecularAssembly.getForceField();
100
101 nMolecules = countMolecules();
102 initialCOM = new double[3][nMolecules];
103 currentCOM = new double[3][nMolecules];
104
105 lambdaTerm = forceField.getBoolean("LAMBDATERM", false);
106 if (lambdaTerm) {
107 lambdaGradient = new double[nAtoms * 3];
108 } else {
109 lambdaGradient = null;
110 lambda = 1.0;
111 lambdaPow = 1.0;
112 dLambdaPow = 0.0;
113 d2LambdaPow = 0.0;
114 }
115 dcomdx = new double[nAtoms];
116 forceConstant = forceField.getDouble("RESTRAINT_K", 10.0);
117
118 computeCOM(initialCOM, nMolecules);
119
120 logger.info("\n COM restraint initialized");
121 }
122
123
124 @Override
125 public double getLambda() {
126 return lambda;
127 }
128
129
130 @Override
131 public void setLambda(double lambda) {
132 if (lambdaTerm) {
133 this.lambda = lambda;
134
135 double lambdaWindow = 1.0;
136
137 if (this.lambda <= lambdaWindow) {
138 double dldgl = 1.0 / lambdaWindow;
139 double l = dldgl * this.lambda;
140 double l2 = l * l;
141 double l3 = l2 * l;
142 double l4 = l2 * l2;
143 double l5 = l4 * l;
144 double c3 = 10.0;
145 double c4 = -15.0;
146 double c5 = 6.0;
147 double threeC3 = 3.0 * c3;
148 double sixC3 = 6.0 * c3;
149 double fourC4 = 4.0 * c4;
150 double twelveC4 = 12.0 * c4;
151 double fiveC5 = 5.0 * c5;
152 double twentyC5 = 20.0 * c5;
153 lambdaPow = c3 * l3 + c4 * l4 + c5 * l5;
154 dLambdaPow = (threeC3 * l2 + fourC4 * l3 + fiveC5 * l4) * dldgl;
155 d2LambdaPow = (sixC3 * l + twelveC4 * l2 + twentyC5 * l3) * dldgl * dldgl;
156 } else {
157 lambdaPow = 1.0;
158 dLambdaPow = 0.0;
159 d2LambdaPow = 0.0;
160 }
161 } else {
162 this.lambda = 1.0;
163 lambdaPow = 1.0;
164 dLambdaPow = 0.0;
165 d2LambdaPow = 0.0;
166 }
167 }
168
169
170 @Override
171 public double getd2EdL2() {
172 if (lambdaTerm) {
173 return d2EdL2;
174 } else {
175 return 0.0;
176 }
177 }
178
179
180 @Override
181 public double getdEdL() {
182 if (lambdaTerm) {
183 return dEdL;
184 } else {
185 return 0.0;
186 }
187 }
188
189
190 @Override
191 public void getdEdXdL(double[] gradient) {
192 if (lambdaTerm) {
193 int n3 = nAtoms * 3;
194 for (int i = 0; i < n3; i++) {
195 gradient[i] += lambdaGradient[i];
196 }
197 }
198 }
199
200
201
202
203
204
205
206
207 public double residual(boolean gradient, boolean print) {
208 if (lambdaTerm) {
209 dEdL = 0.0;
210 d2EdL2 = 0.0;
211 fill(lambdaGradient, 0.0);
212 }
213 double residual = 0.0;
214 double fx2 = forceConstant * 2.0;
215 computeCOM(currentCOM, nMolecules);
216 computedcomdx();
217 for (int i = 0; i < nMolecules; i++) {
218 dx[0] = currentCOM[0][i] - initialCOM[0][i];
219 dx[1] = currentCOM[1][i] - initialCOM[1][i];
220 dx[2] = currentCOM[2][i] - initialCOM[2][i];
221
222 double r2 = length2(dx);
223 residual += r2;
224 for (int j = 0; j < nAtoms; j++) {
225 if (gradient || lambdaTerm) {
226 final double dedx = dx[0] * fx2 * dcomdx[j];
227 final double dedy = dx[1] * fx2 * dcomdx[j];
228 final double dedz = dx[2] * fx2 * dcomdx[j];
229
230 Atom atom = atoms[j];
231 if (gradient) {
232 atom.addToXYZGradient(lambdaPow * dedx, lambdaPow * dedy, lambdaPow * dedz);
233 }
234 if (lambdaTerm) {
235 int j3 = i * 3;
236 lambdaGradient[j3] = dLambdaPow * dedx;
237 lambdaGradient[j3 + 1] = dLambdaPow * dedy;
238 lambdaGradient[j3 + 2] = dLambdaPow * dedz;
239 }
240 }
241 }
242 }
243 if (lambdaTerm) {
244 dEdL = dLambdaPow * forceConstant * residual;
245 d2EdL2 = d2LambdaPow * forceConstant * residual;
246 }
247 return forceConstant * residual * lambdaPow;
248 }
249
250
251
252
253
254
255 public void setLambdaTerm(boolean lambdaTerm) {
256 this.lambdaTerm = lambdaTerm;
257 setLambda(lambda);
258 }
259
260 private void computeCOM(double[][] com, int nMolecules) {
261 int i = 0;
262 while (i < nMolecules) {
263 if (polymers != null && polymers.length > 0) {
264
265 for (Polymer polymer : polymers) {
266 List<Atom> list = polymer.getAtomList();
267 com[0][i] = 0.0;
268 com[1][i] = 0.0;
269 com[2][i] = 0.0;
270 double totalMass = 0.0;
271 for (Atom atom : list) {
272 double m = atom.getMass();
273 com[0][i] += atom.getX() * m;
274 com[1][i] += atom.getY() * m;
275 com[2][i] += atom.getZ() * m;
276 totalMass += m;
277 }
278 com[0][i] /= totalMass;
279 com[1][i] /= totalMass;
280 com[2][i] /= totalMass;
281 i++;
282 }
283 }
284
285
286 for (MSNode molecule : molecules) {
287 List<Atom> list = molecule.getAtomList();
288
289 com[0][i] = 0.0;
290 com[1][i] = 0.0;
291 com[2][i] = 0.0;
292 double totalMass = 0.0;
293 for (Atom atom : list) {
294 double m = atom.getMass();
295 com[0][i] += atom.getX() * m;
296 com[1][i] += atom.getY() * m;
297 com[2][i] += atom.getZ() * m;
298 totalMass += m;
299 }
300 com[0][i] /= totalMass;
301 com[1][i] /= totalMass;
302 com[2][i] /= totalMass;
303 i++;
304 }
305
306
307 for (MSNode water : water) {
308 List<Atom> list = water.getAtomList();
309
310 com[0][i] = 0.0;
311 com[1][i] = 0.0;
312 com[2][i] = 0.0;
313 double totalMass = 0.0;
314 for (Atom atom : list) {
315 double m = atom.getMass();
316 com[0][i] += atom.getX() * m;
317 com[1][i] += atom.getY() * m;
318 com[2][i] += atom.getZ() * m;
319 totalMass += m;
320 }
321 com[0][i] /= totalMass;
322 com[1][i] /= totalMass;
323 com[2][i] /= totalMass;
324 i++;
325 }
326
327
328 for (MSNode ion : ions) {
329 List<Atom> list = ion.getAtomList();
330
331 com[0][i] = 0.0;
332 com[1][i] = 0.0;
333 com[2][i] = 0.0;
334 double totalMass = 0.0;
335 for (Atom atom : list) {
336 double m = atom.getMass();
337 com[0][i] += atom.getX() * m;
338 com[1][i] += atom.getY() * m;
339 com[2][i] += atom.getZ() * m;
340 totalMass += m;
341 }
342 com[0][i] /= totalMass;
343 com[1][i] /= totalMass;
344 com[2][i] /= totalMass;
345 i++;
346 }
347 }
348 }
349
350 private void computedcomdx() {
351
352 int i = 0;
353 while (i < nAtoms) {
354 if (polymers != null && polymers.length > 0) {
355 for (Polymer polymer : polymers) {
356 List<Atom> list = polymer.getAtomList();
357 double totalMass = 0.0;
358 for (Atom atom : list) {
359 double m = atom.getMass();
360 totalMass += m;
361 }
362 for (Atom atom : list) {
363 dcomdx[i] = atom.getMass();
364 dcomdx[i] /= totalMass;
365 i++;
366 }
367 }
368 }
369
370
371 for (MSNode molecule : molecules) {
372 List<Atom> list = molecule.getAtomList();
373 double totalMass = 0.0;
374 for (Atom atom : list) {
375 double m = atom.getMass();
376 totalMass += m;
377 }
378 for (Atom atom : list) {
379 dcomdx[i] = atom.getMass();
380 dcomdx[i] /= totalMass;
381 i++;
382 }
383 }
384
385
386 for (MSNode water : water) {
387 List<Atom> list = water.getAtomList();
388 double totalMass = 0.0;
389 for (Atom atom : list) {
390 double m = atom.getMass();
391 totalMass += m;
392 }
393 for (Atom atom : list) {
394 dcomdx[i] = atom.getMass();
395 dcomdx[i] /= totalMass;
396 i++;
397 }
398 }
399
400
401 for (MSNode ion : ions) {
402 List<Atom> list = ion.getAtomList();
403 double totalMass = 0.0;
404 for (Atom atom : list) {
405 double m = atom.getMass();
406 totalMass += m;
407 }
408 for (Atom atom : list) {
409 dcomdx[i] = atom.getMass();
410 dcomdx[i] /= totalMass;
411 i++;
412 }
413 }
414 }
415
416
417
418
419
420 }
421
422 private int countMolecules() {
423 int count = 0;
424
425 if (polymers != null && polymers.length > 0) {
426 count += polymers.length;
427 }
428 if (molecules != null) {
429 count += molecules.size();
430 }
431 if (water != null) {
432 count += water.size();
433 }
434 if (ions != null) {
435 count += ions.size();
436 }
437 return count;
438 }
439 }