# Benzene phasing through PBC

**URL:** <https://gromacs.bioexcel.eu/t/benzene-phasing-through-pbc/8495>\
**Category:** User discussions\
**Created:** [March 3, 2024, 1:55pm UTC](https://gromacs.bioexcel.eu/t/benzene-phasing-through-pbc/8495 "2024-03-03T13:55:53Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lucag](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/lucag/32/2583_2.png) [@Lucag](https://gromacs.bioexcel.eu/u/Lucag)\
**Post date:** [March 3, 2024, 1:55pm UTC](https://gromacs.bioexcel.eu/t/benzene-phasing-through-pbc/8495/1 "2024-03-03T13:55:53Z")

</div>

When I visualise my molecule trajectory, the benzene molecule seems to split across the pbc.  
WIll this effect my calculation of the diffusivity using MDAnalysis?

Thank you in advance.

 ![image](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/8/8a24fe9c24b25faf6136554ffa61777940c6ed5d.png)

---

<div class="post-metadata">

**Author:** ![Lucag](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/lucag/32/2583_2.png) [@Lucag](https://gromacs.bioexcel.eu/u/Lucag)\
**Post date:** [March 3, 2024, 2:06pm UTC](https://gromacs.bioexcel.eu/t/benzene-phasing-through-pbc/8495/2 "2024-03-03T14:06:01Z")

</div>

Could this be solved by perhaps defining the COM of benzene?

---

<div class="post-metadata">

**Author:** ![obZehn](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/obzehn/32/2629_2.png) [@obZehn](https://gromacs.bioexcel.eu/u/obZehn)\
**Post date:** [March 3, 2024, 8:56pm UTC](https://gromacs.bioexcel.eu/t/benzene-phasing-through-pbc/8495/3 "2024-03-03T20:56:26Z")

</div>

This is a problem of PBC representation, it is not really physically happening in the box. You can fix it by making molecules whole during the trajectory (gmx trjconv -pbc mol). I think MDAnalysis can handle the periodicity of systems but you should double check, eventually just preprocess the trajectory with GROMACS. You can calculate the diffusivity with GROMACS as well, if needed

---

<div class="post-metadata">

**Author:** ![Lucag](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/lucag/32/2583_2.png) [@Lucag](https://gromacs.bioexcel.eu/u/Lucag)\
**Post date:** [March 3, 2024, 9:23pm UTC](https://gromacs.bioexcel.eu/t/benzene-phasing-through-pbc/8495/4 "2024-03-03T21:23:07Z")

</div>

What I am trying to do is get a diffusion coefficient which correspons roughly to the experimental results.  
I am using charmm and tip3p. I minimise, equilibrate and then do the production run. Then I convert trajectories using the -pbc nojump flag (as per [MDAnalysis documentation](https://docs.mdanalysis.org/stable/documentation_pages/analysis/msd.html)). Finally, I calculated the diffusion coefficient using this code:

```auto
import MDAnalysis as mda
import MDAnalysis.analysis.msd as msd
import matplotlib.pyplot as plt
from scipy.stats import linregress
import numpy as np

Diffusion_Coefficient = ""

def Diffusion_Error(tpr, xtc, selection):
    global msd

    u = mda.Universe(str(tpr), str(xtc))
    MSD = msd.EinsteinMSD(u, str(selection), msd_type='xyz', fft=True)
    MSD.run()

    msd = MSD.results.timeseries

    nframes = MSD.n_frames
    
    timestep = 10
    lagtimes = np.arange(nframes)*timestep
	
    start_time = 0
    end_time = 400   

    start_index = int(start_time/timestep)
    end_index = int(end_time/timestep)

    linear_model = linregress(lagtimes[start_index:end_index],
      	                               		 msd[start_index:end_index])

    slope = linear_model.slope            
    Diffusion_Coefficient = ((slope * 1/(2*MSD.dim_fac) ))

    print(Diffusion_Coefficient)
    
a = input("tpr file ")
b = input("xtc file ")
c = input("selection ")
Diffusion_Error(a,b,c)

```

But the diffusion coefficients differ from experimental results. **Furthermore, when I change the box size, the values are different!!!**

What am I doing wrong?

Thank you in advance. Any help is much appreciated.