3. `yaff.sampling` – Phase-space sampling¶

3.1. `yaff.sampling.dof` – Abstraction layer for degrees of freedom¶

Abstraction layer for degrees of freedom

All these classes are called DOF classes, because they specify a set of degrees of freedom. These DOF classes are used for geometry/cell optimization and harmonic approximations.

class yaff.sampling.dof.DOF(ff)¶

Bases: object

Arguments:

ff: A force field object

ndof¶

reset()¶

check_delta(x=None, eps=0.0001, zero=None)¶: Test the analytical derivatives

log()¶

class yaff.sampling.dof.CartesianDOF(ff, gpos_rms=1e-05, dpos_rms=0.001, select=None)¶

Bases: yaff.sampling.dof.DOF

Cartesian degrees of freedom

This DOF is also applicable to periodic systems. Cell parameters are not modified when this DOF is used.

Arguments:

ff: A force field object.

Optional arguments:

gpos_rms, dpos_rms

Thresholds that define the convergence. If all of the actual values drop below these thresholds, the minimizer stops.

For each rms threshold, a corresponding max threshold is included automatically. The maximum of the absolute value of a component should be smaller than 3/sqrt(N) times the rms threshold, where N is the number of degrees of freedom.

select

A selection of atoms whose degrees of freedom are included. If not list is given, all atomic coordinates are included.

Convergence conditions:

gpos_rms: The root-mean-square of the norm of the gradients of the atoms.
dpos_rms: The root-mean-square of the norm of the displacements of the atoms.

fun(x, do_gradient=False)¶

Computes the energy and optionally the gradient.

Arguments:

x: The degrees of freedom

Optional arguments:

do_gradient: When True, the gradient is also returned.

check_convergence()¶

class yaff.sampling.dof.BaseCellDOF(ff, gpos_rms=1e-05, dpos_rms=0.001, grvecs_rms=1e-05, drvecs_rms=0.001, do_frozen=False, freemask=None)¶

Bases: yaff.sampling.dof.DOF

Fractional coordinates and cell parameters

Several subclasses of BaseCellDOF are implemented below. Each one considers a specific representation and subset of the cell parameters.

The following variable names are consistently used (also in subclasses):

cellvars: An array with all variables for the cell (specific for ja BaseCellDOF subclass).
ncellvar: The number of cellvars (at most 9).
celldofs: A selection of the elements in cellvars, based on freemask.
ncelldof: The number of celldofs (less than or equal to ncellvar).
frac: Fractional coordinates.
x: All degrees of freedom, i.e. celldofs and frac (in that order, frac is optional).
The suffix 0: Used for initial values of something.

Arguments:

ff: A force field object.

Optional arguments:

gpos_rms, dpos_rms, grvecs_rms, drvecs_rms

Thresholds that define the convergence. If all of the actual values drop below these thresholds, the minimizer stops.

For each rms threshold, a corresponding max threshold is included automatically. The maximum of the absolute value of a component should be smaller than 3/sqrt(N) times the rms threshold, where N is the number of degrees of freedom.

do_frozen

When True, the fractional coordinates of the atoms are kept fixed.

freemask

When given, this must be an array of booleans indicating which cellvars are free. At least one cellvar must be free.

Convergence conditions:

gpos_rms: The root-mean-square of the norm of the gradients of the atoms.
dpos_rms: The root-mean-square of the norm of the displacements of the atoms.
grvecs_rms: The root-mean-square of the norm of the gradients of the cell vectors.
drvecs_rms: The root-mean-square of the norm of the displacements of the cell vectors.

ncellvar¶: The number of cellvars

ncelldof¶: The number of celldofs (free cellvars)

fun(x, do_gradient=False)¶

Computes the energy and optionally the gradient.

Arguments:

x: All degrees of freedom

Optional arguments:

do_gradient: When True, the gradient is also returned.

check_convergence()¶

log()¶

class yaff.sampling.dof.FullCellDOF(ff, gpos_rms=1e-05, dpos_rms=0.001, grvecs_rms=1e-05, drvecs_rms=0.001, do_frozen=False, freemask=None)¶

Bases: yaff.sampling.dof.BaseCellDOF

DOF that includes all 9 components of the cell vectors

The degrees of freedom are rescaled cell vectors ordered in one row:

3D periodic: [a_x/s, a_y/s, a_z/s, b_x/s, b_y/s, b_z/s, c_x/s, c_y/s, c_z/s] where s is the cube root of the initial cell volume such that the cell DOFs become dimensionless.
2D periodic: [a_x/s, a_y/s, a_z/s, b_x/s, b_y/s, b_z/s] where s is the square root of the initial cell surface such that the cell DOFs become dimensionless.
1D periodic: [a_x/s, a_y/s, a_z/s] where s is the length of the initial cell vector such that the cell DOFs become dimensionless.

Arguments:

ff: A force field object.

Optional arguments:

gpos_rms, dpos_rms, grvecs_rms, drvecs_rms

Thresholds that define the convergence. If all of the actual values drop below these thresholds, the minimizer stops.

For each rms threshold, a corresponding max threshold is included automatically. The maximum of the absolute value of a component should be smaller than 3/sqrt(N) times the rms threshold, where N is the number of degrees of freedom.

do_frozen

When True, the fractional coordinates of the atoms are kept fixed.

freemask

When given, this must be an array of booleans indicating which cellvars are free. At least one cellvar must be free.

Convergence conditions:

gpos_rms: The root-mean-square of the norm of the gradients of the atoms.
dpos_rms: The root-mean-square of the norm of the displacements of the atoms.
grvecs_rms: The root-mean-square of the norm of the gradients of the cell vectors.
drvecs_rms: The root-mean-square of the norm of the displacements of the cell vectors.

class yaff.sampling.dof.StrainCellDOF(ff, gpos_rms=1e-05, dpos_rms=0.001, grvecs_rms=1e-05, drvecs_rms=0.001, do_frozen=False, freemask=None)¶

Bases: yaff.sampling.dof.BaseCellDOF

