PlumedMetadynamics

class PlumedMetadynamics(configuration, timestep, plumed_commands, logfile=None, temporary_plumed_file=None, charges=None)

Constructor for the PlumedMetadynamics hook class.

Parameters:
  • configuration (BulkConfiguration | MoleculeConfiguration) – The configuration on which the metadynamics MD simulation should be run.

  • timestep (PhysicalQuantity of type time) – The timestep used in the MD simulation.

  • plumed_commands (str) – A string containing the native PLUMED commands, that should be executed by PLUMED during the simulation.

  • logfile (str) – The name of the file, plumed writes its log information to.
    Default: ‘plumed.log’

  • temporary_plumed_file – The name of the temporary file, that is used to write the plumed commands to.
    Default: ‘plumed_tmp.dat’

  • charges (PhysicalQuantity of type charge) – The partial charges on the atoms.
    Default: 0.0*elementary_charge.

Usage Example

This example script runs a metadynamics simulation of a germanium interstitial in a silicon crystal. The bias potential acts on the cartesian X-coordinate of the germanium atom, gradually pushing it out of its original interstitial site and facilitating its diffusion in X-direction.

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------

initial_velocity = MaxwellBoltzmannDistribution(
    temperature=300.0*Kelvin,
    remove_center_of_mass_momentum=True
)

method = Langevin(
    time_step=1*femtoSecond,
    reservoir_temperature=300*Kelvin,
    friction=0.01*femtoSecond**-1,
    initial_velocity=initial_velocity,
    heating_rate=0*Kelvin/picoSecond,
)

# The bottom layer of atoms will be fixed to provide an absolute reference frame.
fix_atom_indices_0 = [0, 1, 4, 5, 32, 33, 36, 37]
constraints = [FixAtomConstraints(fix_atom_indices_0)]

# Set up the PLUMED command
plumed_command = """\
# Change the units to Ang, fs, eV
UNITS LENGTH=A TIME=fs ENERGY=96.48533645956869
# Use the cartesian position of the Ge interstitial atom as CV.
p: POSITION ATOM=65
# Put a metadynamics bias on the x component of the position CV p.
METAD ARG=p.x SIGMA=0.2 HEIGHT=0.05 PACE=50 LABEL=restraint
# Print the current position and bias potential value every 10 steps.
PRINT ARG=p.x,p.y,p.z,restraint.bias STRIDE=10 FILE=COLVAR
"""

# Set up the Plumed hook-class.
plumed_hook = PlumedMetadynamics(
    configuration=bulk_configuration,
    timestep=1.0*fs,
    plumed_commands=plumed_command,
    logfile='plumed_example.log',
)

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=constraints,
    trajectory_filename='si_ge_metadynamics_md.hdf5',
    steps=30000,
    log_interval=200,
    post_step_hook=plumed_hook,
    method=method
)

bulk_configuration = md_trajectory.lastImage()

Here, only the Molecular Dynamics block of the script is shown. The full script can be found found in the file script_si_ge_metadynamics.py

Notes

  • This is a hook class around the PLUMED code to run Metadynamics simulations [1]. As it modifies the forces and stress, it must be invoked as a post_step_hook in a MolecularDynamics simulation (more information about hook functions can be found in the MolecularDynamics manual).

  • Metadynamics is a powerful technique for enhancing sampling in molecular dynamics simulations and reconstructing the free-energy surface as a function of few selected degrees of freedom, often referred to as collective variables (CVs) [2].

  • Metadynamics adds an adaptive, history-dependent bias potential, which gradually fills up free energy minima until the system can effortlessly escape these minima and continue exploring the remaining phase space. Once converged, the negative bias potential can be used as an estimate of the free energy profile along the chosen CVs.

  • To run a Metadynamics simulation you need to specify the native PLUMED command as a string in the python script, which is then passed to the PlumedMetadynamics hook object. See the PLUMED manual for a list of all avaliable CVs and commands.

  • Note that PLUMED internally uses the units nm, ps, kJ/mol, meaning that all parameters in the plumed command need to be given in those units. By adding the line.

    UNITS LENGTH=A TIME=fs ENERGY=96.48533645956869
    

    to the plumed commands, you can switch to units Å, fs, and eV, which may be convenient for QuantumATK.

  • Moreover, note that the atom indices in PLUMED start at 1, whereas in QuantumATK they start at 0.

  • The bias potential is stored in a file (default filename HILLS). The free energy profile can be extracted and plotted from that file by invoking the script analyze_hills.py. The script may need to be adapted for the specific application at hand.

  • It is also possible to run an equilibrium simulation without metadynamics bias and just monitor the values of the selected CVs. This can be achieved by invoking the PRINT command.

  • You should cite Ref. [1] in all publications using PLUMED.