Different RDF Plot for DFT and DPMD

[AnuragKr]

Hi,

Objective -- Comparing DPMD accuracy with DFT
Dataset -- examples/water/data
Training Script -- examples/water/se_e2_a_mixed_prec
LAMMPS Code --

units metal
atom_style atomic
boundary p p p

read_data water.lmp
mass 1 16
mass 2 2

pair_style deepmd ../se_e2_a_mixed_prec/graph-compress.pb
pair_coeff * *

timestep 0.0005

group Oatoms type 1
group Hatoms type 2

fix 1 all npt temp 300 300 0.05 iso 1.0 1.0 0.5

thermo_style custom step temp ke pe etotal press vol
thermo 100
thermo_modify format line "%d %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e %.6e"

dump my_dump all custom 100 h2o.xyz id type x y z vx vy vz
dump_modify my_dump sort id
dump my_dump2 all atom 100 water.lammps-dump-text

compute myRDF all rdf 100 1 1 1 2 2 2 
fix 2 all ave/time 100 1 100 c_myRDF[*] file h2o.rdf mode vector ave running

compute msd_O Oatoms msd
compute msd_H Hatoms msd
fix store_msd_O all vector 10 c_msd_O[4]
fix store_msd_H all vector 10 c_msd_H[4]
variable fitslope_O equal slope(f_store_msd_O)/6/(10*dt)
variable fitslope_H equal slope(f_store_msd_H)/6/(10*dt)
fix 3 all ave/time 100 1 100 c_msd_O[4] c_msd_H[4] v_fitslope_O v_fitslope_H file diffusion.txt

run 100000

DFT Code --

import numpy as np
import os

def rdf(sys,
        sel_type = [None, None], 
        max_r = 5, 
        nbins = 100) :
    """
    compute the rdf of a system
    Parameters
    ----------
    sys : System or LabeledSystem
        The dpdata system
    sel_type: list
        List of size 2. The first element specifies the type of the first atom, 
        while the second element specifies the type of the second atom. 
        Both elements can be ints or list of ints. 
        If the element is None, all types are specified.
        Examples are sel_type = [0, 0], sel_type = [0, [0, 1]] or sel_type = [0, None]
    max_r: float
        Maximal range of rdf calculation
    nbins: int
        Number of bins for rdf calculation
    Returns
    -------
    xx: np.array
        The lattice of r
    rdf: np.array
        The value of rdf at r
    coord: np.array
        The coordination number up to r
    """
    return compute_rdf(sys['cells'], sys['coords'], sys['atom_types'], 
                       sel_type = sel_type, 
                       max_r = max_r, 
                       nbins = nbins)

def compute_rdf(box,
                posis, 
                atype, 
                sel_type = [None, None], 
                max_r = 5, 
                nbins = 100) :
    nframes = box.shape[0]
    xx = None
    all_rdf = []
    all_cod = []
    for ii in range(nframes):
        xx, rdf, cod = _compute_rdf_1frame(box[ii], posis[ii], atype, sel_type, max_r, nbins)
        all_rdf.append(rdf)
        all_cod.append(cod)
    all_rdf = np.array(all_rdf).reshape([nframes, -1])
    all_cod = np.array(all_cod).reshape([nframes, -1])
    all_rdf = np.average(all_rdf, axis = 0)
    all_cod = np.average(all_cod, axis = 0)
    return xx, all_rdf, all_cod

