Negative Viscosity from periodic perturbation method

User discussions
1

GROMACS version:
GROMACS modification: Yes/No
Here post your question
Hello, I want to calculate viscosity of water from periodic perturbation method. I am using the simulation protocol same as that given in this paper Cookie Absent. However, the amplitude of cos velocity profile and hence viscosity I am getting from gmx energy are negative. Also when I calculate x-direction velocity profile as a function of box height from gromacs trajectory is not a cosine function of (2piz/lz). I have attached the my .mdp file for your reference. What is exactly wrong with this simulation? Thanks in advance
title = water
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 2500000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 5000 ; save coordinates every 1.0 ps
nstvout = 5000 ; save velocities every 1.0 ps
nstenergy = 5000 ; save energies every 1.0 ps
nstlog = 5000 ; 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 = 300 ; 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 = no ; Velocity generation is off
;cos-acceleration
cos-acceleration = 0.005

2

Hi lakshmi,

The first thing I am noticing is that you simulate 50000 steps (100ps) and compute energies every 5000 steps: it means gmx energy only has 10 data points for the viscosity estimate, most probably too few.

I would say try to simulate at least up to 10ns and compute energy every 1000 steps, or possibly even more frequently.

If the flow profile does not look periodic, most probably it means it hasn’t converged to the steady shape yet; mind that for cos-acceleration=0.005 the signal-to-noise ratio may be small.

Also, I would avoid using Berendsen, try C-rescale.

3

P.S. what water model are you simulating (I assume it’s water from your .mdp file)?

4

Hi Michele,
I am using SPC/E water model. I generated topology from pdb2gmx. I used packmol to generate the initial pdb file with 6912 water molecules. Attaching topology file for your reference
;
; File ‘topol.top’ was generated
; By user: sravya (1000)
; On host: sravya-HP
; At date: Fri Nov 25 14:13:30 2022
;
; This is a standalone topology file
;
; Created by:
; :-) GROMACS - gmx pdb2gmx, 2022.3 (-:
;
; Executable: /usr/local/gromacs/bin/gmx
; Data prefix: /usr/local/gromacs
; Working dir: /home/sravya/png/water
; Command line:
; gmx pdb2gmx -f simbox.pdb -p topol.top -o inp.gro
; Force field was read from the standard GROMACS share directory.
;

; Include forcefield parameters
#include “oplsaa.ff/forcefield.itp”

; Include water topology
#include “oplsaa.ff/spce.itp”

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif

; Include topology for ions
#include “oplsaa.ff/ions.itp”

[ system ]
; Name
Built with Packmol

[ molecules ]
; Compound #mols
SOL 6912