[gmx-users] Viscosity calculation using cos_acceleration
James Cannon
jamesresearching at gmail.com
Thu Apr 11 09:55:37 CEST 2013
Dear Gromacs users,
This question seems to come up periodically in the mailing list, but none
of the previous answers seem helpful in my case.
I'm trying to reproduce the viscosity calculation of SPC water by Berk Hess
(JCP 116, 2002) using cos_acceleration in Gromacs. The answer I get is 2
orders of magnitude out.
My topology file and parameter file is appended at the bottom of this email.
I run g_energy and get
Energy Average Err.Est. RMSD Tot-Drift
-------------------------------------------------------------------------------
1/Viscosity 23.4689 0.13 2.39126 -0.255371 (m
s/kg)
which gives a viscosity of 0.04 kg/(m s), or 40 mPa.s
The value quoted in the paper is about 0.4 mPa.s which is around the
correct value for water, give or take a bit.
So my question is, where is my missing factor of 100?
Any advice is much appreciated.
Thank you.
James
--------------------------
Topology file:
#include "ffgmx.itp"
#include "spc.itp"
[ system ]
Pure Water
[ molecules ]
SOL 3456
--------------------------
Parameter file for system (3456 SPC water molecules, 3.75x3.75x7.5 nm3 box):
; Generic mdp file for SPC water equilibration
; Gromacs 4.3.x
;
; T = 300 K
;
; NVT 1.2 ns
;
define = ; define here posres etc., e.g. -DPOSRES
integrator = md
tinit = 0
dt = 0.001
nsteps = 1200000
; Bond constraints
continuation = no ; switch to 'yes' if need to read in velocities etc.
constraints = none ; constrain all bond lengths
constraint_algorithm = lincs ; default
lincs_order = 4 ; default
nstxout = 1000
nstvout = 1000
nstfout = 0
nstlog = 1000
nstenergy = 1000
; Output frequency and precision for xtc file
nstxtcout = 1000
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps = System
; Neighbor list
ns_type = grid ; neighlist type
nstlist = 5 ; Freq. to update neighbour list
rlist = 0.9 ; nm (cutoff for short-range NL)
pbc = xyz
periodic_molecules = no
; Non-equilibrium MD
;acc_grps = SYSTEM
cos_acceleration = 0.025 ; PPM option for viscosity calculation (nm/ps²)
coulombtype = PME
rcoulomb = 0.9
optimize_fft = yes ; affects only PME calculations
; if you use PME, set also rcoulomb = rlist
; van der Waals interactions
vdwtype = Cut-off ; Van der Waals interactions
rvdw = 0.9 ; nm (LJ cut-off)
DispCorr = EnerPres ; long-range dispersion correction to energy and
pressure
Tcoupl = berendsen
tc-grps = System
tau_t = 2.5
ref_t = 300.0
;Pressure coupling
Pcoupl = no
gen_vel = no
More information about the gromacs.org_gmx-users
mailing list