Eliminates rotations of the unit cell. thus six cell parameters are free.

The degrees of freedom are coefficients in symmetrix matrix transformation, A, that is applied to the initial cell vectors.

3D periodic: [A_00, A_11, A_22, 2*A_12, 2*A_20, 2*A_01]
2D periodic: [A_00, A_11, 2*A_01]
1D periodic: [A_00]

Why does this work? Let R be the array with cell vectors as rows. It can always be written as a product,

R = R_0.F,

where F is an arbitrary 3x3 matrix. Application of SVD to the matrix F :Yields: R = R_0.U.S.V^T = R_0.U.V^T.V.S.V^T

Then W=U.V^T is a orthonormal matrix and A=V.S.V^T is a symmetric matrix. The orthonormal matrix W is merely a rotation of the cell vectors, which can be omitted as the internal energy is invariant under such rotations. The symmetric matrix actually deforms the cell and is the part of interest.

Arguments:

ff: A force field object.

Optional arguments:

gpos_rms, dpos_rms, grvecs_rms, drvecs_rms

Thresholds that define the convergence. If all of the actual values drop below these thresholds, the minimizer stops.

For each rms threshold, a corresponding max threshold is included automatically. The maximum of the absolute value of a component should be smaller than 3/sqrt(N) times the rms threshold, where N is the number of degrees of freedom.

do_frozen

When True, the fractional coordinates of the atoms are kept fixed.

freemask

When given, this must be an array of booleans indicating which cellvars are free. At least one cellvar must be free.

Convergence conditions:

gpos_rms: The root-mean-square of the norm of the gradients of the atoms.
dpos_rms: The root-mean-square of the norm of the displacements of the atoms.
grvecs_rms: The root-mean-square of the norm of the gradients of the cell vectors.
drvecs_rms: The root-mean-square of the norm of the displacements of the cell vectors.

class yaff.sampling.dof.AnisoCellDOF(ff, gpos_rms=1e-05, dpos_rms=0.001, grvecs_rms=1e-05, drvecs_rms=0.001, do_frozen=False, freemask=None)¶

Bases: yaff.sampling.dof.BaseCellDOF

Only the lengths of the cell vectors are free. angles are fixed.

The degrees of freedom are dimensionless scale factors for the cell lengths, using the initial cell vectors as the reference point. (This is one DOF per periodic dimension.)

Arguments:

ff: A force field object.

Optional arguments:

gpos_rms, dpos_rms, grvecs_rms, drvecs_rms

Thresholds that define the convergence. If all of the actual values drop below these thresholds, the minimizer stops.

For each rms threshold, a corresponding max threshold is included automatically. The maximum of the absolute value of a component should be smaller than 3/sqrt(N) times the rms threshold, where N is the number of degrees of freedom.

do_frozen

When True, the fractional coordinates of the atoms are kept fixed.

freemask

When given, this must be an array of booleans indicating which cellvars are free. At least one cellvar must be free.

Convergence conditions:

gpos_rms: The root-mean-square of the norm of the gradients of the atoms.
dpos_rms: The root-mean-square of the norm of the displacements of the atoms.
grvecs_rms: The root-mean-square of the norm of the gradients of the cell vectors.
drvecs_rms: The root-mean-square of the norm of the displacements of the cell vectors.

class yaff.sampling.dof.IsoCellDOF(ff, gpos_rms=1e-05, dpos_rms=0.001, grvecs_rms=1e-05, drvecs_rms=0.001, do_frozen=False, freemask=None)¶

Bases: yaff.sampling.dof.BaseCellDOF

The cell is only allowed to undergo isotropic scaling

The only degree of freedom is an isotropic scaling factor, using the initial cell vectors as a reference.

Arguments:

ff: A force field object.

Optional arguments:

gpos_rms, dpos_rms, grvecs_rms, drvecs_rms

Thresholds that define the convergence. If all of the actual values drop below these thresholds, the minimizer stops.

For each rms threshold, a corresponding max threshold is included automatically. The maximum of the absolute value of a component should be smaller than 3/sqrt(N) times the rms threshold, where N is the number of degrees of freedom.

do_frozen

When True, the fractional coordinates of the atoms are kept fixed.

freemask

When given, this must be an array of booleans indicating which cellvars are free. At least one cellvar must be free.

Convergence conditions:

gpos_rms: The root-mean-square of the norm of the gradients of the atoms.
dpos_rms: The root-mean-square of the norm of the displacements of the atoms.
grvecs_rms: The root-mean-square of the norm of the gradients of the cell vectors.
drvecs_rms: The root-mean-square of the norm of the displacements of the cell vectors.

class yaff.sampling.dof.FixedBCDOF(ff, gpos_rms=1e-05, dpos_rms=0.001, grvecs_rms=1e-05, drvecs_rms=0.001, do_frozen=False, freemask=None)¶

Bases: yaff.sampling.dof.BaseCellDOF

A rectangular cell that can only stretch along one axis

This cell optimization constrains the cell in the y and z direction to the original values, but allows expansion and contraction in the x direction. The system should be rotated such that the initial cell vectors look like:

a = ( ax , 0  , 0  )
b = ( 0  , by , bz )
c = ( 0  , cy , cz )

During optimization, only ax will be allowed to change.

This type of constraint can be used when looking at a structure that is periodic only in one dimension, but you have to fake a 3D structure to be able to use Ewald summation

Arguments:

ff: A force field object.

Optional arguments:

gpos_rms, dpos_rms, grvecs_rms, drvecs_rms

Thresholds that define the convergence. If all of the actual values drop below these thresholds, the minimizer stops.

For each rms threshold, a corresponding max threshold is included automatically. The maximum of the absolute value of a component should be smaller than 3/sqrt(N) times the rms threshold, where N is the number of degrees of freedom.

do_frozen

When True, the fractional coordinates of the atoms are kept fixed.

freemask

When given, this must be an array of booleans indicating which cellvars are free. At least one cellvar must be free.

Convergence conditions:

