Issue with time-dependent diffusivity within a cylinder mesh

[fragrussu]

Hi J.R.,

I am trying to simulate diffusion within a meshed cylinder of diameter 3um (content of PLY file stored here newcyl.txt).

The intrinsic diffusivity is 2.2 um2/ms and I simulate 80ms - it should be long enough to see the effect of restriction and observe diffusivity dropping as a function of simulation time, but all I can see is a fluctuation of the diffusivity about 2.2 um2/ms - essentially free diffusion .

I have loaded the spin trajectories in python as:

Nspins = 500
Ntime = 1000
Nsteps = Ntime + 1
Tdur = 0.080  # total simulation duration in sec
dt = Tdur / float(Ntime)  # duration of each step
T = np.linspace(0,Tdur,Nsteps)
data = np.array(np.fromfile("sim_0.traj", dtype="float32"))

and calculated the X, Y and Z displacements as

x = data[0:np.copy(data).size:3]
X = np.zeros((Nspins,Nsteps))
# Store all spins position in a matrix of size Nspins x Nsteps
for tt in range(0,Nsteps):
    X[:,tt] = x[tt:x.size:Nsteps]
# Calculate displacements from the initial position    
for nn in range(0,Nspins):
    X[nn,:] = X[nn,:] - X[nn,0]
# Calculate diffusion coefficient
Dx = np.mean(X*X,axis=0)/(2*T)

(similarly for Y and Z).

The simulation configuration file is this conf.txt.

Any suggestion on why this is happening would be more than welcome! Thanks in advance!

Francesco

[jonhrafe]

Hello again :). Ok, so, several notes in here.

This case is a little bit trickier, it has to do with the faces' normal direction in the PLY model (see the image below); to correctly orient your mesh, the convention I follow --which is the same as Blender and other software-- is that the face normal should always point outwards:

As you can see, first, the normal vectors of the faces are pointing "inwards", with the exception of the bottom "cap" of the cylinder, which is actually well oriented. See below another mesh I made using Blender with all the normal vectors facing outwards (here attached so you can use it cylinder_3um.txt , it's on the micrometres scale ).

Lastly, checking the PLY model I realized that it's bigger than you described, it actually has a radius of 3 um (so diameter of 6 um), see below (I changed the scale to be able to see it)

So what is happening then is that, first, the cylinder is bigger than you probably thought, and, since the normals are not well oriented the parameter: _ini_walkerspos intra will fail to initialize correctly only inside the mesh (it will actually initialize mostly outside).

Final notes.

So maybe you would like to retry your experiments with the new mesh and check if all is in order. A couple of comments, first, I guess this is a test to try mesh models, but if you will only use arbitrarily placed cylinders, then you could use the _cylinders_filelist type of obstacle, it's much faster and accurate.

Second, there's an option to automatically store the mean squared displacement of all the particles, over several directions, and at several points of time (which I haven't documented, sorry), this way you don't have to store all the trajectory file; to do it basically you have add the following lines to the .conf file:

<log>
<propagator>
directions ./normalize_directions.txt
1 100 1000
</propagator>
</log>

where the list 1,100, 1000, indicates the steps at which the simulator will compute and store the mean squared displacement for each one of the directions in _normalizedirections.txt which is just a file with the list of directions to measure the mean squared displacement, like below.

1 0 0 
0 1 0
0 0 1

Hopefully, this helps, cheers!

[fragrussu]

Hi J.R.,

this indeed helped a lot, thank you so much! As you suspected I am not interested in cylinders, I am just playing around with the simulator before I use it for my study.

I think I have now found my way and I am now getting meaningful simulations.

However, it seems to me that the diffusivity and the simulation duration in the configuration file need to be entered respectively in mm2/ms and ms, rather than m2/s and s as stated in the tutorials. Those are the units that are printed on the command line when MCDC runs:

Here is the example I could successfully run. I simulated diffusion within a cube of sides (x,y,z) = (1um,6um,10um) with intrinsic diffusivity of 2 um2/ms for 300ms. I used this configuration file, this scheme file and this mesh.

Here is an image of the mesh, with normals correctly pointing outwards as per your suggestion:

Here are the x positions of all spins over time:

Here are the y positions of all spins over time:

