Dear all,
I am doing a coarse-grained simulation on a small RNA molecule and I’m interested in finding out how the pairwise interaction energies are calculated. So I printed the pairwise atomic indices and contributions from individual pairs using the following input file:
#LAMMPS Script for CG Model Minimize
units real
dimension 3
boundary p p p
atom_style molecular
special_bonds lj/coul 0.0 0.0 0.0 angle no dihedral no
neighbor 2.0 bin
neigh_modify every 1 delay 1
bond_style harmonic/gauss
angle_style harmonic/gauss
dihedral_style fourier
pair_style morse/gauss2/lj9 13.5
read_data ./data.min.1
log ./log.min.1
thermo 1
thermo_style custom temp etotal epair emol evdwl ecoul ebond eangle edihed eimp
compute p_pair all property/local patom1 patom2
compute e_pair all pair/local dist eng
dump 2 all local 1 pa.txt index c_p_pair[*] c_e_pair[*]
dump 1 all custom 1000 ./dump.min.1.lammpstrj id mol type x y z ix iy iz
min_style sd
minimize 1.0e-4 1.0e-6 500 1000000
The files data.min.1 and pa.txt are provided in the zip file (I cannot attach files to the topic, so here is the Google Drive link).
To understand how pairwise potential is calculated, I wrote a python script to calculate it myself. The script is given by pair_pot.py in the attached zip file. When writing the script, I’ve been speculating what makes some pairs excluded from the pairwise interaction. It seems that all pairs that are already connected either in a bond, angle, and dihedral will be excluded. So I copied the bonds, angles and dihedrals parts from the data.min.1 file into another file called ‘topo.dat’ which is used as input for the Python script, which excludes the pairs accordingly.
The script also needs the coordinates of the last frame, which I extracted them from the dump file, sorted according to atomic index, and put them into the coordinate.dat file. Each line reads (atomic type, x, y, z).
But when I ran my script, although results from all other pairs are consistent with LAMMPS, I found there are 22 pair interactions that are present in my script, but absent from LAMMPS results.
For example, the pair atom 3 and atom 91. Their coordinates are:
2.849060000000000059e+01 -5.511980000000000324e+00 -2.552530000000000143e+01
and
2.272390000000000043e+01 -8.867570000000000618e+00 -1.998349999999999937e+01
respectively. Atom 3 is of type 6, and atom 91 of type 2.
The cutoff radius for pairwise interaction between types 2 and 6 is 10.0, while the distance between these two is 8.643572992. And I also checked there is no bonds, angles or dihedrals involving these two atoms.
A list of all 22 pairs is given below, each line reads (atom index 1, atom index 2, distance):
3 91 8.673313114
4 91 7.989210859
6 88 9.050619365
6 89 8.139045615
13 81 8.902213817
14 81 8.133412624
18 76 9.159147725
19 76 8.078987365
21 73 8.715978913
21 74 8.13284914
26 68 8.820594031
26 69 8.077874046
38 166 8.951728709
39 166 8.214527579
41 163 8.794487254
41 164 7.999854994
46 159 8.331220015
51 153 8.902382068
51 154 8.069052042
56 148 8.699062043
56 149 8.034954938
64 141 8.189484482
I’m wondering what makes these pairs excluded from the calculation of pairwise potentials. I’d be grateful if I can get help on that!
Best,
Wenfei