gpos_rms: The root-mean-square of the norm of the gradients of the atoms.
dpos_rms: The root-mean-square of the norm of the displacements of the atoms.
grvecs_rms: The root-mean-square of the norm of the gradients of the cell vectors.
drvecs_rms: The root-mean-square of the norm of the displacements of the cell vectors.

class yaff.sampling.dof.FixedVolOrthoCellDOF(ff, volume=None, gpos_rms=1e-05, dpos_rms=0.001, grvecs_rms=1e-05, drvecs_rms=0.001, do_frozen=False, freemask=None)¶

Bases: yaff.sampling.dof.BaseCellDOF

Orthorombic cell optimizer with a fixed volume.

These constraints are implemented by using the following cell vectors:

a = (  s*a0*la  ,     0     ,      0         )
b = (     0     ,  s*b0*lb  ,      0         )
c = (     0     ,     0     ,  s*c0/(la*lb)  )

with s = (V/V0)^(1/3)

Optional arguments (in addition to those of BaseCellDOF):

volume: The desired volume of the cell. (When not given, the current volume of the system is not altered.)

3.2. `yaff.sampling.harmonic` – Harmonic models¶

Harmonic models

yaff.sampling.harmonic.estimate_hessian(dof, eps=0.0001)¶

Estimate the Hessian using the symmetric finite difference approximation.

Arguments:

dof: A DOF object

Optional arguments:

eps: The magnitude of the displacements

yaff.sampling.harmonic.estimate_cart_hessian(ff, eps=0.0001, select=None)¶

Estimate the Cartesian Hessian with symmetric finite differences.

Arguments:

ff: A force field object

Optional arguments:

eps: The magnitude of the Cartesian displacements
select: A selection of atoms for which the hessian must be computed. If not given, the entire hessian is computed.

yaff.sampling.harmonic.estimate_elastic(ff, eps=0.0001, do_frozen=False, ridge=0.0001)¶

Estimate the elastic constants using the symmetric finite difference approximation.

Arguments:

ff: A force field object

Optional arguments:

eps: The magnitude of the Cartesian displacements
do_frozen: By default this is False, which means that the changes in fractional atomic coordinates due to cell deformations are properly taken into account. When this is set to True, such displacements (other than uniform scaling) are ignored. The latter is much faster, but only correct for the simplest materials.
ridge: Threshold for the eigenvalues of the Cartesian Hessian. This only matters if do_frozen==False.

The elastic constants are second order derivatives of the strain energy density with respect to uniform deformations. At the molecular scale, uniform deformations can be describe by a linear transformation of the atoms and the cell parameters:

x' = e . x0

When x0 corresponds to the reference point for the second order expansion, small deviations from that reference point correspond to small deviations of e from the unit matrix. The strain energy density is nothing but the energy of the system divided by its volume, minus a the energy density of the fully relaxed system. The adjective strain refers to deviations from the relaxed reference. When computing the second order derivatives, this reference point can be ignored. Assuming the reference point is the relaxed system, the strain energy density can be approximated to second order as:

u = 1/2 \sum_ijkl e_ij c_ijkl e_kl

In this equation the matrix c_ijkl contains all the elastic constants. In principle, it contains 81 values, but due to symmetry considerations, only 21 of these values are independent. [The symmetry considerations are as follows: (i) it is sufficient to consider symmetric deformations e when excluding cell rotations, (ii) the index pairs (ij) and (kl) are interchangeable.] These 21 numbers are typically represented using the (inclusive) lower diagonal of a 6x6 matrix, using the following convention to map pair indexes (ij) and (kl) onto single indexes n and m:

(ij) or (kl)	n or m
11	1
22	2
22	3
23	4
31	5
12	6

(This table is also known as the compressed Voight notation.)

This routine returns the elastic constants in a symmetric 6x6 matrix, using the index conventions by Voight.

3.3. `yaff.sampling.io` – Trajectory writers¶

Trajectory writers

class yaff.sampling.io.HDF5Writer(f, start=0, step=1)¶

Bases: yaff.sampling.iterative.Hook

Argument:

f: A h5.File object to write the trajectory to.

Optional arguments:

start: The first iteration at which this hook should be called.
step: The hook will be called every step iterations.

__call__(iterative)¶

dump_system(system)¶

init_trajectory(iterative)¶

class yaff.sampling.io.XYZWriter(fn_xyz, select=None, start=0, step=1)¶

Bases: yaff.sampling.iterative.Hook

Argument:

fn_xyz: A filename to write the XYZ trajectory too.

Optional arguments:

select: A list of atom indexes that should be written to the trajectory output. If not given, all atoms are included.
start: The first iteration at which this hook should be called.
step: The hook will be called every step iterations.

__call__(iterative)¶

class yaff.sampling.io.RestartWriter(f, start=0, step=1000)¶

Bases: yaff.sampling.iterative.Hook

Argument:

f: A h5.File object to write the restart information to.

Optional arguments:

start: The first iteration at which this hook should be called.
step: The hook will be called every step iterations.

init_state(iterative)¶

__call__(iterative)¶

dump_system(system)¶

init_trajectory(iterative)¶

dump_restart(hook)¶

3.4. `yaff.sampling.iterative` – Base class for iterative algorithms¶

Base class for iterative algorithms

class yaff.sampling.iterative.Iterative(ff, state=None, hooks=None, counter0=0)¶

Bases: object

Arguments:

ff: The ForceField instance used in the iterative algorithm

Optional arguments:

state: A list with state items. State items are simple objects that take or derive a property from the current state of the iterative algorithm.
hooks: A function (or a list of functions) that is called after every iterative.
counter0: The counter value associated with the initial state.

default_state = []¶

log_name = 'ITER'¶

initialize()¶

call_hooks()¶

run(nstep=None)¶

propagate()¶

finalize()¶

class yaff.sampling.iterative.StateItem(key)¶

Bases: object

update(iterative)¶

get_value(iterative)¶

iter_attrs(iterative)¶

copy()¶

