GPUMD icon indicating copy to clipboard operation
GPUMD copied to clipboard

Possible BUG of enhance sampling perfomed by gpumd+plumed

Open XIX-YANG opened this issue 9 months ago • 4 comments

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.

XIX-YANG avatar May 10 '24 08:05 XIX-YANG

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

brucefan1983 avatar May 10 '24 18:05 brucefan1983

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

initqp avatar May 10 '24 22:05 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):

1

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.

initqp avatar May 11 '24 23:05 initqp

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. Dear Yang: I am using gpumd patched with plumed, too. Can you give a working example to make it more clearifier? I mainly using it to do metad recently, and at least the traj seems well and works well, but I haven't compare the result with other md software due to knowledge limit. Thanks! And also thanks the conmment made by the contributor who patched plumed, your work had underreated.

heroiciota avatar May 27 '24 10:05 heroiciota