Dear Vitaly and other Gromacs users,
Thanks for your reply. That's a good suggestion, and actually one of the
first things I tried, though without success. Up to around cos_acceleration
= 0.025 the viscosity is constant around 40 mPa.s, and then for increasing
cos_acceleration the viscosity decreases to around 14 mPa.s before
increasing again and then crashing when it's too high.
I've tried writing my own topology files instead of relying on those
provided by Gromacs (especially since the charges are zero in
ffnonbonded.itp!?) but I still get the same result. I enclose the new files
and grompp command below.
Have you explicitly tried SPC water? Can you see any significant difference
between your input file and mine? Since you say you've had good experience
with this utility, I guess I'm missing something.
Does anyone have any working example that they might kindly be able to
forward? Has anyone else had trouble with this in version 4.5.5?
Best regards,
James
---
SPC.top
---
#include "FF.itp"
#include "mol.itp"
[ system ]
; Name
SPCWater
[ molecules ]
; Compound #mols
SPCWater 3456
---
FF.itp
---
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 3 yes 1 1
[ atomtypes ]
; name bond_type mass charge ptype sigma epsilon
OW OW 15.9994 -0.820 A 0.3166 0.650
HW HW 1.008 0.4100 A 0.000000 0.000000
---
mol.itp
---
[ moleculetype ]
; molname nrexcl
SPCWater 2
[ atoms ]
; id at type res nr residu name at name cg nr charge
1 OW 1 SOL OW 1 -0.8200
2 HW 1 SOL HW 1 0.41
3 HW 1 SOL HW 1 0.41
[ settles ]
; OW funct doh dhh
1 1 0.10000 0.1632991
[ exclusions ]
1 2 3
2 1 3
3 1 2
---
parameters.md
---
integrator = md
tinit = 0
dt = 0.001
nsteps = 1000000
; Bond constraints
continuation = no ; switch to 'yes' if need to read in velocities etc.
constraints = all-angles ; constrain all bond lengths
constraint_algorithm = shake ; default
; X/V/F/E outputs ---- change these according to your system ----
nstxout = 100
nstvout = 100
nstfout = 0
nstlog = 100
nstenergy = 100
; Output frequency and precision for xtc file
nstxtcout = 100
xtc-precision = 10000
; 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 ; xyz(default), no, full(infinite systems)
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 ; berendsen, tau_p = 1.0 for faster equilibration ;Pcoupltype =
isotropic ; semi-isotropic: xy and z separately (CNT)
;tau_p = 1.0 ; ps
;compressibility = 4.5e-5 ; 1/bar (water @ 1 atm, 300 K)
;ref_p = 1.0 ; bar0 4000 1.034798 0.001949
gen_vel = no
;gen_temp = 300.0
;gen_seed = 17352
---
SPC.pdb
---
HEADER 500 MUST_PBC
CRYST1 37.500 37.500 75.000 90.00 90.00 90.00 P 1 1
TITLE SPC water
REMARK
ATOM 1 OW LIG A 1 7.895 13.619 5.312 1.00 0.00
OW
ATOM 2 HW LIG A 1 7.518 12.909 5.907 1.00 0.00
HW
ATOM 3 HW LIG A 1 8.801 13.886 5.645 1.00 0.00
HW
ATOM 4 OW LIG A 2 9.106 26.048 21.988 1.00 0.00
OW
ATOM 5 HW LIG A 2 9.678 25.506 21.372 1.00 0.00
HW
ATOM 6 HW LIG A 2 9.683 26.675 22.512 1.00 0.00
HW
ATOM 7 OW LIG A 3 13.633 12.792 16.603 1.00 0.00
OW
ATOM 8 HW LIG A 3 12.823 13.376 16.662 1.00 0.00
HW
ATOM 9 HW LIG A 3 13.380 11.850 16.829 1.00 0.00
HW
ATOM 10 OW LIG A 4 31.045 18.132 18.190 1.00 0.00
OW
ATOM 11 HW LIG A 4 30.079 17.879 18.240 1.00 0.00
HW
... etc ...
PREP COMMAND: grompp -f parameters.mdp -c SPC.pdb -p SPC.top -maxwarn 2
On 16 April 2013 16:36, Dr. Vitaly Chaban <vvcha...@gmail.com> wrote:
> >Alternatively, has anyone else reproduced the viscosity calculation, or
>
> > tried to?
> > If anyone would be so kind as to forward their input files, that might
> help
> > me narrow down the problem with my own input files.
> >
> > Thank you.
> >
> > James
> >
> > On 11 April 2013 16:55, James Cannon <jamesresearch...@gmail.com> wrote:
> >
> > > 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?
> >
>
>
>
> Non-equilibrium method for viscosity heavily depends on the acceleration
> value. Your selection looks reasonable for me, but try to decrease the
> value. "Ideal" acceleration depends on the particular collective dynamics
> in each system, of course.
>
> I can soothe you ability a workability of the utility. It has been tried
> many times and works really perfectly.
>
>
> Dr. Vitaly Chaban
> --
> gmx-users mailing list gmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> * Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-requ...@gromacs.org.
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