class yaff.sampling.iterative.AttributeStateItem(key)¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

copy()¶

class yaff.sampling.iterative.PosStateItem¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

class yaff.sampling.iterative.DipoleStateItem¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

class yaff.sampling.iterative.DipoleVelStateItem¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

class yaff.sampling.iterative.VolumeStateItem¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

class yaff.sampling.iterative.CellStateItem¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

class yaff.sampling.iterative.EPotContribStateItem¶

Bases: yaff.sampling.iterative.StateItem

Keeps track of all the contributions to the potential energy.

get_value(iterative)¶

iter_attrs(iterative)¶

class yaff.sampling.iterative.EpotBondsStateItem(do_ei=False)¶

Bases: yaff.sampling.iterative.StateItem

Keeps track of all the Valence Bond contributions to the potential energy.

Optional arguments:

do_ei: If True, the electrostatic contributions of the bonded atom pair will also be tracked.

get_value(iterative)¶

iter_atts(iterative)¶

class yaff.sampling.iterative.EpotBendsStateItem(do_ei=False)¶

Bases: yaff.sampling.iterative.StateItem

Keeps track of all the Valence Bend contributions to the potential energy

Optional arguments:

do_ei: If True, the electrostatic contributions of the 1-3 non-bonded atom pair will also be tracked.

get_value(iterative)¶

iter_atts(iterative)¶

class yaff.sampling.iterative.EpotDihedsStateItem(do_ei=False)¶

Bases: yaff.sampling.iterative.StateItem

Keeps track of all the Valence Dihedral contributions to the potential energy

Optional arguments:

do_ei: If True, the electrostatic contributions of the 1-4 non-bonded atom pair will also be tracked.

get_value(iterative)¶

iter_atts(iterative)¶

class yaff.sampling.iterative.Hook(start=0, step=1)¶

Bases: object

Optional arguments:

start: The first iteration at which this hook should be called.
step: The hook will be called every step iterations.

name = None¶

kind = None¶

method = None¶

expects_call(counter)¶

__call__(iterative)¶

3.5. `yaff.sampling.npt` – Barostats¶

Barostats

class yaff.sampling.npt.TBCombination(thermostat, barostat, start=0)¶

Bases: yaff.sampling.verlet.VerletHook

VerletHook combining an arbitrary Thermostat and Barostat instance, which ensures these instances are called in the correct succession, and possible coupling between both is handled correctly.

Arguments:

thermostat: A Thermostat instance
barostat: A Barostat instance

name = 'TBCombination'¶

init(iterative)¶

pre(iterative)¶

post(iterative)¶

expectscall(iterative, kind)¶

verify()¶

class yaff.sampling.npt.McDonaldBarostat(temp, press, start=0, step=1, amp=0.001)¶

Bases: yaff.sampling.verlet.VerletHook

Warning: this code is not fully tested yet!

Arguments:

temp: The average temperature of the NpT ensemble
press: The external pressure of the NpT ensemble

Optional arguments:

start: The first iteration at which this hook is called
step: The number of iterations between two subsequent calls to this hook.
amp: The amplitude of the changes in the logarithm of the volume.

name = 'McDonald'¶

kind = 'stochastic'¶

method = 'barostat'¶

init(iterative)¶

pre(iterative, chainvel0=None)¶

post(iterative, chainvel0=None)¶

class yaff.sampling.npt.BerendsenBarostat(ff, temp, press, start=0, step=1, timecon=41341.37333664682, beta=13445.401558486123, anisotropic=True, vol_constraint=False, restart=False)¶

Bases: yaff.sampling.verlet.VerletHook

This hook implements the Berendsen barostat. The equations are derived in:

Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690

Arguments:

ff: A ForceField instance.
temp: The temperature of thermostat.
press: The applied pressure for the barostat.

Optional arguments:

start: The step at which the thermostat becomes active.
timecon: The time constant of the Berendsen barostat.
beta: The isothermal compressibility, conventionally the compressibility of liquid water
anisotropic: Defines whether anisotropic cell fluctuations are allowed.
vol_constraint: Defines whether the volume is allowed to fluctuate.

name = 'Berendsen'¶

kind = 'deterministic'¶

method = 'barostat'¶

init(iterative)¶

pre(iterative, chainvel0=None)¶

post(iterative, chainvel0=None)¶

class yaff.sampling.npt.LangevinBarostat(ff, temp, press, start=0, step=1, timecon=41341.37333664682, anisotropic=True, vol_constraint=False)¶

Bases: yaff.sampling.verlet.VerletHook

This hook implements the Langevin barostat. The equations are derived in:

Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103, 4613-4621

Arguments:

ff: A ForceField instance.
temp: The temperature of thermostat.
press: The applied pressure for the barostat.

Optional arguments:

start: The step at which the barostat becomes active.
timecon: The time constant of the Langevin barostat.
anisotropic: Defines whether anisotropic cell fluctuations are allowed.
vol_constraint: Defines whether the volume is allowed to fluctuate.

name = 'Langevin'¶

kind = 'stochastic'¶

method = 'barostat'¶

init(iterative)¶

pre(iterative, chainvel0=None)¶

post(iterative, chainvel0=None)¶

baro(iterative, chainvel0)¶

getR()¶

class yaff.sampling.npt.MTKBarostat(ff, temp, press, start=0, step=1, timecon=41341.37333664682, anisotropic=True, vol_constraint=False, baro_thermo=None, vel_press0=None, restart=False)¶

Bases: yaff.sampling.verlet.VerletHook

This hook implements the Martyna-Tobias-Klein barostat. The equations are derived in:

Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177-4189.

The implementation (used here) of a symplectic integrator of this barostat is discussed in

Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117-1157.

Arguments:

ff: A ForceField instance.
temp: The temperature of thermostat.
press: The applied pressure for the barostat.

Optional arguments:

