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
[gmx-users] viscosity calculation using g_energy (3.3.3)
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? --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___ 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
