Monte Carlo Moves¶
A simulation can have an arbitrary number of MC moves operating on molecules, atoms, the volume, or
any other parameter affecting the system energy. Moves are specified in the moves
section
at the top level input. For example:
moves:
- moltransrot: { molecule: water, dp: 2.0, repeat: N,
dprot: 1.0, dir: [1,1,0] }
- volume: { dV: 0.01 }
- ...
random: { seed: hardware }
The pseudo-random number engine used for MC moves can be seeded in three ways,
seed |
Description |
---|---|
fixed |
Deterministic (default if random is absent) |
hardware |
Non-deterministric seed |
engine state | A previously saved stae |
The last option is used to restore the state of the engine as saved along with normal simulation
output as a string containing a lenghty list of numbers.
If initialization from a previously saved state fails – this may happen if
generated on another operating system – a warning is issued and the seed
falls back to fixed
.
Translation and Rotation¶
The following moves are for translation and rotation of atoms, molecules, or clusters.
The dir
keyword restricts translational directions which by default is set to
[1,1,1]
, meaning translation by a unit vector, randomly picked on a sphere, and scaled by
a random number in the interval [0, dp]
.
If dir=[1,1,0]
the unit vector is instead picked on a circle (here x, y)
and if dir=[0,0,1]
on a line (here $z$).
Molecular¶
moltransrot |
Description |
---|---|
molecule |
Molecule name to operate on |
dir=[1,1,1] |
Translational directions |
dirrot=[0,0,0] |
Predefined axis of rotation |
dp |
Translational displacement parameter |
dprot |
Rotational displacement parameter (radians) |
repeat=N |
Number of repeats per MC sweep. N equals $N_{molid}$ times. |
This will simultaneously translate and rotate a molecular group by the following operation
$$ \textbf{r}^N_{trial} = \mbox{Rot}(\textbf{r}^N) + \delta $$
where $\mbox{Rot}$ rotates dprot
$\cdot \left (\zeta-\frac{1}{2} \right )$ radians around a random unit vector
emanating from the mass center,
$\zeta$ is a random number in the interval $[0,1[$, and
$\delta$ is a random unit vector scaled by a random number in the interval [0,dp]
.
A predefined axis of rotation can be specified as dirrot
. For example, setting dirrot
to [1,0,0], [0,1,0] or [0,0,1]
results in rotations about the $x-$, $y-$, and $z-$axis, respectively.
Upon MC movement, the mean squared displacement will be tracked.
Atomic¶
transrot |
Description |
---|---|
molecule |
Molecule name to operate on |
dir=[1,1,1] |
Translational directions |
energy_resolution |
If set to a non-zero value (kT), an energy histogram will be generated. |
As moltransrot
but instead of operating on the molecular mass center, this translates
and rotates individual atoms in the group. The repeat is set to the number of atoms in the specified group and the
displacement parameters dp
and dprot
for the individual atoms are taken from
the atom properties defined in the topology.
Atomic rotation affects only anisotropic particles such as dipoles, spherocylinders, quadrupoles etc.
An energy histogram of each participating species will be written to disk if the energy_resolution
keyword is set. The value (in kT) specifies the resolution of the histogram binning. The analysis is
essentially for free as the energies are already known from the move.
Cluster Move¶
cluster |
Description |
---|---|
molecules |
Array of molecule names; [*] selects all |
threshold |
Mass-center threshold for forming a cluster (number or object) |
dir=[1,1,1] |
Directions to translate |
dirrot=[0,0,0] |
Predefined axis of rotation. If zero, a random unit vector is generated for each move event |
dprot |
Rotational displacement (radians) |
dp |
Translational displacement (Å) |
single_layer=false |
If true , stop cluster-growth after one layer around centered molecule (experimental) |
satellites |
Subset of molecules that cannot be cluster centers |
com_shape=true |
Use mass centers for shape analysis instead of particle positions (affects analysis only) |
analysis |
See below |
This will attempt to rotate and translate clusters of molecular molecules
defined by a distance threshold
between mass centers.
The threshold
can be specified as a single distance or as a complete list of combinations, see example below.
For simulations where small molecules cluster around large macro-molecules, it can be useful to use the satellites
keyword which denotes a list of molecules that can be part of a cluster, but cannot be the cluster nucleus or
starting point. All molecules listed in satellites
must be part of molecules
.
A predefined axis of rotation can be specified as dirrot
. For example, setting dirrot
to [1,0,0], [0,1,0] or [0,0,1]
results in rotations about the $x-$, $y-$, or $z-$axis, respectively.
The move is associated with bias, such that the cluster size and composition remain unaltered. If a cluster is larger than half the simulation box length, only translation will be attempted.
Example:
cluster:
molecules: [protein, cations]
satellites: [cations]
threshold:
protein protein: 25
protein cations: 15
cations cations: 0
dp: 3
dprot: 1
analysis: {file: cluster_shape.dat.gz}
Using analysis
the move also reports on the average cluster size; cluster size distribution; and
relative shape anisotropy, $\kappa^2$ as a function of cluster size. If all particles in the cluster
are distributed on a sphere then $\kappa^2=0$, while if on a straight line, $\kappa^2=1$. See the
polymershape
analysis for more information.
analysis |
Description |
---|---|
com=true |
Use molecule mass center instead of particle positions for shape anisotropy |
file |
If given save shape properties for each sample event |
save_pqr |
If set to true, save PQR files of cluster, one for each size |
interval=10 |
Interval between samples |
Smarter Monte Carlo¶
Preferential selction of particles can be enabled via the region
keyword which instructs
some moves to pick particles or groups preferentially from a given region. As described
in doi:10/frvx8j a bias is introduced which is automatically
accounted for. The preference for sampling inside the region is controlled by p
which
can be regarded as an outside update probability.
If $p=1$ no preferential sampling is performed, whereas if
$p<1$, more sampling is done inside the region.
For example:
- moltransrot:
...
...
region:
policy: ellipsoid
parallel_radius: 5.0
perpendicular_radius: 4.0
index1: 10
index2: 12
p: 0.2
The available regions are:
Ellipsoid¶
The connection vector between two (moving) reference particles defines an ellipsoid centered at the midpoint between the reference particles. The reference particle separation is unimportant, only the direction is used.
policy=ellipsoid |
Description |
---|---|
p |
Number (0,1] where a lower number means higher regional sampling |
index1 |
Index of first reference particle |
index2 |
Index of second reference particle |
parallel_radius |
Radius parallel to axis connecting the two references |
perpendicular_radius |
Radius perpendicular to axis connecting the two references |
group_com=false |
Test only mass center of molecular groups |
Within Molecule (experimental)¶
Samples from within a threshold from a molecule type. This can be useful to e.g.
preferentially sample solvent molecules around dilute solute molecules.
The com
keyword is available if the selected molecule
has a well-defined mass-center,
i.e. if is_atomic=false
.
It is also possible to use only the mass center for the moved groups by setting group_com
.
policy=within_molid |
Description |
---|---|
p |
Number (0,1] where a lower number means higher regional sampling |
molecule |
Name of molecule to search around |
threshold |
Distance threshold to any particle in molecule |
com=false |
Use threshold with respect to mass-center of molecule |
group_com=false |
Test with respect to mass center of molecular groups |
Internal Degrees of Freedom¶
Charge Move¶
charge |
Description |
---|---|
index |
Atom index to operate on |
dq |
Charge displacement |
quadratic=true |
Displace linearly along q^2 instead of q |
This performs a fractional charge move on a specific atom. The charge displacement can be performed linerly along $q$ or linearly along $q^2$. For the latter the following bias energy will be added to ensure uniform sampling of $q$,
$$ u = k_BT\ln \left ( \left | q^{\prime} / q \right |\right ) $$
Limitations:
This move changes the particle charge and therefore cannot be used with
splined pair-potentials where the initial charges from are read from atomlist
.
Instead, use a hard-coded variant like nonbonded_coulomblj
etc.
Conformational Swap¶
conformationswap |
Description |
---|---|
molecule |
Molecule name to operate on |
repeat=N |
Number of repeats per MC sweep |
keeppos=False |
Keep original positions of traj |
copy_policy=all |
What to copy from library. See table below. |
This will swap between different molecular conformations
as defined in the Molecule Properties with traj
and trajweight
If defined, the weight
distribution is respected, otherwise all conformations
have equal intrinsic weight. Upon insertion, the new conformation
is randomly oriented and placed on top of the mass-center of
an exising molecule. That is, there is no mass center movement.
If keeppos
is activated, the raw coordinates from the conformation
is used, i.e. no rotation and no mass-center overlay.
By default all information from the conformation is copied (copy_policy=all
), including charges and particle type.
To for example copy only positions, use copy_policy=positions
. This can be useful when using speciation moves.
copy_policy |
What is copied |
---|---|
all |
All particle properties |
positions |
Positions, only |
charges |
Charges, only |
patches |
Spherocylinder patch and length, but keep directions |
Pivot¶
pivot |
Description |
---|---|
molecule |
Molecule name to operate on |
dprot |
Rotational displacement |
repeat=N |
Number of repeats per MC sweep |
skiplarge=true |
Skip too large molecules |
Performs a rotation around a random, harmonic bond vector in molecule
, moving all atoms
either before or after the bond with equal probability. Current implementation assumes
unbranched chains with all atoms as links, i.e., no side chains are present.
For long polymers (compared to the box size), a large displacement parameter may cause
problems with mass center calculation in periodic systems.
This can be caught with the sanity
analysis and should it occur, try one of the following:
- enable
skiplarge
- decrease
dprot
- increase the simulation container.
The first option will simply reject troublesome configurations and the final output contains information of the skipped fraction. Skipping is unphysical so make sure the skipped fraction is small.
The default value of repeat
is the number of harmonic bonds in the molecule
(multiplied by the number of molecules).
Limitations: Chain bonds have to be ordered sequentially in the topology.
Crankshaft¶
crankshaft |
Description |
---|---|
molecule |
Molecule name to operate on |
dprot |
Rotational displacement |
repeat=N |
Number of repeats per MC sweep |
skiplarge=true |
Skip too large molecules |
joint_max= $\infty$ |
Maximum number of bonds between randomly selected atoms |
Performs a rotation of a chain segment between two randomly selected atoms in the molecule
.
The default value of repeat
is the number of atoms in the molecule
minus two
(multiplied by the number of molecules).
Parallel Tempering¶
temper |
Description |
---|---|
format=xyzqi |
Particle properties to copy between replicas (xyzqi , xyzq , xyz ) |
volume_scale=isotropic |
How to apply exchanged volumes: z , xy , isotropic , isochoric |
nstep=1 |
Number of sweeps between samples. |
partner_policy=oddeven |
Policy used to create partner pairs (currently only oddeven ) |
We consider an extended ensemble, consisting of n sub-systems or replicas, each in a distinct thermodynamic state (different Hamiltonians) and with the total energy
$$ U = \sum_i^n\mathcal{H}_i(\mathcal{R}_i) $$
The parallel tempering move performs a swap move where coordinate spaces (positions, volume) between random, neighboring sub-systems, i and j, are exchanged,
$$ \mathcal{R}_i^{\prime} = \mathcal{R}_j \quad \text{and} \quad \mathcal{R}_j^{\prime} = \mathcal{R}_i $$
and the energy change of the extended ensemble, $\Delta U_{i\leftrightarrow j}$, is used in the Metropolis acceptance criteria.
Parallel tempering requires compilation with MPI and the number
of replicas, n, exactly matches the number of processes. Each
replica prefixes input and output files with mpi0.
, mpi1.
,
etc. and only exchange between neighboring processes is performed.
The move is is performed exactly every nstep
Monte Carlo cycle.
By default, particle positions (xyz
), charge (q
), and atom id (i
) are exchanged
between replicas and can be controlled with format
.
Support for fluctuating number of particles, i.e. grand canonical moves is currently untested and should be considered experimental.
Volume Move¶
volume |
Description |
---|---|
dV |
Volume displacement parameter |
repeat=1 |
Number of repeats per MC sweep. |
method=isotropic |
Scaling method: z , xy , isotropic , isochoric |
Performs a random walk in logarithmic volume,
$$ V^{\prime} = e^{\ln V + \left (\zeta-\frac{1}{2} \right )\cdot dV } $$
and scales:
- molecular mass centers
- positions of free atoms (groups with
atomic=true
)
by $(V^{\prime}/V)^{1/3}$.
This is typically used for the $NPT$ ensemble, and for this an additional pressure term should be added to the Hamiltonian.
In the case of isochoric
scaling, the total volume is kept constant and dV
refers to an area change and reported output
statistics on volume should be regarded as area.
The table below explains the scaling behavior in different geometries:
method |
Geometry | Description |
---|---|---|
isotropic |
cuboid |
Scales x, y, z |
isotropic |
cylinder |
Scales radius |
isotropic |
sphere |
Scales radius |
z |
cuboid |
Scales z, xy untouched. |
xy |
cuboid |
Scales xy, z untouched. |
isochoric |
cuboid |
Scales xy/z, const. volume |
For cuboidal geometries, the scaling in each of the specified dimensions is $(V^{\prime}/V)^{1/d}$,
where $d=3$ for isotropic
, $d=2$ for xy
, and $d=1$ for z
.
Warning: Untested for cylinders, slits.
Reactive Canonical Monte Carlo¶
The speciation move handles density fluctuations and particle transformations and is the main move
for particle insertion, deletion, and swapping used in (semi)-grand canonical ensembles.
A reaction from reactionlist
is randomly picked from the topology and is either propagated forward or backward.
In Faunus, the total number of atoms and molecules is constant, but these can be either active or inactive.
Deleting a molecule simply deactivates it, while insertion vice versa activates an inactive molecule.
Thus, it is important that the capacity or reservoir of particles (active plus inactive) is
sufficiently large to allow for fluctuations.
This is ensured using insertmolecules
(see Topology).
A runtime warning will be given, should you run low on particles.
Besides deleting/inserting molecules (mono- or polyatomic), the speciation move performs reactions involving a
single-atom ID transformation (e.g., acid-base reactions).
In this case, an particle of type A (part of a mono- or polyatomic molecule) is randomly picked from the system
and all its properties, except its position, are replaced with those of an atom of type B.
Such ID transormations can also involve the addition/deletion of molecules or implicit atoms.
For a reaction
$$
\sum_i \nu_i M_i = 0
$$
where $M_i$ is the chemical symbol and $\nu_i$ is the stoichiometric coefficient of species $i$ (positive for products and negative for reagents),
the contribution of a speciation move to the energy change is
$$
\beta \Delta U = - \sum_i \ln{ \left ( \frac{ N_i! }{(N_i+\nu_i)!} V^{\nu_i} \right ) } - \ln{ \prod_i a_i^{\nu_i} },
$$
where $N_i$ is the number of particles of species $i$ in the current state and $a_i$ is the activity of species $i$.
For more information, see the Topology section and doi:10/fqcpg3.
rcmc |
Description |
---|---|
repeat=1 |
Average number of moves per sweep |
Replay¶
replay |
Description |
---|---|
file |
Trajectory file to read (xtc) |
Use next frame of the recorded trajectory as a move. The move is always unconditionally accepted, hence it may be used to replay a simulation, e.g., for analysis. Currently only Gromacs compressed trajectory file format (XTC) is supported. Note that total number of steps (macro × micro) should correspond to the number of frames in the trajectory.