start: The step at which the barostat becomes active.
timecon: The time constant of the Martyna-Tobias-Klein barostat.
anisotropic: Defines whether anisotropic cell fluctuations are allowed.
vol_constraint: Defines whether the volume is allowed to fluctuate.
baro_thermo: NHCThermostat instance, coupled directly to the barostat
vel_press0: The initial barostat velocity tensor
restart: If true, the cell is not symmetrized initially

name = 'MTTK'¶

kind = 'deterministic'¶

method = 'barostat'¶

init(iterative)¶

pre(iterative, chainvel0=None)¶

post(iterative, chainvel0=None)¶

baro(iterative, chainvel0)¶

add_press_cont()¶

class yaff.sampling.npt.PRBarostat(ff, temp, press, start=0, step=1, timecon=41341.37333664682, anisotropic=True, vol_constraint=False, baro_thermo=None, vel_press0=None, restart=False)¶

Bases: yaff.sampling.verlet.VerletHook

This hook implements the Parrinello-Rahman barostat for finite strains. The equations are derived in:

Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182-7190.

Arguments:

ff: A ForceField instance.
temp: The temperature of thermostat.
press: The applied pressure for the barostat.

Optional arguments:

start: The step at which the barostat becomes active.
timecon: The time constant of the Tadmor-Miller barostat.
anisotropic: Defines whether anisotropic cell fluctuations are allowed.
vol_constraint: Defines whether the volume is allowed to fluctuate.
baro_thermo: NHCThermostat instance, coupled directly to the barostat
vel_press0: The initial barostat velocity tensor
restart: If true, the cell is not symmetrized initially

name = 'PR'¶

kind = 'deterministic'¶

method = 'barostat'¶

init(iterative)¶

pre(iterative, chainvel0=None)¶

post(iterative, chainvel0=None)¶

baro(iterative, chainvel0)¶

add_press_cont()¶

class yaff.sampling.npt.TadmorBarostat(ff, temp, press, start=0, step=1, timecon=41341.37333664682, anisotropic=True, vol_constraint=False, baro_thermo=None, vel_press0=None, restart=False)¶

Bases: yaff.sampling.verlet.VerletHook

This hook implements the Tadmor-Miller barostat for finite strains. The equations are derived in:

Tadmor, E. B.; Miller, R. E. Modeling Materials: Continuum, Atomistic and Multiscale Techniques 2011, 520-527.

Arguments:

ff: A ForceField instance.
temp: The temperature of thermostat.
press: The applied second Piola-Kirchhoff tensor for the barostat.

Optional arguments:

start: The step at which the barostat becomes active.
timecon: The time constant of the Tadmor-Miller barostat.
anisotropic: Defines whether anisotropic cell fluctuations are allowed.
vol_constraint: Defines whether the volume is allowed to fluctuate.
baro_thermo: NHCThermostat instance, coupled directly to the barostat
vel_press0: The initial barostat velocity tensor
restart: If true, the cell is not symmetrized initially

name = 'Tadmor'¶

kind = 'deterministic'¶

method = 'barostat'¶

init(iterative)¶

pre(iterative, chainvel0=None)¶

post(iterative, chainvel0=None)¶

baro(iterative, chainvel0)¶

add_press_cont()¶

class yaff.sampling.npt.MTKAttributeStateItem(attr)¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

copy()¶

3.6. `yaff.sampling.nvt` – Thermostats¶

Thermostats

class yaff.sampling.nvt.AndersenThermostat(temp, start=0, step=1, select=None, annealing=1.0)¶

Bases: yaff.sampling.verlet.VerletHook

This is an implementation of the Andersen thermostat. The method is described in:

Andersen, H. C. J. Chem. Phys. 1980, 72, 2384-2393.

Arguments:

temp: The average temperature of the NVT ensemble

Optional arguments:

start: The first iteration at which this hook is called
step: The number of iterations between two subsequent calls to this hook.
select: An array of atom indexes to indicate which atoms controlled by the thermostat.
annealing: After every call to this hook, the temperature is multiplied with this annealing factor. This effectively cools down the system.

name = 'Andersen'¶

kind = 'stochastic'¶

method = 'thermostat'¶

init(iterative)¶

pre(iterative, G1_add=None)¶

post(iterative, G1_add=None)¶

class yaff.sampling.nvt.BerendsenThermostat(temp, start=0, timecon=4134.137333664683, restart=False)¶

Bases: yaff.sampling.verlet.VerletHook

This is an implementation of the Berendsen thermostat. The algorithm is described in:

Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690

Arguments:

temp: The temperature of thermostat.

Optional arguments:

start

The step at which the thermostat becomes active.

timecon

The time constant of the Berendsen thermostat.

restart: Indicates whether the initalisation should be carried out.

name = 'Berendsen'¶

kind = 'deterministic'¶

method = 'thermostat'¶

init(iterative)¶

pre(iterative, G1_add=None)¶

post(iterative, G1_add=None)¶

class yaff.sampling.nvt.LangevinThermostat(temp, start=0, timecon=4134.137333664683)¶

Bases: yaff.sampling.verlet.VerletHook

This is an implementation of the Langevin thermostat. The algorithm is described in:

Bussi, G.; Parrinello, M. Phys. Rev. E 2007, 75, 056707

Arguments:

temp: The temperature of thermostat.

Optional arguments:

start: The step at which the thermostat becomes active.
timecon: The time constant of the Langevin thermostat.

name = 'Langevin'¶

kind = 'stochastic'¶

method = 'thermostat'¶

init(iterative)¶

pre(iterative, G1_add=None)¶

post(iterative, G1_add=None)¶

thermo(iterative)¶

class yaff.sampling.nvt.CSVRThermostat(temp, start=0, timecon=4134.137333664683)¶

Bases: yaff.sampling.verlet.VerletHook

This is an implementation of the CSVR thermostat. The equations are derived in:

Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126, 014101

The implementation (used here) is derived in

Bussi, G.; Parrinello, M. Comput. Phys. Commun. 2008, 179, 26-29

Arguments:

temp: The temperature of thermostat.

Optional arguments:

start: The step at which the thermostat becomes active.
timecon: The time constant of the CSVR thermostat.

