tools/i-pi/ipi/engine/barostats.py

   1 """Contains the classes that deal with constant pressure dynamics.
   2
   3 Copyright (C) 2013, Joshua More and Michele Ceriotti
   4
   5 This program is free software: you can redistribute it and/or modify
   6 it under the terms of the GNU General Public License as published by
   7 the Free Software Foundation, either version 3 of the License, or
   8 (at your option) any later version.
   9
  10 This program is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this program. If not, see <http.//www.gnu.org/licenses/>.
  17
  18
  19 Contains the algorithms which propagate the position and momenta steps in the
  20 constant pressure ensemble. Holds the properties directly related to
  21 these ensembles, such as the internal and external pressure and stress.
  22
  23 Classes:
  24    Barostat: Base barostat class with the generic methods and attributes.
  25    BaroBZP: Generates dynamics with a stochastic barostat -- see
  26             Ceriotti, More, Manolopoulos, Comp. Phys. Comm. 2013 for
  27             implementation details.
  28 """
  29
  30 # NB: this file also contains a 'BaroMHT' class, that follows more closely the
  31 # Martyna, Hughes, Tuckerman implementation of a PIMD barostat. However it is so
  32 # close to the BZP implementation that we disabled it for the sake of simplicity
  33 # BaroMHT: Generates dynamics according to the method of G. Martyna, A.
  34 # Hughes and M. Tuckerman, J. Chem. Phys., 110, 3275.
  35
  36 __all__ = ['Barostat', 'BaroBZP']
  37
  38 import numpy as np
  39 from ipi.utils.depend import *
  40 from ipi.utils.units import *
  41 from ipi.utils.mathtools import eigensystem_ut3x3, invert_ut3x3, exp_ut3x3, det_ut3x3
  42 from ipi.inputs.thermostats import InputThermo
  43 from ipi.engine.thermostats import Thermostat
  44
  45 class Barostat(dobject):
  46    """Base barostat class.
  47
  48    Gives the standard methods and attributes needed in all the barostat classes.
  49
  50    Attributes:
  51       beads: A beads object giving the atoms positions
  52       cell: A cell object giving the system box.
  53       forces: A forces object giving the virial and the forces acting on
  54          each bead.
  55       nm: An object to do the normal mode transformation.
  56       thermostat: A thermostat coupled to the barostat degrees of freedom.
  57       mdof: The number of atomic degrees of freedom
  58
  59    Depend objects:
  60       dt: The time step used in the algorithms. Depends on the simulation dt.
  61       temp: The (classical) simulation temperature. Higher than the physical
  62          temperature by a factor of the number of beads.
  63       tau: The timescale associated with the piston
  64       pext: The external pressure
  65       ebaro: The conserved quantity associated with the barostat.
  66       pot: The potential energy associated with the barostat.
  67       kstress: The system kinetic stress tensor.
  68       stress: The system stress tensor.
  69       press: The system pressure.
  70    """
  71
  72    def __init__(self, dt=None, temp=None, pext=None, tau=None, ebaro=None, thermostat=None):
  73       """Initialises base barostat class.
  74
  75       Note that the external stress and the external pressure are synchronized.
  76       This makes most sense going from the stress to the pressure, but if you
  77       must go in the other direction the stress is assumed to be isotropic.
  78
  79       Args:
  80          dt: Optional float giving the time step for the algorithms. Defaults
  81             to the simulation dt.
  82          temp: Optional float giving the temperature for the thermostat.
  83             Defaults to the simulation temp.
  84          pext: Optional float giving the external pressure.
  85          tau: Optional float giving the time scale associated with the barostat.
  86          ebaro: Optional float giving the conserved quantity already stored
  87             in the barostat initially. Used on restart.
  88          thermostat: The thermostat connected to the barostat degree of freedom.
  89       """
  90
  91       dset(self,"dt",depend_value(name='dt'))
  92       if not dt is None:
  93          self.dt = dt
  94       else: self.dt = 1.0
  95
  96       dset(self, "temp", depend_value(name="temp"))
  97       if not temp is None:
  98          self.temp = temp
  99       else: self.temp = 1.0
 100
 101       dset(self,"tau",depend_value(name='tau'))
 102       if not tau is None:
 103          self.tau = tau
 104       else: self.tau = 1.0
 105
 106       dset(self,"pext",depend_value(name='pext'))
 107       if not pext is None:
 108          self.pext = pext
 109       else: self.pext = 0.0
 110
 111       dset(self,"ebaro",depend_value(name='ebaro'))
 112       if not ebaro is None:
 113          self.ebaro = ebaro
 114       else: self.ebaro = 0.0
 115
 116       if thermostat is None:
 117          thermostat = Thermostat()
 118       self.thermostat = thermostat
 119
 120       # pipes timestep and temperature to the thermostat
 121       deppipe(self,"dt", self.thermostat, "dt")
 122       deppipe(self, "temp", self.thermostat,"temp")
 123
 124
 125    def bind(self, beads, nm, cell, forces, prng=None, fixdof=None):
 126       """Binds beads, cell and forces to the barostat.
 127
 128       This takes a beads object, a cell object and a forcefield object and
 129       makes them members of the barostat. It also then creates the objects that
 130       will hold the data needed in the barostat algorithms and the dependency
 131       network.
 132
 133       Args:
 134          beads: The beads object from which the bead positions are taken.
 135          nm: The normal modes propagator object
 136          cell: The cell object from which the system box is taken.
 137          forces: The forcefield object from which the force and virial are
 138             taken.
 139          prng: The parent PRNG to bind the thermostat to
 140          fixdof: The number of blocked degrees of freedom.
 141       """
 142
 143       self.beads = beads
 144       self.cell = cell
 145       self.forces = forces
 146       self.nm = nm
 147
 148       dset(self,"pot",
 149          depend_value(name='pot', func=self.get_pot,
 150             dependencies=[ dget(cell,"V"), dget(self,"pext") ]))
 151       dset(self,"kstress",
 152          depend_value(name='kstress', func=self.get_kstress,
 153             dependencies=[ dget(beads,"q"), dget(beads,"qc"), dget(beads,"pc"), dget(forces,"f") ]))
 154       dset(self,"stress",
 155          depend_value(name='stress', func=self.get_stress,
 156             dependencies=[ dget(self,"kstress"), dget(cell,"V"), dget(forces,"vir") ]))
 157       dset(self,"press",
 158          depend_value(name='press', func=self.get_press,
 159             dependencies=[ dget(self,"stress") ]))
 160
 161       if fixdof is None:
 162          self.mdof = float(self.beads.natoms)*3.0
 163       else:
 164          self.mdof = float(self.beads.natoms)*3.0 - float(fixdof)
 165
 166    def get_pot(self):
 167       """Calculates the elastic strain energy of the cell."""
 168
 169       # NOTE: since there are nbeads replicas of the unit cell, the enthalpy contains a nbeads factor
 170       return self.cell.V*self.pext*self.beads.nbeads
 171
 172    def get_kstress(self):
 173       """Calculates the quantum centroid virial kinetic stress tensor
 174       estimator.
 175       """
 176
 177       kst = np.zeros((3,3),float)
 178       q = depstrip(self.beads.q)
 179       qc = depstrip(self.beads.qc)
 180       pc = depstrip(self.beads.pc)
 181       m = depstrip(self.beads.m)
 182       na3 = 3*self.beads.natoms
 183       fall = depstrip(self.forces.f)
 184
 185       for b in range(self.beads.nbeads):
 186          for i in range(3):
 187             for j in range(i,3):
 188                kst[i,j] -= np.dot(q[b,i:na3:3] - qc[i:na3:3],
 189                   fall[b,j:na3:3])
 190
 191       # NOTE: In order to have a well-defined conserved quantity, the Nf kT term in the
 192       # diagonal stress estimator must be taken from the centroid kinetic energy.
 193       for i in range(3):
 194          kst[i,i] += np.dot(pc[i:na3:3],pc[i:na3:3]/m) *self.beads.nbeads
 195
 196       return kst
 197
 198    def get_stress(self):
 199       """Calculates the internal stress tensor."""
 200
 201       return (self.kstress + self.forces.vir)/self.cell.V
 202
 203    def get_press(self):
 204       """Calculates the internal pressure."""
 205
 206       return np.trace(self.stress)/3.0
 207
 208    def pstep(self):
 209       """Dummy momenta propagator step."""
 210
 211       pass
 212
 213    def qcstep(self):
 214       """Dummy centroid position propagator step."""
 215
 216       pass
 217
 218
 219 class BaroBZP(Barostat):
 220    """Bussi-Zykova-Parrinello barostat class.
 221
 222    Just extends the standard class adding finite-dt propagators for the barostat
 223    velocities, positions, piston.
 224
 225    Depend objects:
 226       p: The momentum associated with the volume degree of freedom.
 227       m: The mass associated with the volume degree of freedom.
 228    """
 229
 230    def __init__(self, dt=None, temp=None, pext=None, tau=None, ebaro=None, thermostat=None, p=None):
 231       """Initializes BZP barostat.
 232
 233       Args:
 234          dt: Optional float giving the time step for the algorithms. Defaults
 235             to the simulation dt.
 236          temp: Optional float giving the temperature for the thermostat.
 237             Defaults to the simulation temp.
 238          pext: Optional float giving the external pressure.
 239          tau: Optional float giving the time scale associated with the barostat.
 240          ebaro: Optional float giving the conserved quantity already stored
 241             in the barostat initially. Used on restart.
 242          thermostat: The thermostat connected to the barostat degree of freedom.
 243          p: Optional initial volume conjugate momentum. Defaults to 0.
 244       """
 245
 246
 247       super(BaroBZP, self).__init__(dt, temp, pext, tau, ebaro, thermostat)
 248
 249       dset(self,"p", depend_array(name='p', value=np.atleast_1d(0.0)))
 250
 251       if not p is None:
 252          self.p = np.asarray([p])
 253       else:
 254          self.p = 0.0
 255
 256    def bind(self, beads, nm, cell, forces, prng=None, fixdof=None):
 257       """Binds beads, cell and forces to the barostat.
 258
 259       This takes a beads object, a cell object and a forcefield object and
 260       makes them members of the barostat. It also then creates the objects that
 261       will hold the data needed in the barostat algorithms and the dependency
 262       network.
 263
 264       Args:
 265          beads: The beads object from which the bead positions are taken.
 266          nm: The normal modes propagator object
 267          cell: The cell object from which the system box is taken.
 268          forces: The forcefield object from which the force and virial are
 269             taken.
 270          prng: The parent PRNG to bind the thermostat to
 271          fixdof: The number of blocked degrees of freedom.
 272       """
 273
 274       super(BaroBZP, self).bind(beads, nm, cell, forces, prng, fixdof)
 275
 276       # obtain the thermostat mass from the given time constant
 277       # note that the barostat temperature is nbeads times the physical T
 278       dset(self,"m", depend_array(name='m', value=np.atleast_1d(0.0),
 279          func=(lambda:np.asarray([self.tau**2*3*self.beads.natoms*Constants.kb*self.temp])),
 280             dependencies=[ dget(self,"tau"), dget(self,"temp") ] ))
 281
 282       # binds the thermostat to the piston degrees of freedom
 283       self.thermostat.bind(pm=[ self.p, self.m ], prng=prng)
 284
 285       dset(self,"kin",depend_value(name='kin',
 286          func=(lambda:0.5*self.p[0]**2/self.m[0]),
 287             dependencies= [dget(self,"p"), dget(self,"m")] ) )
 288
 289       # the barostat energy must be computed from bits & pieces (overwrite the default)
 290       dset(self, "ebaro", depend_value(name='ebaro', func=self.get_ebaro,
 291          dependencies=[ dget(self, "kin"), dget(self, "pot"),
 292             dget(self.cell, "V"), dget(self, "temp"),
 293                dget(self.thermostat,"ethermo")] ))
 294
 295    def get_ebaro(self):
 296       """Calculates the barostat conserved quantity."""
 297
 298       return self.thermostat.ethermo + self.kin + self.pot - np.log(self.cell.V)*Constants.kb*self.temp
 299
 300
 301    def pstep(self):
 302       """Propagates the momenta for half a time step."""
 303
 304       dthalf = self.dt*0.5
 305       dthalf2 = dthalf**2
 306       dthalf3 = dthalf**3/3.0
 307
 308       # This differs from the BZP thermostat in that it uses just one kT in the propagator.
 309       # This leads to an ensemble equaivalent to Martyna-Hughes-Tuckermann for both fixed and moving COM
 310       # Anyway, it is a small correction so whatever.
 311       self.p += dthalf*3.0*( self.cell.V* ( self.press - self.beads.nbeads*self.pext ) +
 312                 Constants.kb*self.temp )
 313
 314       fc = np.sum(depstrip(self.forces.f),0)/self.beads.nbeads
 315       m = depstrip(self.beads.m3)[0]
 316       pc = depstrip(self.beads.pc)
 317
 318       # I am not 100% sure, but these higher-order terms come from integrating the pressure virial term,
 319       # so they should need to be multiplied by nbeads to be consistent with the equations of motion in the PI context
 320       # again, these are tiny tiny terms so whatever.
 321       self.p += (dthalf2*np.dot(pc,fc/m) + dthalf3*np.dot(fc,fc/m)) * self.beads.nbeads
 322
 323       self.beads.p += depstrip(self.forces.f)*dthalf
 324
 325    def qcstep(self):
 326       """Propagates the centroid position and momentum and the volume."""
 327
 328       v = self.p[0]/self.m[0]
 329       expq, expp = (np.exp(v*self.dt), np.exp(-v*self.dt))
 330
 331       m = depstrip(self.beads.m3)[0]
 332
 333       self.nm.qnm[0,:] *= expq
 334       self.nm.qnm[0,:] += ((expq-expp)/(2.0*v))* (depstrip(self.nm.pnm)[0,:]/m)
 335       self.nm.pnm[0,:] *= expp
 336
 337       self.cell.h *= expq
 338
 339
 340 class BaroMHT(Barostat):
 341    """Martyna-Hughes-Tuckerman barostat class.
 342
 343    Just extends the standard class adding finite-dt propagators for the barostat
 344    velocities, positions, piston.
 345
 346    Depend objects:
 347       p: The momentum associated with the volume degree of freedom.
 348       m: The mass associated with the volume degree of freedom.
 349    """
 350
 351    def __init__(self, dt=None, temp=None, pext=None, tau=None, ebaro=None, thermostat=None, p=None):
 352       """Initializes MHT barostat.
 353
 354       Args:
 355          dt: Optional float giving the time step for the algorithms. Defaults
 356             to the simulation dt.
 357          temp: Optional float giving the temperature for the thermostat.
 358             Defaults to the simulation temp.
 359          pext: Optional float giving the external pressure.
 360          tau: Optional float giving the time scale associated with the barostat.
 361          ebaro: Optional float giving the conserved quantity already stored
 362             in the barostat initially. Used on restart.
 363          thermostat: The thermostat connected to the barostat degree of freedom.
 364          p: Optional initial volume conjugate momentum. Defaults to 0.
 365       """
 366
 367       super(BaroMHT, self).__init__(dt, temp, pext, tau, ebaro, thermostat)
 368
 369       dset(self,"p", depend_array(name='p', value=np.atleast_1d(0.0)))
 370
 371       if not p is None:
 372          self.p = np.asarray([p])
 373       else:
 374          self.p = 0.0
 375
 376    def bind(self, beads, nm, cell, forces, prng=None, fixdof=None):
 377       """Binds beads, cell and forces to the barostat.
 378
 379       This takes a beads object, a cell object and a forcefield object and
 380       makes them members of the barostat. It also then creates the objects that
 381       will hold the data needed in the barostat algorithms and the dependency
 382       network.
 383
 384       Args:
 385          beads: The beads object from which the bead positions are taken.
 386          nm: The normal modes propagator object
 387          cell: The cell object from which the system box is taken.
 388          forces: The forcefield object from which the force and virial are
 389             taken.
 390          prng: The parent PRNG to bind the thermostat to
 391          fixdof: The number of blocked degrees of freedom.
 392       """
 393
 394       super(BaroMHT, self).bind(beads, nm, cell, forces, prng, fixdof)
 395
 396       # obtain the thermostat mass from the given time constant
 397       # note that the barostat temperature is nbeads times the physical T
 398       dset(self,"m", depend_array(name='m', value=np.atleast_1d(0.0),
 399          func=(lambda:np.asarray([self.tau**2*3*self.beads.natoms*Constants.kb*self.temp])),
 400             dependencies=[ dget(self,"tau"), dget(self,"temp") ] ))
 401
 402       # binds the thermostat to the piston degrees of freedom
 403       self.thermostat.bind(pm=[ self.p, self.m ], prng=prng)
 404
 405       dset(self,"kin",depend_value(name='kin',
 406          func=(lambda:0.5*self.p[0]**2/self.m[0]),
 407             dependencies=[dget(self,"p"), dget(self,"m")] ) )
 408
 409       # the barostat energy must be computed from bits & pieces (overwrite the default)
 410       dset(self, "ebaro", depend_value(name='ebaro', func=self.get_ebaro,
 411          dependencies=[ dget(self, "kin"), dget(self, "pot"),
 412             dget(self.cell, "V"), dget(self, "temp"),
 413                dget(self.thermostat,"ethermo")]))
 414
 415    def get_ebaro(self):
 416       """Calculates the barostat conserved quantity."""
 417
 418       return self.thermostat.ethermo + self.kin + self.pot
 419
 420    def pstep(self):
 421       """Propagates the momenta for half a time step."""
 422
 423       dthalf = self.dt*0.5
 424       dthalf2 = dthalf**2
 425       dthalf3 = dthalf**3/3.0
 426
 427       fc = np.sum(depstrip(self.forces.f),0)/float(self.beads.nbeads)
 428       m = depstrip(self.beads.m3)[0]
 429       pc = depstrip(self.beads.pc)
 430
 431       self.p += dthalf*3.0*( self.cell.V* ( self.press - self.beads.nbeads*self.pext ) +
 432                 float(self.beads.nbeads)/self.mdof*np.dot(pc,pc/m) )
 433
 434       self.beads.p += depstrip(self.forces.f)*dthalf
 435
 436    def qcstep(self):
 437       """Propagates the centroid position and momentum and the volume."""
 438
 439       v = self.p[0]/self.m[0]
 440       adof = (1 + 3.0/self.mdof)
 441       expq, expp = (np.exp(v*self.dt), np.exp( -v*self.dt * adof  ) )
 442
 443       m = depstrip(self.beads.m3)[0]
 444
 445       self.nm.qnm[0,:] *= expq
 446       self.nm.qnm[0,:] += ((expq-expp)/(v*(1+adof)) *
 447                     (depstrip(self.nm.pnm)[0,:])/m)
 448       self.nm.pnm[0,:] *= expp
 449
 450       self.cell.h *= expq