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

   1 """Contains the classes that are used to initialize data in the simulation.
   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 These classes can either be used to restart a simulation with some different
  20 data or used to start a calculation. Any data given in these classes will
  21 overwrite data given elsewhere.
  22
  23 Classes:
  24    Initializer: Holds the functions that are required to initialize objects in
  25       the code. Data can be initialized from a file, or according to a
  26       particular parameter. An example of the former would be initializing
  27       the configurations from a xyz file, an example of the latter would be
  28       initializing the velocities according to the physical temperature.
  29    InitBase: Simple class that reads data from a string or file.
  30    InitIndexed: The same as init base, but can also optionally hold
  31       information about which atom or bead to initialize from.
  32
  33 Functions:
  34    init_xyz: Reads beads data from a xyz file.
  35    init_pdb: Reads beads and cell data from a pdb file.
  36    init_chk: Reads beads, cell and thermostat data from a checkpoint file.
  37    init_beads: Initializes a beads object from an Initializer object.
  38    init_vector: Initializes a vector from an Initializer object.
  39    set_vector: Initializes a vector from another vector.
  40 """
  41
  42 import numpy as np
  43
  44 from ipi.engine.beads import Beads
  45 from ipi.engine.cell import Cell
  46 from ipi.engine.normalmodes import NormalModes
  47 from ipi.engine.ensembles import Ensemble
  48 from ipi.utils.io.io_xyz import read_xyz
  49 from ipi.utils.io.io_pdb import read_pdb
  50 from ipi.utils.io.io_xml import xml_parse_file
  51 from ipi.utils.depend import dobject
  52 from ipi.utils.units import Constants, unit_to_internal
  53 from ipi.utils.nmtransform import nm_rescale
  54 from ipi.utils.messages import verbosity, warning, info
  55
  56 __all__ = ['Initializer', 'InitBase', 'InitIndexed']
  57
  58 class InitBase(dobject):
  59    """Base class for initializer objects.
  60
  61    Attributes:
  62       value: A duck-typed stored value.
  63       mode: A string that determines how the value is to be interpreted.
  64       units: A string giving which unit the value is in.
  65    """
  66
  67    def __init__(self, value="", mode="", units="", **others):
  68       """Initializes InitFile.
  69
  70       Args:
  71          value: A string which specifies what value to initialize the
  72             simulation property to.
  73          mode: A string specifiying what style of initialization should be
  74             used to read the data.
  75          units: A string giving which unit the value is in.
  76       """
  77
  78       self.value = value
  79       self.mode = mode
  80       self.units = units
  81
  82       for (o, v) in others.items():
  83          self.__dict__[o] = v
  84
  85
  86 class InitIndexed(InitBase):
  87    """Class to initialize objects which can be set for a particular bead.
  88
  89    Attributes:
  90       index: Which atom to initialize the value of.
  91       bead: Which bead to initialize the value of.
  92    """
  93
  94    def __init__(self, value="", mode="", units="", index=-1, bead=-1):
  95       """Initializes InitFile.
  96
  97       Args:
  98          value: A string which specifies what value to initialize the
  99             simulation property to.
 100          mode: A string specifiying what style of initialization should be
 101             used to read the data.
 102          units: A string giving which unit the value is in.
 103          index: Which atom to initialize the value of.
 104          bead: Which bead to initialize the value of.
 105       """
 106
 107       super(InitIndexed,self).__init__(value=value, mode=mode, units=units, index=index, bead=bead)
 108
 109
 110 def init_xyz(filename):
 111    """Reads an xyz file and returns the data contained in it.
 112
 113    Args:
 114       filename: A string giving the name of the xyz file to be read from.
 115
 116    Returns:
 117       A list of Atoms objects as read from each frame of the xyz file.
 118    """
 119
 120    rfile = open(filename,"r")
 121    ratoms = []
 122    while True:
 123    #while loop, so that more than one configuration can be given
 124    #so multiple beads can be initialized at once.
 125       try:
 126          myatoms = read_xyz(rfile)
 127       except EOFError:
 128          break
 129       ratoms.append(myatoms)
 130    return ratoms
 131
 132 def init_pdb(filename):
 133    """Reads an pdb file and returns the data contained in it.
 134
 135    Args:
 136       filename: A string giving the name of the pdb file to be read from.
 137
 138    Returns:
 139       A list of Atoms objects as read from each frame of the pdb file, and
 140       a Cell object as read from the final pdb frame.
 141    """
 142
 143    rfile = open(filename,"r")
 144    ratoms = []
 145    while True:
 146    #while loop, so that more than one configuration can be given
 147    #so multiple beads can be initialized at once.
 148       try:
 149          myatoms, rcell  = read_pdb(rfile)
 150       except EOFError:
 151          break
 152       ratoms.append(myatoms)
 153    return ( ratoms, rcell ) # if multiple frames, the last cell is returned
 154
 155 def init_chk(filename):
 156    """Reads an checkpoint file and returns the data contained in it.
 157
 158    Args:
 159       filename: A string giving the name of the checkpoint file to be read from.
 160
 161    Returns:
 162       A Beads object, Cell object and Thermostat object as read from the
 163       checkpoint file.
 164    """
 165
 166    # reads configuration from a checkpoint file
 167    rfile = open(filename,"r")
 168    xmlchk = xml_parse_file(rfile) # Parses the file.
 169
 170    from ipi.inputs.simulation import InputSimulation
 171    simchk = InputSimulation()
 172    simchk.parse(xmlchk.fields[0][1])
 173    rcell = simchk.cell.fetch()
 174    rbeads = simchk.beads.fetch()
 175    rthermo = simchk.ensemble.thermostat.fetch()
 176
 177    return (rbeads, rcell, rthermo)
 178
 179 def init_beads(iif, nbeads):
 180    """A file to initialize a beads object from an appropriate initializer
 181    object.
 182
 183    Args:
 184       iif: An Initializer object which has information on the bead positions.
 185       nbeads: The number of beads.
 186
 187    Raises:
 188       ValueError: If called using an Initializer object with a 'manual' mode.
 189    """
 190
 191    mode = iif.mode; value = iif.value
 192    if mode == "xyz" or mode == "pdb":
 193       if mode == "xyz": ratoms = init_xyz(value)
 194       if mode == "pdb": ratoms = init_pdb(value)[0]
 195       rbeads = Beads(ratoms[0].natoms,len(ratoms))
 196       for i in range(len(ratoms)): rbeads[i] = ratoms[i]
 197    elif mode == "chk":
 198       rbeads = init_chk(value)[0]
 199    elif mode == "manual":
 200       raise ValueError("Cannot initialize manually a whole beads object.")
 201
 202    return rbeads
 203
 204 def init_vector(iif, nbeads, momenta=False):
 205    """A file to initialize a vector from an appropriate initializer
 206    object.
 207
 208    Args:
 209       iif: An Initializer object specifying the value of a vector.
 210       nbeads: The number of beads.
 211       momenta: If bead momenta rather than positions are being initialized
 212          from a checkpoint file, this is set to True.
 213    """
 214
 215    mode = iif.mode; value = iif.value
 216    if mode == "xyz" or mode == "pdb":
 217       rq = init_beads(iif, nbeads).q
 218    elif mode == "chk":
 219       if momenta: rq = init_beads(iif, nbeads).p
 220       else:       rq = init_beads(iif, nbeads).q
 221    elif mode == "manual":
 222       rq = value
 223
 224    # determines the size of the input data
 225    if mode == "manual":
 226       if iif.bead >= 0: # if there is a bead specifier then we return a single bead slice
 227          nbeads = 1
 228       natoms = len(rq)/nbeads/3
 229       rq.shape = (nbeads,3*natoms)
 230
 231    return rq
 232
 233 def set_vector(iif, dq, rq):
 234    """A file to initialize a vector from an another vector.
 235
 236    If the first dimension is different, i.e. the two vectors correspond
 237    to a different number of beads, then the ring polymer contraction/expansion
 238    is used to rescale the original vector to the one used in the simulation,
 239    as described in the paper T. E. Markland and D. E. Manolopoulos, J. Chem.
 240    Phys. 129, 024105, (2008).
 241
 242    Args:
 243       iif: An Initializer object specifying the value of a vector.
 244       dq: The vector to be initialized.
 245       rq: The vector to initialize from.
 246    """
 247
 248    (nbeads, natoms) = rq.shape; natoms /= 3
 249    (dbeads, datoms) = dq.shape; datoms /= 3
 250
 251    # Check that indices make sense
 252    if iif.index < 0 and natoms != datoms:
 253       raise ValueError("Initialization tries to mix up structures with different atom numbers.")
 254    if iif.index >= datoms:
 255       raise ValueError("Cannot initialize single atom as atom index %d is larger than the number of atoms" % iif.index)
 256    if iif.bead >= dbeads:
 257       raise ValueError("Cannot initialize single bead as bead index %d is larger than the number of beads" % iif.bead)
 258
 259    if iif.bead < 0:   # we are initializing the path
 260       res = nm_rescale(nbeads,dbeads)  # path rescaler
 261       if nbeads != dbeads:
 262          info(" # Initialize is rescaling from %5d beads to %5d beads" % (nbeads, dbeads), verbosity.low)
 263       if iif.index < 0:
 264          dq[:] = res.b1tob2(rq)
 265       else: # we are initializing a specific atom
 266          dq[:,3*iif.index:3*(iif.index+1)] = res.b1tob2(rq)
 267    else:  # we are initializing a specific bead
 268       if iif.index < 0:
 269          dq[iif.bead] = rq
 270       else:
 271          dq[iif.bead,3*iif.index:3*(iif.index+1)] = rq
 272
 273 class Initializer(dobject):
 274    """Class that deals with the initialization of data.
 275
 276    This can either be used to initialize the atom positions and the cell data
 277    from a file, or to initialize them from a beads, atoms or cell object.
 278
 279    Currently, we use a ring polymer contraction scheme to create a new beads
 280    object from one given in initialize if they have different numbers of beads,
 281    as described in the paper T. E. Markland and D. E. Manolopoulos, J. Chem.
 282    Phys. 129, 024105, (2008). If the new beads object has more beads than
 283    the beads object it was initialized from, we set the higher ring polymer
 284    normal modes to zero.
 285
 286    Attributes:
 287       queue: A list of things to initialize. Each member of the list is a tuple
 288          of the form ('type', 'object'), where 'type' specifies what kind of
 289          initialization is being done, and 'object' gives the data to
 290          initialize it from.
 291    """
 292
 293    def __init__(self, nbeads=0, queue=None):
 294       """Initializes Initializer.
 295
 296       Arguments:
 297          nbeads: The number of beads that we need in the simulation. Not
 298             necessarily the same as the number of beads of the objects we are
 299             initializing the data from.
 300          queue: A list of things to initialize. Each member of the list is a
 301             tuple of the form ('type', 'object'), where 'type' specifies what
 302             kind of initialization is being done, and 'object' gives the data to
 303             initialize it from.
 304       """
 305
 306       self.nbeads = nbeads
 307
 308       if queue is None:
 309          self.queue = []
 310       else:
 311          self.queue = queue
 312
 313    def init_stage1(self, simul):
 314       """Initializes the simulation -- first stage.
 315
 316       Takes a simulation object, and uses all the data in the initialization
 317       queue to fill up the beads and cell data needed to run the simulation.
 318
 319       Args:
 320          simul: A simulation object to be initialized.
 321
 322       Raises:
 323          ValueError: Raised if there is a problem with the initialization,
 324             if something that should have been has not been, or if the objects
 325             that have been specified are not compatible with each other.
 326       """
 327
 328       if simul.beads.nbeads == 0:
 329          fpos = fmom = fmass = flab = fcell = False   # we don't have an explicitly defined beads object yet
 330       else:
 331          fpos = fmom = fmass = flab = fcell = True
 332
 333       for (k,v) in self.queue:
 334          info(" # Initializer (stage 1) parsing " + str(k) + " object.", verbosity.high)
 335
 336          if k == "cell":
 337             if fcell :
 338                warning("Overwriting previous cell parameters", verbosity.medium)
 339             if v.mode == "pdb":
 340                rh = init_pdb(v.value)[1].h
 341             elif v.mode == "chk":
 342                rh = init_chk(v.value)[1].h
 343             else:
 344                rh = v.value.reshape((3,3))
 345             rh *= unit_to_internal("length",v.units,1.0)
 346
 347             simul.cell.h = rh
 348             if simul.cell.V == 0.0:
 349                ValueError("Cell provided has zero volume")
 350
 351             fcell = True
 352          elif k == "masses":
 353             if simul.beads.nbeads == 0:
 354                raise ValueError("Cannot initialize the masses before the size of the system is known")
 355             if fmass:
 356                warning("Overwriting previous atomic masses", verbosity.medium)
 357             if v.mode == "manual":
 358                rm = v.value
 359             else:
 360                rm = init_beads(v, self.nbeads).m
 361             rm *= unit_to_internal("mass",v.units,1.0)
 362
 363             if v.bead < 0:   # we are initializing the path
 364                if (fmom and fmass):
 365                   warning("Rescaling momenta to make up for changed mass", verbosity.medium)
 366                   simul.beads.p /= simul.beads.sm3   # go to mass-scaled momenta, that are mass-invariant
 367                if v.index < 0:
 368                   simul.beads.m = rm
 369                else: # we are initializing a specific atom
 370                   simul.beads.m[v.index:v.index+1] = rm
 371                if (fmom and fmass):  # finishes correcting the momenta
 372                   simul.beads.p *= simul.beads.sm3  # back to normal momenta
 373             else:
 374                raise ValueError("Cannot change the mass of a single bead")
 375             fmass = True
 376
 377          elif k == "labels":
 378             if simul.beads.nbeads == 0:
 379                raise ValueError("Cannot initialize the labels before the size of the system is known")
 380             if flab:
 381                warning("Overwriting previous atomic labels", verbosity.medium)
 382             if v.mode == "manual":
 383                rn = v.value
 384             else:
 385                rn = init_beads(v, self.nbeads).names
 386
 387             if v.bead < 0:   # we are initializing the path
 388                if v.index < 0:
 389                   simul.beads.names = rn
 390                else: # we are initializing a specific atom
 391                   simul.beads.names[v.index:v.index+1] = rn
 392             else:
 393                raise ValueError("Cannot change the label of a single bead")
 394             flab = True
 395
 396          elif k == "positions":
 397             if fpos:
 398                warning("Overwriting previous atomic positions", verbosity.medium)
 399             # read the atomic positions as a vector
 400             rq = init_vector(v, self.nbeads)
 401             rq *= unit_to_internal("length",v.units,1.0)
 402             (nbeads, natoms) = rq.shape;   natoms /= 3
 403
 404             # check if we must initialize the simulation beads
 405             if simul.beads.nbeads == 0:
 406                if v.index >= 0:
 407                   raise ValueError("Cannot initialize single atoms before the size of the system is known")
 408                simul.beads.resize(natoms,self.nbeads)
 409
 410             set_vector(v, simul.beads.q, rq)
 411             fpos = True
 412
 413          elif (k == "velocities" or k == "momenta") and v.mode == "thermal" :   # intercept here thermal initialization, so we don't need to check further down
 414             if fmom:
 415                warning("Overwriting previous atomic momenta", verbosity.medium)
 416             if simul.beads.natoms == 0:
 417                raise ValueError("Cannot initialize momenta before the size of the system is known.")
 418             if not fmass:
 419                raise ValueError("Trying to resample velocities before having masses.")
 420
 421             rtemp = v.value * unit_to_internal("temperature",v.units,1.0)
 422             if rtemp <= 0:
 423                warning("Using the simulation temperature to resample velocities", verbosity.low)
 424                rtemp = simul.ensemble.temp
 425             else:
 426                info(" # Resampling velocities at temperature %s %s" % (v.value, v.units), verbosity.low)
 427
 428             # pull together a mock initialization to get NM masses right
 429             #without too much code duplication
 430             if v.bead >= 0:
 431                raise ValueError("Cannot thermalize a single bead")
 432             if v.index >= 0:
 433                rnatoms = 1
 434             else:
 435                rnatoms = simul.beads.natoms
 436             rbeads = Beads(rnatoms, simul.beads.nbeads)
 437             if v.index < 0:
 438                rbeads.m[:] = simul.beads.m
 439             else:
 440                rbeads.m[:] = simul.beads.m[v.index]
 441             rnm = NormalModes(mode=simul.nm.mode, transform_method=simul.nm.transform_method, freqs=simul.nm.nm_freqs)
 442             rens = Ensemble(dt=simul.ensemble.dt, temp=simul.ensemble.temp)
 443             rnm.bind(rbeads,rens)
 444             # then we exploit the sync magic to do a complicated initialization
 445             # in the NM representation
 446             # with (possibly) shifted-frequencies NM
 447             rnm.pnm = simul.prng.gvec((rbeads.nbeads,3*rbeads.natoms))*np.sqrt(rnm.dynm3)*np.sqrt(rbeads.nbeads*rtemp*Constants.kb)
 448
 449             if v.index < 0:
 450                simul.beads.p = rbeads.p
 451             else:
 452                simul.beads.p[:,3*v.index:3*(v.index+1)] = rbeads.p
 453             fmom = True
 454
 455          elif k == "momenta":
 456             if fmom:
 457                warning("Overwriting previous atomic momenta", verbosity.medium)
 458             # read the atomic momenta as a vector
 459             rp = init_vector(v, self.nbeads, momenta = True)
 460             rp *= unit_to_internal("momentum",v.units,1.0)
 461             (nbeads, natoms) = rp.shape;   natoms /= 3
 462
 463             # checks if we must initialize the simulation beads
 464             if simul.beads.nbeads == 0:
 465                if v.index >= 0 :
 466                   raise ValueError("Cannot initialize single atoms before the size of the system is known")
 467                simul.beads.resize(natoms,self.nbeads)
 468
 469             rp *= np.sqrt(self.nbeads/nbeads)
 470             set_vector(v, simul.beads.p, rp)
 471             fmom = True
 472
 473          elif k == "velocities":
 474             if fmom:
 475                warning("Overwriting previous atomic momenta", verbosity.medium)
 476             # read the atomic velocities as a vector
 477             rv = init_vector(v, self.nbeads)
 478             rv *= unit_to_internal("velocity",v.units,1.0)
 479             (nbeads, natoms) = rv.shape;   natoms /= 3
 480
 481             # checks if we must initialize the simulation beads
 482             if simul.beads.nbeads == 0 or not fmass:
 483                ValueError("Cannot initialize velocities before the masses of the atoms are known")
 484                simul.beads.resize(natoms,self.nbeads)
 485
 486             warning("Initializing from velocities uses the previously defined masses -- not the masses inferred from the file -- to build momenta", verbosity.low)
 487             if v.index >= 0:
 488                rv *= simul.beads.m[v.index]
 489             elif v.bead >= 0:
 490                rv *= simul.beads.m3[0]
 491             else:
 492                rv *= simul.beads.m3
 493             rv *= np.sqrt(self.nbeads/nbeads)
 494             set_vector(v, simul.beads.p, rv)
 495             fmom = True
 496          elif k == "thermostat": pass   # thermostats must be initialised in a second stage
 497
 498       if simul.beads.natoms == 0:
 499          raise ValueError("Initializer could not initialize the atomic positions")
 500       if simul.cell.V == 0:
 501          raise ValueError("Initializer could not initialize the cell")
 502       for i in range(simul.beads.natoms):
 503          if simul.beads.m[i] <= 0:
 504             raise ValueError("Initializer could not initialize the masses")
 505          if simul.beads.names[i] == "":
 506             raise ValueError("Initializer could not initialize the atom labels")
 507       if not fmom:
 508          warning("Momenta not specified in initialize. Will start with zero velocity if they are not specified in beads.", verbosity.low)
 509
 510    def init_stage2(self, simul):
 511       """Initializes the simulation -- second stage.
 512
 513       Takes a simulation object which has been fully generated,
 514       and restarts additional information such as the thermostat internal state.
 515
 516       Args:
 517          simul: A simulation object to be initialized.
 518
 519       Raises:
 520          ValueError: Raised if there is a problem with the initialization,
 521             if something that should have been has not been, or if the objects
 522             that have been specified are not compatible with each other.
 523       """
 524
 525       for (k,v) in self.queue:
 526          info(" # Initializer (stage 2) parsing " + str(k) + " object.", verbosity.high)
 527
 528          if k == "gle":
 529             # read thermostat parameters from file
 530             if not ( hasattr(simul.ensemble, "thermostat") ):
 531                raise ValueError("Ensemble does not have a thermostat to initialize")
 532             if not ( hasattr(simul.ensemble.thermostat, "s") ):
 533                raise ValueError("There is nothing to initialize in non-GLE thermostats")
 534             ssimul = simul.ensemble.thermostat.s
 535             if v.mode == "manual":
 536                sinput = v.value.copy()
 537                if (sinput.size() != ssimul.size() ):
 538                   raise ValueError("Size mismatch in thermostat initialization data")
 539                sinput.shape = ssimul.shape
 540             elif v.mode == "chk":
 541                rthermo = init_chk(v.value)[2]
 542                if not hasattr(rthermo,"s"):
 543                   raise ValueError("Checkpoint file does not contain usable thermostat data")
 544                sinput = rthermo.s.copy()
 545                if sinput.shape != ssimul.shape :
 546                   raise ValueError("Shape mismatch in thermostat initialization data")
 547
 548             # if all the preliminary checks are good, we can initialize the s's
 549             ssimul[:] = sinput