name = 'CSVR'¶

kind = 'stochastic'¶

method = 'thermostat'¶

init(iterative)¶

pre(iterative, G1_add=None)¶

post(iterative, G1_add=None)¶

class yaff.sampling.nvt.GLEThermostat(temp, a_p, c_p=None, start=0)¶

Bases: yaff.sampling.verlet.VerletHook

This hook implements the coloured noise thermostat. The equations are derived in:

Ceriotti, M.; Bussi, G.; Parrinello, M J. Chem. Theory Comput. 2010, 6, 1170-1180.

Arguments:

temp: The temperature of thermostat.
a_p: Square drift matrix, with elements fitted to the specific problem.

Optional arguments:

c_p: Square static covariance matrix. In equilibrium, its elements are fixed. For non-equilibrium dynamics, its elements should be fitted.
start: The step at which the thermostat becomes active.

name = 'GLE'¶

kind = 'stochastic'¶

method = 'thermostat'¶

init(iterative)¶

pre(iterative, G1_add=None)¶

post(iterative, G1_add=None)¶

thermo(iterative)¶

class yaff.sampling.nvt.NHCThermostat(temp, start=0, timecon=4134.137333664683, chainlength=3, chain_pos0=None, chain_vel0=None, restart=False)¶

Bases: yaff.sampling.verlet.VerletHook

This hook implements the Nose-Hoover chain thermostat. The equations are derived in:

Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635-2643.

The implementation (used here) of a symplectic integrator of the Nose-Hoover chain thermostat is discussed in:

Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Mol. Phys. 1996, 87, 1117-1157.

Arguments:

temp: The temperature of thermostat.

Optional arguments:

start: The step at which the thermostat becomes active.
timecon: The time constant of the Nose-Hoover thermostat.
chainlength: The number of beads in the Nose-Hoover chain.
chain_pos0: The initial thermostat chain positions
chain_vel0: The initial thermostat chain velocities
restart: Indicates whether the initalisation should be carried out

name = 'NHC'¶

kind = 'deterministic'¶

method = 'thermostat'¶

init(iterative)¶

pre(iterative, G1_add=None)¶

post(iterative, G1_add=None)¶

class yaff.sampling.nvt.NHCAttributeStateItem(attr)¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

copy()¶

3.7. `yaff.sampling.opt` – Geometry/Cell optimization¶

Geometry/Cell optimization

class yaff.sampling.opt.OptScreenLog(start=0, step=1)¶

Bases: yaff.sampling.iterative.Hook

__call__(iterative)¶

class yaff.sampling.opt.BaseOptimizer(dof, state=None, hooks=None, counter0=0)¶

Bases: yaff.sampling.iterative.Iterative

Arguments:

dof: A specification of the degrees of freedom. The convergence criteria are also part of this argument. This must be a DOF instance.

Optional arguments:

state: A list with state items. State items are simple objects that take or derive a property from the current state of the iterative algorithm.
hooks: A function (or a list of functions) that is called after every iterative.
counter0: The counter value associated with the initial state.

default_state = [<yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.PosStateItem object>, <yaff.sampling.iterative.DipoleStateItem object>, <yaff.sampling.iterative.VolumeStateItem object>, <yaff.sampling.iterative.CellStateItem object>, <yaff.sampling.iterative.EPotContribStateItem object>]¶

log_name = 'XXOPT'¶

fun(x, do_gradient=False)¶

initialize()¶

propagate()¶

finalize()¶

class yaff.sampling.opt.CGOptimizer(dof, state=None, hooks=None, counter0=0)¶

Bases: yaff.sampling.opt.BaseOptimizer

log_name = 'CGOPT'¶

initialize()¶

propagate()¶

class yaff.sampling.opt.BFGSHessianModel(ndof, hessian0=None)¶

Bases: yaff.sampling.opt.HessianModel

Arguments:

ndof: The number of degrees of freedom

Optional arguments:

hessian0: An initial guess for the hessian

update(dx, dg)¶

class yaff.sampling.opt.SR1HessianModel(ndof, hessian0=None)¶

Bases: yaff.sampling.opt.HessianModel

Arguments:

ndof: The number of degrees of freedom

Optional arguments:

hessian0: An initial guess for the hessian

update(dx, dg)¶

class yaff.sampling.opt.QNOptimizer(dof, state=None, hooks=None, counter0=0, trust_radius=1.0, small_radius=1e-05, too_small_radius=1e-10, hessian0=None)¶

Bases: yaff.sampling.opt.BaseOptimizer

A Quasi-newton optimizer

This is just a basic implementation of the algorithm, but it has the potential to become more advanced and efficient. The following improvements will be made when time permits:

Support for non-linear constraints. This should be relatively easy. We need a routine that can bring the unknowns back to the constraints, and a routine to solve a constrained second order problem with linear equality/inequality constraints. These should be methods of an object that is an attribute of the dof object, which is need to give the constraint code access to the Cartesian coordinates. In the code below, some comments are added to mark where the constraint methods should be called.
The Hessian updates and the diagonalization are currently very slow for big systems. This can be fixed with a rank-1 update algorithm for the spectral decomposition.
The optimizer would become much more efficient if redundant coordinates were used. This can be implemented efficiently by using the same machinery as the constraint code, but using the dlist and iclist concepts for the sake of efficiency.
It is in practice not needed to keep track of the full Hessian. The L-BFGS algorithm is a nice method to obtain a linear memory usage and computational cost. However, L-BFGS is not compatible with the trust radius used in this class, while we want to keep the trust radius for the sake of efficiency, robustness and support for constraints. Using the rank-1 updates mentioned above, it should be relatively easy to keep track of the decomposition of a subspace of the Hessian. This subspace can be defined as the basis of the last N rank-1 updates. Simple assumptions about the remainder of the spectrum should be sufficient to keep the algorithm efficient.

Arguments:

dof: A specification of the degrees of freedom. The convergence criteria are also part of this argument. This must be a DOF instance.

Optional arguments:

