[I-PI] RDF Graph Plot Error

[AnuragKr]

Hello,

I have tried to draw a rdf plot using the md_analysis package from the I-PI simulation of liquid_water. But the Rdf graph not coming accurate. Please guide me on what I am missing if someone tried before to draw rdf graph following the same procedure. Or If any other better way to draw an RDF plot Please suggest.

Simulation has been performed using water example data.

RDF graph --

Code Snipped for input.xml

<simulation verbosity='high'>
  <output prefix='simulation'>
    <properties stride='50' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, pressure_cv{megapascal}, pressure_md{megapascal}, density{g/cm3} ] </properties>
    <!--trajectory filename='pos' stride='1'> positions </trajectory-->
    <trajectory stride='50' format='xyz' filename='pos'> positions{angstrom} </trajectory>
    <checkpoint stride='200'/>
  </output>
  <total_steps>625000</total_steps>
  <prng>
    <seed>31415</seed>
  </prng>
  <ffsocket name='lammps' mode='unix'>
    <address>localhost</address>
    <latency>0.001</latency>
    <port>8080</port>
  </ffsocket>
  <system>
    <initialize nbeads='8'>
      <file mode='xyz'> conf.xyz </file>
      <velocities mode='thermal' units='kelvin'> 300 </velocities>
    </initialize>
    <forces>
      <force forcefield='lammps'> </force>
    </forces>
    <motion mode='dynamics'>
      <dynamics mode='npt'>
       <timestep units='femtosecond'> 0.48 </timestep>
       <barostat mode='isotropic'>
        <thermostat mode='gle'>
          <A shape='(5,5)'>
            [   1.119252684345e-2,    8.975945537096e-6,   -2.728650447177e-5,   -7.011051989908e-4,    8.414373543550e-3,
          1.194065829660e-7,    1.071876106695e-6,   -9.439643019253e-6,   -1.329262116026e-5,   -2.667325152958e-6,
          2.161979961890e-4,    9.439643019253e-6,    4.639122293442e-5,    4.329184279724e-6,    8.132076333400e-5,
         -7.322257663569e-4,    1.329262116026e-5,   -4.329184279724e-6,    5.173717780694e-4,    1.590872642196e-5,
          8.299189140989e-3,    2.667325152958e-6,   -8.132076333400e-5,   -1.590872642196e-5,    6.992095202254e-3
           ]
          </A>
        </thermostat>
        <tau units='femtosecond'> 1000 </tau>
       </barostat>
 <thermostat mode='gle'>
          <A shape='(5,5)'>
            [   1.119252684345e-2,    8.975945537096e-6,   -2.728650447177e-5,   -7.011051989908e-4,    8.414373543550e-3,
          1.194065829660e-7,    1.071876106695e-6,   -9.439643019253e-6,   -1.329262116026e-5,   -2.667325152958e-6,
          2.161979961890e-4,    9.439643019253e-6,    4.639122293442e-5,    4.329184279724e-6,    8.132076333400e-5,
         -7.322257663569e-4,    1.329262116026e-5,   -4.329184279724e-6,    5.173717780694e-4,    1.590872642196e-5,
          8.299189140989e-3,    2.667325152958e-6,   -8.132076333400e-5,   -1.590872642196e-5,    6.992095202254e-3
           ]
          </A>
       </thermostat>
      </dynamics>
    </motion>
    <ensemble>
      <temperature units='kelvin'> 300 </temperature>
      <pressure units='megapascal'> 0.1 </pressure>
    </ensemble>
  </system>
</simulation>

Code snipped for ---water.json

{
    "verbose":          false,
    "use_unix":         true,
    "port":             8080,
    "host":             "localhost",
    "graph_file":       "graph.pb",
    "coord_file":       "conf.xyz",
    "atom_type" : {
        "OW":           0,
        "HW1":          1,
        "HW2":          1
    }
}

Code snippet for -- rdf.py

import MDAnalysis as mda
from MDAnalysis.analysis import rdf
from matplotlib import pyplot as plt

# Load simulation results with a single line
u = mda.Universe('conf.xyz','simulation.pos_0.xyz','simulation.pos_1.xyz','simulation.pos_2.xyz','simulation.pos_3.xyz','simulation.pos_4.xyz','simulation.pos_5.xyz','simulation.pos_6.xyz','simulation.pos_7.xyz')

print(f'No. Of Frames : {u.trajectory.n_frames}')
print('------------------------')
print(f'Atoms --- {u.atoms}')

def reset_volume(ts, dimensions=[30,30,30,90,90,90]):
   ts.dimensions = dimensions
   return ts

u.trajectory.add_transformations(reset_volume)
V = [ts.dimensions.copy() for ts in u.trajectory]

print(f'Dimensions -- {u.trajectory.ts.dimensions}')
print('--------------------------')

g1 = u.select_atoms('type H')
g2 = u.select_atoms("type O")

irdf1 = rdf.InterRDF(g1, g1, nbins=100, range=(0.0, 5), verbose=True)
irdf1.run()
irdf2 = rdf.InterRDF(g2, g2, nbins=100, range=(0.0, 5), verbose=True)
irdf2.run()

x1=irdf1.bins
y1=irdf1.rdf
x2=irdf2.bins
y2=irdf2.rdf

#plt.plot(x1,y1)
plt.plot(x2,y2)
plt.ylabel('g(r)')
plt.xlabel('r(A)')
plt.savefig('rdf.png',bbox_inches = 'tight')