Page MenuHomec4science
Files F70488044

beads.py
No OneTemporary
Actions

File Metadata

Created
Sun, Jul 7, 03:11

beads.py
View Options

"""Contains the classes which deal with all the beads.
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/>.
Used for holding information about the beads, including their positions, masses
momenta and kinetic energy. Has different objects for the position and normal
mode representations, and has a special centroid atoms object for when the
centroid coordinate is required.
Classes:
Beads: Class with methods dealing with all the beads.
"""
__all__ = ['Beads']
import numpy as np
from ipi.utils.depend import *
from ipi.engine.atoms import Atoms
from ipi.utils import units
class Beads(dobject):
"""Storage for the beads positions and velocities.
Everything is stored as (nbeads,3*natoms) sized contiguous arrays,
and a convenience-access to each replica of the system is provided through a
list of Atoms objects. Contains arrays of both the normal mode representation
and the position representation, and various sized arrays for the atom
labels and masses. Also contains the potential and force between
neighbouring replicas.
Attributes:
natoms: The number of atoms.
nbeads: The number of beads.
_blist: A list of Atoms objects for each replica of the system. Each
replica is assumed to have the same mass and atom label.
centroid: An atoms object giving the centroid coordinate of the beads.
Depend objects:
names: An array giving the atom names.
m: An array giving the atom masses.
m3: An array giving the mass associated with each degree of freedom.
sm3: An array giving the square root of m3.
q: An array giving all the bead positions.
p: An array giving all the bead momenta.
qc: An array giving the centroid positions. Depends on qnm.
pc: An array giving the centroid momenta. Depends on pnm.
vpath: The spring potential between the beads, divided by omegan**2.
Depends on q.
fpath: The spring force between the beads, divided by omegan**2.
Depends on q.
kins: A list of the kinetic energy of each replica.
kin: The total kinetic energy of the system. Note that this is not the
same as the estimate of the kinetic energy of the system, which is
contained in the properties module.
kstress: The total kinetic stress tensor for the system.
rg: An array giving the radius of gyration of each atom.
"""
def __init__(self, natoms, nbeads):
"""Initialises Beads.
Args:
natoms: Number of atoms.
nbeads: Number of beads.
"""
self.resize(natoms, nbeads)
def resize(self, natoms, nbeads):
"""Creates all the data arrays needed in the simulation.
Effectively initializes the whole Beads object, according to the
specified number of atoms and beads. Is also used, as the name suggests,
to resize the data to a new number of beads when this is necessary, for
example in initialization from a simulation with a different number of
beads.
Also creates, or recreates, the dependency network, as this requires
the data arrays to be created for it to work.
Args:
natoms: The number of atoms.
nbeads: The number of beads.
"""
self.natoms = natoms
self.nbeads = nbeads
dset(self,"names",
depend_array(name="names",value=np.zeros(natoms, np.dtype('|S6'))) )
# atom masses, and mass-related arrays
dset(self,"m",depend_array(name="m",value=np.zeros(natoms, float)) ) # this is the prototype mass array (just one independent of bead n)
dset(self,"m3",
depend_array(name="m3",value=np.zeros((nbeads,3*natoms), float), # this is m conveniently replicated to be (nb,3*nat)
func=self.mtom3, dependencies=[dget(self,"m")]))
dset(self,"sm3",
depend_array(name="sm3",value=np.zeros((nbeads,3*natoms), float), # this is just the square root of m3
func=self.m3tosm3, dependencies=[dget(self,"m3")]))
# positions and momenta. bead representation, base storage used everywhere
dset(self,"q",
depend_array(name="q",value=np.zeros((nbeads,3*natoms), float)) )
dset(self,"p",
depend_array(name="p",value=np.zeros((nbeads,3*natoms), float)) )
# position and momentum of the centroid
dset(self,"qc",
depend_array(name="qc",value=np.zeros(3*natoms, float),
func=self.get_qc, dependencies=[dget(self,"q")] ) )
dset(self,"pc",
depend_array(name="pc",value=np.zeros(3*natoms, float),
func=self.get_pc, dependencies=[dget(self,"p")] ) )
# create proxies to access the centroid and the individual beads as Atoms objects
self.centroid = Atoms(natoms, _prebind=(self.qc, self.pc, self.m, self.names))
self._blist = [Atoms(natoms, _prebind=( self.q[i,:], self.p[i,:], self.m, self.names )) for i in range(nbeads) ]
# path springs potential and force
dset(self,"vpath",
depend_value(name="vpath", func=self.get_vpath,
dependencies=[dget(self,"q")]))
dset(self,"fpath",
depend_array(name="fpath", value=np.zeros((nbeads,3*natoms), float),
func=self.get_fpath, dependencies=[dget(self,"q")]))
# kinetic energies of thhe beads, and total (classical) kinetic stress tensor
dset(self,"kins",
depend_array(name="kins",value=np.zeros(nbeads, float),
func=self.kin_gather,
dependencies=[dget(b,"kin") for b in self._blist]))
dset(self,"kin",
depend_value(name="kin", func=self.get_kin,
dependencies=[dget(self,"kins")]))
dset(self,"kstress",
depend_array(name="kstress",value=np.zeros((3,3), float),
func=self.get_kstress,
dependencies=[dget(b,"kstress") for b in self._blist]))
def copy(self):
"""Creates a new beads object from the original.
Returns:
A Beads object with the same q, p, m and names arrays as the original.
"""
newbd = Beads(self.natoms, self.nbeads)
newbd.q[:] = self.q
newbd.p[:] = self.p
newbd.m[:] = self.m
newbd.names[:] = self.names
return newbd
def m3tosm3(self):
"""Takes the mass array and returns the square rooted mass array."""
return np.sqrt(depstrip(self.m3))
def mtom3(self):
"""Takes the mass array for each bead and returns one with an element
for each degree of freedom.
Returns:
An array of size (nbeads,3*natoms), with each element corresponding
to the mass associated with the appropriate degree of freedom in q.
"""
m3 = np.zeros((self.nbeads,3*self.natoms),float)
m3[:,0:3*self.natoms:3] = self.m
m3[:,1:3*self.natoms:3] = m3[:,0:3*self.natoms:3]
m3[:,2:3*self.natoms:3] = m3[:,0:3*self.natoms:3]
return m3
def get_qc(self):
"""Gets the centroid coordinates."""
return np.dot(np.ones(self.nbeads,float),depstrip(self.q))/float(self.nbeads)
def get_pc(self):
"""Gets the centroid momenta."""
return np.dot(np.ones(self.nbeads,float),depstrip(self.p))/float(self.nbeads)
def kin_gather(self):
"""Gets the kinetic energy for all the replicas.
Returns:
A list of the kinetic energy for each system.
"""
return np.array([b.kin for b in self._blist])
def get_kin(self):
"""Gets the total kinetic energy of all the replicas.
Note that this does not correspond to the total kinetic energy estimate
for the system.
Returns:
The sum of the kinetic energy of each replica.
"""
return self.kins.sum()
def get_kstress(self):
"""Calculates the total kinetic stress tensor of all the replicas.
Note that this does not correspond to the quantum kinetic stress tensor
estimate for the system.
Returns:
The sum of the kinetic stress tensor of each replica.
"""
ks = np.zeros((3,3),float)
for b in range(self.nbeads):
ks += self[b].kstress
return ks
def get_vpath(self):
"""Calculates the spring potential between the replicas.
Note that this is actually the harmonic potential without being
multiplied by the factor omegan**2, which is only available in the
ensemble as the temperature is required to calculate it.
"""
epath = 0.0
q = depstrip(self.q)
m = depstrip(self.m3)[0]
for b in range(self.nbeads):
if b > 0:
dq = q[b,:] - q[b-1,:]
else:
dq = q[b,:] - q[self.nbeads-1,:]
epath += np.dot(dq, m*dq)
return epath*0.5
def get_fpath(self):
"""Calculates the spring force between the replicas.
Note that this is actually the harmonic force without being
multiplied by the factor omegan**2, which is only available in the
ensemble as the temperature is required to calculate it.
"""
nbeads = self.nbeads
natoms = self.natoms
f = np.zeros((nbeads,3*natoms),float)
q = depstrip(self.q)
m = depstrip(self.m3)[0]
for b in range(nbeads):
if b > 0:
dq = q[b,:] - q[b-1,:]
else:
dq = q[b,:] - q[self.nbeads-1,:]
dq *= m
f[b] -= dq
if b > 0:
f[b-1] += dq
else:
f[nbeads-1] += dq
return f
# A set of functions to access individual beads as Atoms objects
def __len__(self):
"""Length function.
This is called whenever the standard function len(beads) is used.
Returns:
The number of beads.
"""
return self.nbeads
def __getitem__(self,index):
"""Overwrites standard getting function.
This is called whenever the standard function beads[index] is used.
Returns an Atoms object with the appropriate position and momenta arrays.
Args:
index: The index of the replica of the system to be accessed.
Returns:
The replica of the system given by the index.
"""
return self._blist[index]
def __setitem__(self,index,value):
"""Overwrites standard setting function.
This is called whenever the standard function beads[index]=value is used.
Changes the position and momenta of the appropriate slice of the global
position and momentum arrays to those given by value.
Args:
index: The replica of the system to be changed.
value: The Atoms object that holds the new values.
"""
self._blist[index].p[:] = value.p
self._blist[index].q[:] = value.q
self._blist[index].m[:] = value.m
self._blist[index].names[:] = value.names

Event Timeline

c4science · Help