[gmx-users] Using lennard jones and buckingham terms simultaneously
Matthew Watkins
matthew.watkins at ucl.ac.uk
Tue Mar 30 17:09:10 CEST 2010
Hi Gareth,
Gareth Tribello wrote:
> Hello again
>
> Still trying to get this carbonate forcefield to work with gromacs. I
> now know that the tables and so on are being read in correctly as I can
> get gromacs to reproduce the energies that I get for the various terms
> when I calculate them by hand/with a separate program. There is stilla
> a problem with the charge groups however. As when I attempt to split up
> the atoms in the water molecule the simulation fails. Again checking
> that I am doing this correctly I would replace:
>
> ; at type res nr res name at name cg nr
> charge mass
> 1 amber99_61 1 SOL OW 1 0
> 16.00000
> 2 amber99_60 1 SOL HW2 1 0.52422
> 1.00800
> 3 amber99_60 1 SOL HW3 1 0.52422
> 1.00800
> 4 MW 1 SOL MW4 1
> -1.04844 0.00000
>
> with
>
> ; at type res nr res name at name cg nr
> charge mass
> 1 amber99_61 1 SOL OW 1 0
> 16.00000
> 2 amber99_60 1 SOL HW2 2 0.52422
> 1.00800
> 3 amber99_60 1 SOL HW3 3 0.52422
> 1.00800
> 4 MW 1 SOL MW4 4
> -1.04844 0.00000
>
> The problem I get (even if I just run water without any
> carbonate/tabulated potentials) if I do the above is that the settles
> algorithm fails.
extract from *.top file for TIP4P_2005 water, below. I think you can
use two charge groups - one with OW and one with HW + MW - both groups
are then neutral, which might be beneficial.
You then need an index group with a HW and MW group (or whatever you
call them). Something like below in you *.mdp.
"
energygrps = F Ca OW HW_AND_T_MW
energygrp_table = OW OW F F F Ca Ca OW F OW F HW_AND_T_MW
"
I am not sure that settles likes the split charge group, I'm using lincs
at the moment.
Cheers,
Matt
"
[moleculetype]
; name nrexcl
water 1
[atoms]
; nr type resnr residu atom cgnr charge
1 OW 1 OW OW 1 0 15.994
2 HW 1 HW HW 2 0.5564 1.008
3 HW 1 HW HW 2 0.5564 1.008
4 MW 1 MW MW 2 -1.1128 0.0
[constraints]
;i j funct doh dhh
1 2 1 0.09572
1 3 1 0.09572
2 3 1 0.15139
[exclusions]
1 2 3 4
2 1 2 3
3 1 2 4
4 1 2 3
....
"
>
> I'm obviously missing something fundamental - I'm not even sure that cg
> nr stands for the charge group. Any help would be greatly appreciated.
>
> Gareth
>
> On Thu, Mar 25, 2010 at 11:07 PM, Mark Abraham <Mark.Abraham at anu.edu.au
> <mailto:Mark.Abraham at anu.edu.au>> wrote:
>
> On 26/03/2010 7:03 AM, Gareth Tribello wrote:
>
> Hello again
>
> I have tried to do as you suggest and use tables but I have a new
> problem. First let me describe my process and then you can let
> me know
> if there is anything wrong in the stages:
>
> OK so first you include the following directives into the mdp file:
>
> coulombtype = pme (or whatever sort of coulomb interaction you
> are using)
> vdw-type = user
>
> energygrps = Ca CCA OCA OW HW
> energygrp_table = Ca OCA Ca CCA OCA OCA OCA OW OCA HW
>
> Gromacs is then (at some stage) going to look for a series of
> files called
>
> table.xvg - which is the default 6-12 Lennard Jones that will
> be used
> for most of the atoms
> table_Ca_OCA.xvg - which are the Buckingham interactions
> between your
> various atom types.
> table_Ca_CCA.xvg
> table_OCA_OCA.xvg
> table_OCA_OW.xvg
> table_OCA_HW.xvg
>
> These files have the format (and contents) described in section
> 6.7 of
> the manual. Finally, you define the various energy groups Ca,
> OCA and
> so on in your index.ndx file.
>
> The problem is that grompp gives me the following error:
>
> "atoms 1 and 2 in charge group 1 of molecule type 'SOL' are in
> different
> energy groups"
>
> (incidentally these atoms 1 and 2 are OW and HW)
>
> Does this mean that I cannot use different tabulated potentials for
> different atoms in a molecules? By which I mean that I can't use
> different tabulated potentials for the OW Ca and HW Ca
> interactions for
> example?
>
>
> Charge groups are the fundamental unit GROMACS uses in constructing
> a simulation. Energy groups are the next "higher" layer in the data
> structures, and these must be sets of whole charge groups. With some
> electrostatics models, looping over charge groups whose charge is
> preferably an integer is essential for modelling correct behaviour.
> GROMACS does a complex sorting of all the interactions between
> charge groups into lists that allow it to iterate over charge groups
> and energy groups. A user table then gets applied to a whole intra-
> or inter- energy-group loop. Thus your attempt violates this
> precondition.
>
> However, PME does not require the use of charge groups for accurate
> results, since all inter-atomic electrostatic interactions get
> treated, regardless of distance. So you could decompose your water
> molecules into two charge groups, O and Hs. (Caveat, a near-brokenly
> bad PME approximation might get a little worse with arbitrary charge
> groups)
>
>
> Final question, as its not clear to me from the manual, if you use a
> tabulated potential for Lennard Jones and you use mix type 2 (so
> are you
> are providing epsilon and sigma in the input rather than A and
> B) does
> gromacs still know that it has to manipulate the input parameters in
> order to get the coefficients of the (tabulated) g(r) and h(r)
> dispersion and repulsion functions (I mean the g(r) and h(r)
> defined in
> section 6.7 of the manual here)? At the same time does it also
> know not
> to do anything to the parameters you input for the (tabulated)
> buckingham potentials (as for a buckingham you are providing A
> and C)?
>
>
> I expect the point of the tables is that GROMACS just uses them per
> equation 6.23. Thus I'd expect C6 and C12 in that equation to be
> constructed according to whatever combination rule is in force. If
> you've specified them explicitly in the topology (see chapter 5),
> then they will not be constructed.
>
> You should be very careful to test your assumptions and deductions
> about how all this machinery is working. Do take the time to set up
> a tabulated version of a quick non-tabulated calculation to at least
> make sure you've done the simple things right! To test the function
> of eq 6.23, do a non-table constructed-parameter calculation, a
> non-table topology-specified-parameter calculation, a table
> constructed-parameter calculation, etc. to be sure you understand
> what is going on - and that the code works right! Please consider
> contributing any insights to a wiki page on the GROMACS webpage.
>
> Before you go to all this work please consider my advice of last
> email. Unless you have an existing reason to expect some ad-hoc
> combination of different interaction functions to work well
> together, you may find that your best possible result is a
> correctly-functioning random number generator.
>
> Mark
>
> On Wed, Mar 24, 2010 at 2:25 AM, Mark Abraham
> <mark.abraham at anu.edu.au <mailto:mark.abraham at anu.edu.au>
> <mailto:mark.abraham at anu.edu.au
> <mailto:mark.abraham at anu.edu.au>>> wrote:
>
>
>
> ----- Original Message -----
> From: Matthew Watkins <matthew.watkins at ucl.ac.uk
> <mailto:matthew.watkins at ucl.ac.uk>
> <mailto:matthew.watkins at ucl.ac.uk
> <mailto:matthew.watkins at ucl.ac.uk>>>
> Date: Wednesday, March 24, 2010 2:59
> Subject: Re: [gmx-users] Using lennard jones and buckingham terms
> simultaneously
> To: Discussion list for GROMACS users <gmx-users at gromacs.org
> <mailto:gmx-users at gromacs.org>
> <mailto:gmx-users at gromacs.org <mailto:gmx-users at gromacs.org>>>
>
> > Hi Gareth,
> >
> > as Vitaly suggested tabulated potentials seem to be the
> only way
> > to go, it took me a while to get up to speed on the
> Gromacs way
> > of doing this, so get in touch if you wish.
> >
> > The tables for buck potentials need to include the
> standard 1/r6
> > term whilst what would be the 1/r12 term needs to contain
> exp(-
> > Bx.rho), the C6 and C12 coefficients can then be put in a
> > standard nonbonded section. You'll need a separated table
> > for each pair of interactions that interact with buckingham
> > potential. Each pair must be an energy group as well.
> >
> > If there is a simpler method I'd love to hear it.
>
> There's probably not a simpler method because it's not a
> widely-used
> procedure. It shouldn't be used at all unless you have
> established
> that the combination of functional forms is effective...
>
> Mark
> --
> gmx-users mailing list gmx-users at gromacs.org
> <mailto:gmx-users at gromacs.org>
> <mailto:gmx-users at gromacs.org <mailto:gmx-users at gromacs.org>>
>
> http://lists.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-request at gromacs.org
> <mailto:gmx-users-request at gromacs.org>
> <mailto:gmx-users-request at gromacs.org
> <mailto:gmx-users-request at gromacs.org>>.
>
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
>
> --
> gmx-users mailing list gmx-users at gromacs.org
> <mailto:gmx-users at gromacs.org>
> http://lists.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-request at gromacs.org
> <mailto:gmx-users-request at gromacs.org>.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>
>
More information about the gromacs.org_gmx-users
mailing list