Here are the z positions of all spins over time:

Here is the time-dipendent diffusivity and the visulation of a random walk of a spin:

This is the code used to analyse the MCDC results:

# coding: utf-8

# In[1]:

import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

# In[2]:

Nspins = 750
Ntime = 1000
Nsteps = Ntime + 1
Tdur = 0.3  # total simulation duration in sec
dt = Tdur / float(Ntime)  # duration of each step
T = np.linspace(0,Tdur,Nsteps)

# In[3]:

data = np.array(np.fromfile("aniso/sim_0.traj", dtype="float32")) # Positions are provided in mm!
data = 0.001*data # Scale positions to meters

# In[4]:

###  Calculate x diffusivity
x = data[0:np.copy(data).size:3]
X = np.zeros((Nspins,Nsteps))

# Store all spins position in a matrix of size Nspins x Nsteps
for tt in range(0,Nsteps):
    X[:,tt] = x[tt:x.size:Nsteps]

# Calculate displacements from the initial position   
dX = np.zeros((Nspins,Nsteps))
for nn in range(0,Nspins):
    dX[nn,:] = X[nn,:] - X[nn,0]

# Calculate diffusion coefficient
Dx = np.mean(dX*dX,axis=0)/(2*T)

# In[5]:

###  Calculate y diffusivity
y = data[1:np.copy(data).size:3]
Y = np.zeros((Nspins,Nsteps))

# Store all spins position in a matrix of size Nspins x Nsteps
for tt in range(0,Nsteps):
    Y[:,tt] = y[tt:y.size:Nsteps]

# Calculate displacements from the initial position    
dY = np.zeros((Nspins,Nsteps))
for nn in range(0,Nspins):
    dY[nn,:] = Y[nn,:] - Y[nn,0]

# Calculate diffusion coefficient
Dy = np.mean(dY*dY,axis=0)/(2*T)

# In[6]:

### Calculate z diffusivity
z = data[2:np.copy(data).size:3]
Z = np.zeros((Nspins,Nsteps))

# Store all spins position in a matrix of size Nspins x Nsteps
for tt in range(0,Nsteps):
    Z[:,tt] = z[tt:z.size:Nsteps]

# Calculate displacements from the initial position
dZ = np.zeros((Nspins,Nsteps))
for nn in range(0,Nspins):
    dZ[nn,:] = Z[nn,:] - Z[nn,0]

# Calculate diffusion coefficient
Dz = np.mean(dZ*dZ,axis=0)/(2*T)

# In[7]:

### Calculate overall diffusivity
D = np.mean(dX*dX + dY*dY + dZ*dZ,axis=0)/(6*T)

# In[8]:

plt.plot(T,np.transpose(X[0:1000,:]))
plt.title('All spins')
plt.xlabel('time [sec]')
plt.ylabel('x(t) [m]')
plt.show()

plt.plot(T,np.transpose(Y[0:1000,:]))
plt.title('All spins')
plt.xlabel('time [sec]')
plt.ylabel('y(t) [m]')
plt.show()

plt.plot(T,np.transpose(Z[0:1000,:]))
plt.title('All spins')
plt.xlabel('time [sec]')
plt.ylabel('z(t) [m]')
plt.show()

# In[17]:

plt.plot(T,Dx,label='Dx')
plt.plot(T,Dy,label='Dy')
plt.plot(T,Dz,label='Dz')
plt.plot(T,D,label='D')
plt.title('Diffusion-time dependent diffusivity')
plt.xlabel('time [sec]')
plt.ylabel('D [m2/s]')
plt.legend()
plt.show()

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(X[0,:],Y[0,:],Z[0,:])
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_zlabel('z [m]')
plt.title('Random walk of one spin')
plt.show()

[jonhrafe]

Great catch! Indeed, I forgot to add an important line to the free diffusion tutorial .conf file:

_scale_fromstu 1

If this option is set to 0 (default), then the parameters for the duration, diffusivity and the scheme_file, are expected to be in mm, ms, and mm^s/ms. I have corrected the example/tutorial.

Here I added the .conf file with this option added and the aforementioned command in case you want to save the propagators.

corrected_conf.zip

Thanks!

[fragrussu]

Amazing, thanks a lot for this.

Cheers, Francesco