state: A list with state items. State items are simple objects that take or derive a property from the current state of the iterative algorithm.
hooks: A function (or a list of functions) that is called after every iterative.
counter0: The counter value associated with the initial state.
trust_radius: The initial value for the trust radius. It is adapted by the algorithm after every step. The adapted trust radius is never allowed to increase above this initial value.
small_radius: If the trust radius goes below this limit, the decrease in energy is no longer essential. Instead a decrease in the norm of the gradient is used to accept/reject a step.
too_small_radius: If the trust radius becomes smaller than this parameter, the optimizer gives up. Insanely small trust radii are typical for potential energy surfaces that are not entirely smooth.
hessian0: An initial guess for the Hessian

log_name = 'QNOPT'¶

initialize()¶

propagate()¶

make_step()¶

yaff.sampling.opt.solve_trust_radius(grad, evals, radius, threshold=1e-05)¶: Find a step in eigen space with the given radius

3.8. `yaff.sampling.trajectory` – Computations on a reference trajectory¶

Computations on a reference trajectory

class yaff.sampling.trajectory.TrajScreenLog(start=0, step=1)¶

Bases: yaff.sampling.iterative.Hook

__call__(iterative)¶

class yaff.sampling.trajectory.RefTrajectory(ff, fn_traj, state=None, hooks=None, counter0=0)¶

Bases: yaff.sampling.iterative.Iterative

Arguments:

ff: A ForceField instance

fn_traj

A hdf5 file name containing the trajectory

Optional arguments:

state: A list with state items. State items are simple objects that take or derive a property from the current state of the iterative algorithm.
hooks: A function (or a list of functions) that is called after every iterative.
counter0: The counter value associated with the initial state.

default_state = [<yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.PosStateItem object>, <yaff.sampling.iterative.DipoleStateItem object>, <yaff.sampling.iterative.VolumeStateItem object>, <yaff.sampling.iterative.CellStateItem object>, <yaff.sampling.iterative.EPotContribStateItem object>]¶

log_name = 'TRAJEC'¶

initialize()¶

propagate()¶

finalize()¶

3.9. `yaff.sampling.utils` – Auxiliary routines for initial velocities¶

Auxiliary routines for initial velocities

yaff.sampling.utils.get_random_vel(temp0, scalevel0, masses, select=None)¶

Generate random atomic velocities using a Maxwell-Boltzmann distribution

Arguments:

temp0: The temperature for the Maxwell-Boltzmann distribution.
scalevel0: When set to True, the velocities are rescaled such that the instantaneous temperature coincides with temp0.
masses: An (N,) array with atomic masses.

Optional arguments:

select: When given, this must be an array of integer indexes. Only for these atoms (masses) initial velocities will be generated.

Returns: An (N, 3) array with random velocities. When the select option is used, the shape of the results is (M, 3), where M is the length of the select array.

yaff.sampling.utils.remove_com_moment(vel, masses)¶

Zero the linear center-of-mass momentum.

Arguments:

vel: An (N, 3) array with atomic velocities. This array is modified in-place.
masses: An (N,) array with atomic masses

The zero linear COM momentum is achieved by subtracting translational rigid body motion from the atomic velocities.

yaff.sampling.utils.remove_angular_moment(pos, vel, masses)¶

Zero the global angular momentum.

Arguments:

pos: An (N, 3) array with atomic positions. This array is not modified.
vel: An (N, 3) array with atomic velocities. This array is modified in-place.
masses: An (N,) array with atomic masses

The zero angular momentum is achieved by subtracting angular rigid body motion from the atomic velocities. (The angular momentum is measured with respect to the center of mass to avoid that this routine reintroduces a linear COM velocity. This is also beneficial for the numerical stability.)

yaff.sampling.utils.clean_momenta(pos, vel, masses, cell)¶

Remove any relevant external momenta

Arguments:

pos: An (N, 3) array with atomic positions. This array is not modified.
vel: An (N, 3) array with atomic velocities. This array is modified in-place.
masses: An (N,) array with atomic masses
cell: A Cell instance describing the periodic boundary conditions.

yaff.sampling.utils.angular_moment(pos, vel, masses)¶

Compute the angular moment of a set of point particles

Arguments:

pos: An (N, 3) array with atomic positions.
vel: An (N, 3) array with atomic velocities.
masses: An (N,) array with atomic masses.

Returns: a (3,) array with the angular momentum vector.

yaff.sampling.utils.get_ndof_internal_md(natom, nper)¶

Return the effective number of internal degrees of freedom for MD simulations

Arguments:

natom: The number of atoms
nper: The number of periodic boundary conditions (0 for isolated systems)

yaff.sampling.utils.cell_symmetrize(ff, vector_list=None, tensor_list=None)¶

Symmetrizes the unit cell tensor, and updates the position vectors

Arguments:

ff: A ForceField instance

Optional arguments:

vector_list: A list of numpy vectors which should be transformed under the symmetrization. Note that the positions are already transformed automatically
tensor_list: A list of numpy tensors of rank 2 which should be transformed under the symmetrization.

yaff.sampling.utils.get_random_vel_press(mass, temp)¶

Generates symmetric tensor of barostat velocities

Arguments:

mass: The Barostat mass
temp: The temperature at which the velocities are selected

yaff.sampling.utils.get_ndof_baro(dim, anisotropic, vol_constraint)¶

Calculates the number of degrees of freedom associated with the cell fluctuation

Arguments:

dim: The dimension of the system
anisotropic: Boolean value determining whether anisotropic cell fluctuations are allowed
vol_constraint: Boolean value determining whether the cell volume can change

yaff.sampling.utils.stabilized_cholesky_decomp(mat)¶: Do LDL^T and transform to MM^T with negative diagonal entries of D put equal to zero Assume mat is square and symmetric (but not necessarily positive definite).

3.10. `yaff.sampling.verlet` – Generic Verlet integrator¶

Generic Verlet integrator

class yaff.sampling.verlet.TemperatureStateItem¶

Bases: yaff.sampling.iterative.StateItem

