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

   1 """Contains the class that deals with the running of the simulation and
   2 outputting the results.
   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 The root class for the whole simulation. Contains references to all the top
  21 level objects used in the simulation, and controls all the steps that are
  22 not inherently system dependent, like the running of each time step,
  23 choosing which properties to initialise, and which properties to output.
  24
  25 Classes:
  26    Simulation: Deals with running the simulation and outputting the results.
  27 """
  28
  29 __all__ = ['Simulation']
  30
  31 import numpy as np
  32 import os.path, sys, time
  33 from ipi.utils.depend import *
  34 from ipi.utils.units  import *
  35 from ipi.utils.prng   import *
  36 from ipi.utils.io     import *
  37 from ipi.utils.io.io_xml import *
  38 from ipi.utils.messages import verbosity, info
  39 from ipi.utils.softexit import softexit
  40 from ipi.engine.atoms import *
  41 from ipi.engine.cell import *
  42 from ipi.engine.forces import Forces
  43 from ipi.engine.beads import Beads
  44 from ipi.engine.normalmodes import NormalModes
  45 from ipi.engine.properties import Properties, Trajectories
  46 from ipi.engine.outputs import CheckpointOutput
  47
  48 class Simulation(dobject):
  49    """Main simulation object.
  50
  51    Contains all the references and the main dynamics loop. Also handles the
  52    initialisation and output.
  53
  54    Attributes:
  55       beads: A beads object giving the atom positions.
  56       cell: A cell object giving the system box.
  57       prng: A random number generator object.
  58       flist: A list of forcefield objects giving different ways to partially
  59          calculate the forces.
  60       forces: A Forces object for calculating the total force for all the
  61          replicas.
  62       ensemble: An ensemble object giving the objects necessary for producing
  63          the correct ensemble.
  64       tsteps: The total number of steps.
  65       ttime: The wall clock time (in seconds).
  66       format: A string specifying both the format and the extension of traj
  67          output.
  68       outputs: A list of output objects that should be printed during the run
  69       nm:  A helper object dealing with normal modes transformation
  70       properties: A property object for dealing with property output.
  71       trajs: A trajectory object for dealing with trajectory output.
  72       chk: A checkpoint object for dealing with checkpoint output.
  73       rollback: If set to true, the state of the simulation at the start
  74          of the step will be output to a restart file rather than
  75          the current state of the simulation. This is because we cannot
  76          restart from half way through a step, only from the beginning of a
  77          step, so this is necessary for the trajectory to be continuous.
  78
  79    Depend objects:
  80       step: The current simulation step.
  81    """
  82
  83    def __init__(self, beads, cell, forces, ensemble, prng, outputs, nm, init, step=0, tsteps=1000, ttime=0):
  84       """Initialises Simulation class.
  85
  86       Args:
  87          beads: A beads object giving the atom positions.
  88          cell: A cell object giving the system box.
  89          forces: A forcefield object giving the force calculator for each
  90             replica of the system.
  91          ensemble: An ensemble object giving the objects necessary for
  92             producing the correct ensemble.
  93          prng: A random number object.
  94          outputs: A list of output objects.
  95          nm: A class dealing with path NM operations.
  96          init: A class to deal with initializing the simulation object.
  97          step: An optional integer giving the current simulation time step.
  98             Defaults to 0.
  99          tsteps: An optional integer giving the total number of steps. Defaults
 100             to 1000.
 101          ttime: The simulation running time. Used on restart, to keep a
 102             cumulative total.
 103       """
 104
 105       info(" # Initializing simulation object ", verbosity.low )
 106       self.prng = prng
 107       self.ensemble = ensemble
 108       self.beads = beads
 109       self.cell = cell
 110       self.nm = nm
 111
 112       # initialize the configuration of the system
 113       self.init = init
 114       init.init_stage1(self)
 115
 116       self.flist = forces
 117       self.forces = Forces()
 118       self.outputs = outputs
 119
 120       dset(self, "step", depend_value(name="step", value=step))
 121       self.tsteps = tsteps
 122       self.ttime = ttime
 123
 124       self.properties = Properties()
 125       self.trajs = Trajectories()
 126       self.chk = None
 127       self.rollback = True
 128
 129    def bind(self):
 130       """Calls the bind routines for all the objects in the simulation."""
 131
 132       # binds important computation engines
 133       self.nm.bind(self.beads, self.ensemble)
 134       self.forces.bind(self.beads, self.cell, self.flist)
 135       self.ensemble.bind(self.beads, self.nm, self.cell, self.forces, self.prng)
 136       self.init.init_stage2(self)
 137
 138       # binds output management objects
 139       self.properties.bind(self)
 140       self.trajs.bind(self)
 141       for o in self.outputs:
 142          o.bind(self)
 143
 144       self.chk = CheckpointOutput("RESTART", 1, True, 0)
 145       self.chk.bind(self)
 146
 147       # registers the softexit routine
 148       softexit.register(self.softexit)
 149
 150    def softexit(self):
 151       """Deals with a soft exit request.
 152
 153       Tries to ensure that a consistent restart checkpoint is
 154       written out.
 155       """
 156
 157       if self.step < self.tsteps:
 158          self.step += 1
 159       if not self.rollback:
 160          self.chk.store()
 161       self.chk.write(store=False)
 162
 163       self.forces.stop()
 164
 165    def run(self):
 166       """Runs the simulation.
 167
 168       Does all the simulation steps, and outputs data to the appropriate files
 169       when necessary. Also deals with starting and cleaning up the threads used
 170       in the communication between the driver and the PIMD code.
 171       """
 172
 173       self.forces.run()
 174
 175       # prints inital configuration -- only if we are not restarting
 176       if (self.step == 0):
 177          self.step = -1
 178          for o in self.outputs:
 179             o.write()
 180          self.step = 0
 181
 182       steptime = 0.0
 183       simtime =  time.time()
 184
 185       cstep = 0
 186       tptime = 0.0
 187       tqtime = 0.0
 188       tttime = 0.0
 189       ttot = 0.0
 190       # main MD loop
 191       for self.step in range(self.step,self.tsteps):
 192          # stores the state before doing a step.
 193          # this is a bit time-consuming but makes sure that we can honor soft
 194          # exit requests without screwing the trajectory
 195
 196          steptime = -time.time()
 197          self.chk.store()
 198
 199          self.ensemble.step()
 200
 201          for o in self.outputs:
 202             o.write()
 203
 204          if os.path.exists("EXIT"): # soft-exit
 205             self.rollback = False
 206             softexit.trigger()
 207
 208          steptime += time.time()
 209          ttot += steptime
 210          tptime += self.ensemble.ptime
 211          tqtime += self.ensemble.qtime
 212          tttime += self.ensemble.ttime
 213          cstep += 1
 214
 215          if verbosity.high or (verbosity.medium and self.step%100 == 0) or (verbosity.low and self.step%1000 == 0):
 216             info(" # Average timings at MD step % 7d. t/step: %10.5e [p: %10.5e  q: %10.5e  t: %10.5e]" %
 217                ( self.step, ttot/cstep, tptime/cstep, tqtime/cstep, tttime/cstep ) )
 218             cstep = 0
 219             tptime = 0.0
 220             tqtime = 0.0
 221             tttime = 0.0
 222             ttot = 0.0
 223             info(" # MD diagnostics: V: %10.5e    Kcv: %10.5e   Ecns: %10.5e" %
 224                (self.properties["potential"], self.properties["kinetic_cv"], self.properties["conserved"] ) )
 225
 226          if (self.ttime > 0 and time.time() - simtime > self.ttime):
 227             info(" # Wall clock time expired! Bye bye!", verbosity.low )
 228             break
 229
 230       info(" # Simulation ran successfully for the prescribed total_step! Bye bye!", verbosity.low )
 231       self.rollback = False
 232       softexit.trigger()