Viscosity Calculation

[molecular dynamics]

Hello Hoomd users, 

I am working on a simulation of a gas with the end goal of calculating its viscosity. I have completed the simulation. I was wondering what is the best way to obtain the viscosity? Does hoomd have way to calculate viscosity maybe using the velocity correlation function? If the best way is to use hoomd.md.update.mueller_plathe_flow then how does this command work and how would I use it to get a viscosity value?

Thank you in advance 

[mohamad nabizadeh]

Hi,

This is one way that I know of:

- Save the pressure quantities that I need. for instance, if there is shearing in xy:

hoomd.analyze.log(filename='Name_Of_File.log',

                      quantities=['pressure_xy],
                      period=log_dump_period(Preferrarably smaller than 10, set to 1 if your number of time steps is not too large),
                      overwrite=True);

- P = -\tau, Therefore, you are simply saving the negative value of the stress.

And

\mu = -\tau / \Dot{\gamma};

Use Matlab, Python etc. to post-process the output file. I recommend to average over each 500-1000 time steps to cancel out the noises. 

P.S. 

If you need to modify the pressure equation according to your shearing system etc. you can simply modify the ComputeThermo.cc file (/hoomd-Directory/hoomd/ComputeThermo.cc) and recompile the package. That's where the pressure components are calculated.

Best,

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/98cd37aa-7eae-432a-855e-6895e82fed50%40googlegroups.com.

--

Mohammad Nabizadeh

Office : Richards Hall 254

Research Assistant

Department of Mechanical and Industrial Engineering

Northeastern University, Boston, MA

Phone: +18573336911

[Michael Howard]

This method only works if you have Lees-Edwards BCs because this puts the system under a constant shear stress. But, HOOMD currently only deforms the simulation cell and does not have the sllod equations of motion, so you your mileage may vary trying to do this. It will not give you the viscosity using Muller-Plathe's method, since the average shear stress the method produces in the simulation cell is actually 0.

To get the viscosity using an equilibrium measurement, you would use the Green-Kubo relation that computes the viscosity from the stress relaxation function. Note that because the stress fluctuates quite a bit, you probably need to sample it for a very long time with high fidelity to get a reliable value. This can become pretty costly.

To get the viscosity using Muller-Plathe's nonequilibrium method, you should go read about how the method works first:

https://doi.org/10.1103/PhysRevE.59.4894

along with some caveats:

https://doi.org/10.1103/PhysRevFluids.4.043905

https://arxiv.org/abs/1907.07086

You need to extract the generated flow field from the simulation. Comparison of the measured shear rate with the imposed shear stress in the regime of linear response gives the viscosity.

Regards,

Mike

On Tuesday, October 8, 2019 at 1:18:51 PM UTC-5, mohamad nabizadeh wrote:

Hi,

This is one way that I know of:

- Save the pressure quantities that I need. for instance, if there is shearing in xy:

hoomd.analyze.log(filename='Name_Of_File.log',

                      quantities=['pressure_xy],
                      period=log_dump_period(Preferrarably smaller than 10, set to 1 if your number of time steps is not too large),
                      overwrite=True);

- P = -\tau, Therefore, you are simply saving the negative value of the stress.

And

\mu = -\tau / \Dot{\gamma};

Use Matlab, Python etc. to post-process the output file. I recommend to average over each 500-1000 time steps to cancel out the noises. 

P.S. 

If you need to modify the pressure equation according to your shearing system etc. you can simply modify the ComputeThermo.cc file (/hoomd-Directory/hoomd/ComputeThermo.cc) and recompile the package. That's where the pressure components are calculated.

Best,

On Tue, Oct 8, 2019 at 1:15 PM molecular dynamics <marm...@gmail.com> wrote:

Hello Hoomd users, 

I am working on a simulation of a gas with the end goal of calculating its viscosity. I have completed the simulation. I was wondering what is the best way to obtain the viscosity? Does hoomd have way to calculate viscosity maybe using the velocity correlation function? If the best way is to use hoomd.md.update.mueller_plathe_flow then how does this command work and how would I use it to get a viscosity value?

Thank you in advance 

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/98cd37aa-7eae-432a-855e-6895e82fed50%40googlegroups.com.

[molecular dynamics]

Thank you for the suggestion. I have found a way to calculate viscosity using the green-kubo relation method with a python script. Unfortunately, this script (seen below) does not use hoomd and requires paid software. Is there any way to calculate the viscosity using the green-kubo relation in hoomd?

import pylab

import scipy

data_file = 'md_methanol_data.hdf5'
time = nlread(data_file, object_id="Time")[-1]
pressure_tensor = nlread(data_file, object_id="Pressure")[-1]
volume_avg = nlread(data_file, object_id="Volume_Average")[-1]
temperature = nlread(data_file, object_id="Temperature")[-1]
time_step = nlread(data_file, object_id="Time_Step")[-1]

# Reduce down the data by the factor given
# 1 uses all of the gathered data, 10 uses one 10th
# This can be used to speed getting the autocorrelation function
factor = 1
pressure_tensor = pressure_tensor[::factor,:,:]
time = time[::factor]

# Calculate the shear components of the pressure tensor and place them in a convenient array
N = pressure_tensor.shape[0]
P_shear = numpy.zeros((5,N), dtype=float) * Pa
P_shear[0] = pressure_tensor[:N,0,1]
P_shear[1] = pressure_tensor[:N,0,2]
P_shear[2] = pressure_tensor[:N,1,2]
P_shear[3] = (pressure_tensor[:N,0,0] - pressure_tensor[:N,1,1]) / 2
P_shear[4] = (pressure_tensor[:N,1,1] - pressure_tensor[:N,2,2]) / 2

# Calculate the 5 ACF together then average them
# Note that they are only N/2 the length of the data
# because beyond that there are not enough available time windows
# to be accurate
size = int(N/2)
ACF = numpy.zeros((5,size), dtype=float) * Pa**2
for t in range(size):
	ACF[:,t] = numpy.mean(P_shear[:,:N-t] * P_shear[:,t:])
ACF_avg = ACF.mean(axis=0)

# Plotting the ACF
pylab.figure()
pylab.plot(time[:ACF_avg.shape[0]].inUnitsOf(ps), ACF_avg.inUnitsOf(GPa**2), label='ACF')
pylab.xlabel('Time (ps)')
pylab.ylabel('ACF (GPa)')
pylab.legend()
pylab.show()

# Integrate the ACF to get the viscosity
# Note that units need to be put back in because the 
# scipy function call strips them
kbT = boltzmann_constant * temperature
ACF_avg *= volume_avg / (kbT)
gk_raw = scipy.integrate.cumtrapz(
 	y=ACF_avg.inUnitsOf(Pa),
	dx=time_step.inUnitsOf(Second)
)
viscosity = gk_raw * Pa * Second

# Plotting the time evolution of the viscosity estimate
pylab.figure()
pylab.plot(time[:viscosity.shape[0]].inUnitsOf(ps), viscosity.inUnitsOf(Pa*millisecond), label='Viscosity')
pylab.xlabel('Time (ps)')
pylab.ylabel('Viscosity (cP)')
pylab.legend()
pylab.show()

[Michael Howard]

There is nothing magic or proprietary about most of this script--it is just reading in simulation data, autocorrelating the stress tensor (i.e., computing the stress relaxation function), and then integrating it to get the viscosity. So, use HOOMD to run your simulation and log the pressure tensor with an appropriate period to sample the data, like:

hoomd.analyze.log(filename='pressure.log', quantities=('pressure_xx','pressure_xy','pressure_xz','pressure_yy','pressure_yz','pressure_zz'), period=10)

Then, you could borrow the numpy/scipy parts of this script to do the autocorrelation, or write your own code for it (it isn't very difficult to implement).

The same caveats I gave before still apply: sampling requirements are system specific, so you will need to play around with the sampling period, the length of the autocorrelation window you compute (so that the integral converges), and how long you need to sample for to get reliable statistics.

Regards,

Mike

[Ludwig Schneider]

Mike is right. If you want to use Green-Kubo, do as he suggests.

If you want to use the MuellerPlathe method, I would like to point you
to these two papers:

http://dx.doi.org/10.1103/PhysRevE.59.4894

Explains the original method, while in

https://doi.org/10.1016/j.commatsci.2019.109107

I explain the rationale and an application of the implementation in
HOOMD-blue.
I hope, after reading both papers, it becomes a bit clearer how to use
the implementation in HOOMD-blue.

Mike already mentioned, that there are some box geometries that lead to
unexpected flow patterns, it is relatively easy to avoid these.