[gmx-users] RE: Important: Bugs in NEMD calculation
Xiaohu Li
xiaohuli914 at gmail.com
Wed Jan 19 23:41:25 CET 2011
On Wed, Jan 19, 2011 at 4:40 PM, Xiaohu Li <xiaohuli914 at gmail.com> wrote:
>
> Message: 2
>> Date: Wed, 19 Jan 2011 23:23:12 +0100
>> From: Berk Hess <gmx3 at hotmail.com>
>> Subject: RE: [gmx-users] RE: Important: Bugs in NEMD calculation
>> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>> Message-ID: <COL113-W651FF8B51B6E10B91E0E4F8EF60 at phx.gbl>
>> Content-Type: text/plain; charset="iso-8859-1"
>>
>>
>>
>>
>> It seems you are right.
>> The density is calculated a few lines before 1/visc is stored, but it is
>> not being used.
>> Strange that this always went unnoticed, since it already seems to have
>> been wrong in 4.0.
>> The density of water is so close to 1000 that you wouldn't notice the
>> difference
>>
>> But in a quick test I don't get the correct viscosity either way.
>> I will have to find time for a thorough check.
>> Do you get the right answer by dividing by the missing density term in
>> src/mdlib/mdebin.c?
>>
>> Berk
>>
> Well, my force field can not reproduce the experiment anyway.
> However, the wrong way of calculating in NEMD(vol instead of dens) is
> getting visocisty much much lower than what I got from using g_energy -vis.
> If I'm understanding correctly, the one with visco.xvg(column 1 vs column
> 2) should be the one I'm supposed to look at, even with 12ns simulation,
> this does not converge, but roughly speaking the NEMD results with the vol
> replaced by rho(I wrote a outside script to calculate this) is on the same
> order of magnitude as the g_energy -vis using green-kudo(is that right?)
> expression.
>
> Xiaohu
>
>
>>
>> Date: Wed, 19 Jan 2011 14:48:32 -0600
>> From: xiaohuli914 at gmail.com
>> To: gmx-users at gromacs.org
>> Subject: [gmx-users] RE: Important: Bugs in NEMD calculation
>>
>>
>>
>> Date: Wed, 19 Jan 2011 20:43:20 +0100
>>
>> From: Berk Hess <gmx3 at hotmail.com>
>>
>> Subject: RE: [gmx-users] Important: Bugs in NEMD calculation
>>
>> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>>
>> Message-ID: <COL113-W2019E7174AFFBAEF819A1A8EF60 at phx.gbl>
>>
>> Content-Type: text/plain; charset="iso-8859-1"
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> > Date: Wed, 19 Jan 2011 19:13:12 +0100
>>
>> > From: spoel at xray.bmc.uu.se
>>
>> > To: gmx-users at gromacs.org
>>
>> > Subject: Re: [gmx-users] Important: Bugs in NEMD calculation
>>
>> >
>>
>> > On 2011-01-19 18.36, Xiaohu Li wrote:
>>
>> > > Hi, All,
>>
>> > > I've found a bug in the NEMD calculation for viscosity. This has
>>
>> > > been reviewed in /*Hess's paper at JCP 116 209 2002.*/
>>
>> > > The version of gromacs I'm using is the development version.
>>
>> > > Notice that this version correct a previous but of 4.5.3, where you
>> uses
>>
>> > > NEMD, both the term V(eq. 21 on Hess's paper) and 1/eta(shear
>> viscosity
>>
>> > > inverse) are supposed to be written to the *.edr file, however,
>>
>> > > the 4.5.3 versions does not have this. This version can be retrieved
>> at
>>
>> > >
>> http://repo.or.cz/w/gromacs.git/commit/c83de86d65ce7135be6cef15e9100d7516e6d9a7
>>
>> > > *However, even this version is buggy since eta=A*rho/(V*k**2)(eq. 20
>>
>> > > Hess's paper)*. *I have performed simulation and has found out that
>> the
>>
>> > > V and eta which are written in *edr file does not match according to
>> the
>>
>> > > formula, a little further check on the source code mdebin.c under the
>>
>> > > directory src/mdlib shows that it is actually calculating
>>
>> > > eta=A*Volume/(V*k**2) where density of rho should have been used.
>> (This
>>
>> > > is at line 755 of mdebin.c ).
>>
>> > > I hope everyone who is using this can be aware of this, if you
>> ever
>>
>> > > used this code to produce data, the V is correct from *edr, however,
>> you
>>
>> > > need to manuelly get your eta using the above formula.
>>
>> > > For the GMX developers, I hope anyone of you can correct this
>> bug.
>>
>> > >
>>
>>
>>
>> I fixed this bug recently is the git release-4-5-patches branch.
>>
>> You can get the fix from git and it will be in the 4.5.4 release
>>
>> (no date yet).
>> Which one are you referring to? The one I got is the one you uploaded that
>> fixed the zero viscosity and 2*cosZ*vel-x in edr file.
>> This bug refers the wrong calculation of 1/eta.
>>
>>
>>
>>
>>
>> > Thanks for pointing that out. There also is a small issue in that the
>>
>> > volume is computed for a rectangular box
>>
>> > vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
>>
>> > dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
>>
>> >
>>
>> > which would be incorrect for a non-rectangular box. You should however
>>
>> > use a rectangular box for this kind of calculations, although this is
>>
>> > not enforced by grompp.
>>
>>
>>
>> No.
>>
>> That formula is correct for any triclinic box!
>>
>>
>>
>> Berk
>>
>>
>>
>> >
>>
>> > >
>>
>> > >
>>
>> > > Xiaohu
>>
>> > > *
>>
>> > >
>>
>> >
>>
>> >
>>
>> > --
>>
>> > David van der Spoel, Ph.D., Professor of Biology
>>
>> > Dept. of Cell & Molec. Biol., Uppsala University.
>>
>> > Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
>>
>> > spoel at xray.bmc.uu.se http://folding.bmc.uu.se
>>
>> > --
>>
>> > gmx-users mailing list gmx-users at 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-request at gromacs.org.
>>
>> > Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>>
>>
>> -------------- next part --------------
>>
>> An HTML attachment was scrubbed...
>>
>> URL:
>> http://lists.gromacs.org/pipermail/gmx-users/attachments/20110119/1ffc9f59/attachment.html
>>
>>
>>
>>
>> ------------------------------
>>
>>
>>
>> --
>>
>> gmx-users mailing list
>>
>> gmx-users at 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!
>>
>>
>>
>> End of gmx-users Digest, Vol 81, Issue 127
>>
>> ******************************************
>>
>>
>>
>>
>>
>> --
>> gmx-users mailing list gmx-users at 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-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>> -------------- next part --------------
>> An HTML attachment was scrubbed...
>> URL:
>> http://lists.gromacs.org/pipermail/gmx-users/attachments/20110119/78703d32/attachment-0001.html
>>
>> ------------------------------
>>
>> Message: 3
>> Date: Wed, 19 Jan 2011 17:29:02 -0500
>> From: "Justin A. Lemkul" <jalemkul at vt.edu>
>> Subject: Re: [gmx-users] CNT
>> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>> Message-ID: <4D3765AE.6080307 at vt.edu>
>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>
>>
>>
>> trevor brown wrote:
>> > Dear Justin,
>> > Regarding our previous discussion under "CNT" topic, you want to know
>> > exactly what I did to define C-C interaction. I have summarized it
>> below.
>> >
>> > 1. I copied oplsaa.ff folder in my working directory
>> > 2. I added following lines to atomname2type.n2t
>> > C opls_995 0 12.011 2 C 0.142 C 0.142
>> > C opls_996 0 12.011 3 C 0.142 C 0.142 C 0.142
>> > C opls_997 0 12.011 4 C 0.142 C 0.142 C 0.142 C 0.142
>> > C opls_998 0 12.011 5 C 0.142 C 0.142 C 0.142 C 0.142
>> > C 0.142
>> > 3. I added those to atomtypes.atp
>> > opls_995 12.01100 ;
>> > opls_996 12.01100 ;
>> > opls_997 12.01100 ;
>> > opls_998 12.01100 ;
>> >
>> > 4. I added those to ffbonded.itp
>> > [ bondtypes ]
>> > C C 1 0.14210 478900
>> > [ angletypes ]
>> > C C C 1 120.000 562.2
>> > [ dihedraltypes ]
>> > 5. I constructed .top file for CNT with the following, then I converted
>> > it into .itp by vi editor
>> > g_x2top -f cnt.gro -o cnt.top -nopairs -nexcl 5
>> > 6. By pdb2gmx, I obtained .top for peptide then I converted it into .itp
>> > 7. I combined CNT and peptide's pdbs in pymol and saved.
>> > 8. I wrote a .top file given below,
>> > ;
>> > ; File 'topol.top' was generated
>> > ; By user: onbekend (0)
>> > ; On host: onbekend
>> > ; At date: Mon Jan 10 02:51:19 2011
>> > ;
>> > ; This is a standalone topology file
>> > ;
>> > ; It was generated using program:
>> > ; pdb2gmx - VERSION 4.5.3
>> > ;
>> > ; Command line was:
>> > ; pdb2gmx -ignh -f vpgvg10.pdb
>> > ;
>> > ; Force field was read from the standard Gromacs share directory.
>> > ;
>> > ; Include forcefield parameters
>> > #include "forcefield.itp"
>> >
>>
>> Be careful with this - if you move this topology between directories, it
>> will
>> call any forcefield.itp file that may be in the working directory. What
>> did you
>> remove the "oplsaa.ff" prefix?
>>
>> > ; Include topology for CNT
>> > #include "cnt.itp"
>> > ; Include topology for UW1
>> > #include "uw1.itp"
>> > ; Include water topology
>> > #include "spc.itp"
>> > [ system ]
>> > ; Name
>> > Protein and CNT in water
>> > [ molecules ]
>> > ; Compound #mols
>> > Protein 1
>> > CNT 1
>> > SOL 4019
>> > 9. When I make grompp for minimization it givea me the following error
>> > WARNING 1 [file ffbonded.itp, line 2703]:
>> > Overriding Bond parameters.
>> > old: 0.151 292880 0.151 292880
>> > new: C C 1 0.14210 478900
>> >
>> > ERROR 1 [file ffbonded.itp, line 2707]:
>> > Not enough parameters
>> > Generated 334153 of the 334153 non-bonded parameter combinations
>> > Generating 1-4 interactions: fudge = 0.5
>> > Generated 334153 of the 334153 1-4 parameter combinations
>> > -------------------------------------------------------
>> > Program grompp, VERSION 4.5.3
>> > Source code file: toppush.c, line: 1166
>> > Fatal error:
>> > Atomtype opls_996 not found
>> >
>> > That is all what I exactly did. Is anything wrong or missing?
>> >
>>
>> You never defined nonbonded parameters for your new atom type. You've
>> declared
>> that it exists in the .atp file, but then you didn't add that atom type in
>> ffnonbonded.itp.
>>
>> -Justin
>>
>> > best wishes
>> > trevor
>> >
>> >
>> >
>> >
>>
>> --
>> ========================================
>>
>> 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 list
>> gmx-users at 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!
>>
>> End of gmx-users Digest, Vol 81, Issue 130
>> ******************************************
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20110119/5b5139da/attachment.html>
More information about the gromacs.org_gmx-users
mailing list