Shear viscosity calculation

GROMACS version:2024.1
GROMACS modification: Yes/No
Here post your question
Is there any smaple mdp file to calculate shear viscosity of drug molecule in water? I already obtained diffusion coefficient of drug molecule in water. Do the files obtained during msd run work for calculating viscosity? I will be much thankful if any one could help me to know the steps to be followed to calculate shear viscosity.

Hi,

I am not sure if it makes much sense to calculate the viscosity of ‘a molecule’. Do you mean the viscosity of an acqueous solution of drug molecules?

You can use gmx energy -evisco to calculate the viscosity coefficient using Einstein’s formula. Otherwise you can look into cos-acceleration or deform (.mdp options), which also offer a way to compute viscosity.

@MichelePellegrino
yes i want to calculate viscosity of a drug equilibrated in water I have calculated diffusion using equilibration md but I think to calculate viscosity I need non eq md but I cant find any reliable mdp file to get reference how should I handle it can you suggest me any mdp file what should be it like

There is no such thing as “the viscosity of a molecule in water”. Why do you want to compute this?

One could compute the viscosity of a certain concentration of drug molecules solvated in water. If this concentration is low, the viscosity will be very close to that of pure water.

hello i calculated viscosityusing gmx energy -evisco after equlibration run with . but I got 4 different viscosity with velocity rescaling thermostat for temp coupling at 310 K and pressure was regulated using the Berendsen barostat. upon plotting I get following graph


now how can I calculate viscosity?? i am confused by these graph please help

This is only apparently confusing. The erratic behaviour of the observable has to do with the long-time diffusive behaviour of the integral (squared) of pressure. Either average the results from several runs, or take a long run but only consider the behaviour for short times (until 20000 all the graphs seem to more or less agree). It also seem that the output is very coarse: consider reducing nstcalcenergy (and possibly nstenergy), and increasing the number of samples with (-einstein_restarts in gmx energy, from Gromacs 2024).

hi used cos acceleration to calculate shear viscosity after a bit of study in forums using this mdp file
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 25000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 50 ; save coordinates every 1.0 ps
nstvout = 50 ; save velocities every 1.0 ps
nstenergy = 50 ; save energies every 1.0 ps
nstlog = 50 ; update log file every 1.0 ps
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 0.9 ; short-range electrostatic cutoff (in nm)
rvdw = 0.9 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
lj-pme-comb-rule = Lorentz-Berthelot ;
; Temperature coupling is on
tcoupl = Berendsen ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 310 ; reference temperature, one for each group, in K
Pressure coupling is on
pcoupl = no ; Pressure coupling on in NPT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes
;cos-acceleration
cos-acceleration = 0.005
and it give this type of graph


now this is good right on averaging value it nearly coincides with experimental range for SPC/E
water now why there are negative values at the end

Hi, I cannot say on top of my mind if cos-acceleration=0.005 is too large, but 10 ps is for sure too short. Also, what you are showing is not the viscosity extracted from the periodic perturbation method, but again the Einstein one. You should select the 1/viscosity term of gmx energy, no need to provide the -evisco flag if you are running the cosine perturbation.

plz can you provide full code for gmx energy viscosity calculation using periodic perturabation method and one query we have to calculate gmx energy from equilibration or prooduction plz help its taking much more time than i anticipated.

hi i read your paper on Determining the shear viscosity of model liquids from molecular dynamics simulations but i am not able to figure out which ensemble should i use to calculate shear viscosity i am using production md to calculate viscosity of your fifth type simulation of PME is it correct or not and also while using periodic perturbation method whiich ensemble should be used and how much cos acceleration value will be feasible??

I think I used NVT for all simulations.

But nowadays I would suggest to use a simple Couette flow setup, using the deform option, instead of the cosine acceleration method. You should then use GROMACS version 2024 as that has some major improvements to the implementation of the deform option.

rather using cosine function i used equilibration md as used in your paper I obtained this graph for spce water after five 10 ns simulation

then I averaged all 5 simulation and averaged that also and I get this graph


is this sufficient as I get 5.823 e-4 kg m-1 s-1 after averaging which is near to experimental