2-RDM AIMD (cumulant Δ)

NVE molecular dynamics for water at 300 K from δγ and the cumulant Δ = Γ − γ ∧ γ.

Two QMModels are trained: one for delta_gamma and one for gamma2cum. QMLCalculator is called with method='gamma2cum' and qmmodel2 holding the cumulant model.

Point dbfile at an HDF5 database that stores both delta_gamma and gamma2cum. Production examples used 20 000 steps of 0.5 fs.

import numpy as np
from ase.build import molecule
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, force_temperature
from ase.md.verlet import VelocityVerlet
from ase import units
from sklearn.kernel_ridge import KernelRidge

from qmlearn.io.model import db2qmmodel
from qmlearn.api.api4ase import QMLCalculator

dbfile = "train.hdf5"

T = 300
nsteps = 200
timestep = 0.5 * units.fs
np.random.seed(8888)

models = {
    "gamma2cum": KernelRidge(alpha=0.0, kernel="rbf"),
    "delta_gamma": KernelRidge(alpha=0.0, kernel="rbf"),
}

qmmodel = db2qmmodel(
    dbfile,
    names="*",
    mmodels=models,
    target="delta_gamma",
    method="delta_gamma",
    purify_gamma=False,
)
qmmodel2 = db2qmmodel(
    dbfile,
    names="*",
    mmodels=models,
    target="gamma2cum",
    method="gamma2cum",
    purify_gamma=False,
)

atoms = molecule("H2O")
atoms.calc = QMLCalculator(
    qmmodel=qmmodel,
    qmmodel2=qmmodel2,
    method="gamma2cum",
    properties=("energy",),
)

MaxwellBoltzmannDistribution(atoms, temperature_K=T, force_temp=True)
p = atoms.get_momenta()
p -= p.sum(axis=0) / len(atoms)
atoms.set_momenta(p)
force_temperature(atoms, T)

dyn = VelocityVerlet(
    atoms,
    timestep=timestep,
    trajectory="md_nve_cumulant.traj",
    logfile="md_nve_cumulant.log",
)
dyn.run(nsteps)
print("Finished", nsteps, "steps; energy =", atoms.get_potential_energy())