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.algorithms.dynamics;
39
40 import ffx.crystal.Crystal;
41 import ffx.crystal.CrystalPotential;
42 import ffx.crystal.SpaceGroup;
43 import ffx.numerics.Potential;
44 import ffx.numerics.math.RunningStatistics;
45 import ffx.potential.MolecularAssembly;
46 import ffx.potential.bonded.Atom;
47 import ffx.potential.parsers.XYZFilter;
48 import org.apache.commons.io.FilenameUtils;
49
50 import java.io.File;
51 import java.util.ArrayList;
52 import java.util.List;
53 import java.util.logging.Level;
54 import java.util.logging.Logger;
55
56 import static ffx.crystal.LatticeSystem.check;
57 import static ffx.numerics.math.ScalarMath.mirrorDegrees;
58 import static ffx.utilities.Constants.AVOGADRO;
59 import static ffx.utilities.Constants.PRESCON;
60 import static ffx.utilities.Constants.R;
61 import static java.lang.String.format;
62 import static org.apache.commons.math3.util.FastMath.exp;
63 import static org.apache.commons.math3.util.FastMath.floor;
64 import static org.apache.commons.math3.util.FastMath.log;
65 import static org.apache.commons.math3.util.FastMath.min;
66 import static org.apache.commons.math3.util.FastMath.random;
67
68
69
70
71
72
73
74
75
76
77 public class Barostat implements CrystalPotential {
78
79 private static final Logger logger = Logger.getLogger(Barostat.class.getName());
80
81
82
83
84 private final MolecularAssembly molecularAssembly;
85
86
87
88 private Crystal crystal;
89
90
91
92 private Crystal unitCell;
93
94
95
96 private final double mass;
97
98
99
100 private final CrystalPotential potential;
101
102
103
104 private final double[] x;
105
106
107
108 private final int nSymm;
109
110
111
112 private final SpaceGroup spaceGroup;
113
114
115
116 private final int nMolecules;
117
118
119
120 private double kT;
121
122
123
124 private double pressure = 1.0;
125
126
127
128 private boolean isotropic = false;
129
130
131
132 private boolean active = true;
133
134
135
136 private double maxAngleMove = 0.5;
137
138
139
140 private double maxVolumeMove = 1.0;
141
142
143
144 private double minDensity = 0.75;
145
146
147
148 private double maxDensity = 1.60;
149
150
151
152 private int meanBarostatInterval = 1;
153
154
155
156 private int barostatCount = 0;
157
158
159
160 private final RunningStatistics unitCellMoves = new RunningStatistics();
161
162
163
164 private final RunningStatistics sideMoves = new RunningStatistics();
165
166
167
168 private final RunningStatistics angleMoves = new RunningStatistics();
169
170
171
172 private STATE state = STATE.BOTH;
173
174
175
176 private MoveType moveType = MoveType.SIDE;
177
178
179
180 private boolean moveAccepted = false;
181
182
183
184 private double currentDensity = 0;
185
186
187
188 private double a = 0;
189
190
191
192 private double b = 0;
193
194
195
196 private double c = 0;
197
198
199
200 private double alpha = 0;
201
202
203
204 private double beta = 0;
205
206
207
208 private double gamma = 0;
209
210
211
212 private final RunningStatistics aStats = new RunningStatistics();
213
214
215
216 private final RunningStatistics bStats = new RunningStatistics();
217
218
219
220 private final RunningStatistics cStats = new RunningStatistics();
221
222
223
224 private final RunningStatistics alphaStats = new RunningStatistics();
225
226
227
228 private final RunningStatistics betaStats = new RunningStatistics();
229
230
231
232 private final RunningStatistics gammaStats = new RunningStatistics();
233
234
235
236 private final RunningStatistics densityStats = new RunningStatistics();
237
238
239
240 private int printFrequency = 1000;
241
242
243
244
245
246
247
248 public Barostat(MolecularAssembly molecularAssembly, CrystalPotential potential) {
249 this(molecularAssembly, potential, 298.15);
250 }
251
252
253
254
255
256
257
258
259 public Barostat(MolecularAssembly molecularAssembly, CrystalPotential potential,
260 double temperature) {
261
262 this.molecularAssembly = molecularAssembly;
263 this.potential = potential;
264 this.kT = temperature * R;
265
266 crystal = potential.getCrystal();
267 unitCell = crystal.getUnitCell();
268 spaceGroup = unitCell.spaceGroup;
269
270 Atom[] atoms = molecularAssembly.getAtomArray();
271
272 int nAtoms = atoms.length;
273 nSymm = spaceGroup.getNumberOfSymOps();
274 mass = molecularAssembly.getMass();
275 x = new double[3 * nAtoms];
276 nMolecules = molecularAssembly.fractionalCount();
277 }
278
279
280
281
282
283
284 public double density() {
285 return (mass * nSymm / AVOGADRO) * (1.0e24 / unitCell.volume);
286 }
287
288
289
290
291 @Override
292 public boolean destroy() {
293
294 return potential.destroy();
295 }
296
297
298
299
300 @Override
301 public double energy(double[] x) {
302
303 return potential.energy(x);
304 }
305
306
307
308
309 @Override
310 public double energyAndGradient(double[] x, double[] g) {
311
312
313 double energy = potential.energyAndGradient(x, g);
314
315
316 if (active && state != STATE.FAST) {
317 if (random() < (1.0 / meanBarostatInterval)) {
318
319
320 moveAccepted = false;
321
322 applyBarostat(energy);
323
324
325 collectStats();
326
327
328
329 if (moveAccepted) {
330 energy = potential.energyAndGradient(x, g);
331 }
332 }
333 }
334 return energy;
335 }
336
337
338
339
340
341
342
343 public boolean isIsotropic() {
344 return isotropic;
345 }
346
347
348
349
350
351
352
353 public void setIsotropic(boolean isotropic) {
354 this.isotropic = isotropic;
355 }
356
357
358
359
360 @Override
361 public double[] getAcceleration(double[] acceleration) {
362 return potential.getAcceleration(acceleration);
363 }
364
365
366
367
368 @Override
369 public double[] getCoordinates(double[] parameters) {
370 return potential.getCoordinates(parameters);
371 }
372
373
374
375
376 @Override
377 public void setCoordinates(double[] parameters) {
378 potential.setCoordinates(parameters);
379 }
380
381
382
383
384 @Override
385 public Crystal getCrystal() {
386 return potential.getCrystal();
387 }
388
389
390
391
392 @Override
393 public void setCrystal(Crystal crystal) {
394 potential.setCrystal(crystal);
395 }
396
397
398
399
400 @Override
401 public STATE getEnergyTermState() {
402 return state;
403 }
404
405
406
407
408 @Override
409 public void setEnergyTermState(STATE state) {
410 this.state = state;
411 potential.setEnergyTermState(state);
412 }
413
414
415
416
417 @Override
418 public double[] getMass() {
419 return potential.getMass();
420 }
421
422
423
424
425
426
427 public int getMeanBarostatInterval() {
428 return meanBarostatInterval;
429 }
430
431
432
433
434
435
436 public void setMeanBarostatInterval(int meanBarostatInterval) {
437 this.meanBarostatInterval = meanBarostatInterval;
438 }
439
440
441
442
443 @Override
444 public int getNumberOfVariables() {
445 return potential.getNumberOfVariables();
446 }
447
448
449
450
451
452
453 public double getPressure() {
454 return pressure;
455 }
456
457
458
459
460
461
462 public void setPressure(double pressure) {
463 this.pressure = pressure;
464 }
465
466
467
468
469 @Override
470 public double[] getPreviousAcceleration(double[] previousAcceleration) {
471 return potential.getPreviousAcceleration(previousAcceleration);
472 }
473
474
475
476
477 @Override
478 public double[] getScaling() {
479 return potential.getScaling();
480 }
481
482
483
484
485 @Override
486 public void setScaling(double[] scaling) {
487 potential.setScaling(scaling);
488 }
489
490
491
492
493 @Override
494 public double getTotalEnergy() {
495 return potential.getTotalEnergy();
496 }
497
498 @Override
499 public List<Potential> getUnderlyingPotentials() {
500 List<Potential> underlying = new ArrayList<>();
501 underlying.add(potential);
502 underlying.addAll(potential.getUnderlyingPotentials());
503 return underlying;
504 }
505
506
507
508
509
510
511 public CrystalPotential getCrystalPotential() {
512 return potential;
513 }
514
515
516
517
518 @Override
519 public VARIABLE_TYPE[] getVariableTypes() {
520 return potential.getVariableTypes();
521 }
522
523
524
525
526 @Override
527 public double[] getVelocity(double[] velocity) {
528 return potential.getVelocity(velocity);
529 }
530
531 public boolean isActive() {
532 return active;
533 }
534
535
536
537
538
539
540 public void setActive(boolean active) {
541 this.active = active;
542 }
543
544
545
546
547 @Override
548 public void setAcceleration(double[] acceleration) {
549 potential.setAcceleration(acceleration);
550 }
551
552
553
554
555
556
557 public void setDensity(double density) {
558 molecularAssembly.computeFractionalCoordinates();
559 crystal.setDensity(density, mass);
560 potential.setCrystal(crystal);
561 molecularAssembly.moveToFractionalCoordinates();
562 }
563
564
565
566
567
568
569 public void setMaxAngleMove(double maxAngleMove) {
570 this.maxAngleMove = maxAngleMove;
571 }
572
573
574
575
576
577
578 public void setMaxVolumeMove(double maxVolumeMove) {
579 this.maxVolumeMove = maxVolumeMove;
580 }
581
582
583
584
585
586
587 public void setMaxDensity(double maxDensity) {
588 this.maxDensity = maxDensity;
589 }
590
591
592
593
594
595
596 public void setMinDensity(double minDensity) {
597 this.minDensity = minDensity;
598 }
599
600
601
602
603
604
605 public void setBarostatPrintFrequency(int frequency) {
606 this.printFrequency = frequency;
607 }
608
609
610
611
612 @Override
613 public void setPreviousAcceleration(double[] previousAcceleration) {
614 potential.setPreviousAcceleration(previousAcceleration);
615 }
616
617
618
619
620
621
622 public void setTemperature(double temperature) {
623 this.kT = temperature * R;
624 }
625
626
627
628
629 @Override
630 public void setVelocity(double[] velocity) {
631 potential.setVelocity(velocity);
632 }
633
634 private double mcStep(double currentE, double currentV) {
635
636
637 double den = density();
638 if(Double.isNaN(den) || Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(c) || Double.isNaN(alpha) || Double.isNaN(beta) || Double.isNaN(gamma)){
639 logger.warning(format(" Found value of NaN: density: %7.4f a: %7.4f b: %7.4f c: %7.4f alpha: %7.4f beta: %7.4f gamma: %7.4f",
640 den, a, b, c, alpha, beta, gamma));
641 File errorFile = new File(FilenameUtils.removeExtension(molecularAssembly.getFile().getName()) + "_err.xyz");
642 XYZFilter.version(errorFile);
643 XYZFilter writeFilter = new XYZFilter(errorFile, molecularAssembly, molecularAssembly.getForceField(), molecularAssembly.getProperties());
644 writeFilter.writeFile(errorFile, true, null);
645 }
646 if (den < minDensity || den > maxDensity) {
647 if (logger.isLoggable(Level.FINE)) {
648 logger.fine(
649 format(" MC Density %10.6f is outside the range %10.6f - %10.6f.", den, minDensity,
650 maxDensity));
651 }
652
653 crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
654 return currentE;
655 }
656
657
658
659 double minInterfacialRadius = 1.2;
660 double currentMinInterfacialRadius = min(min(unitCell.interfacialRadiusA, unitCell.interfacialRadiusB),
661 unitCell.interfacialRadiusC);
662
663 if (currentMinInterfacialRadius < minInterfacialRadius) {
664 if (logger.isLoggable(Level.FINE)) {
665 logger.fine(format(" MC An interfacial radius (%10.6f,%10.6f,%10.6f) is below the minimum %10.6f",
666 unitCell.interfacialRadiusA, unitCell.interfacialRadiusB,
667 unitCell.interfacialRadiusC, minInterfacialRadius));
668 }
669
670 crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
671 return currentE;
672 }
673
674
675
676 potential.setCrystal(crystal);
677
678
679 molecularAssembly.moveToFractionalCoordinates();
680
681
682 double newV = unitCell.volume / nSymm;
683
684 potential.getCoordinates(x);
685
686
687 double newE = potential.energy(x);
688
689
690 double dPE = newE - currentE;
691
692
693 double dEV = pressure * (newV - currentV) / PRESCON;
694
695
696 double dES = -nMolecules * kT * log(newV / currentV);
697
698
699 double dE = dPE + dEV + dES;
700
701
702 if (dE < 0.0) {
703 if (logger.isLoggable(Level.FINE)) {
704 logger.fine(format(" MC Accept: %12.6f (dPE) + %12.6f (dEV) + %12.6f (dES) = %12.6f with (V=%12.6f, E=%12.6f)",
705 dPE, dEV, dES, dE, newV, newE));
706 }
707 moveAccepted = true;
708 switch (moveType) {
709 case SIDE -> sideMoves.addValue(1.0);
710 case ANGLE -> angleMoves.addValue(1.0);
711 case UNIT -> unitCellMoves.addValue(1.0);
712 }
713 return newE;
714 }
715
716
717 double acceptanceProbability = exp(-dE / kT);
718
719
720
721 double metropolis = random();
722
723
724 if (metropolis > acceptanceProbability) {
725 rejectMove();
726 if (logger.isLoggable(Level.FINE)) {
727 logger.fine(format(" MC Reject: %12.6f (dPE) + %12.6f (dEV) + %12.6f (dES) = %12.6f with (V=%12.6f, E=%12.6f)",
728 dPE, dEV, dES, dE, currentV, currentE));
729 }
730 switch (moveType) {
731 case SIDE -> sideMoves.addValue(0.0);
732 case ANGLE -> angleMoves.addValue(0.0);
733 case UNIT -> unitCellMoves.addValue(0.0);
734 }
735 return currentE;
736 }
737
738
739 if (logger.isLoggable(Level.FINE)) {
740 logger.fine(format(" MC Accept: %12.6f (dPE) + %12.6f (dEV) + %12.6f (dES) = %12.6f with (V=%12.6f, E=%12.6f)",
741 dPE, dEV, dES, dE, newV, newE));
742 }
743
744 moveAccepted = true;
745 switch (moveType) {
746 case SIDE -> sideMoves.addValue(1.0);
747 case ANGLE -> angleMoves.addValue(1.0);
748 case UNIT -> unitCellMoves.addValue(1.0);
749 }
750
751 return newE;
752 }
753
754
755
756
757 private void rejectMove() {
758
759 crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
760
761
762 potential.setCrystal(crystal);
763
764
765 molecularAssembly.moveToFractionalCoordinates();
766 }
767
768
769
770
771
772
773
774 private double mcA(double currentE) {
775 moveType = MoveType.SIDE;
776 double currentAUV = unitCell.volume / nSymm;
777 double dAUVolume = maxVolumeMove * (2.0 * random() - 1.0);
778 double dVdA = unitCell.dVdA / nSymm;
779 double dA = dAUVolume / dVdA;
780 boolean succeed = crystal.changeUnitCellParameters(a + dA, b, c, alpha, beta, gamma);
781 if (succeed) {
782 if (logger.isLoggable(Level.FINE)) {
783 logger.fine(format(" Proposing MC change of the a-axis to (%6.3f) with volume change %6.3f (A^3)", a, dAUVolume));
784 }
785 return mcStep(currentE, currentAUV);
786 }
787 return currentE;
788 }
789
790
791
792
793
794
795
796 private double mcB(double currentE) {
797 moveType = MoveType.SIDE;
798 double currentAUV = unitCell.volume / nSymm;
799 double dAUVolume = maxVolumeMove * (2.0 * random() - 1.0);
800 double dVdB = unitCell.dVdB / nSymm;
801 double dB = dAUVolume / dVdB;
802 boolean succeed = crystal.changeUnitCellParameters(a, b + dB, c, alpha, beta, gamma);
803 if (succeed) {
804 if (logger.isLoggable(Level.FINE)) {
805 logger.fine(format(" Proposing MC change of the b-axis to (%6.3f) with volume change %6.3f (A^3)", b, dAUVolume));
806 }
807 return mcStep(currentE, currentAUV);
808 }
809 return currentE;
810 }
811
812
813
814
815
816
817
818 private double mcC(double currentE) {
819 moveType = MoveType.SIDE;
820 double currentAUV = unitCell.volume / nSymm;
821 double dAUVolume = maxVolumeMove * (2.0 * random() - 1.0);
822 double dVdC = unitCell.dVdC / nSymm;
823 double dC = dAUVolume / dVdC;
824 boolean succeed = crystal.changeUnitCellParameters(a, b, c + dC, alpha, beta, gamma);
825 if (succeed) {
826 if (logger.isLoggable(Level.FINE)) {
827 logger.fine(format(" Proposing MC change to the c-axis to (%6.3f) with volume change %6.3f (A^3)", c, dAUVolume));
828 }
829 return mcStep(currentE, currentAUV);
830 }
831 return currentE;
832 }
833
834
835
836
837
838
839
840 private double mcAB(double currentE) {
841 moveType = MoveType.SIDE;
842 double currentAUV = unitCell.volume / nSymm;
843 double dAUVolume = maxVolumeMove * (2.0 * random() - 1.0);
844 double dVdAB = (unitCell.dVdA + unitCell.dVdB) / nSymm;
845 double dAB = dAUVolume / dVdAB;
846 boolean succeed = crystal.changeUnitCellParametersAndVolume(a + dAB, b + dAB, c, alpha, beta,
847 gamma, currentAUV + dAUVolume);
848 if (succeed) {
849 if (logger.isLoggable(Level.FINE)) {
850 logger.fine(format(" Proposing MC change to the a,b-axes to (%6.3f) with volume change %6.3f (A^3)", a, dAUVolume));
851 }
852 return mcStep(currentE, currentAUV);
853 }
854 return currentE;
855 }
856
857 private double mcABC(double currentE) {
858 moveType = MoveType.SIDE;
859 double currentAUV = unitCell.volume / nSymm;
860 double dAUVolume = maxVolumeMove * (2.0 * random() - 1.0);
861 double dVdABC = (unitCell.dVdA + unitCell.dVdB + unitCell.dVdC) / nSymm;
862 double dABC = dAUVolume / dVdABC;
863 boolean succeed = crystal.changeUnitCellParametersAndVolume(a + dABC, b + dABC, c + dABC, alpha,
864 beta, gamma, currentAUV + dAUVolume);
865 if (succeed) {
866 if (logger.isLoggable(Level.FINE)) {
867 logger.fine(format(" Proposing MC change to the a,b,c-axes to (%6.3f) with volume change %6.3f (A^3)", a, dAUVolume));
868 }
869 return mcStep(currentE, currentAUV);
870 }
871 return currentE;
872 }
873
874
875
876
877
878
879
880 private double mcAlpha(double currentE) {
881 moveType = MoveType.ANGLE;
882 double move = maxAngleMove * (2.0 * random() - 1.0);
883 double currentAUV = unitCell.volume / nSymm;
884 double newAlpha = mirrorDegrees(alpha + move);
885 boolean succeed = crystal.changeUnitCellParametersAndVolume(a, b, c, newAlpha, beta, gamma, currentAUV);
886 if (succeed) {
887 if (logger.isLoggable(Level.FINE)) {
888 logger.fine(format(" Proposing MC change of alpha to %6.3f.", newAlpha));
889 }
890 return mcStep(currentE, currentAUV);
891 }
892 return currentE;
893 }
894
895
896
897
898
899
900
901 private double mcBeta(double currentE) {
902 moveType = MoveType.ANGLE;
903 double move = maxAngleMove * (2.0 * random() - 1.0);
904 double currentAUV = unitCell.volume / nSymm;
905 double newBeta = mirrorDegrees(beta + move);
906 boolean succeed = crystal.changeUnitCellParametersAndVolume(a, b, c, alpha, newBeta, gamma, currentAUV);
907 if (succeed) {
908 if (logger.isLoggable(Level.FINE)) {
909 logger.fine(format(" Proposing MC change of beta to %6.3f.", newBeta));
910 }
911 return mcStep(currentE, currentAUV);
912 }
913 return currentE;
914 }
915
916
917
918
919
920
921
922 private double mcGamma(double currentE) {
923 moveType = MoveType.ANGLE;
924 double move = maxAngleMove * (2.0 * random() - 1.0);
925 double currentAUV = unitCell.volume / nSymm;
926 double newGamma = mirrorDegrees(gamma + move);
927 boolean succeed = crystal.changeUnitCellParametersAndVolume(a, b, c, alpha, beta, newGamma, currentAUV);
928 if (succeed) {
929 if (logger.isLoggable(Level.FINE)) {
930 logger.fine(format(" Proposing MC change of gamma to %6.3f.", newGamma));
931 }
932 return mcStep(currentE, currentAUV);
933 }
934 return currentE;
935 }
936
937
938
939
940
941
942
943 private double mcABG(double currentE) {
944 moveType = MoveType.ANGLE;
945 double move = maxAngleMove * (2.0 * random() - 1.0);
946 double currentAUV = unitCell.volume / nSymm;
947 double newAlpha = mirrorDegrees(alpha + move);
948 double newBeta = mirrorDegrees(beta + move);
949 double newGamma = mirrorDegrees(gamma + move);
950 boolean succeed = crystal.changeUnitCellParametersAndVolume(a, b, c, newAlpha, newBeta, newGamma, currentAUV);
951 if (succeed) {
952 if (logger.isLoggable(Level.FINE)) {
953 logger.fine(format(" Proposing MC change to all angles to %6.3f.", newAlpha));
954 }
955 return mcStep(currentE, currentAUV);
956 }
957 return currentE;
958 }
959
960
961
962
963
964
965
966 private double applyBarostat(double currentE) {
967
968 molecularAssembly.computeFractionalCoordinates();
969
970
971 crystal = potential.getCrystal();
972 unitCell = crystal.getUnitCell();
973 a = unitCell.a;
974 b = unitCell.b;
975 c = unitCell.c;
976 alpha = unitCell.alpha;
977 beta = unitCell.beta;
978 gamma = unitCell.gamma;
979
980 if (isotropic) {
981 currentE = mcABC(currentE);
982 } else {
983 switch (spaceGroup.latticeSystem) {
984 case MONOCLINIC_LATTICE -> {
985
986 int move = (int) floor(random() * 4.0);
987 switch (move) {
988 case 0 -> currentE = mcA(currentE);
989 case 1 -> currentE = mcB(currentE);
990 case 2 -> currentE = mcC(currentE);
991 case 3 -> currentE = mcBeta(currentE);
992 default -> logger.severe(" Barostat programming error.");
993 }
994 }
995 case ORTHORHOMBIC_LATTICE -> {
996
997 int move = (int) floor(random() * 3.0);
998 switch (move) {
999 case 0 -> currentE = mcA(currentE);
1000 case 1 -> currentE = mcB(currentE);
1001 case 2 -> currentE = mcC(currentE);
1002 default -> logger.severe(" Barostat programming error.");
1003 }
1004 }
1005 case TETRAGONAL_LATTICE -> {
1006
1007 int move = (int) floor(random() * 2.0);
1008 switch (move) {
1009 case 0 -> currentE = mcAB(currentE);
1010 case 1 -> currentE = mcC(currentE);
1011 default -> logger.severe(" Barostat programming error.");
1012 }
1013 }
1014 case RHOMBOHEDRAL_LATTICE -> {
1015
1016 int move = (int) floor(random() * 2.0);
1017 switch (move) {
1018 case 0 -> currentE = mcABC(currentE);
1019 case 1 -> currentE = mcABG(currentE);
1020 default -> logger.severe(" Barostat programming error.");
1021 }
1022 }
1023 case HEXAGONAL_LATTICE -> {
1024
1025 int move = (int) floor(random() * 2.0);
1026 switch (move) {
1027 case 0 -> currentE = mcAB(currentE);
1028 case 1 -> currentE = mcC(currentE);
1029 default -> logger.severe(" Barostat programming error.");
1030 }
1031 }
1032 case CUBIC_LATTICE ->
1033
1034 currentE = mcABC(currentE);
1035 case TRICLINIC_LATTICE -> {
1036 if (check(a, b) && check(b, c) && check(alpha, 90.0) && check(beta, 90.0) && check(gamma, 90.0)) {
1037 currentE = mcABC(currentE);
1038 } else {
1039 int move = (int) floor(random() * 6.0);
1040 switch (move) {
1041 case 0 -> currentE = mcA(currentE);
1042 case 1 -> currentE = mcB(currentE);
1043 case 2 -> currentE = mcC(currentE);
1044 case 3 -> currentE = mcAlpha(currentE);
1045 case 4 -> currentE = mcBeta(currentE);
1046 case 5 -> currentE = mcGamma(currentE);
1047 default -> logger.severe(" Barostat programming error.");
1048 }
1049 }
1050 }
1051 }
1052 }
1053
1054 currentDensity = density();
1055 if (moveAccepted) {
1056 if (logger.isLoggable(Level.FINE)) {
1057 StringBuilder sb = new StringBuilder(" MC Barostat Acceptance:");
1058 if (sideMoves.getCount() > 0) {
1059 sb.append(format(" Side %5.1f%%", sideMoves.getMean() * 100.0));
1060 }
1061 if (angleMoves.getCount() > 0) {
1062 sb.append(format(" Angle %5.1f%%", angleMoves.getMean() * 100.0));
1063 }
1064 if (unitCellMoves.getCount() > 0) {
1065 sb.append(format(" UC %5.1f%%", unitCellMoves.getMean() * 100.0));
1066 }
1067 sb.append(format("\n Density: %5.3f UC: %s", currentDensity, unitCell.toShortString()));
1068 logger.fine(sb.toString());
1069 }
1070 } else {
1071
1072 if (unitCell.a != a || unitCell.b != b || unitCell.c != c || unitCell.alpha != alpha
1073 || unitCell.beta != beta || unitCell.gamma != gamma) {
1074 logger.severe(" Reversion of unit cell parameters did not succeed after failed Barostat MC move.");
1075 }
1076 }
1077
1078 return currentE;
1079 }
1080
1081
1082
1083
1084 private void collectStats() {
1085
1086 barostatCount++;
1087
1088
1089 if(Double.isNaN(currentDensity)||Double.isNaN(a)||Double.isNaN(b)||Double.isNaN(c)||Double.isNaN(alpha)||Double.isNaN(beta)||Double.isNaN(gamma)){
1090 logger.warning(format(" Statistic Value was NaN: Density: %5.3f A: %5.3f B: %5.3f C: %5.3f Alpha: %5.3f Beta: %5.3f Gamma: %5.3f",
1091 currentDensity, a, b, c, alpha, beta, gamma));
1092 }
1093 densityStats.addValue(currentDensity);
1094 aStats.addValue(a);
1095 bStats.addValue(b);
1096 cStats.addValue(c);
1097 alphaStats.addValue(alpha);
1098 betaStats.addValue(beta);
1099 gammaStats.addValue(gamma);
1100 if (barostatCount % printFrequency == 0) {
1101 logger.info(format(" Barostat statistics for the last %d moves:", printFrequency));
1102
1103 logger.info(format(" Density: %5.3f±%5.3f with range %5.3f .. %5.3f (g/cc)", densityStats.getMean(),
1104 densityStats.getStandardDeviation(), densityStats.getMin(), densityStats.getMax()));
1105 densityStats.reset();
1106
1107 logger.info(format(" Lattice a-axis: %5.2f±%3.2f b-axis: %5.2f±%3.2f c-axis: %5.2f±%3.2f",
1108 aStats.getMean(), aStats.getStandardDeviation(), bStats.getMean(), bStats.getStandardDeviation(), cStats.getMean(),
1109 cStats.getStandardDeviation()));
1110 logger.info(format(" alpha: %5.2f±%3.2f beta: %5.2f±%3.2f gamma: %5.2f±%3.2f",
1111 alphaStats.getMean(), alphaStats.getStandardDeviation(), betaStats.getMean(),
1112 betaStats.getStandardDeviation(), gammaStats.getMean(), gammaStats.getStandardDeviation()));
1113 aStats.reset();
1114 bStats.reset();
1115 cStats.reset();
1116 alphaStats.reset();
1117 betaStats.reset();
1118 gammaStats.reset();
1119
1120 if (sideMoves.getCount() > 0) {
1121 logger.info(format(" Axis length moves: %4d (%5.2f%%)", sideMoves.getCount(), sideMoves.getMean() * 100.0));
1122 sideMoves.reset();
1123 }
1124 if (angleMoves.getCount() > 0) {
1125 logger.info(format(" Angle moves: %4d (%5.2f%%)", angleMoves.getCount(), angleMoves.getMean() * 100.0));
1126 angleMoves.reset();
1127 }
1128 if (unitCellMoves.getCount() > 0) {
1129 logger.info(format(" Unit cell moves: %4d (%5.2f%%)", unitCellMoves.getCount(), unitCellMoves.getMean() * 100.0));
1130 unitCellMoves.reset();
1131 }
1132 }
1133 }
1134
1135
1136
1137
1138 private enum MoveType {
1139 SIDE, ANGLE, UNIT
1140 }
1141 }