[gmx-users] Viscosity calculation using cos_acceleration

[James Cannon]

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