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 }