def _compute_rdf_1frame(box, 
                        posis, 
                        atype, 
                        sel_type = [None, None], 
                        max_r = 5, 
                        nbins = 100) :
    all_types = list(set(list(np.sort(atype))))
    if sel_type[0] is None:
        sel_type[0] = all_types
    if sel_type[1] is None:
        sel_type[1] = all_types
    if type(sel_type[0]) is not list:
        sel_type[0] = [sel_type[0]]
    if type(sel_type[1]) is not list:
        sel_type[1] = [sel_type[1]]
    natoms = len(posis)
    from ase import Atoms
    import ase.neighborlist    
    atoms = Atoms(positions=posis, cell=box, pbc=[1,1,1])
    nlist = ase.neighborlist.NeighborList(max_r, self_interaction=False, bothways=True, primitive=ase.neighborlist.NewPrimitiveNeighborList)
    nlist.update(atoms)
    stat = np.zeros(nbins)
    stat_acc = np.zeros(nbins)
    hh = max_r / float(nbins)
    for ii in range(natoms) :
        # atom "0"
        if atype[ii] in sel_type[0]:
            indices, offsets = nlist.get_neighbors(ii)
            for jj, os in zip(indices, offsets):
                # atom "1"
                if atype[jj] in sel_type[1]:
                    posi_jj = atoms.positions[jj] + np.dot(os, atoms.get_cell())
                    diff = posi_jj - atoms.positions[ii]
                    dr = np.linalg.norm(diff)                    
                    # if (np.linalg.norm(diff- diff_1)) > 1e-12 :
                    #     raise RuntimeError
                    si = int(dr / hh)
                    if si < nbins:
                        stat[si] += 1
    # count the number of atom1
    c0 = 0
    for ii in sel_type[0]:
        c0 += np.sum(atype == ii)
    # count the number of atom1
    c1 = 0
    for ii in sel_type[1]:
        c1 += np.sum(atype == ii)
    rho1 = c1 / np.linalg.det(box)
    # compute coordination number
    for ii in range(1, nbins):
        stat_acc[ii] = stat_acc[ii-1] + stat[ii-1]
    stat_acc = stat_acc / c0
    # compute rdf
    for ii in range(nbins):
        vol = 4./3. * np.pi * ( ((ii+1)*hh) ** 3 - ((ii)*hh) ** 3 )
        rho = stat[ii] / vol
        stat[ii] = rho / rho1 / c0
    xx = np.arange(0, max_r-1e-12, hh)
    return xx, stat, stat_acc

def get_systems(path):
    return [dpdata.LabeledSystem(os.path.join(path, s), fmt='deepmd/npy') for s in os.listdir(path)]

if __name__ == '__main__':
    import dpdata
    systems = get_systems('../../data')
    weights = []
    rdf_list = []
    for sys in systems:
        coord_dim = sys['coords'].shape
        frame_num = coord_dim[0]
        atom_num = coord_dim[1]
        weights.append(frame_num*atom_num)
        xx,rdf_values,_ = rdf(sys, sel_type = [1,1], max_r = 6, nbins = 100)
        rdf_list.append(rdf_values)
    rdf = np.average(rdf_list, axis=0, weights=weights)
    res = np.concatenate([xx, rdf]).reshape([2, -1])
    np.savetxt('rdf_HH.out', res.T)

RDF Plot --

Graph Vizualization Code -- link

Another RDF plot follows the same above process with the descriptor se_e2_a.
Dataset -- H2O, PBE0TS
Graph --

Is the DFT code above enough for creating accurate RDF compared to AIMD software?

Please let me know where I am going wrong.

[wanghan-iapcm]

It is warned that the example water data (in folder examples/water/data) is of very limited amount, is provided only for testing purposes, and should not be used to train a production model (i.e. the rdf may off!).

The rdf should converge to 1 as r goes to infinity, but it is not the case in your plot of deepmd rdfs.

[AnuragKr]

@benshpiro I have used a 95%-5% data split for training and validation/testing. In some discussions on GitHub author suggested trying batch size 1 for accuracy improvement, You also can try different batch sizes to see when you are getting better results.

[benshpiro]

I seem to have the same problem you had at some point for the RDF of O-O, for some reason, the DeepMD has non-physical peaks at very close distances. How did you manage to solve this problem? When you are running the MD run, what time step did you use?

[AnuragKr]

For Water --
Timestep - 0.48fs
Run - 625000
Simulation - 300ps

Few Checkpoints --

1. Model should have good accuracy (Not much difference between train and test data).
2. Run the MD simulation in precisely the same thermodynamic condition at which data is created.
3. Have the same number of MD simulation input molecules as it is used for training.
4. Check the structure of MD simulation input data.
5. Also How do you create RDF for DFT data and code for creating RDF graphs.

[benshpiro]

Hi Anurag,
Thanks for the info. I tried training with the same parameters and I'm still getting it wrong. Regarding your points:

1. I have good accuracy
2. simulation is exactly the same conditions to the training data
3. exact same number of molecules
4. input structure is correct
5. For the RDF I am using the "compute rdf" command in lammps (at this point I am training the deepmd model on data frames of the system that I produce in a Lammps simulation, so the RDF plot is produced the same for both methods in the comparison)

How many frames do you use for training? I tried, 6000.

[AnuragKr]

I have used out of 10000, either 9500 for training or 9000 for training and the remaining for validation/test. Iteration will be 1000000.
For water lammps simulation use the deepmd example input data that would be enough for RDF graph creation.
Also, could you try batch sizes 1 or 4?
Please share the graph and code to understand your issue better.