RDF calculation for MultiSystem

[AnuragKr]

Hello,
I hope you are doing well. In a deepmd-kit repo, data is given with the same water system but in different folders. I want to calculate RDF which will consider data_0,data_1,data_2.

Code I am referring to for calculating RDF -- https://github.com/deepmodeling/dpdata/blob/master/dpdata/md/rdf.py

As we can see in the above code it accepts System or LabeledSystem.What are the changes need to make in the above if I give MultiSystem?

I used the following code to create MultiSystem from this link

`
def get_system_cls(jdata):
if jdata.get("labeled", False):
return dpdata.LabeledSystem
return dpdata.System

def get_systems(path, jdata):
system = get_system_cls(jdata)
systems = dpdata.MultiSystems(
*[system(os.path.join(path, s), fmt='deepmd/npy') for s in os.listdir(path)])
return systems
`
I tried to pass MultiSystem into rdf function but it gave error - Unable to access cells.

I also tried to calculate iteratively for every data folder -- xx coordinates came the same and RDF values were coming different.

Please suggest what changes do I need to make or Do I need to use data from all the folders for RDF generation ?

Thank You

[wanghan-iapcm]

Weighted average over the rdfs calculated from different System or LabeledSystem should work.

[AnuragKr]

Thanks, @wanghan-iapcm for the response.

I have modified the code based on your suggestion.

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 = [1/len(systems) for _ in range(len(systems))]
    rdf_list = []
    for sys in systems:
        xx,rdf_values,_ = rdf(sys, sel_type = [0,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_OH.out', res.T)

Currently I have given equal weight to every system,if there is some other criteria for allocating weight Please let me know.
Please suggest if I did anything wrong in above code.

[wanghan-iapcm]

I would suggest using the weight: number of frames time number of atoms (properly normalized to guarantee that the summation is equal to 1)

[AnuragKr]

Thanks once again, @wanghan-iapcm for the response.

I have modified the code based on your suggestion.

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 = [0,0], 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_OO.out', res.T)

Output for example/water/se_e2_a. --

Please suggest if I did anything wrong in the above code. I calculated it is normalizing to 1.

[wanghan-iapcm]

Please make sure you are not using example data for productive purpose.
https://github.com/deepmodeling/deepmd-kit/blob/master/doc/train/training.md#warning

you may find productive datasets from https://dplibrary.deepmd.net/

[AnuragKr]

Definitely @wanghan-iapcm I will keep it in mind.
Thanks for your help.

[AnuragKr]

Final Solution --
We need to take the weighted average. And weight will be - frame_no * atom_no
Code --

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 = [0,0], 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_OO.out', res.T)