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;
39
40 import static org.apache.commons.math3.util.FastMath.toDegrees;
41
42 import ffx.crystal.Crystal;
43 import ffx.crystal.SpaceGroup;
44 import ffx.numerics.Potential;
45 import ffx.potential.MolecularAssembly.FractionalMode;
46 import ffx.potential.bonded.Atom;
47
48
49
50
51
52
53
54
55 public class XtalEnergy implements Potential {
56
57
58
59
60 private final MolecularAssembly molecularAssembly;
61
62
63
64 private final ForceFieldEnergy forceFieldEnergy;
65
66
67
68 private final Atom[] activeAtoms;
69
70
71
72 private final int nActive;
73
74
75
76 private final double[] xyz;
77 private final double[] gr;
78 private final int nParams;
79
80 private final Crystal crystal;
81 private final VARIABLE_TYPE[] type;
82 private final double[] mass;
83 private final Crystal unitCell;
84 private double[] scaling;
85 private double totalEnergy;
86
87 private FractionalMode fractionalMode = FractionalMode.OFF;
88
89
90
91
92
93
94
95 public XtalEnergy(ForceFieldEnergy forceFieldEnergy, MolecularAssembly molecularAssembly) {
96 this.forceFieldEnergy = forceFieldEnergy;
97 this.molecularAssembly = molecularAssembly;
98 Atom[] atoms = molecularAssembly.getAtomArray();
99
100 int n = 0;
101 for (Atom a : atoms) {
102 if (a.isActive()) {
103 n++;
104 }
105 }
106 nActive = n;
107
108 activeAtoms = new Atom[nActive];
109 int index = 0;
110 for (Atom a : atoms) {
111 if (a.isActive()) {
112 activeAtoms[index++] = a;
113 }
114 }
115
116 nParams = 3 * nActive + 6;
117 crystal = forceFieldEnergy.getCrystal();
118 unitCell = crystal.getUnitCell();
119 xyz = new double[3 * nActive];
120 gr = new double[3 * nActive];
121 type = new VARIABLE_TYPE[nParams];
122 mass = new double[nParams];
123
124 index = 0;
125 for (int i = 0; i < nActive; i++) {
126 double m = activeAtoms[i].getMass();
127 mass[index] = m;
128 mass[index + 1] = m;
129 mass[index + 2] = m;
130 type[index] = VARIABLE_TYPE.X;
131 type[index + 1] = VARIABLE_TYPE.Y;
132 type[index + 2] = VARIABLE_TYPE.Z;
133 index += 3;
134 }
135 for (int i = nActive * 3; i < nActive * 3 + 6; i++) {
136 mass[i] = 1.0;
137 type[i] = VARIABLE_TYPE.OTHER;
138 }
139 }
140
141
142
143
144 @Override
145 public boolean destroy() {
146 return forceFieldEnergy.destroy();
147 }
148
149
150
151
152 @Override
153 public double energy(double[] x) {
154
155
156 unscaleCoordinates(x);
157
158
159 setCoordinates(x);
160
161 totalEnergy = forceFieldEnergy.energy(false, false);
162
163
164 scaleCoordinates(x);
165
166 return totalEnergy;
167 }
168
169
170
171
172 @Override
173 public double energyAndGradient(double[] x, double[] g) {
174
175 unscaleCoordinates(x);
176
177
178 setCoordinates(x);
179
180
181 double e = forceFieldEnergy.energyAndGradient(xyz, gr);
182
183
184 packGradient(x, g);
185
186
187 unitCellParameterDerivatives(x, g);
188
189 totalEnergy = e;
190
191 return totalEnergy;
192 }
193
194
195
196
197 @Override
198 public double[] getAcceleration(double[] acceleration) {
199 throw new UnsupportedOperationException("Not supported yet.");
200 }
201
202
203
204
205 @Override
206 public double[] getCoordinates(double[] x) {
207 int n = getNumberOfVariables();
208 if (x == null || x.length < n) {
209 x = new double[n];
210 }
211 int index = 0;
212 for (int i = 0; i < nActive; i++) {
213 Atom a = activeAtoms[i];
214 x[index] = a.getX();
215 index++;
216 x[index] = a.getY();
217 index++;
218 x[index] = a.getZ();
219 index++;
220 }
221 x[index] = unitCell.a;
222 index++;
223 x[index] = unitCell.b;
224 index++;
225 x[index] = unitCell.c;
226 index++;
227 x[index] = unitCell.alpha;
228 index++;
229 x[index] = unitCell.beta;
230 index++;
231 x[index] = unitCell.gamma;
232 return x;
233 }
234
235
236
237
238 @Override
239 public STATE getEnergyTermState() {
240 return forceFieldEnergy.getEnergyTermState();
241 }
242
243
244
245
246 @Override
247 public void setEnergyTermState(STATE state) {
248 forceFieldEnergy.setEnergyTermState(state);
249 }
250
251
252
253
254 @Override
255 public double[] getMass() {
256 return mass;
257 }
258
259
260
261
262 @Override
263 public int getNumberOfVariables() {
264 return nParams;
265 }
266
267
268
269
270 @Override
271 public double[] getPreviousAcceleration(double[] previousAcceleration) {
272 throw new UnsupportedOperationException("Not supported yet.");
273 }
274
275
276
277
278 @Override
279 public double[] getScaling() {
280 return scaling;
281 }
282
283
284
285
286 @Override
287 public void setScaling(double[] scaling) {
288 this.scaling = scaling;
289 }
290
291
292
293
294 @Override
295 public double getTotalEnergy() {
296 return totalEnergy;
297 }
298
299
300
301
302 @Override
303 public VARIABLE_TYPE[] getVariableTypes() {
304 return type;
305 }
306
307
308
309
310 @Override
311 public double[] getVelocity(double[] velocity) {
312 throw new UnsupportedOperationException("Not supported yet.");
313 }
314
315
316
317
318 @Override
319 public void setAcceleration(double[] acceleration) {
320 throw new UnsupportedOperationException("Not supported yet.");
321 }
322
323
324
325
326
327
328 public void setFractionalCoordinateMode(FractionalMode fractionalMode) {
329 this.fractionalMode = fractionalMode;
330 molecularAssembly.setFractionalMode(fractionalMode);
331 }
332
333
334
335
336 @Override
337 public void setPreviousAcceleration(double[] previousAcceleration) {
338 throw new UnsupportedOperationException("Not supported yet.");
339 }
340
341
342
343
344 @Override
345 public void setVelocity(double[] velocity) {
346 throw new UnsupportedOperationException("Not supported yet.");
347 }
348
349
350
351
352
353
354
355 private void unitCellParameterDerivatives(double[] x, double[] g) {
356
357 double eps = 1.0e-5;
358 double deps = toDegrees(eps);
359 SpaceGroup spaceGroup = crystal.spaceGroup;
360
361 int index = 3 * nActive;
362 switch (spaceGroup.latticeSystem) {
363 case TRICLINIC_LATTICE -> {
364 g[index] = finiteDifference(x, index, eps);
365 index++;
366 g[index] = finiteDifference(x, index, eps);
367 index++;
368 g[index] = finiteDifference(x, index, eps);
369 index++;
370 g[index] = finiteDifference(x, index, deps);
371 index++;
372 g[index] = finiteDifference(x, index, deps);
373 index++;
374 g[index] = finiteDifference(x, index, deps);
375 }
376 case MONOCLINIC_LATTICE -> {
377
378 g[index] = finiteDifference(x, index, eps);
379 index++;
380 g[index] = finiteDifference(x, index, eps);
381 index++;
382 g[index] = finiteDifference(x, index, eps);
383 index++;
384 g[index] = 0.0;
385 index++;
386 g[index] = finiteDifference(x, index, deps);
387 index++;
388 g[index] = 0.0;
389 }
390 case ORTHORHOMBIC_LATTICE -> {
391
392 g[index] = finiteDifference(x, index, eps);
393 index++;
394 g[index] = finiteDifference(x, index, eps);
395 index++;
396 g[index] = finiteDifference(x, index, eps);
397 index++;
398 g[index] = 0.0;
399 index++;
400 g[index] = 0.0;
401 index++;
402 g[index] = 0.0;
403 }
404 case TETRAGONAL_LATTICE,
405 HEXAGONAL_LATTICE
406 -> {
407 g[index] = finiteDifference2(x, index, index + 1, eps);
408 index++;
409 g[index] = g[index - 1];
410 index++;
411 g[index] = finiteDifference(x, index, eps);
412 index++;
413 g[index] = 0.0;
414 index++;
415 g[index] = 0.0;
416 index++;
417 g[index] = 0.0;
418 }
419 case RHOMBOHEDRAL_LATTICE -> {
420
421 g[index] = finiteDifference3(x, index, index + 1, index + 2, eps);
422 index++;
423 g[index] = g[index - 1];
424 index++;
425 g[index] = g[index - 2];
426 index++;
427 g[index] = finiteDifference3(x, index, index + 1, index + 2, deps);
428 index++;
429 g[index] = g[index - 1];
430 index++;
431 g[index] = g[index - 2];
432 }
433
434 case CUBIC_LATTICE -> {
435
436 g[index] = finiteDifference3(x, index, index + 1, index + 2, eps);
437 index++;
438 g[index] = g[index - 1];
439 index++;
440 g[index] = g[index - 2];
441 index++;
442 g[index] = 0.0;
443 index++;
444 g[index] = 0.0;
445 index++;
446 g[index] = 0.0;
447 }
448 }
449
450
451 if (scaling != null) {
452 index = 3 * nActive;
453 g[index] /= scaling[index];
454 index++;
455 g[index] /= scaling[index];
456 index++;
457 g[index] /= scaling[index];
458 index++;
459 g[index] /= scaling[index];
460 index++;
461 g[index] /= scaling[index];
462 index++;
463 g[index] /= scaling[index];
464 }
465 }
466
467
468
469
470
471
472
473
474
475 private double finiteDifference(double[] x, int index, double eps) {
476 double scale = 1.0;
477 if (scaling != null) {
478 scale = scaling[index];
479 }
480 double x1 = x[index];
481 final double param = x1 / scale;
482
483 x[index] = (param + eps / 2.0) * scale;
484 double ePlus = energy(x);
485 x[index] = (param - eps / 2.0) * scale;
486 double eMinus = energy(x);
487
488 x[index] = x1;
489
490 return (ePlus - eMinus) / eps;
491 }
492
493
494
495
496
497
498
499
500 private double finiteDifference2(double[] x, int index1, int index2, double eps) {
501 double scale1 = 1.0;
502 double scale2 = 1.0;
503
504 if (scaling != null) {
505 scale1 = scaling[index1];
506 scale2 = scaling[index2];
507 }
508
509 final double x1 = x[index1];
510 final double x2 = x[index2];
511 final double param1 = x1 / scale1;
512 final double param2 = x2 / scale2;
513
514 x[index1] = (param1 + eps / 2.0) * scale1;
515 x[index2] = (param2 + eps / 2.0) * scale2;
516 double ePlus = energy(x);
517 x[index1] = (param1 - eps / 2.0) * scale1;
518 x[index2] = (param2 - eps / 2.0) * scale2;
519 double eMinus = energy(x);
520
521 x[index1] = x1;
522 x[index2] = x2;
523
524 return (ePlus - eMinus) / eps;
525 }
526
527
528
529
530
531
532
533
534
535 private double finiteDifference3(double[] x, int index1, int index2, int index3, double eps) {
536 double scale1 = 1.0;
537 double scale2 = 1.0;
538 double scale3 = 1.0;
539 if (scaling != null) {
540 scale1 = scaling[index1];
541 scale2 = scaling[index2];
542 scale3 = scaling[index3];
543 }
544
545 final double x1 = x[index1];
546 final double x2 = x[index2];
547 final double x3 = x[index3];
548
549 final double param1 = x1 / scale1;
550 final double param2 = x2 / scale2;
551 final double param3 = x3 / scale3;
552 x[index1] = (param1 + eps / 2.0) * scale1;
553 x[index2] = (param2 + eps / 2.0) * scale2;
554 x[index3] = (param3 + eps / 2.0) * scale3;
555 final double ePlus = energy(x);
556 x[index1] = (param1 - eps / 2.0) * scale1;
557 x[index2] = (param2 - eps / 2.0) * scale2;
558 x[index3] = (param3 - eps / 2.0) * scale3;
559 final double eMinus = energy(x);
560
561 x[index1] = x1;
562 x[index2] = x2;
563 x[index3] = x3;
564
565 return (ePlus - eMinus) / eps;
566 }
567
568
569
570
571
572
573
574 private void packGradient(double[] x, double[] g) {
575
576 if (scaling != null) {
577 int len = x.length;
578 for (int i = 0; i < len; i++) {
579 if (fractionalMode == FractionalMode.OFF) {
580 g[i] /= scaling[i];
581 } else {
582
583
584 g[i] = 0.0;
585 }
586 x[i] *= scaling[i];
587 }
588 }
589 }
590
591
592
593
594
595
596 private void setCoordinates(double[] x) {
597 assert (x != null);
598
599
600 if (fractionalMode != FractionalMode.OFF) {
601 molecularAssembly.computeFractionalCoordinates();
602 }
603
604 int index = nActive * 3;
605 double a = x[index];
606 double b = x[index + 1];
607 double c = x[index + 2];
608 double alpha = x[index + 3];
609 double beta = x[index + 4];
610 double gamma = x[index + 5];
611
612 SpaceGroup spaceGroup = crystal.spaceGroup;
613
614
615 switch (spaceGroup.latticeSystem) {
616 case TRICLINIC_LATTICE -> {
617 }
618 case MONOCLINIC_LATTICE -> {
619
620 alpha = 90.0;
621 gamma = 90.0;
622 }
623 case ORTHORHOMBIC_LATTICE -> {
624
625 alpha = 90.0;
626 beta = 90.0;
627 gamma = 90.0;
628 }
629 case TETRAGONAL_LATTICE -> {
630
631 double ab = 0.5 * (a + b);
632 a = ab;
633 b = ab;
634 alpha = 90.0;
635 beta = 90.0;
636 gamma = 90.0;
637 }
638 case RHOMBOHEDRAL_LATTICE -> {
639
640 double abc = (a + b + c) / 3.0;
641 a = abc;
642 b = abc;
643 c = abc;
644 double angles = (alpha + beta + gamma) / 3.0;
645 alpha = angles;
646 beta = angles;
647 gamma = angles;
648 }
649 case HEXAGONAL_LATTICE -> {
650
651 double ab = 0.5 * (a + b);
652 a = ab;
653 b = ab;
654 alpha = 90.0;
655 beta = 90.0;
656 gamma = 120.0;
657 }
658 case CUBIC_LATTICE -> {
659
660 double abc = (a + b + c) / 3.0;
661 a = abc;
662 b = abc;
663 c = abc;
664 alpha = 90.0;
665 beta = 90.0;
666 gamma = 90.0;
667 }
668 }
669
670 crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
671 forceFieldEnergy.setCrystal(crystal);
672
673
674 if (fractionalMode == FractionalMode.OFF) {
675 index = 0;
676 for (int i = 0; i < nActive; i++) {
677 Atom atom = activeAtoms[i];
678 double xx = x[index];
679 double yy = x[index + 1];
680 double zz = x[index + 2];
681 xyz[index] = xx;
682 xyz[index + 1] = yy;
683 xyz[index + 2] = zz;
684 index += 3;
685 atom.moveTo(xx, yy, zz);
686 }
687 } else {
688
689 molecularAssembly.moveToFractionalCoordinates();
690 index = 0;
691 for (int i = 0; i < nActive; i++) {
692 Atom atom = activeAtoms[i];
693 double xx = atom.getX();
694 double yy = atom.getY();
695 double zz = atom.getZ();
696 xyz[index] = xx;
697 xyz[index + 1] = yy;
698 xyz[index + 2] = zz;
699 index += 3;
700 }
701 }
702 }
703 }