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

   1 """Contains the classes that deal with the different dynamics required in
   2 different types of ensembles.
   3
   4 Copyright (C) 2013, Joshua More and Michele Ceriotti
   5
   6 This program is free software: you can redistribute it and/or modify
   7 it under the terms of the GNU General Public License as published by
   8 the Free Software Foundation, either version 3 of the License, or
   9 (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  14 GNU General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program. If not, see <http.//www.gnu.org/licenses/>.
  18
  19
  20 Holds the algorithms required for normal mode propagators, and the objects to
  21 do the constant temperature and pressure algorithms. Also calculates the
  22 appropriate conserved energy quantity for the ensemble of choice.
  23
  24 Classes:
  25    Ensemble: Base ensemble class with generic methods and attributes.
  26    NVEEnsemble: Deals with constant energy dynamics.
  27    NVTEnsemble: Deals with constant temperature dynamics.
  28    NPTEnsemble: Deals with constant pressure dynamics.
  29    ReplayEnsemble: Takes a trajectory, and simply sets the atom positions to
  30       match it, rather than doing dynamics. In this way new properties can
  31       be calculated on an old simulation, without having to rerun it from
  32       scratch.
  33 """
  34
  35 __all__ = ['Ensemble', 'NVEEnsemble', 'NVTEnsemble', 'NPTEnsemble', 'ReplayEnsemble']
  36
  37 import numpy as np
  38 import time
  39
  40 from ipi.utils.depend import *
  41 from ipi.utils import units
  42 from ipi.utils.softexit import softexit
  43 from ipi.utils.io.io_xyz import read_xyz
  44 from ipi.utils.io.io_pdb import read_pdb
  45 from ipi.utils.io.io_xml import xml_parse_file
  46 from ipi.utils.units import Constants, unit_to_internal
  47 from ipi.inputs.thermostats import InputThermo
  48 from ipi.inputs.barostats import InputBaro
  49 from ipi.engine.thermostats import *
  50 from ipi.engine.barostats import *
  51
  52
  53 class Ensemble(dobject):
  54    """Base ensemble class.
  55
  56    Gives the standard methods and attributes needed in all the
  57    ensemble classes.
  58
  59    Attributes:
  60       beads: A beads object giving the atoms positions.
  61       cell: A cell object giving the system box.
  62       forces: A forces object giving the virial and the forces acting on
  63          each bead.
  64       prng: A random number generator object.
  65       nm: An object which does the normal modes transformation.
  66       fixcom: A boolean which decides whether the centre of mass
  67          motion will be constrained or not.
  68
  69    Depend objects:
  70       econs: The conserved energy quantity appropriate to the given
  71          ensemble. Depends on the various energy terms which make it up,
  72          which are different depending on the ensemble.
  73       temp: The system temperature.
  74       dt: The timestep for the algorithms.
  75       ntemp: The simulation temperature. Will be nbeads times higher than
  76          the system temperature as PIMD calculations are done at this
  77          effective classical temperature.
  78    """
  79
  80    def __init__(self, dt, temp, fixcom=False):
  81       """Initialises Ensemble.
  82
  83       Args:
  84          dt: The timestep of the simulation algorithms.
  85          temp: The temperature.
  86          fixcom: An optional boolean which decides whether the centre of mass
  87             motion will be constrained or not. Defaults to False.
  88       """
  89
  90       dset(self, "econs", depend_value(name='econs', func=self.get_econs))
  91       dset(self, "temp",  depend_value(name='temp',  value=temp))
  92       dset(self, "dt",    depend_value(name='dt',    value=dt))
  93       self.fixcom = fixcom
  94
  95
  96    def bind(self, beads, nm, cell, bforce, prng):
  97       """Binds beads, cell, bforce and prng to the ensemble.
  98
  99       This takes a beads object, a cell object, a forcefield object and a
 100       random number generator object and makes them members of the ensemble.
 101       It also then creates the objects that will hold the data needed in the
 102       ensemble algorithms and the dependency network. Note that the conserved
 103       quantity is defined in the init, but as each ensemble has a different
 104       conserved quantity the dependencies are defined in bind.
 105
 106       Args:
 107          beads: The beads object from whcih the bead positions are taken.
 108          nm: A normal modes object used to do the normal modes transformation.
 109          cell: The cell object from which the system box is taken.
 110          bforce: The forcefield object from which the force and virial are
 111             taken.
 112          prng: The random number generator object which controls random number
 113             generation.
 114       """
 115
 116       # store local references to the different bits of the simulation
 117       self.beads = beads
 118       self.cell = cell
 119       self.forces = bforce
 120       self.prng = prng
 121       self.nm = nm
 122
 123       # n times the temperature
 124       dset(self,"ntemp", depend_value(name='ntemp',func=self.get_ntemp,
 125          dependencies=[dget(self,"temp")]))
 126
 127       # dependencies of the conserved quantity
 128       dget(self,"econs").add_dependency(dget(self.beads, "kin"))
 129       dget(self,"econs").add_dependency(dget(self.forces, "pot"))
 130       dget(self,"econs").add_dependency(dget(self.beads, "vpath"))
 131
 132
 133    def get_ntemp(self):
 134       """Returns the PI simulation temperature (P times the physical T)."""
 135
 136       return self.temp*self.beads.nbeads
 137
 138
 139    def pstep(self):
 140       """Dummy momenta propagator which does nothing."""
 141
 142       pass
 143
 144    def qcstep(self):
 145       """Dummy centroid position propagator which does nothing."""
 146
 147       pass
 148
 149    def step(self):
 150       """Dummy simulation time step which does nothing."""
 151
 152       pass
 153
 154    def get_econs(self):
 155       """Calculates the conserved energy quantity for constant energy
 156       ensembles.
 157       """
 158
 159       return self.beads.vpath*self.nm.omegan2 + self.nm.kin + self.forces.pot
 160
 161
 162 class NVEEnsemble(Ensemble):
 163    """Ensemble object for constant energy simulations.
 164
 165    Has the relevant conserved quantity and normal mode propagator for the
 166    constant energy ensemble. Note that a temperature of some kind must be
 167    defined so that the spring potential can be calculated.
 168
 169    Attributes:
 170       ptime: The time taken in updating the velocities.
 171       qtime: The time taken in updating the positions.
 172       ttime: The time taken in applying the thermostat steps.
 173
 174    Depend objects:
 175       econs: Conserved energy quantity. Depends on the bead kinetic and
 176          potential energy, and the spring potential energy.
 177    """
 178
 179    def __init__(self, dt, temp, fixcom=False):
 180       """Initialises NVEEnsemble.
 181
 182       Args:
 183          dt: The simulation timestep.
 184          temp: The system temperature.
 185          fixcom: An optional boolean which decides whether the centre of mass
 186             motion will be constrained or not. Defaults to False.
 187       """
 188
 189       super(NVEEnsemble,self).__init__(dt=dt,temp=temp, fixcom=fixcom)
 190
 191    def rmcom(self):
 192       """This removes the centre of mass contribution to the kinetic energy.
 193
 194       Calculates the centre of mass momenta, then removes the mass weighted
 195       contribution from each atom. If the ensemble defines a thermostat, then
 196       the contribution to the conserved quantity due to this subtraction is
 197       added to the thermostat heat energy, as it is assumed that the centre of
 198       mass motion is due to the thermostat.
 199
 200       If there is a choice of thermostats, the thermostat
 201       connected to the centroid is chosen.
 202       """
 203
 204       if (self.fixcom):
 205          pcom = np.zeros(3,float);
 206
 207          na3 = self.beads.natoms*3
 208          nb = self.beads.nbeads
 209          p = depstrip(self.beads.p)
 210          m = depstrip(self.beads.m3)[:,0:na3:3]
 211          M = self.beads[0].M
 212
 213          for i in range(3):
 214             pcom[i] = p[:,i:na3:3].sum()
 215
 216          if hasattr(self,"thermostat"):
 217             if hasattr(self.thermostat, "_thermos"):
 218                self.thermostat._thermos[0].ethermo += np.dot(pcom,pcom)/(2.0*M*nb)
 219             else:
 220                self.thermostat.ethermo += np.dot(pcom,pcom)/(2.0*M*nb)
 221
 222          # subtracts COM _velocity_
 223          pcom *= 1.0/(nb*M)
 224          for i in range(3):
 225             self.beads.p[:,i:na3:3] -= m*pcom[i]
 226
 227    def pstep(self):
 228       """Velocity Verlet momenta propagator."""
 229
 230       self.beads.p += depstrip(self.forces.f)*(self.dt*0.5)
 231
 232    def qcstep(self):
 233       """Velocity Verlet centroid position propagator."""
 234
 235       self.nm.qnm[0,:] += depstrip(self.nm.pnm)[0,:]/depstrip(self.beads.m3)[0]*self.dt
 236
 237    def step(self):
 238       """Does one simulation time step."""
 239
 240       self.ptime = -time.time()
 241       self.pstep()
 242       self.ptime += time.time()
 243
 244       self.qtime = -time.time()
 245       self.qcstep()
 246
 247       self.nm.free_qstep()
 248       self.qtime += time.time()
 249
 250       self.ptime -= time.time()
 251       self.pstep()
 252       self.ptime += time.time()
 253
 254       self.ttime = -time.time()
 255       self.rmcom()
 256       self.ttime += time.time()
 257
 258
 259 class NVTEnsemble(NVEEnsemble):
 260    """Ensemble object for constant temperature simulations.
 261
 262    Has the relevant conserved quantity and normal mode propagator for the
 263    constant temperature ensemble. Contains a thermostat object containing the
 264    algorithms to keep the temperature constant.
 265
 266    Attributes:
 267       thermostat: A thermostat object to keep the temperature constant.
 268
 269    Depend objects:
 270       econs: Conserved energy quantity. Depends on the bead kinetic and
 271          potential energy, the spring potential energy and the heat
 272          transferred to the thermostat.
 273    """
 274
 275    def __init__(self, dt, temp, thermostat=None, fixcom=False):
 276       """Initialises NVTEnsemble.
 277
 278       Args:
 279          dt: The simulation timestep.
 280          temp: The system temperature.
 281          thermostat: A thermostat object to keep the temperature constant.
 282             Defaults to Thermostat()
 283          fixcom: An optional boolean which decides whether the centre of mass
 284             motion will be constrained or not. Defaults to False.
 285       """
 286
 287       super(NVTEnsemble,self).__init__(dt=dt,temp=temp, fixcom=fixcom)
 288
 289       if thermostat is None:
 290          self.thermostat = Thermostat()
 291       else:
 292          self.thermostat = thermostat
 293
 294    def bind(self, beads, nm, cell, bforce, prng):
 295       """Binds beads, cell, bforce and prng to the ensemble.
 296
 297       This takes a beads object, a cell object, a forcefield object and a
 298       random number generator object and makes them members of the ensemble.
 299       It also then creates the objects that will hold the data needed in the
 300       ensemble algorithms and the dependency network. Also note that the
 301       thermostat timestep and temperature are defined relative to the system
 302       temperature, and the the thermostat temperature is held at the
 303       higher simulation temperature, as is appropriate.
 304
 305       Args:
 306          beads: The beads object from whcih the bead positions are taken.
 307          nm: A normal modes object used to do the normal modes transformation.
 308          cell: The cell object from which the system box is taken.
 309          bforce: The forcefield object from which the force and virial are
 310             taken.
 311          prng: The random number generator object which controls random number
 312             generation.
 313       """
 314
 315       super(NVTEnsemble,self).bind(beads, nm, cell, bforce, prng)
 316       fixdof = None
 317       if self.fixcom:
 318          fixdof = 3
 319
 320       # first makes sure that the thermostat has the correct temperature, then proceed with binding it.
 321       deppipe(self,"ntemp", self.thermostat,"temp")
 322       deppipe(self,"dt", self.thermostat, "dt")
 323
 324       #decides whether the thermostat will work in the normal mode or
 325       #the bead representation.
 326       if isinstance(self.thermostat,ThermoNMGLE) or isinstance(self.thermostat,ThermoNMGLEG) or isinstance(self.thermostat,ThermoPILE_L) or isinstance(self.thermostat,ThermoPILE_G):
 327          self.thermostat.bind(nm=self.nm,prng=prng,fixdof=fixdof )
 328       else:
 329          self.thermostat.bind(beads=self.beads,prng=prng, fixdof=fixdof)
 330
 331       dget(self,"econs").add_dependency(dget(self.thermostat, "ethermo"))
 332
 333    def step(self):
 334       """Does one simulation time step."""
 335
 336       self.ttime = -time.time()
 337       self.thermostat.step()
 338       self.rmcom()
 339       self.ttime += time.time()
 340
 341       self.ptime = -time.time()
 342       self.pstep()
 343       self.ptime += time.time()
 344
 345       self.qtime = -time.time()
 346       self.qcstep()
 347       self.nm.free_qstep()
 348       self.qtime += time.time()
 349
 350       self.ptime -= time.time()
 351       self.pstep()
 352       self.ptime += time.time()
 353
 354       self.ttime -= time.time()
 355       self.thermostat.step()
 356       self.rmcom()
 357       self.ttime += time.time()
 358
 359    def get_econs(self):
 360       """Calculates the conserved energy quantity for constant temperature
 361       ensemble.
 362       """
 363
 364       return NVEEnsemble.get_econs(self) + self.thermostat.ethermo
 365
 366
 367 class NPTEnsemble(NVTEnsemble):
 368    """Ensemble object for constant pressure simulations.
 369
 370    Has the relevant conserved quantity and normal mode propagator for the
 371    constant pressure ensemble. Contains a thermostat object containing the
 372    algorithms to keep the temperature constant, and a barostat to keep the
 373    pressure constant.
 374
 375    Attributes:
 376       barostat: A barostat object to keep the pressure constant.
 377
 378    Depend objects:
 379       econs: Conserved energy quantity. Depends on the bead and cell kinetic
 380          and potential energy, the spring potential energy, the heat
 381          transferred to the beads and cell thermostat, the temperature and
 382          the cell volume.
 383       pext: External pressure.
 384    """
 385
 386    def __init__(self, dt, temp, pext, thermostat=None, barostat=None, fixcom=False):
 387       """Initialises NPTEnsemble.
 388
 389       Args:
 390          dt: The simulation timestep.
 391          temp: The system temperature.
 392          pext: The external pressure.
 393          thermostat: A thermostat object to keep the temperature constant.
 394             Defaults to Thermostat().
 395          barostat: A barostat object to keep the pressure constant.
 396             Defaults to Barostat().
 397          fixcom: An optional boolean which decides whether the centre of mass
 398             motion will be constrained or not. Defaults to False.
 399       """
 400
 401       super(NPTEnsemble,self).__init__(dt, temp, thermostat, fixcom=fixcom)
 402       if barostat == None:
 403          self.barostat = Barostat()
 404       else:
 405          self.barostat = barostat
 406
 407       dset(self,"pext",depend_value(name='pext'))
 408       if not pext is None:
 409          self.pext = pext
 410       else: self.pext = 0.0
 411
 412
 413    def bind(self, beads, nm, cell, bforce, prng):
 414       """Binds beads, cell, bforce and prng to the ensemble.
 415
 416       This takes a beads object, a cell object, a forcefield object and a
 417       random number generator object and makes them members of the ensemble.
 418       It also then creates the objects that will hold the data needed in the
 419       ensemble algorithms and the dependency network. Also note that the cell
 420       thermostat timesteps and temperatures are defined relative to the system
 421       temperature, and the the thermostat temperatures are held at the
 422       higher simulation temperature, as is appropriate.
 423
 424       Args:
 425          beads: The beads object from whcih the bead positions are taken.
 426          nm: A normal modes object used to do the normal modes transformation.
 427          cell: The cell object from which the system box is taken.
 428          bforce: The forcefield object from which the force and virial are
 429             taken.
 430          prng: The random number generator object which controls random number
 431             generation.
 432       """
 433
 434
 435       fixdof = None
 436       if self.fixcom:
 437          fixdof = 3
 438
 439       super(NPTEnsemble,self).bind(beads, nm, cell, bforce, prng)
 440       self.barostat.bind(beads, nm, cell, bforce, prng=prng, fixdof=fixdof)
 441
 442
 443       deppipe(self,"ntemp", self.barostat, "temp")
 444       deppipe(self,"dt", self.barostat, "dt")
 445       deppipe(self,"pext", self.barostat, "pext")
 446       dget(self,"econs").add_dependency(dget(self.barostat, "ebaro"))
 447
 448    def get_econs(self):
 449       """Calculates the conserved energy quantity for the constant pressure
 450       ensemble.
 451       """
 452
 453       return NVTEnsemble.get_econs(self) + self.barostat.ebaro
 454
 455    def step(self):
 456       """NPT time step.
 457
 458       Note that the barostat only propagates the centroid coordinates. If this
 459       approximation is made a centroid virial pressure and stress estimator can
 460       be defined, so this gives the best statistical convergence. This is
 461       allowed as the normal mode propagation is approximately unaffected
 462       by volume fluctuations as long as the system box is much larger than
 463       the radius of gyration of the ring polymers.
 464       """
 465
 466       self.ttime = -time.time()
 467       self.thermostat.step()
 468       self.barostat.thermostat.step()
 469       self.rmcom()
 470       self.ttime += time.time()
 471
 472       self.ptime = -time.time()
 473       self.barostat.pstep()
 474       self.ptime += time.time()
 475
 476       self.qtime = -time.time()
 477       self.barostat.qcstep()
 478       self.nm.free_qstep()
 479       self.qtime += time.time()
 480
 481       self.ptime -= time.time()
 482       self.barostat.pstep()
 483       self.ptime += time.time()
 484
 485       self.ttime -= time.time()
 486       self.barostat.thermostat.step()
 487       self.thermostat.step()
 488       self.rmcom()
 489       self.ttime += time.time()
 490
 491
 492 class ReplayEnsemble(Ensemble):
 493    """Ensemble object that just loads snapshots from an external file in sequence.
 494
 495    Has the relevant conserved quantity and normal mode propagator for the
 496    constant energy ensemble. Note that a temperature of some kind must be
 497    defined so that the spring potential can be calculated.
 498
 499    Attributes:
 500       intraj: The input trajectory file.
 501       ptime: The time taken in updating the velocities.
 502       qtime: The time taken in updating the positions.
 503       ttime: The time taken in applying the thermostat steps.
 504
 505    Depend objects:
 506       econs: Conserved energy quantity. Depends on the bead kinetic and
 507          potential energy, and the spring potential energy.
 508    """
 509
 510    def __init__(self, dt, temp, fixcom=False, intraj=None):
 511       """Initialises ReplayEnsemble.
 512
 513       Args:
 514          dt: The simulation timestep.
 515          temp: The system temperature.
 516          fixcom: An optional boolean which decides whether the centre of mass
 517             motion will be constrained or not. Defaults to False.
 518          intraj: The input trajectory file.
 519       """
 520
 521       super(ReplayEnsemble,self).__init__(dt=dt,temp=temp,fixcom=fixcom)
 522       if intraj == None:
 523          raise ValueError("Must provide an initialized InitFile object to read trajectory from")
 524       self.intraj = intraj
 525       if intraj.mode == "manual":
 526          raise ValueError("Replay can only read from PDB or XYZ files -- or a single frame from a CHK file")
 527       self.rfile = open(self.intraj.value,"r")
 528
 529    def step(self):
 530       """Does one simulation time step."""
 531
 532       self.ptime = self.ttime = 0
 533       self.qtime = -time.time()
 534
 535       try:
 536          if (self.intraj.mode == "xyz"):
 537             for b in self.beads:
 538                myatoms = read_xyz(self.rfile)
 539                myatoms.q *= unit_to_internal("length",self.intraj.units,1.0)
 540                b.q[:] = myatoms.q
 541          elif (self.intraj.mode == "pdb"):
 542             for b in self.beads:
 543                myatoms, mycell = read_pdb(self.rfile)
 544                myatoms.q *= unit_to_internal("length",self.intraj.units,1.0)
 545                mycell.h  *= unit_to_internal("length",self.intraj.units,1.0)
 546                b.q[:] = myatoms.q
 547             self.cell.h[:] = mycell.h
 548          elif (self.intraj.mode == "chk" or self.intraj.mode == "checkpoint"):
 549             # reads configuration from a checkpoint file
 550             xmlchk = xml_parse_file(self.rfile) # Parses the file.
 551
 552             from ipi.inputs.simulation import InputSimulation
 553             simchk = InputSimulation()
 554             simchk.parse(xmlchk.fields[0][1])
 555             mycell = simchk.cell.fetch()
 556             mybeads = simchk.beads.fetch()
 557             self.cell.h[:] = mycell.h
 558             self.beads.q[:] = mybeads.q
 559             softexit.trigger(" # Read single checkpoint")
 560       except EOFError:
 561          softexit.trigger(" # Finished reading re-run trajectory")
 562
 563       self.qtime += time.time()
 564
 565