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 <[email protected] <mailto:[email protected]>> 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
        <[email protected] <mailto:[email protected]>
        <mailto:[email protected]
        <mailto:[email protected]>>> wrote:



           ----- Original Message -----
           From: Matthew Watkins <[email protected]
        <mailto:[email protected]>
           <mailto:[email protected]
        <mailto:[email protected]>>>
           Date: Wednesday, March 24, 2010 2:59
           Subject: Re: [gmx-users] Using lennard jones and buckingham terms
                simultaneously
           To: Discussion list for GROMACS users <[email protected]
        <mailto:[email protected]>
           <mailto:[email protected] <mailto:[email protected]>>>

            > 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 [email protected]
        <mailto:[email protected]>
           <mailto:[email protected] <mailto:[email protected]>>

           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 [email protected]
        <mailto:[email protected]>
           <mailto:[email protected]
        <mailto:[email protected]>>.

           Can't post? Read http://www.gromacs.org/mailing_lists/users.php


-- gmx-users mailing list [email protected]
    <mailto:[email protected]>
    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 [email protected]
    <mailto:[email protected]>.
    Can't post? Read http://www.gromacs.org/mailing_lists/users.php



--
gmx-users mailing list    [email protected]
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 [email protected].
Can't post? Read http://www.gromacs.org/mailing_lists/users.php

Reply via email to