get_value(iterative)¶

iter_attrs(iterative)¶

class yaff.sampling.verlet.VerletIntegrator(ff, timestep=None, state=None, hooks=None, vel0=None, temp0=300, scalevel0=True, time0=None, ndof=None, counter0=None, restart_h5=None)¶

Bases: yaff.sampling.iterative.Iterative

Arguments:

ff: A ForceField instance

Optional arguments:

timestep: The integration time step (in atomic units)
state: A list with state items. State items are simple objects that take or derive a property from the current state of the iterative algorithm.
hooks: A function (or a list of functions) that is called after every iterative.
vel0: An array with initial velocities. If not given, random velocities are sampled from the Maxwell-Boltzmann distribution corresponding to the optional arguments temp0 and scalevel0
temp0: The (initial) temperature for the random initial velocities
scalevel0: If True (the default), the random velocities are rescaled such that the instantaneous temperature coincides with temp0.
time0: The time associated with the initial state.
ndof: When given, this option overrides the number of degrees of freedom determined from internal heuristics. When ndof is not given, its default value depends on the thermostat used. In most cases it is 3*natom, except for the NHC thermostat where the number if internal degrees of freedom is counted. The ndof attribute is used to derive the temperature from the kinetic energy.
counter0: The counter value associated with the initial state.
restart_h5: HDF5 object containing the restart information

default_state = [<yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.PosStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.verlet.TemperatureStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.AttributeStateItem object>, <yaff.sampling.iterative.DipoleStateItem object>, <yaff.sampling.iterative.DipoleVelStateItem object>, <yaff.sampling.iterative.VolumeStateItem object>, <yaff.sampling.iterative.CellStateItem object>, <yaff.sampling.iterative.EPotContribStateItem object>]¶

log_name = 'VERLET'¶

initialize()¶

propagate()¶

compute_properties(restart_h5=None)¶

finalize()¶

call_verlet_hooks(kind)¶

class yaff.sampling.verlet.VerletHook(start=0, step=1)¶

Bases: yaff.sampling.iterative.Hook

Specialized Verlet hook.

This is mainly used for the implementation of thermostats and barostats.

Optional arguments:

start: The first iteration at which this hook should be called.
step: The hook will be called every step iterations.

__call__(iterative)¶

init(iterative)¶

pre(iterative)¶

post(iterative)¶

class yaff.sampling.verlet.VerletScreenLog(start=0, step=1)¶

Bases: yaff.sampling.iterative.Hook

A screen logger for the Verlet algorithm

__call__(iterative)¶

class yaff.sampling.verlet.ConsErrTracker(restart_h5=None)¶

Bases: object

A class that tracks the errors on the conserved quantity. Given its superior numerical accuracy, the algorithm below is used to calculate the running average. Its properties are discussed in Donald Knuth’s Art of Computer Programming, vol. 2, p. 232, 3rd edition.

update(ekin, econs)¶

get()¶

class yaff.sampling.verlet.KineticAnnealing(annealing=0.99999, select=None, start=0, step=1)¶

Bases: yaff.sampling.verlet.VerletHook

This annealing hook is designed to be used with a plain Verlet integrator. At every call, the velocities are rescaled with the annealing parameter.

Arguments:

annealing: After every call to this hook, the temperature is multiplied with this annealing factor. This effectively cools down the system.
select: An array mask or a list of indexes to indicate which atomic velocities should be annealed.
start: The first iteration at which this hook is called
step: The number of iterations between two subsequent calls to this hook.

init(iterative)¶

pre(iterative)¶

post(iterative)¶

3. `yaff.sampling` – Phase-space sampling¶

3.1. `yaff.sampling.dof` – Abstraction layer for degrees of freedom¶

3.2. `yaff.sampling.harmonic` – Harmonic models¶

3.3. `yaff.sampling.io` – Trajectory writers¶

3.4. `yaff.sampling.iterative` – Base class for iterative algorithms¶

3.5. `yaff.sampling.npt` – Barostats¶

3.6. `yaff.sampling.nvt` – Thermostats¶

3.7. `yaff.sampling.opt` – Geometry/Cell optimization¶

3.8. `yaff.sampling.trajectory` – Computations on a reference trajectory¶

3.9. `yaff.sampling.utils` – Auxiliary routines for initial velocities¶

3.10. `yaff.sampling.verlet` – Generic Verlet integrator¶

Table Of Contents

Previous topic

Next topic

This Page

3. yaff.sampling – Phase-space sampling¶

3.1. yaff.sampling.dof – Abstraction layer for degrees of freedom¶

3.2. yaff.sampling.harmonic – Harmonic models¶

3.3. yaff.sampling.io – Trajectory writers¶

3.4. yaff.sampling.iterative – Base class for iterative algorithms¶

3.5. yaff.sampling.npt – Barostats¶

3.6. yaff.sampling.nvt – Thermostats¶

3.7. yaff.sampling.opt – Geometry/Cell optimization¶

3.8. yaff.sampling.trajectory – Computations on a reference trajectory¶

3.9. yaff.sampling.utils – Auxiliary routines for initial velocities¶

3.10. yaff.sampling.verlet – Generic Verlet integrator¶

3. `yaff.sampling` – Phase-space sampling¶

3.1. `yaff.sampling.dof` – Abstraction layer for degrees of freedom¶

3.2. `yaff.sampling.harmonic` – Harmonic models¶

3.3. `yaff.sampling.io` – Trajectory writers¶

3.4. `yaff.sampling.iterative` – Base class for iterative algorithms¶

3.5. `yaff.sampling.npt` – Barostats¶

3.6. `yaff.sampling.nvt` – Thermostats¶

3.7. `yaff.sampling.opt` – Geometry/Cell optimization¶

3.8. `yaff.sampling.trajectory` – Computations on a reference trajectory¶

3.9. `yaff.sampling.utils` – Auxiliary routines for initial velocities¶

3.10. `yaff.sampling.verlet` – Generic Verlet integrator¶