Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80937225
initializer.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 3, 08:41
Size
21 KB
Mime Type
text/x-python
Expires
Thu, Sep 5, 08:41 (2 d)
Engine
blob
Format
Raw Data
Handle
20482590
Attached To
rLAMMPS lammps
initializer.py
View Options
"""Contains the classes that are used to initialize data in the simulation.
Copyright (C) 2013, Joshua More and Michele Ceriotti
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http.//www.gnu.org/licenses/>.
These classes can either be used to restart a simulation with some different
data or used to start a calculation. Any data given in these classes will
overwrite data given elsewhere.
Classes:
Initializer: Holds the functions that are required to initialize objects in
the code. Data can be initialized from a file, or according to a
particular parameter. An example of the former would be initializing
the configurations from a xyz file, an example of the latter would be
initializing the velocities according to the physical temperature.
InitBase: Simple class that reads data from a string or file.
InitIndexed: The same as init base, but can also optionally hold
information about which atom or bead to initialize from.
Functions:
init_xyz: Reads beads data from a xyz file.
init_pdb: Reads beads and cell data from a pdb file.
init_chk: Reads beads, cell and thermostat data from a checkpoint file.
init_beads: Initializes a beads object from an Initializer object.
init_vector: Initializes a vector from an Initializer object.
set_vector: Initializes a vector from another vector.
"""
import
numpy
as
np
from
ipi.engine.beads
import
Beads
from
ipi.engine.cell
import
Cell
from
ipi.engine.normalmodes
import
NormalModes
from
ipi.engine.ensembles
import
Ensemble
from
ipi.utils.io.io_xyz
import
read_xyz
from
ipi.utils.io.io_pdb
import
read_pdb
from
ipi.utils.io.io_xml
import
xml_parse_file
from
ipi.utils.depend
import
dobject
from
ipi.utils.units
import
Constants
,
unit_to_internal
from
ipi.utils.nmtransform
import
nm_rescale
from
ipi.utils.messages
import
verbosity
,
warning
,
info
__all__
=
[
'Initializer'
,
'InitBase'
,
'InitIndexed'
]
class
InitBase
(
dobject
):
"""Base class for initializer objects.
Attributes:
value: A duck-typed stored value.
mode: A string that determines how the value is to be interpreted.
units: A string giving which unit the value is in.
"""
def
__init__
(
self
,
value
=
""
,
mode
=
""
,
units
=
""
,
**
others
):
"""Initializes InitFile.
Args:
value: A string which specifies what value to initialize the
simulation property to.
mode: A string specifiying what style of initialization should be
used to read the data.
units: A string giving which unit the value is in.
"""
self
.
value
=
value
self
.
mode
=
mode
self
.
units
=
units
for
(
o
,
v
)
in
others
.
items
():
self
.
__dict__
[
o
]
=
v
class
InitIndexed
(
InitBase
):
"""Class to initialize objects which can be set for a particular bead.
Attributes:
index: Which atom to initialize the value of.
bead: Which bead to initialize the value of.
"""
def
__init__
(
self
,
value
=
""
,
mode
=
""
,
units
=
""
,
index
=-
1
,
bead
=-
1
):
"""Initializes InitFile.
Args:
value: A string which specifies what value to initialize the
simulation property to.
mode: A string specifiying what style of initialization should be
used to read the data.
units: A string giving which unit the value is in.
index: Which atom to initialize the value of.
bead: Which bead to initialize the value of.
"""
super
(
InitIndexed
,
self
)
.
__init__
(
value
=
value
,
mode
=
mode
,
units
=
units
,
index
=
index
,
bead
=
bead
)
def
init_xyz
(
filename
):
"""Reads an xyz file and returns the data contained in it.
Args:
filename: A string giving the name of the xyz file to be read from.
Returns:
A list of Atoms objects as read from each frame of the xyz file.
"""
rfile
=
open
(
filename
,
"r"
)
ratoms
=
[]
while
True
:
#while loop, so that more than one configuration can be given
#so multiple beads can be initialized at once.
try
:
myatoms
=
read_xyz
(
rfile
)
except
EOFError
:
break
ratoms
.
append
(
myatoms
)
return
ratoms
def
init_pdb
(
filename
):
"""Reads an pdb file and returns the data contained in it.
Args:
filename: A string giving the name of the pdb file to be read from.
Returns:
A list of Atoms objects as read from each frame of the pdb file, and
a Cell object as read from the final pdb frame.
"""
rfile
=
open
(
filename
,
"r"
)
ratoms
=
[]
while
True
:
#while loop, so that more than one configuration can be given
#so multiple beads can be initialized at once.
try
:
myatoms
,
rcell
=
read_pdb
(
rfile
)
except
EOFError
:
break
ratoms
.
append
(
myatoms
)
return
(
ratoms
,
rcell
)
# if multiple frames, the last cell is returned
def
init_chk
(
filename
):
"""Reads an checkpoint file and returns the data contained in it.
Args:
filename: A string giving the name of the checkpoint file to be read from.
Returns:
A Beads object, Cell object and Thermostat object as read from the
checkpoint file.
"""
# reads configuration from a checkpoint file
rfile
=
open
(
filename
,
"r"
)
xmlchk
=
xml_parse_file
(
rfile
)
# Parses the file.
from
ipi.inputs.simulation
import
InputSimulation
simchk
=
InputSimulation
()
simchk
.
parse
(
xmlchk
.
fields
[
0
][
1
])
rcell
=
simchk
.
cell
.
fetch
()
rbeads
=
simchk
.
beads
.
fetch
()
rthermo
=
simchk
.
ensemble
.
thermostat
.
fetch
()
return
(
rbeads
,
rcell
,
rthermo
)
def
init_beads
(
iif
,
nbeads
):
"""A file to initialize a beads object from an appropriate initializer
object.
Args:
iif: An Initializer object which has information on the bead positions.
nbeads: The number of beads.
Raises:
ValueError: If called using an Initializer object with a 'manual' mode.
"""
mode
=
iif
.
mode
;
value
=
iif
.
value
if
mode
==
"xyz"
or
mode
==
"pdb"
:
if
mode
==
"xyz"
:
ratoms
=
init_xyz
(
value
)
if
mode
==
"pdb"
:
ratoms
=
init_pdb
(
value
)[
0
]
rbeads
=
Beads
(
ratoms
[
0
]
.
natoms
,
len
(
ratoms
))
for
i
in
range
(
len
(
ratoms
)):
rbeads
[
i
]
=
ratoms
[
i
]
elif
mode
==
"chk"
:
rbeads
=
init_chk
(
value
)[
0
]
elif
mode
==
"manual"
:
raise
ValueError
(
"Cannot initialize manually a whole beads object."
)
return
rbeads
def
init_vector
(
iif
,
nbeads
,
momenta
=
False
):
"""A file to initialize a vector from an appropriate initializer
object.
Args:
iif: An Initializer object specifying the value of a vector.
nbeads: The number of beads.
momenta: If bead momenta rather than positions are being initialized
from a checkpoint file, this is set to True.
"""
mode
=
iif
.
mode
;
value
=
iif
.
value
if
mode
==
"xyz"
or
mode
==
"pdb"
:
rq
=
init_beads
(
iif
,
nbeads
)
.
q
elif
mode
==
"chk"
:
if
momenta
:
rq
=
init_beads
(
iif
,
nbeads
)
.
p
else
:
rq
=
init_beads
(
iif
,
nbeads
)
.
q
elif
mode
==
"manual"
:
rq
=
value
# determines the size of the input data
if
mode
==
"manual"
:
if
iif
.
bead
>=
0
:
# if there is a bead specifier then we return a single bead slice
nbeads
=
1
natoms
=
len
(
rq
)
/
nbeads
/
3
rq
.
shape
=
(
nbeads
,
3
*
natoms
)
return
rq
def
set_vector
(
iif
,
dq
,
rq
):
"""A file to initialize a vector from an another vector.
If the first dimension is different, i.e. the two vectors correspond
to a different number of beads, then the ring polymer contraction/expansion
is used to rescale the original vector to the one used in the simulation,
as described in the paper T. E. Markland and D. E. Manolopoulos, J. Chem.
Phys. 129, 024105, (2008).
Args:
iif: An Initializer object specifying the value of a vector.
dq: The vector to be initialized.
rq: The vector to initialize from.
"""
(
nbeads
,
natoms
)
=
rq
.
shape
;
natoms
/=
3
(
dbeads
,
datoms
)
=
dq
.
shape
;
datoms
/=
3
# Check that indices make sense
if
iif
.
index
<
0
and
natoms
!=
datoms
:
raise
ValueError
(
"Initialization tries to mix up structures with different atom numbers."
)
if
iif
.
index
>=
datoms
:
raise
ValueError
(
"Cannot initialize single atom as atom index
%d
is larger than the number of atoms"
%
iif
.
index
)
if
iif
.
bead
>=
dbeads
:
raise
ValueError
(
"Cannot initialize single bead as bead index
%d
is larger than the number of beads"
%
iif
.
bead
)
if
iif
.
bead
<
0
:
# we are initializing the path
res
=
nm_rescale
(
nbeads
,
dbeads
)
# path rescaler
if
nbeads
!=
dbeads
:
info
(
" # Initialize is rescaling from
%5d
beads to
%5d
beads"
%
(
nbeads
,
dbeads
),
verbosity
.
low
)
if
iif
.
index
<
0
:
dq
[:]
=
res
.
b1tob2
(
rq
)
else
:
# we are initializing a specific atom
dq
[:,
3
*
iif
.
index
:
3
*
(
iif
.
index
+
1
)]
=
res
.
b1tob2
(
rq
)
else
:
# we are initializing a specific bead
if
iif
.
index
<
0
:
dq
[
iif
.
bead
]
=
rq
else
:
dq
[
iif
.
bead
,
3
*
iif
.
index
:
3
*
(
iif
.
index
+
1
)]
=
rq
class
Initializer
(
dobject
):
"""Class that deals with the initialization of data.
This can either be used to initialize the atom positions and the cell data
from a file, or to initialize them from a beads, atoms or cell object.
Currently, we use a ring polymer contraction scheme to create a new beads
object from one given in initialize if they have different numbers of beads,
as described in the paper T. E. Markland and D. E. Manolopoulos, J. Chem.
Phys. 129, 024105, (2008). If the new beads object has more beads than
the beads object it was initialized from, we set the higher ring polymer
normal modes to zero.
Attributes:
queue: A list of things to initialize. Each member of the list is a tuple
of the form ('type', 'object'), where 'type' specifies what kind of
initialization is being done, and 'object' gives the data to
initialize it from.
"""
def
__init__
(
self
,
nbeads
=
0
,
queue
=
None
):
"""Initializes Initializer.
Arguments:
nbeads: The number of beads that we need in the simulation. Not
necessarily the same as the number of beads of the objects we are
initializing the data from.
queue: A list of things to initialize. Each member of the list is a
tuple of the form ('type', 'object'), where 'type' specifies what
kind of initialization is being done, and 'object' gives the data to
initialize it from.
"""
self
.
nbeads
=
nbeads
if
queue
is
None
:
self
.
queue
=
[]
else
:
self
.
queue
=
queue
def
init_stage1
(
self
,
simul
):
"""Initializes the simulation -- first stage.
Takes a simulation object, and uses all the data in the initialization
queue to fill up the beads and cell data needed to run the simulation.
Args:
simul: A simulation object to be initialized.
Raises:
ValueError: Raised if there is a problem with the initialization,
if something that should have been has not been, or if the objects
that have been specified are not compatible with each other.
"""
if
simul
.
beads
.
nbeads
==
0
:
fpos
=
fmom
=
fmass
=
flab
=
fcell
=
False
# we don't have an explicitly defined beads object yet
else
:
fpos
=
fmom
=
fmass
=
flab
=
fcell
=
True
for
(
k
,
v
)
in
self
.
queue
:
info
(
" # Initializer (stage 1) parsing "
+
str
(
k
)
+
" object."
,
verbosity
.
high
)
if
k
==
"cell"
:
if
fcell
:
warning
(
"Overwriting previous cell parameters"
,
verbosity
.
medium
)
if
v
.
mode
==
"pdb"
:
rh
=
init_pdb
(
v
.
value
)[
1
]
.
h
elif
v
.
mode
==
"chk"
:
rh
=
init_chk
(
v
.
value
)[
1
]
.
h
else
:
rh
=
v
.
value
.
reshape
((
3
,
3
))
rh
*=
unit_to_internal
(
"length"
,
v
.
units
,
1.0
)
simul
.
cell
.
h
=
rh
if
simul
.
cell
.
V
==
0.0
:
ValueError
(
"Cell provided has zero volume"
)
fcell
=
True
elif
k
==
"masses"
:
if
simul
.
beads
.
nbeads
==
0
:
raise
ValueError
(
"Cannot initialize the masses before the size of the system is known"
)
if
fmass
:
warning
(
"Overwriting previous atomic masses"
,
verbosity
.
medium
)
if
v
.
mode
==
"manual"
:
rm
=
v
.
value
else
:
rm
=
init_beads
(
v
,
self
.
nbeads
)
.
m
rm
*=
unit_to_internal
(
"mass"
,
v
.
units
,
1.0
)
if
v
.
bead
<
0
:
# we are initializing the path
if
(
fmom
and
fmass
):
warning
(
"Rescaling momenta to make up for changed mass"
,
verbosity
.
medium
)
simul
.
beads
.
p
/=
simul
.
beads
.
sm3
# go to mass-scaled momenta, that are mass-invariant
if
v
.
index
<
0
:
simul
.
beads
.
m
=
rm
else
:
# we are initializing a specific atom
simul
.
beads
.
m
[
v
.
index
:
v
.
index
+
1
]
=
rm
if
(
fmom
and
fmass
):
# finishes correcting the momenta
simul
.
beads
.
p
*=
simul
.
beads
.
sm3
# back to normal momenta
else
:
raise
ValueError
(
"Cannot change the mass of a single bead"
)
fmass
=
True
elif
k
==
"labels"
:
if
simul
.
beads
.
nbeads
==
0
:
raise
ValueError
(
"Cannot initialize the labels before the size of the system is known"
)
if
flab
:
warning
(
"Overwriting previous atomic labels"
,
verbosity
.
medium
)
if
v
.
mode
==
"manual"
:
rn
=
v
.
value
else
:
rn
=
init_beads
(
v
,
self
.
nbeads
)
.
names
if
v
.
bead
<
0
:
# we are initializing the path
if
v
.
index
<
0
:
simul
.
beads
.
names
=
rn
else
:
# we are initializing a specific atom
simul
.
beads
.
names
[
v
.
index
:
v
.
index
+
1
]
=
rn
else
:
raise
ValueError
(
"Cannot change the label of a single bead"
)
flab
=
True
elif
k
==
"positions"
:
if
fpos
:
warning
(
"Overwriting previous atomic positions"
,
verbosity
.
medium
)
# read the atomic positions as a vector
rq
=
init_vector
(
v
,
self
.
nbeads
)
rq
*=
unit_to_internal
(
"length"
,
v
.
units
,
1.0
)
(
nbeads
,
natoms
)
=
rq
.
shape
;
natoms
/=
3
# check if we must initialize the simulation beads
if
simul
.
beads
.
nbeads
==
0
:
if
v
.
index
>=
0
:
raise
ValueError
(
"Cannot initialize single atoms before the size of the system is known"
)
simul
.
beads
.
resize
(
natoms
,
self
.
nbeads
)
set_vector
(
v
,
simul
.
beads
.
q
,
rq
)
fpos
=
True
elif
(
k
==
"velocities"
or
k
==
"momenta"
)
and
v
.
mode
==
"thermal"
:
# intercept here thermal initialization, so we don't need to check further down
if
fmom
:
warning
(
"Overwriting previous atomic momenta"
,
verbosity
.
medium
)
if
simul
.
beads
.
natoms
==
0
:
raise
ValueError
(
"Cannot initialize momenta before the size of the system is known."
)
if
not
fmass
:
raise
ValueError
(
"Trying to resample velocities before having masses."
)
rtemp
=
v
.
value
*
unit_to_internal
(
"temperature"
,
v
.
units
,
1.0
)
if
rtemp
<=
0
:
warning
(
"Using the simulation temperature to resample velocities"
,
verbosity
.
low
)
rtemp
=
simul
.
ensemble
.
temp
else
:
info
(
" # Resampling velocities at temperature
%s
%s
"
%
(
v
.
value
,
v
.
units
),
verbosity
.
low
)
# pull together a mock initialization to get NM masses right
#without too much code duplication
if
v
.
bead
>=
0
:
raise
ValueError
(
"Cannot thermalize a single bead"
)
if
v
.
index
>=
0
:
rnatoms
=
1
else
:
rnatoms
=
simul
.
beads
.
natoms
rbeads
=
Beads
(
rnatoms
,
simul
.
beads
.
nbeads
)
if
v
.
index
<
0
:
rbeads
.
m
[:]
=
simul
.
beads
.
m
else
:
rbeads
.
m
[:]
=
simul
.
beads
.
m
[
v
.
index
]
rnm
=
NormalModes
(
mode
=
simul
.
nm
.
mode
,
transform_method
=
simul
.
nm
.
transform_method
,
freqs
=
simul
.
nm
.
nm_freqs
)
rens
=
Ensemble
(
dt
=
simul
.
ensemble
.
dt
,
temp
=
simul
.
ensemble
.
temp
)
rnm
.
bind
(
rbeads
,
rens
)
# then we exploit the sync magic to do a complicated initialization
# in the NM representation
# with (possibly) shifted-frequencies NM
rnm
.
pnm
=
simul
.
prng
.
gvec
((
rbeads
.
nbeads
,
3
*
rbeads
.
natoms
))
*
np
.
sqrt
(
rnm
.
dynm3
)
*
np
.
sqrt
(
rbeads
.
nbeads
*
rtemp
*
Constants
.
kb
)
if
v
.
index
<
0
:
simul
.
beads
.
p
=
rbeads
.
p
else
:
simul
.
beads
.
p
[:,
3
*
v
.
index
:
3
*
(
v
.
index
+
1
)]
=
rbeads
.
p
fmom
=
True
elif
k
==
"momenta"
:
if
fmom
:
warning
(
"Overwriting previous atomic momenta"
,
verbosity
.
medium
)
# read the atomic momenta as a vector
rp
=
init_vector
(
v
,
self
.
nbeads
,
momenta
=
True
)
rp
*=
unit_to_internal
(
"momentum"
,
v
.
units
,
1.0
)
(
nbeads
,
natoms
)
=
rp
.
shape
;
natoms
/=
3
# checks if we must initialize the simulation beads
if
simul
.
beads
.
nbeads
==
0
:
if
v
.
index
>=
0
:
raise
ValueError
(
"Cannot initialize single atoms before the size of the system is known"
)
simul
.
beads
.
resize
(
natoms
,
self
.
nbeads
)
rp
*=
np
.
sqrt
(
self
.
nbeads
/
nbeads
)
set_vector
(
v
,
simul
.
beads
.
p
,
rp
)
fmom
=
True
elif
k
==
"velocities"
:
if
fmom
:
warning
(
"Overwriting previous atomic momenta"
,
verbosity
.
medium
)
# read the atomic velocities as a vector
rv
=
init_vector
(
v
,
self
.
nbeads
)
rv
*=
unit_to_internal
(
"velocity"
,
v
.
units
,
1.0
)
(
nbeads
,
natoms
)
=
rv
.
shape
;
natoms
/=
3
# checks if we must initialize the simulation beads
if
simul
.
beads
.
nbeads
==
0
or
not
fmass
:
ValueError
(
"Cannot initialize velocities before the masses of the atoms are known"
)
simul
.
beads
.
resize
(
natoms
,
self
.
nbeads
)
warning
(
"Initializing from velocities uses the previously defined masses -- not the masses inferred from the file -- to build momenta"
,
verbosity
.
low
)
if
v
.
index
>=
0
:
rv
*=
simul
.
beads
.
m
[
v
.
index
]
elif
v
.
bead
>=
0
:
rv
*=
simul
.
beads
.
m3
[
0
]
else
:
rv
*=
simul
.
beads
.
m3
rv
*=
np
.
sqrt
(
self
.
nbeads
/
nbeads
)
set_vector
(
v
,
simul
.
beads
.
p
,
rv
)
fmom
=
True
elif
k
==
"thermostat"
:
pass
# thermostats must be initialised in a second stage
if
simul
.
beads
.
natoms
==
0
:
raise
ValueError
(
"Initializer could not initialize the atomic positions"
)
if
simul
.
cell
.
V
==
0
:
raise
ValueError
(
"Initializer could not initialize the cell"
)
for
i
in
range
(
simul
.
beads
.
natoms
):
if
simul
.
beads
.
m
[
i
]
<=
0
:
raise
ValueError
(
"Initializer could not initialize the masses"
)
if
simul
.
beads
.
names
[
i
]
==
""
:
raise
ValueError
(
"Initializer could not initialize the atom labels"
)
if
not
fmom
:
warning
(
"Momenta not specified in initialize. Will start with zero velocity if they are not specified in beads."
,
verbosity
.
low
)
def
init_stage2
(
self
,
simul
):
"""Initializes the simulation -- second stage.
Takes a simulation object which has been fully generated,
and restarts additional information such as the thermostat internal state.
Args:
simul: A simulation object to be initialized.
Raises:
ValueError: Raised if there is a problem with the initialization,
if something that should have been has not been, or if the objects
that have been specified are not compatible with each other.
"""
for
(
k
,
v
)
in
self
.
queue
:
info
(
" # Initializer (stage 2) parsing "
+
str
(
k
)
+
" object."
,
verbosity
.
high
)
if
k
==
"gle"
:
# read thermostat parameters from file
if
not
(
hasattr
(
simul
.
ensemble
,
"thermostat"
)
):
raise
ValueError
(
"Ensemble does not have a thermostat to initialize"
)
if
not
(
hasattr
(
simul
.
ensemble
.
thermostat
,
"s"
)
):
raise
ValueError
(
"There is nothing to initialize in non-GLE thermostats"
)
ssimul
=
simul
.
ensemble
.
thermostat
.
s
if
v
.
mode
==
"manual"
:
sinput
=
v
.
value
.
copy
()
if
(
sinput
.
size
()
!=
ssimul
.
size
()
):
raise
ValueError
(
"Size mismatch in thermostat initialization data"
)
sinput
.
shape
=
ssimul
.
shape
elif
v
.
mode
==
"chk"
:
rthermo
=
init_chk
(
v
.
value
)[
2
]
if
not
hasattr
(
rthermo
,
"s"
):
raise
ValueError
(
"Checkpoint file does not contain usable thermostat data"
)
sinput
=
rthermo
.
s
.
copy
()
if
sinput
.
shape
!=
ssimul
.
shape
:
raise
ValueError
(
"Shape mismatch in thermostat initialization data"
)
# if all the preliminary checks are good, we can initialize the s's
ssimul
[:]
=
sinput
Event Timeline
Log In to Comment