The BUG of enhance sampling perfomed by gpumd+plumed

[XIX-YANG]

The results of lammps +plumed cannot be replicated with gpumd+plumed.
This may be due to the incompleteness of the gpumd interface with plumed.
One possible reason is that, the thermostat and barostat need to be modified, as the bias is applied to the system.

[brucefan1983]

@initqp Do you think there is a potential bug for the PLUMED interface in GPUMD?

[initqp]

Maybe. However, without a working example, I will not be able to find out the true reason.

[initqp]

As far as I could tell, the “force” part of the interface should be bug free. Here is an example of calculating the rotational barrier of an ethane molecule in a vacuum, using ASE and GPUMD with PLUMED (you can download the example files from here):

run.in:

potential     ../nep.txt
velocity      300

plumed        ../plumed.inp 1 0
ensemble      nvt_lan 300 300 200
time_step     0.5
dump_thermo   10000
dump_position 5000
run           5000000

run.py:

import sys
import ase
import numpy as np
from ase import units
from ase.md import MDLogger
from ase.md.langevin import Langevin
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.io import read, write, Trajectory, NetCDFTrajectory
from ase.calculators.plumed import Plumed

import pynep.calculate


atoms = read('./topo.pdb')
calc = pynep.calculate.NEP('../nep.txt')
atoms.set_calculator(calc)

setup = open('../plumed.inp', 'r').read().splitlines()

plumed = Plumed(
    calc=atoms.calc,
    input=setup,
    timestep=0.5 * units.fs,
    atoms=atoms,
    kT=300 * units.kB,
    log='plumed.log'
)

MaxwellBoltzmannDistribution(atoms, temperature_K=360)

dyn = Langevin(
    atoms,
    timestep=0.5 * units.fs,
    temperature_K=300,
    friction=0.01 / ase.units.fs
)

dyn.attach(
    MDLogger(
        dyn,
        atoms,
        'ener.log',
        header=True,
        stress=False,
        peratom=False,
        mode='a'
    ),
    interval=10000
)

traj = NetCDFTrajectory('traj.nc', mode='a')
def write_frame():
    traj.write(dyn.atoms)
    print(dyn.get_number_of_steps(), flush=True)
dyn.attach(write_frame, interval=5000)

dyn.run(5000000)

Despite the sampling error, which can not be avoided, results from the two packages agree very well.

For the virial part, I am not 100 percent sure, since I do not usually work with isobaric processes.