RE: [gmx-users] Viscosity calculation using cos_acceleration
Hi, I fixed an issue with the 1/visc output in 4.5.4, but incorretly. I filed a bug report and a fix, you can easily correct your numbers: http://redmine.gromacs.org/issues/1244 https://gerrit.gromacs.org/#/c/2370/ Thanks for reporting this. Cheers, Berk Date: Thu, 11 Apr 2013 16:55:37 +0900 From: jamesresearch...@gmail.com To: gmx-users@gromacs.org Subject: [gmx-users] Viscosity calculation using cos_acceleration 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.132.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 = 120 ; 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 -- gmx-users mailing listgmx-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 listgmx-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
Re: [gmx-users] Viscosity calculation problems (periodic perturbation method)
There have been several bug fixes related to viscosity calculations since the release of version 4.5.3, of which at least one (on 1/10/2011) specifically mentions zero viscosity in the .edr and .log files. I would suggest pulling the latest release-4-5-patches from the git repository and try again. -Justin Xiang Gu wrote: Hi, all, I was trying to reproduce the viscosity results of SPC water reported in Berk Hess' paper (JCP, vol.116 No.1, 209-217, 2002) by using the periodic perturbation method. After reading some old tips on the maillist exchanged by Song, David, etc., I tried to set up the NVT simulation just according to Hess has done, and sent it to calculation. Probably because not everything is properly understood, I processed the edr file with g_energy_d but found both the columns 2CosZ*Vel-X 1/Viscosity always zero. I played with several options of Tcoupl and tau_t, still the same happened. Could anybody experienced in this tutor me what has been improperly specified in the mdp file (see below; system is 3456 SPC water molecules contained in a 3.75*3.75*7.5 nm^3 box--just as in Hess' paper Table II, already well equilibrated at 300 K in advance; perhaps it is unnecessary to run this double precisely, but I was using double precision version (4.5.3.d) just because I forgot to switch to the single precision one, however this shouldn't cause the problem I had ...)? Thanks and wish to see your suggestions soon! Xiang Gu ;Generic mdp file for SPC water equilibration ;Gromacs 4.3.x ; ;T = 300 K ; ;NVT equilibration run, 1.2 ns ; define = ; define here posres etc., e.g. -DPOSRES integrator = md tinit = 0 dt = 0.001 nsteps = 120 ; 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 ; X/V/F/E outputs change these according to your system 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 ; 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= 173529 -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-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
Re: [gmx-users] viscosity calculation
Payman Pirzadeh wrote: Hello, I used g_energy –f *.edr –vis *.xvg –s *.tpr to calculate the viscosity of my system which is water. Two files are generate: *.xvg and the enecorr.xvg. Now, what should I do to calculate the viscosity of my system with these two files? Sorry for such naïve question. Plot the output .xvg; you will find it should contain two quantities: bulk and shear viscosity. -Justin Regards, Payman -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-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
RE: [gmx-users] viscosity calculation
Well, two questions: 1. My bulk values are fluctuating around 60 (I assume cP). But it looks like very weird to me since experimental value for water viscosity is 0.854 cP. My system is at 265K with 4202 molecules. Is sth strange going on? Or I am missing some unit conversions in the plot. 2. For calculation of protein diffusion constant, should we use the bulk or shear viscosity of the solvent (water)? Payman -Original Message- From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of Justin A. Lemkul Sent: September 29, 2010 12:28 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] viscosity calculation Payman Pirzadeh wrote: Hello, I used g_energy f *.edr vis *.xvg s *.tpr to calculate the viscosity of my system which is water. Two files are generate: *.xvg and the enecorr.xvg. Now, what should I do to calculate the viscosity of my system with these two files? Sorry for such naïve question. Plot the output .xvg; you will find it should contain two quantities: bulk and shear viscosity. -Justin Regards, Payman -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-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 listgmx-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
Re: [gmx-users] viscosity calculation
On 2010-09-29 20.52, Payman Pirzadeh wrote: Well, two questions: 1. My bulk values are fluctuating around 60 (I assume cP). But it looks like very weird to me since experimental value for water viscosity is 0.854 cP. My system is at 265K with 4202 molecules. Is sth strange going on? Or I am missing some unit conversions in the plot. Is the experimental value 0.854 at that temperature? Could be the the unit is different in the output too. Doesn't it say in the xvg file? 2. For calculation of protein diffusion constant, should we use the bulk or shear viscosity of the solvent (water)? Why not use g_msd? Payman -Original Message- From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of Justin A. Lemkul Sent: September 29, 2010 12:28 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] viscosity calculation Payman Pirzadeh wrote: Hello, I used g_energy –f *.edr –vis *.xvg –s *.tpr to calculate the viscosity of my system which is water. Two files are generate: *.xvg and the enecorr.xvg. Now, what should I do to calculate the viscosity of my system with these two files? Sorry for such naïve question. Plot the output .xvg; you will find it should contain two quantities: bulk and shear viscosity. -Justin Regards, Payman -- David van der Spoel, Ph.D., Professor of Biology Dept. of Cell Molec. Biol., Uppsala University. Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. sp...@xray.bmc.uu.sehttp://folding.bmc.uu.se -- gmx-users mailing listgmx-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
Re: [gmx-users] viscosity calculation
265 K is -8C, so ice, not liquid water. That would be more viscous than you might expect. Waren Gallin On 2010-09-29, at 12:52 PM, Payman Pirzadeh wrote: Well, two questions: 1. My bulk values are fluctuating around 60 (I assume cP). But it looks like very weird to me since experimental value for water viscosity is 0.854 cP. My system is at 265K with 4202 molecules. Is sth strange going on? Or I am missing some unit conversions in the plot. 2. For calculation of protein diffusion constant, should we use the bulk or shear viscosity of the solvent (water)? Payman -Original Message- From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of Justin A. Lemkul Sent: September 29, 2010 12:28 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] viscosity calculation Payman Pirzadeh wrote: Hello, I used g_energy –f *.edr –vis *.xvg –s *.tpr to calculate the viscosity of my system which is water. Two files are generate: *.xvg and the enecorr.xvg. Now, what should I do to calculate the viscosity of my system with these two files? Sorry for such naïve question. Plot the output .xvg; you will find it should contain two quantities: bulk and shear viscosity. -Justin Regards, Payman -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-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 listgmx-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 listgmx-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
RE: [gmx-users] viscosity calculation
Apparently, in the file the unit in the file is cP. It might be the case that I plotted the wrong column. I think I should get the third column. Regarding the diffusion, I read in papers about rotational diffusion of the proteins which could be calculated having the viscosity of the liquid and hydrodynamic volume of the protein. So, volume could be calculated with radius of gyration + 0.3nm of hydration shell and viscosity is given by g_energy. So, a rotational diffusion constant could be theoretically calculated. Payman -Original Message- From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of David van der Spoel Sent: September 29, 2010 1:03 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] viscosity calculation On 2010-09-29 20.52, Payman Pirzadeh wrote: Well, two questions: 1. My bulk values are fluctuating around 60 (I assume cP). But it looks like very weird to me since experimental value for water viscosity is 0.854 cP. My system is at 265K with 4202 molecules. Is sth strange going on? Or I am missing some unit conversions in the plot. Is the experimental value 0.854 at that temperature? Could be the the unit is different in the output too. Doesn't it say in the xvg file? 2. For calculation of protein diffusion constant, should we use the bulk or shear viscosity of the solvent (water)? Why not use g_msd? Payman -Original Message- From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of Justin A. Lemkul Sent: September 29, 2010 12:28 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] viscosity calculation Payman Pirzadeh wrote: Hello, I used g_energy f *.edr vis *.xvg s *.tpr to calculate the viscosity of my system which is water. Two files are generate: *.xvg and the enecorr.xvg. Now, what should I do to calculate the viscosity of my system with these two files? Sorry for such naïve question. Plot the output .xvg; you will find it should contain two quantities: bulk and shear viscosity. -Justin Regards, Payman -- David van der Spoel, Ph.D., Professor of Biology Dept. of Cell Molec. Biol., Uppsala University. Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. sp...@xray.bmc.uu.sehttp://folding.bmc.uu.se -- gmx-users mailing listgmx-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 listgmx-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
Re: RE: [gmx-users] viscosity calculation using g_energy (3.3.3)
Dear Berk, Thanks for the reply. I have read through your paper. I have still some doubts. When used tcaf#347; I got files called tcaf_all.xvg tcaf_fit.xvg tcaf.xvg visc_k.xvg ( all default names). I think, the file which I should use for fitting, using the formula eta(k) = eta (0) (1-ak²) + O(k^4) the viscoty ( viscosity - k vector plot) is visc_k.xvg. Am I right? IF yes, what all other files ( tcf_all.xvg, tcaf_fit.xvg and tcaf.xvg) do? expecting your reply, Jestin On Fri, 13 Mar 2009 Berk Hess wrote : Hi, I don't understand what you are actually doing now. You seem to be mixing multiple methods. First off all, I would use NPT for all methods, except the one that uses the pressure fluctuation. The pressure will have a large effect on the viscosity and if you run NVT you need to have exactly the right volume. If you use the cosine acceleration method, the 1/viscosity is printed in the energy file, g_energy will plot it for you. g_tcaf is only for use with an equilibrium simulation. If you read the paper, you will have seen an expression to extrapolate the k=0. Berk Date: Fri, 13 Mar 2009 07:33:30 + From: jesb...@rediffmail.com To: gmx-users@gromacs.org Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3) CC: Dear Berk and David, Thank you very much for your appropriate and informative replies. I tried another method (traverse current method) to calculate the shear viscosity ( a non equilibrium method, which has been described in Berk#347; paper : Journal of Chemical Physics, 116, page 209 ( Determining the shear viscosity of model liquids from molecular dynamics simulations)), I used the g_tcaf utility (ie g_tcaf -f traj1.trr -s binary.tpr -oc test.xvg) . As suggested by David, I increased the system size ( from 500 to 2048 TIP4P molecules). I ran in NVT ensemble which allows the pressure to fluctuate. Apart from that I added following options to my mdp file, where accelaration of 1A/ps² was given to the system. ;NON EQUILIBRIUM STUFF acc_grps = system accelerate= 0.1 0.0 0.0 cos_acceleration = 0.1 Moreover, I saved the trajectory in every 1ps ( so total 500 frames for a 500ps simulation) then, I got the following output: k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) - Which shows a strong k dependence over the property: shorter k, better the viscosity, as pointed out in the paper. However, the value obtained is around 0.01 times less than the experimental value (1pa-second). Adding to that, the results obtained by this method seems to be very convincing unlike the g_energy that shows a great divergence!! So the situation is getting better now. Now, I would like to know whether this can be improved if I save the trajectories more frequently ( 500 fs) and run for longer, say 2ns or change value of accelaration . Any thoughts ? regards, Jes. On Thu, 12 Mar 2009 Berk Hess wrote : Hi, This is a very inefficient method for determining the viscosity. Also you need really perfect pressure fluctuations: NVT, shifted potentials, probably even double precision. There was a mail about this recently. There are better methods, have a look at: http://dx.doi.org/10.1063/1.1421362 Berk Date: Thu, 12 Mar 2009 07:39:52 + From: jesb...@rediffmail.com To: gmx-users@gromacs.org Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3) CC: David, Thanks for the quick reply. Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg The output file created includes three columns. 1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity. It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second). The bulk viscosity of water at 300 K is approximately 0.7 cp. But the value ( Bulk viscosity) I got from the program
Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3)
Dear Berk and David, Thank you very much for your appropriate and informative replies. I tried another method (traverse current method) to calculate the shear viscosity ( a non equilibrium method, which has been described in Berk#347; paper : Journal of Chemical Physics, 116, page 209 ( Determining the shear viscosity of model liquids from molecular dynamics simulations)), I used the g_tcaf utility (ie g_tcaf -f traj1.trr -s binary.tpr -oc test.xvg) . As suggested by David, I increased the system size ( from 500 to 2048 TIP4P molecules). I ran in NVT ensemble which allows the pressure to fluctuate. Apart from that I added following options to my mdp file, where accelaration of 1A/ps² was given to the system. ;NON EQUILIBRIUM STUFF acc_grps = system accelerate= 0.1 0.0 0.0 cos_acceleration = 0.1 Moreover, I saved the trajectory in every 1ps ( so total 500 frames for a 500ps simulation) then, I got the following output: k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) - Which shows a strong k dependence over the property: shorter k, better the viscosity, as pointed out in the paper. However, the value obtained is around 0.01 times less than the experimental value (1pa-second). Adding to that, the results obtained by this method seems to be very convincing unlike the g_energy that shows a great divergence!! So the situation is getting better now. Now, I would like to know whether this can be improved if I save the trajectories more frequently ( 500 fs) and run for longer, say 2ns or change value of accelaration . Any thoughts ? regards, Jes. On Thu, 12 Mar 2009 Berk Hess wrote : Hi, This is a very inefficient method for determining the viscosity. Also you need really perfect pressure fluctuations: NVT, shifted potentials, probably even double precision. There was a mail about this recently. There are better methods, have a look at: http://dx.doi.org/10.1063/1.1421362 Berk Date: Thu, 12 Mar 2009 07:39:52 + From: jesb...@rediffmail.com To: gmx-users@gromacs.org Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3) CC: David, Thanks for the quick reply. Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg The output file created includes three columns. 1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity. It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second). The bulk viscosity of water at 300 K is approximately 0.7 cp. But the value ( Bulk viscosity) I got from the program gives me 100 pa-s, an increase of two order of magnitude. I wonder whether I have done anything wrong while specifying the frequency of saving energy file. I have saved the energy file in every 2ps. Isn´t that enough for a simple system like water? OR should I have to save trajectories in every 5fs as suggested by one in a previous post. I post the first 20 lines of the output file. --- # This file was created Thu Mar 12 16:20:09 2009 # by the following command: # g_energy -f water.edr -vis test.xvg # # g_energy is part of G R O M A C S: # # GROup of MAchos and Cynical Suckers # @title Bulk Viscosity @xaxis label Time (ps) @yaxis label \8h\4 (cp) @TYPE xy @ view 0.15, 0.15, 0.75, 0.85 @ legend on @ legend box on @ legend loctype view @ legend 0.78, 0.8 @ legend length 2 @ s0 legend Shear @ s1 legend Bulk 1.99203 9.6633 96.3893 3.98406 11.1625 98.1365 5.9761 12.6631 99.838 7.96813 13.4652 101.366 9.96016 13.7012 100.249 - ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list
RE: [gmx-users] viscosity calculation using g_energy (3.3.3)
Hi, I don't understand what you are actually doing now. You seem to be mixing multiple methods. First off all, I would use NPT for all methods, except the one that uses the pressure fluctuation. The pressure will have a large effect on the viscosity and if you run NVT you need to have exactly the right volume. If you use the cosine acceleration method, the 1/viscosity is printed in the energy file, g_energy will plot it for you. g_tcaf is only for use with an equilibrium simulation. If you read the paper, you will have seen an expression to extrapolate the k=0. Berk Date: Fri, 13 Mar 2009 07:33:30 + From: jesb...@rediffmail.com To: gmx-users@gromacs.org Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3) CC: Dear Berk and David, Thank you very much for your appropriate and informative replies. I tried another method (traverse current method) to calculate the shear viscosity ( a non equilibrium method, which has been described in Berk#347; paper : Journal of Chemical Physics, 116, page 209 ( Determining the shear viscosity of model liquids from molecular dynamics simulations)), I used the g_tcaf utility (ie g_tcaf -f traj1.trr -s binary.tpr -oc test.xvg) . As suggested by David, I increased the system size ( from 500 to 2048 TIP4P molecules). I ran in NVT ensemble which allows the pressure to fluctuate. Apart from that I added following options to my mdp file, where accelaration of 1A/ps² was given to the system. ;NON EQUILIBRIUM STUFF acc_grps = system accelerate= 0.1 0.0 0.0 cos_acceleration = 0.1 Moreover, I saved the trajectory in every 1ps ( so total 500 frames for a 500ps simulation) then, I got the following output: k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s) - Which shows a strong k dependence over the property: shorter k, better the viscosity, as pointed out in the paper. However, the value obtained is around 0.01 times less than the experimental value (1pa-second). Adding to that, the results obtained by this method seems to be very convincing unlike the g_energy that shows a great divergence!! So the situation is getting better now. Now, I would like to know whether this can be improved if I save the trajectories more frequently ( 500 fs) and run for longer, say 2ns or change value of accelaration . Any thoughts ? regards, Jes. On Thu, 12 Mar 2009 Berk Hess wrote : Hi, This is a very inefficient method for determining the viscosity. Also you need really perfect pressure fluctuations: NVT, shifted potentials, probably even double precision. There was a mail about this recently. There are better methods, have a look at: http://dx.doi.org/10.1063/1.1421362 Berk Date: Thu, 12 Mar 2009 07:39:52 + From: jesb...@rediffmail.com To: gmx-users@gromacs.org Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3) CC: David, Thanks for the quick reply. Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg The output file created includes three columns. 1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity. It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second). The bulk viscosity of water at 300 K is approximately 0.7 cp. But the value ( Bulk viscosity) I got from the program gives me 100 pa-s, an increase of two order of magnitude. I wonder whether I have done anything wrong while specifying the frequency of saving energy file. I have saved the energy file in every 2ps. Isn´t that enough for a simple system like water? OR should I have to save trajectories in every 5fs as suggested by one in a previous post. I post the first 20 lines of the output file. --- # This file was created Thu Mar 12 16:20:09 2009 # by the following command: # g_energy -f water.edr -vis test.xvg
Re: [gmx-users] viscosity calculation using g_energy (3.3.3)
JMandumpal wrote: Dear GROMACS users, As explained in the manual ( page 139, section 6.5/3.3.3) I would like to calculate viscosity of my system ( water) using g_energy. I opted(40 Mu-X ) from the g-energy selection. But the unit written on the Y axis of the corresponding xvg file is (kJ mol\S-1\N). I think GROMACS uses SI unit , which is Pa-Second in the case of viscosity. Then why is this discrepancy.? Or did I make any mistake? Mu is the dipole (in Debye). The units of these things are incorrect for everything that is not an energy. This will be fixed in the next gmx version. g_energy -h tells you what to do: g_energy -f ener -vis viscosity --I give the command on the prompt: g-energy -f ener.edr - o viscosity.xvg ; then chose option 40 ( Mu-X). system details: ** My system consists of 500 TIP4P water molecules, ran for 3.5 nanoseconds at 150K, the ensemble is NPT. The version I am using is 3.3.3 regards, Jes Shopping http://adworks.rediff.com/cgi-bin/AdWorks/click.cgi/www.rediff.com/signature-home.htm/1050715...@middle5/2705639_2677138/2680796/1?PARTNER=3OAS_QUERY=null ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php -- David van der Spoel, Ph.D., Professor of Biology Molec. Biophys. group, Dept. of Cell Molec. Biol., Uppsala University. Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. Fax: +4618511755. sp...@xray.bmc.uu.sesp...@gromacs.org http://folding.bmc.uu.se ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3)
David, Thanks for the quick reply. Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg The output file created includes three columns. 1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity. It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second). The bulk viscosity of water at 300 K is approximately 0.7 cp. But the value ( Bulk viscosity) I got from the program gives me 100 pa-s, an increase of two order of magnitude. I wonder whether I have done anything wrong while specifying the frequency of saving energy file. I have saved the energy file in every 2ps. Isn´t that enough for a simple system like water? OR should I have to save trajectories in every 5fs as suggested by one in a previous post. I post the first 20 lines of the output file. --- # This file was created Thu Mar 12 16:20:09 2009 # by the following command: # g_energy -f water.edr -vis test.xvg # # g_energy is part of G R O M A C S: # # GROup of MAchos and Cynical Suckers # @title Bulk Viscosity @xaxis label Time (ps) @yaxis label \8h\4 (cp) @TYPE xy @ view 0.15, 0.15, 0.75, 0.85 @ legend on @ legend box on @ legend loctype view @ legend 0.78, 0.8 @ legend length 2 @ s0 legend Shear @ s1 legend Bulk 1.99203 9.6633 96.3893 3.98406 11.1625 98.1365 5.9761 12.6631 99.838 7.96813 13.4652 101.366 9.96016 13.7012 100.249 - regards, Jes On Thu, 12 Mar 2009 David van der Spoel wrote : JMandumpal wrote: Dear GROMACS users, As explained in the manual ( page 139, section 6.5/3.3.3) I would like to calculate viscosity of my system ( water) using g_energy. I opted(40 Mu-X ) from the g-energy selection. But the unit written on the Y axis of the corresponding xvg file is (kJ mol\S-1\N). I think GROMACS uses SI unit , which is Pa-Second in the case of viscosity. Then why is this discrepancy.? Or did I make any mistake? Mu is the dipole (in Debye). The units of these things are incorrect for everything that is not an energy. This will be fixed in the next gmx version. g_energy -h tells you what to do: g_energy -f ener -vis viscosity --I give the command on the prompt: g-energy -f ener.edr - o viscosity.xvg ; then chose option 40 ( Mu-X). system details: ** My system consists of 500 TIP4P water molecules, ran for 3.5 nanoseconds at 150K, the ensemble is NPT. The version I am using is 3.3.3 ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
RE: [gmx-users] viscosity calculation using g_energy (3.3.3)
Hi, This is a very inefficient method for determining the viscosity. Also you need really perfect pressure fluctuations: NVT, shifted potentials, probably even double precision. There was a mail about this recently. There are better methods, have a look at: http://dx.doi.org/10.1063/1.1421362 Berk Date: Thu, 12 Mar 2009 07:39:52 + From: jesb...@rediffmail.com To: gmx-users@gromacs.org Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3) CC: David, Thanks for the quick reply. Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg The output file created includes three columns. 1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity. It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second). The bulk viscosity of water at 300 K is approximately 0.7 cp. But the value ( Bulk viscosity) I got from the program gives me 100 pa-s, an increase of two order of magnitude. I wonder whether I have done anything wrong while specifying the frequency of saving energy file. I have saved the energy file in every 2ps. Isn´t that enough for a simple system like water? OR should I have to save trajectories in every 5fs as suggested by one in a previous post. I post the first 20 lines of the output file. --- # This file was created Thu Mar 12 16:20:09 2009 # by the following command: # g_energy -f water.edr -vis test.xvg # # g_energy is part of G R O M A C S: # # GROup of MAchos and Cynical Suckers # @title Bulk Viscosity @xaxis label Time (ps) @yaxis label \8h\4 (cp) @TYPE xy @ view 0.15, 0.15, 0.75, 0.85 @ legend on @ legend box on @ legend loctype view @ legend 0.78, 0.8 @ legend length 2 @ s0 legend Shear @ s1 legend Bulk 1.99203 9.6633 96.3893 3.98406 11.1625 98.1365 5.9761 12.6631 99.838 7.96813 13.4652 101.366 9.96016 13.7012 100.249 - regards, Jes On Thu, 12 Mar 2009 David van der Spoel wrote : JMandumpal wrote: Dear GROMACS users, As explained in the manual ( page 139, section 6.5/3.3.3) I would like to calculate viscosity of my system ( water) using g_energy. I opted(40 Mu-X ) from the g-energy selection. But the unit written on the Y axis of the corresponding xvg file is (kJ mol\S-1\N). I think GROMACS uses SI unit , which is Pa-Second in the case of viscosity. Then why is this discrepancy.? Or did I make any mistake? Mu is the dipole (in Debye). The units of these things are incorrect for everything that is not an energy. This will be fixed in the next gmx version. g_energy -h tells you what to do: g_energy -f ener -vis viscosity --I give the command on the prompt: g-energy -f ener.edr - o viscosity.xvg ; then chose option 40 ( Mu-X). system details: ** My system consists of 500 TIP4P water molecules, ran for 3.5 nanoseconds at 150K, the ensemble is NPT. The version I am using is 3.3.3 _ What can you do with the new Windows Live? Find out http://www.microsoft.com/windows/windowslive/default.aspx___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
