PLUMEDを使ったmetadynamics simulationをMatlantisで実行する方法 – Matlantis株式会社

この記事ではPLUMEDを使ったmetadynamics simulationをMatlantisで実行する方法について解説します。

PLUMEDをMatlantisにインストールする方法についてはこちらの記事を参考にしてください。

2024年7月現在、MatlantisではLAMMPSとASE (=> 3.23)について、PLUMEDの連携が可能です。

Matlantis上でLAMMPSを使ってmetadynamicsを実行する方法

Matlantis上でLAMMPSを使ってmetadynamicsを実行する方法については、こちらのMatlantis Contrib Githubで公開されているnotebookを参考に実行してください。

Matlantis上でASEを使ってmetadynamisを実行する方法

PLUMEDがインストールされているディレクトリにパスを通します。

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します。

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

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)

結果を保存するディレクトリを生成します。日付+時間でフォルダ名を生成しています。

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)

シミュレーションしたい構造を読み込み、viewerで表示します。ここでは各原子のindexも表示されるようにしています。

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

Collective Variables (CV)を設定します。設定可能なCVについてはPLUMEDの公式documentを確認してください。https://www.plumed.org/doc-v2.9/user-doc/html/_colvar.html
例えば、C-N bond形成反応を起こしたい場合は、C原子、N原子それぞれのindexを上記5のviewerで確認してPLUMEDのsettingに `"d1: DISTANCE ATOMS=3,31"` のように入力します。
注意点: aseのatoms objectのindexは0から、PLUMEDのindexは1から始まりますので、PLUMEDのsettingに入力する際は1を足したindex値を入力してください。
DISTANCE以外にもCVはありますのでPLUMEDの公式documentを確認してください。

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",
]

6のsettingをplumed calculatorに渡してmolecular dynamicsを実行することでmetadynamicsの計算が進みます。
```
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)
```
metadynamicsの計算中にHILLSとCOLVARSというファイルが生成し、計算が進むごとに更新されます。HILLSはmetadynamics中に追加されたポテンシャルに関する情報、COLVARSはmetadynamics中のCVと外部バイアスに関する情報が入っており、FESはこれらのファイルから計算することが可能です。その解析についてはこちらのMatlantis Contrib Githubで公開されているnotebookのCalculating FES by sum_hills以降を参考にしてください。

Matlantis上でLAMMPSを使ってmetadynamicsを実行する方法

Matlantis上でASEを使ってmetadynamisを実行する方法

関連記事