# 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)

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:

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)

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.
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:

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)

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.
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)

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.
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)

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.
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)

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.
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)

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.
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)

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)

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)

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)

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)
get_value(iterative)
copy()
class yaff.sampling.iterative.PosStateItem
get_value(iterative)
class yaff.sampling.iterative.DipoleStateItem
get_value(iterative)
class yaff.sampling.iterative.DipoleVelStateItem
get_value(iterative)
class yaff.sampling.iterative.VolumeStateItem
get_value(iterative)
class yaff.sampling.iterative.CellStateItem
get_value(iterative)
class yaff.sampling.iterative.EPotContribStateItem

Keeps track of all the contributions to the potential energy.

get_value(iterative)
iter_attrs(iterative)
class yaff.sampling.iterative.EpotBondsStateItem(do_ei=False)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)
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)

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)

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)

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)

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)

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)

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)
get_value(iterative)
copy()

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

Geometry/Cell optimization

class yaff.sampling.opt.OptScreenLog(start=0, step=1)
__call__(iterative)
class yaff.sampling.opt.BaseOptimizer(dof, state=None, hooks=None, counter0=0)

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)
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)

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:

1. 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.
2. 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.
3. 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.
4. 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.
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.
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.
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)
__call__(iterative)
class yaff.sampling.trajectory.RefTrajectory(ff, fn_traj, state=None, hooks=None, counter0=0)

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
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)

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)

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)

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)

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)