View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * The Barostat class maintains constant pressure using random trial moves in lattice parameters,
70   * which are consistent with the space group.
71   *
72   * <p>Frenkel and Smit, "Understanding Molecular Simulation, 2nd Edition", Academic Press, San
73   * Diego, CA, 2002; Section 5.4
74   *
75   * @author Michael J. Schnieders
76   */
77  public class Barostat implements CrystalPotential {
78  
79    private static final Logger logger = Logger.getLogger(Barostat.class.getName());
80  
81    /**
82     * MolecularAssembly being simulated.
83     */
84    private final MolecularAssembly molecularAssembly;
85    /**
86     * Boundary conditions and symmetry operators (could be a ReplicatedCrystal).
87     */
88    private Crystal crystal;
89    /**
90     * The unit cell.
91     */
92    private Crystal unitCell;
93    /**
94     * Mass of the system.
95     */
96    private final double mass;
97    /**
98     * ForceFieldEnergy that describes the system.
99     */
100   private final CrystalPotential potential;
101   /**
102    * Atomic coordinates.
103    */
104   private final double[] x;
105   /**
106    * The number of space group symmetry operators.
107    */
108   private final int nSymm;
109   /**
110    * Number of independent molecules in the simulation cell.
111    */
112   private final SpaceGroup spaceGroup;
113   /**
114    * Number of molecules in the systems.
115    */
116   private final int nMolecules;
117   /**
118    * Ideal gas constant * temperature (kcal/mol).
119    */
120   private double kT;
121   /**
122    * Sampling pressure (atm)
123    */
124   private double pressure = 1.0;
125   /**
126    * Only isotropic MC moves.
127    */
128   private boolean isotropic = false;
129   /**
130    * Flag to turn the Barostat on or off. If false, MC moves will not be tried.
131    */
132   private boolean active = true;
133   /**
134    * Default angle move (degrees).
135    */
136   private double maxAngleMove = 0.5;
137   /**
138    * Maximum volume move (Angstroms^3).
139    */
140   private double maxVolumeMove = 1.0;
141   /**
142    * Minimum density constraint.
143    */
144   private double minDensity = 0.75;
145   /**
146    * Maximum density constraint.
147    */
148   private double maxDensity = 1.60;
149   /**
150    * Number of energy evaluations between application of MC moves.
151    */
152   private int meanBarostatInterval = 1;
153   /**
154    * A counter for the number of barostat calls.
155    */
156   private int barostatCount = 0;
157   /**
158    * Number of unit cell Monte Carlo moves attempted.
159    */
160   private final RunningStatistics unitCellMoves = new RunningStatistics();
161   /**
162    * Number of lattice axis Monte Carlo moves attempted.
163    */
164   private final RunningStatistics sideMoves = new RunningStatistics();
165   /**
166    * Number of Monte Carlo moves attempted.
167    */
168   private final RunningStatistics angleMoves = new RunningStatistics();
169   /**
170    * Energy STATE.
171    */
172   private STATE state = STATE.BOTH;
173   /**
174    * Barostat move type.
175    */
176   private MoveType moveType = MoveType.SIDE;
177   /**
178    * True when a Monte Carlo move is accepted.
179    */
180   private boolean moveAccepted = false;
181   /**
182    * Current density value.
183    */
184   private double currentDensity = 0;
185   /**
186    * Current A-axis value.
187    */
188   private double a = 0;
189   /**
190    * Current B-axis value.
191    */
192   private double b = 0;
193   /**
194    * Current C-axis value.
195    */
196   private double c = 0;
197   /**
198    * Current alpha value.
199    */
200   private double alpha = 0;
201   /**
202    * Current beta value.
203    */
204   private double beta = 0;
205   /**
206    * Current gamma value.
207    */
208   private double gamma = 0;
209   /**
210    * A-axis statistics.
211    */
212   private final RunningStatistics aStats = new RunningStatistics();
213   /**
214    * B-axis statistics.
215    */
216   private final RunningStatistics bStats = new RunningStatistics();
217   /**
218    * C-axis statistics.
219    */
220   private final RunningStatistics cStats = new RunningStatistics();
221   /**
222    * Alpha statistics.
223    */
224   private final RunningStatistics alphaStats = new RunningStatistics();
225   /**
226    * Beta statistics.
227    */
228   private final RunningStatistics betaStats = new RunningStatistics();
229   /**
230    * Gamma statistics.
231    */
232   private final RunningStatistics gammaStats = new RunningStatistics();
233   /**
234    * Density statistics.
235    */
236   private final RunningStatistics densityStats = new RunningStatistics();
237   /**
238    * Barostat print frequency.
239    */
240   private int printFrequency = 1000;
241 
242   /**
243    * Initialize the Barostat.
244    *
245    * @param molecularAssembly The molecular assembly to apply the MC barostat to.
246    * @param potential         a {@link ffx.crystal.CrystalPotential} object.
247    */
248   public Barostat(MolecularAssembly molecularAssembly, CrystalPotential potential) {
249     this(molecularAssembly, potential, 298.15);
250   }
251 
252   /**
253    * Initialize the Barostat.
254    *
255    * @param molecularAssembly The molecular assembly to apply the MC barostat to.
256    * @param potential         a {@link ffx.crystal.CrystalPotential} object.
257    * @param temperature       The Metropolis Monte Carlo temperature (Kelvin).
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     // Atoms in the system.
270     Atom[] atoms = molecularAssembly.getAtomArray();
271     // Number of atoms.
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    * density.
281    *
282    * @return a double.
283    */
284   public double density() {
285     return (mass * nSymm / AVOGADRO) * (1.0e24 / unitCell.volume);
286   }
287 
288   /**
289    * {@inheritDoc}
290    */
291   @Override
292   public boolean destroy() {
293     // Nothing at this level to destroy.
294     return potential.destroy();
295   }
296 
297   /**
298    * {@inheritDoc}
299    */
300   @Override
301   public double energy(double[] x) {
302     // Do not apply the Barostat for energy only evaluations.
303     return potential.energy(x);
304   }
305 
306   /**
307    * {@inheritDoc}
308    */
309   @Override
310   public double energyAndGradient(double[] x, double[] g) {
311 
312     // Calculate the energy and gradient as usual.
313     double energy = potential.energyAndGradient(x, g);
314 
315     // Apply the barostat during computation of slowly varying forces.
316     if (active && state != STATE.FAST) {
317       if (random() < (1.0 / meanBarostatInterval)) {
318 
319         // Attempt to change the unit cell parameters.
320         moveAccepted = false;
321 
322         applyBarostat(energy);
323 
324         // Collect Statistics.
325         collectStats();
326 
327         // If a move was accepted, then re-calculate the gradient so
328         // that it's consistent with the current unit cell parameters.
329         if (moveAccepted) {
330           energy = potential.energyAndGradient(x, g);
331         }
332       }
333     }
334     return energy;
335   }
336 
337   /**
338    * Restrict the MC Barostat to isotropic moves. The lattice angles are held fixed, and lattice
339    * lengths are scaled equally.
340    *
341    * @return Returns true if only isotropic moves are allowed.
342    */
343   public boolean isIsotropic() {
344     return isotropic;
345   }
346 
347   /**
348    * Restrict the MC Barostat to isotropic moves. The lattice angles are held fixed, and lattice
349    * lengths are scaled equally.
350    *
351    * @param isotropic If true, if only isotropic moves are allowed.
352    */
353   public void setIsotropic(boolean isotropic) {
354     this.isotropic = isotropic;
355   }
356 
357   /**
358    * {@inheritDoc}
359    */
360   @Override
361   public double[] getAcceleration(double[] acceleration) {
362     return potential.getAcceleration(acceleration);
363   }
364 
365   /**
366    * {@inheritDoc}
367    */
368   @Override
369   public double[] getCoordinates(double[] parameters) {
370     return potential.getCoordinates(parameters);
371   }
372 
373   /**
374    * {@inheritDoc}
375    */
376   @Override
377   public void setCoordinates(double[] parameters) {
378     potential.setCoordinates(parameters);
379   }
380 
381   /**
382    * {@inheritDoc}
383    */
384   @Override
385   public Crystal getCrystal() {
386     return potential.getCrystal();
387   }
388 
389   /**
390    * {@inheritDoc}
391    */
392   @Override
393   public void setCrystal(Crystal crystal) {
394     potential.setCrystal(crystal);
395   }
396 
397   /**
398    * {@inheritDoc}
399    */
400   @Override
401   public STATE getEnergyTermState() {
402     return state;
403   }
404 
405   /**
406    * {@inheritDoc}
407    */
408   @Override
409   public void setEnergyTermState(STATE state) {
410     this.state = state;
411     potential.setEnergyTermState(state);
412   }
413 
414   /**
415    * {@inheritDoc}
416    */
417   @Override
418   public double[] getMass() {
419     return potential.getMass();
420   }
421 
422   /**
423    * Returns the mean number of steps between barostat applications.
424    *
425    * @return Mean steps between barostat applications.
426    */
427   public int getMeanBarostatInterval() {
428     return meanBarostatInterval;
429   }
430 
431   /**
432    * Setter for the field <code>meanBarostatInterval</code>.
433    *
434    * @param meanBarostatInterval The mean number of steps between barostat applications.
435    */
436   public void setMeanBarostatInterval(int meanBarostatInterval) {
437     this.meanBarostatInterval = meanBarostatInterval;
438   }
439 
440   /**
441    * {@inheritDoc}
442    */
443   @Override
444   public int getNumberOfVariables() {
445     return potential.getNumberOfVariables();
446   }
447 
448   /**
449    * Gets the pressure of this Barostat in atm.
450    *
451    * @return Pressure in atm.
452    */
453   public double getPressure() {
454     return pressure;
455   }
456 
457   /**
458    * Setter for the field <code>pressure</code>.
459    *
460    * @param pressure a double.
461    */
462   public void setPressure(double pressure) {
463     this.pressure = pressure;
464   }
465 
466   /**
467    * {@inheritDoc}
468    */
469   @Override
470   public double[] getPreviousAcceleration(double[] previousAcceleration) {
471     return potential.getPreviousAcceleration(previousAcceleration);
472   }
473 
474   /**
475    * {@inheritDoc}
476    */
477   @Override
478   public double[] getScaling() {
479     return potential.getScaling();
480   }
481 
482   /**
483    * {@inheritDoc}
484    */
485   @Override
486   public void setScaling(double[] scaling) {
487     potential.setScaling(scaling);
488   }
489 
490   /**
491    * {@inheritDoc}
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    * Get the CrystalPotential that this Barostat is applying to.
508    *
509    * @return The CrystalPotential.
510    */
511   public CrystalPotential getCrystalPotential() {
512     return potential;
513   }
514 
515   /**
516    * {@inheritDoc}
517    */
518   @Override
519   public VARIABLE_TYPE[] getVariableTypes() {
520     return potential.getVariableTypes();
521   }
522 
523   /**
524    * {@inheritDoc}
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    * Setter for the field <code>active</code>.
537    *
538    * @param active a boolean.
539    */
540   public void setActive(boolean active) {
541     this.active = active;
542   }
543 
544   /**
545    * {@inheritDoc}
546    */
547   @Override
548   public void setAcceleration(double[] acceleration) {
549     potential.setAcceleration(acceleration);
550   }
551 
552   /**
553    * setDensity.
554    *
555    * @param density a double.
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    * Setter for the field <code>maxAngleMove</code>.
566    *
567    * @param maxAngleMove a double.
568    */
569   public void setMaxAngleMove(double maxAngleMove) {
570     this.maxAngleMove = maxAngleMove;
571   }
572 
573   /**
574    * Setter for the field <code>maxVolumeMove</code>.
575    *
576    * @param maxVolumeMove a double.
577    */
578   public void setMaxVolumeMove(double maxVolumeMove) {
579     this.maxVolumeMove = maxVolumeMove;
580   }
581 
582   /**
583    * Setter for the field <code>maxDensity</code>.
584    *
585    * @param maxDensity a double.
586    */
587   public void setMaxDensity(double maxDensity) {
588     this.maxDensity = maxDensity;
589   }
590 
591   /**
592    * Setter for the field <code>minDensity</code>.
593    *
594    * @param minDensity a double.
595    */
596   public void setMinDensity(double minDensity) {
597     this.minDensity = minDensity;
598   }
599 
600   /**
601    * Set the Barostat print frequency.
602    *
603    * @param frequency The number of Barostat moves between print statements.
604    */
605   public void setBarostatPrintFrequency(int frequency) {
606     this.printFrequency = frequency;
607   }
608 
609   /**
610    * {@inheritDoc}
611    */
612   @Override
613   public void setPreviousAcceleration(double[] previousAcceleration) {
614     potential.setPreviousAcceleration(previousAcceleration);
615   }
616 
617   /**
618    * Set the Metropolis Monte Carlo temperature.
619    *
620    * @param temperature Temperature (Kelvin).
621    */
622   public void setTemperature(double temperature) {
623     this.kT = temperature * R;
624   }
625 
626   /**
627    * {@inheritDoc}
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     // Enforce minimum & maximum density constraints.
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       // Reject moves outside the specified density range.
653       crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
654       return currentE;
655     }
656 
657     // A carbon atom cannot fit into a unit cell without interfacial radii greater than ~1.2
658     // Angstroms.
659     double minInterfacialRadius = 1.2;
660     double currentMinInterfacialRadius = min(min(unitCell.interfacialRadiusA, unitCell.interfacialRadiusB),
661         unitCell.interfacialRadiusC);
662     // Enforce minimum interfacial radii of 1.2 Angstroms.
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       // Reject due to small interfacial radius.
670       crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
671       return currentE;
672     }
673 
674     // Apply the boundary condition for the proposed move. If the move is
675     // rejected, then the previous boundary conditions must be restored.
676     potential.setCrystal(crystal);
677 
678     // Update atomic coordinates to maintain molecular fractional centers of mass.
679     molecularAssembly.moveToFractionalCoordinates();
680 
681     // Save the new volume
682     double newV = unitCell.volume / nSymm;
683 
684     potential.getCoordinates(x);
685 
686     // Compute the new energy
687     double newE = potential.energy(x);
688 
689     // Compute the change in potential energy
690     double dPE = newE - currentE;
691 
692     // Compute the pressure-volume work for the asymmetric unit.
693     double dEV = pressure * (newV - currentV) / PRESCON;
694 
695     // Compute the volume entropy
696     double dES = -nMolecules * kT * log(newV / currentV);
697 
698     // Add up the contributions
699     double dE = dPE + dEV + dES;
700 
701     // Energy decreased so the move is accepted.
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     // Apply the Metropolis criteria.
717     double acceptanceProbability = exp(-dE / kT);
718 
719     // Draw a pseudorandom double greater than or equal to 0.0 and less than 1.0
720     // from an (approximately) uniform distribution from that range.
721     double metropolis = random();
722 
723     // Energy increase without the Metropolis criteria satisfied.
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     // Energy increase with Metropolis criteria satisfied.
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    * Reset the to state prior to trial move.
756    */
757   private void rejectMove() {
758     // Reset the unit cell parameters
759     crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
760 
761     // Reset the potential PBC.
762     potential.setCrystal(crystal);
763 
764     // Reset the atomic coordinates to maintain molecular fractional centers of mass.
765     molecularAssembly.moveToFractionalCoordinates();
766   }
767 
768   /**
769    * Propose to change the A-axis length.
770    *
771    * @param currentE Current energy.
772    * @return Energy after the MC trial.
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    * Propose to change the B-axis length.
792    *
793    * @param currentE Current energy.
794    * @return Energy after the MC trial.
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    * Propose to change the C-axis length.
814    *
815    * @param currentE Current energy.
816    * @return Energy after the MC trial.
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    * Propose the same change to the A-axis and B-axis lengths.
836    *
837    * @param currentE Current energy.
838    * @return Energy after the MC trial.
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    * Propose to change the Alpha angle.
876    *
877    * @param currentE Current energy.
878    * @return Energy after the MC trial.
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    * Propose to change the Beta angle.
897    *
898    * @param currentE Current energy.
899    * @return Energy after the MC trial.
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    * Propose to change the Gamma angle.
918    *
919    * @param currentE Current energy.
920    * @return Energy after the MC trial.
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    * Propose the same change to all 3 angles.
939    *
940    * @param currentE Current energy.
941    * @return Energy after the MC trial.
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    * Attempt an MC move on a lattice parameter.
962    *
963    * @param currentE The current potential energy.
964    * @return The potential energy after attempting the MC move.
965    */
966   private double applyBarostat(double currentE) {
967     // Determine the current molecular centers of mass in fractional coordinates.
968     molecularAssembly.computeFractionalCoordinates();
969 
970     // Collect the current unit cell parameters.
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           // alpha = gamma = 90
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           // alpha = beta = gamma = 90
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           // a = b, alpha = beta = gamma = 90
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           // a = b = c, alpha = beta = gamma
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           // a = b, alpha = beta = 90, gamma = 120
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             // a = b = c, alpha = beta = gamma = 90
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       // Check that the unit cell parameters have not changed.
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    * Collect statistics on lattice parameters and density.
1083    */
1084   private void collectStats() {
1085     // Collect statistics.
1086     barostatCount++;
1087 
1088     // Sanity check for values.
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       // Density
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       // Lattice Parameters
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       // MC Move Statistics
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    * The type of Barostat move.
1137    */
1138   private enum MoveType {
1139     SIDE, ANGLE, UNIT
1140   }
1141 }