How to Run Metadynamics Simulations Using PLUMED in Matlantis – Matlantis Corporation

In this article, we will explain how to run metadynamics simulations using PLUMED in Matlantis.

Please refer to this article for instructions on how to install PLUMED in Matlantis.

As of July 2024, Matlantis can integrate with PLUMED for LAMMPS and ASE (>= 3.23).

How to Run Metadynamics Using LAMMPS in Matlantis

For instructions on how to run metadynamics using LAMMPS in Matlantis, please refer to the notebook available on the Matlantis Contrib GitHub.

How to Run Metadynamics Using ASE in Matlantis

Add the directory where PLUMED is installed to your PATH.

import os
import sys
os.environ["PLUMED_KERNEL"]="/home/jovyan/plumed_local/lib/libplumedKernel.so"
os.environ["L"]="/home/jovyan/plumed_local/lib/libplumedKernel.so"
os.environ["PLUMED_TYPESAFE_IGNORE"]="yes"
sys.path.append('/home/jovyan/plumed_local/')

Import the necessary modules.

from ase.calculators.plumed import Plumed
from ase import units
from ase.md.langevin import Langevin
from ase.atoms import Atoms
from ase.io import read, write
from ase.optimize import LBFGS, FIRE
from ase.constraints import FixAtoms, FixedPlane, FixBondLength, ExpCellFilter, UnitCellFilter, FixedPlane, Hookean, FixSymmetry
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import math, os, glob, optuna, random, scipy, shutil, datetime

Prepare the PFP calculator.

from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorElementStatus
calc_mode="CRYSTAL_U0_PLUS_D3"
model_version="v7.0.0"
estimator = Estimator(calc_mode=calc_mode, model_version=model_version)
calculator = ASECalculator(estimator)

Create a directory to save the results. The folder name is generated based on the date and time.

now = datetime.datetime.now(datetime.timezone(datetime.timedelta(hours=9)))
yyyymmddhhmm = now.strftime('%Y%m%d%H%M')
outpath = "output_plumed/"+yyyymmddhhmm+"/"
os.makedirs(outpath, exist_ok=True)

Load the structure you want to simulate and display it in the viewer. Here, the index of each atom is also shown.

atoms = read("filename.xyz")
v = view_ngl(atoms, representations=["ball+stick"])
v.view.add_label(
    color="black", labelType="text",
    labelText=[atoms[i].symbol + str(i) for i in range(atoms.get_global_number_of_atoms())],
    zOffset=1.5, attachment='middle_center', radius=1.0)
v

Set the Collective Variables (CV). For the available CV options, please refer to the official PLUMED documentation: https://www.plumed.org/doc-v2.9/user-doc/html/_colvar.html
For example, if you want to induce a C-N bond formation reaction, check the indices of the C and N atoms using the viewer mentioned in step 5, and enter the settings into PLUMED like `"d1: DISTANCE ATOMS=3,31"`.
Note: In an ASE atoms object, the indices start from 0, while in PLUMED, the indices start from 1. Therefore, when entering the indices into the PLUMED settings, add 1 to the indices.
There are other CV types besides DISTANCE, so please refer to the official PLUMED documentation for more information.

setting = [f"UNITS LENGTH=A TIME=ps ENERGY={units.mol/units.kJ}",
　　　　　　　　　　　　　　　　　"d1: DISTANCE ATOMS=3,31",
　　　　　　　　　　　　　　　　　"uwall: UPPER_WALLS ARG=d1 AT=4.0 KAPPA=10.0 EXP=2 EPS=1 OFFSET=0",
　　　　　　　　　　　　　　　　　f"METAD ARG=d1 SIGMA=0.2 HEIGHT=0.05 PACE=200 LABEL=restraint FILE={outpath}HILLS",
　　　　　　　　　　　　　　　　　f"PRINT ARG=d1,restraint.bias STRIDE=10 FILE={outpath}COLVAR",
]

Use the "setting" defined in step 6 as arguments for the PLUMED calculator and run the molecular dynamics. This will advance the metadynamics calculation.
```
timestep = 1*units.fs
ps = 1000 * units.fs
temp = 300

atoms.calc = Plumed(calc=calculator,
            input=setting,
            timestep=timestep,
            atoms=atoms,
            kT=1
)

dyn = Langevin(atoms, 
            timestep, 
            temperature_K=temp, 
            friction=0.002,
            trajectory=outpath+'MTdyn.traj', 
            logfile=outpath+"MTdyn.log",
            loginterval=100,
)
dyn.run(2_000_000)
```
During the metadynamics calculation, files named HILLS and COLVARS are generated and updated as the calculation progresses. HILLS contains information about the potentials added during the metadynamics, while COLVARS contains information about the CVs and external biases. The FES can be calculated from these files. For analysis, please refer to the "Calculating FES by sum_hills" section in the notebook available on the Matlantis Contrib GitHub.

Related to

How to Run Metadynamics Using LAMMPS in Matlantis

How to Run Metadynamics Using ASE in Matlantis

Related articles