1-RDM AIMD

NVE molecular dynamics for water at 300 K using a machine-learned 1-RDM.

The model learns δγ = γ − γ_HF from the external potential, then a second-learn step maps the predicted γ to energy and forces (d_energy, d_forces). The ASE calculator uses method='gamma'.

Point dbfile at an HDF5 database that contains delta_gamma (and energy/force targets for the second-learn models). Production examples used 20 000 steps of 0.5 fs; this notebook defaults to a short run.

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 sklearn.linear_model import LinearRegression

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

# Path to a QMLearn HDF5 training database
dbfile = "train.hdf5"

T = 300
nsteps = 200          # production examples used 20000
timestep = 0.5 * units.fs
np.random.seed(8888)

models = {
    "delta_gamma": KernelRidge(alpha=0.0, kernel="rbf"),
    "d_energy": LinearRegression(),
    "d_forces": KernelRidge(alpha=0.0, kernel="rbf", gamma=1e-1),
}
second_learn = {
    "energy": "d_energy",
    "forces": "d_forces",
}

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

atoms = molecule("H2O")
atoms.calc = QMLCalculator(
    qmmodel=qmmodel,
    second_learn=second_learn,
    method="gamma",
    properties=("energy", "forces"),
)

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_1rdm.traj",
    logfile="md_nve_1rdm.log",
)
dyn.run(nsteps)
print("Finished", nsteps, "steps; energy =", atoms.get_potential_energy())