add csvr barostat

ipsuite/calculators/__init__.pyi

@@ -8,6 +8,7 @@ from .ase_md import (
  NPTThermostat,
  PressureRampModifier,
  RescaleBoxModifier,
+ SVCRBarostat,
  TemperatureOscillatingRampModifier,
  TemperatureRampModifier,
  VelocityVerletDynamic,
@@ -40,4 +41,5 @@ __all__ = [
  "LammpsSimulator",
  "FixedLayerConstraint",
  "MixCalculator",
+ "SVCRBarostat",
 ]

ipsuite/calculators/ase_md.py

@@ -19,6 +19,7 @@
 from tqdm import trange
 
 from ipsuite import base
+from ipsuite.calculators.integrators import StochasticVelocityCellRescaling
 from ipsuite.utils.ase_sim import freeze_copy_atoms, get_box_from_density, get_energy
 
 log = logging.getLogger(__name__)
@@ -359,6 +360,46 @@ def get_thermostat(self, atoms):
  return thermostat
 
 
+class SVCRBarostat(base.IPSNode):
+ """Initialize the CSVR thermostat
+
+ Attributes
+ ----------
+ time_step: float
+ time step of simulation
+
+ temperature: float
+ temperature in K to simulate at
+ betaT: float
+ Very approximate compressibility of the system.
+ pressure_au: float
+ Pressure in atomic units.
+ taut: float
+ Temperature coupling time scale.
+ taup: float
+ Pressure coupling time scale.
+ """
+
+ time_step: int = zntrack.params()
+ temperature: float = zntrack.params()
+ betaT: float = zntrack.params(4.57e-5 / units.bar)
+ pressure_au: float = zntrack.params(1.01325 * units.bar)
+ taut: float = zntrack.params(100 * units.fs)
+ taup: float = zntrack.params(1000 * units.fs)
+
+ def get_thermostat(self, atoms):
+ thermostat = StochasticVelocityCellRescaling(
+ atoms=atoms,
+ timestep=self.time_step * units.fs,
+ temperature_K=self.temperature,
+ betaT=self.betaT,
+ pressure_au=self.pressure_au,
+ taut=self.taut,
+ taup=self.taup,
+ )
+ return thermostat
+
+
 class FixedSphereConstraint(base.IPSNode):
  """Attributes
  ----------

ipsuite/calculators/integrators.py

@@ -0,0 +1,150 @@
+import math
+
+import numpy as np
+from ase import units
+from ase.md.md import MolecularDynamics
+from ase.parallel import world
+
+
+class StochasticVelocityCellRescaling(MolecularDynamics):
+ """Bussi stochastic velocity and cell rescaling (NVT / NPT) molecular dynamics.
+ Based on the paper from Bussi et al. (https://arxiv.org/abs/0803.4060)
+ Thermostat is based on the ASE implementation of SVR.
+
+ Parameters
+ ----------
+ atoms : Atoms
+ The atoms object.
+ timestep : float
+ The time step in ASE time units.
+ temperature_K : float
+ The desired temperature, in Kelvin.
+ taut : float
+ Time constant for Bussi temperature coupling in ASE time units.
+ rng : numpy.random, optional
+ Random number generator.
+ **md_kwargs : dict, optional
+ Additional arguments passed to :class:~ase.md.md.MolecularDynamics
+ base class.
+ """
+
+ def __init__(
+ self,
+ atoms,
+ timestep,
+ temperature_K,
+ betaT=4.57e-5 / units.bar,
+ pressure_au=1.01325 * units.bar,
+ taut=100 * units.fs,
+ taup=1000 * units.fs,
+ rng=np.random,
+ **md_kwargs,
+ ):
+ super().__init__(
+ atoms,
+ timestep,
+ **md_kwargs,
+ )
+
+ self.temp = temperature_K * units.kB
+ self.taut = taut
+ self.taup = taup
+ self.communicator = world
+ self.rng = rng
+ self.betaT = betaT
+ self.pressure = pressure_au
+
+ self.ndof = len(self.atoms) * 3
+
+ self.target_kinetic_energy = 0.5 * self.temp * self.ndof
+
+ if np.isclose(self.atoms.get_kinetic_energy(), 0.0, rtol=0, atol=1e-12):
+ raise ValueError(
+ "Initial kinetic energy is zero. "
+ "Please set the initial velocities before running Bussi NVT."
+ )
+
+ self._exp_term = math.exp(-self.dt / self.taut)
+ self._masses = self.atoms.get_masses()[:, np.newaxis]
+
+ self.transferred_energy = 0.0
+
+ def set_temperature(self, temperature=None, temperature_K=None):
+ self.temp = units.kB * self._process_temperature(temperature, temperature_K, "eV")
+
+ def scale_velocities(self):
+ """Do the NVT Bussi stochastic velocity scaling."""
+ kinetic_energy = self.atoms.get_kinetic_energy()
+ alpha = self.calculate_alpha(kinetic_energy)
+
+ momenta = self.atoms.get_momenta()
+ self.atoms.set_momenta(alpha * momenta)
+
+ self.transferred_energy += (alpha**2 - 1.0) * kinetic_energy
+
+ def calculate_alpha(self, kinetic_energy):
+ """Calculate the scaling factor alpha using equation (A7)
+ from the Bussi paper."""
+
+ energy_scaling_term = (
+ (1 - self._exp_term) * self.target_kinetic_energy / kinetic_energy / self.ndof
+ )
+
+ # R1 in Eq. (A7)
+ normal_noise = self.rng.standard_normal()
+ # \sum_{i=2}^{Nf} R_i^2 in Eq. (A7)
+ # 2 * standard_gamma(n / 2) is equal to chisquare(n)
+ sum_of_noises = 2.0 * self.rng.standard_gamma(0.5 * (self.ndof - 1))
+
+ return math.sqrt(
+ self._exp_term
+ + energy_scaling_term * (sum_of_noises + normal_noise**2)
+ + 2 * normal_noise * math.sqrt(self._exp_term * energy_scaling_term)
+ )
+
+ def scale_positions_and_cell(self):
+ """Do the Berendsen pressure coupling,
+ scale the atom position and the simulation cell."""
+
+ stress = self.atoms.get_stress(voigt=False, include_ideal_gas=True)
+
+ volume = self.atoms.cell.volume
+ pint = -stress.trace() / 3
+ pint += 2.0 / 3.0 * self.atoms.get_kinetic_energy() / volume
+
+ dW = self.rng.standard_normal(size=1)
+ deterministic = -self.betaT / self.taup * (self.pressure - pint) * self.dt
+ stochastic = (
+ np.sqrt(2 * self.temp * self.betaT / volume / self.taup * self.dt) * dW
+ )
+ depsilon = deterministic + stochastic
+
+ scaling = np.exp(depsilon / 3.0)
+
+ cell = self.atoms.get_cell()
+ cell = scaling * cell
+ self.atoms.set_cell(cell, scale_atoms=True)
+
+ velocities = self.atoms.get_velocities()
+ velocities = velocities / scaling
+ self.atoms.set_velocities(velocities)
+
+ def step(self, forces=None):
+ """Move one timestep forward using Bussi NVT molecular dynamics."""
+ if forces is None:
+ forces = self.atoms.get_forces(md=True)
+
+ self.scale_positions_and_cell()
+
+ self.scale_velocities()
+
+ self.atoms.set_momenta(self.atoms.get_momenta() + 0.5 * self.dt * forces)
+ momenta = self.atoms.get_momenta()
+
+ self.atoms.set_positions(self.atoms.positions + self.dt * momenta / self._masses)
+
+ forces = self.atoms.get_forces(md=True)
+
+ self.atoms.set_momenta(momenta + 0.5 * self.dt * forces)
+
+ return forces

ipsuite/nodes.py

@@ -94,6 +94,7 @@ class _Nodes:
  LangevinThermostat = "ipsuite.calculators.LangevinThermostat"
  VelocityVerletDynamic = "ipsuite.calculators.VelocityVerletDynamic"
  NPTThermostat = "ipsuite.calculators.NPTThermostat"
+ CSVR = "ipsuite.calculators.CSVR"
  RescaleBoxModifier = "ipsuite.calculators.RescaleBoxModifier"
  BoxOscillatingRampModifier = "ipsuite.calculators.BoxOscillatingRampModifier"
  TemperatureRampModifier = "ipsuite.calculators.TemperatureRampModifier"

tests/integration/calculators/test_ase_md.py

@@ -8,7 +8,7 @@
 from ipsuite.utils.ase_sim import get_density_from_atoms
 
 
-def test_ase_md(proj_path, cu_box):
+def test_ase_run_md(proj_path, cu_box):
  atoms = []
  for _ in range(5):
  atoms.extend(cu_box)
@@ -20,11 +20,15 @@ def test_ase_md(proj_path, cu_box):
  temperature=1,
  friction=1,
  )
+ barostat = ips.calculators.SVCRBarostat(
+ time_step=1,
+ temperature=1,
+ )
  rescale_box = ips.calculators.RescaleBoxModifier(cell=10)
  temperature_ramp = ips.calculators.TemperatureRampModifier(temperature=300)
  with ips.Project() as project:
  data = ips.AddData(file="cu_box.xyz")
- model = ips.calculators.EMTSinglePoint(data=data.atoms)
+ model = ips.calculators.LJSinglePoint(data=data.atoms)
  md = ips.calculators.ASEMD(
  data=data.atoms,
  model=model,
@@ -46,6 +50,16 @@ def test_ase_md(proj_path, cu_box):
  sampling_rate=1,
  dump_rate=33,
  )
+ md2 = ips.calculators.ASEMD(
+ data=data.atoms,
+ model=model,
+ checks=[checker],
+ modifiers=[rescale_box, temperature_ramp],
+ thermostat=barostat,
+ steps=30,
+ sampling_rate=1,
+ dump_rate=33,
+ )
 
  project.run()