[gmx-users] Tabulated potentials make newbies crazy

[Mark Abraham]

ms wrote:
> Dear Mark,
> 
> First of all, thanks for your patient and thorough reply!
> 
> Mark Abraham ha scritto:
>> Yeah, that's an unfortunate fact of life. The manual can't be organized
>> so that everybody finds all the information they want in one location
>> with the detail they need right now. Using tabulated potentials is an
>> advanced feature and if the user needs to search the manual diligently
>> to find the information, that's probably good for them in the long run.
>> There's a good table of contents, index and some cross references, and
>> text-searching is also available...
> 
> Yes, that's what I am doing now but it seems I am not good enough at it
> :) .Don't get me wrong, the manual is truly, amazingly, incredibly good
> -but probably a wiki page on the subject could help people tie the
> strings together more easily. When I'll get all of this correct I wish
> to draft one so that you can correct it.

Sure - wiki contributions are always welcome. Any registered user can 
put one up, so feel free to have a go and if you want to solicit 
feedback, post here and doubtless some people will edit mercilessly :-) 
That's how wikis work.

>>> Just to try, therefore, I removed all occurences of energy groups in a
>>> .ndx files, except for system and C-alpha (I have no solvent). And there
>>> it complains that atoms belong to different energy groups, but belong to
>>> the same charge group.
>>>
>>> I understand therefore that energy groups and charge groups must be
>>> identically defined.
>> No. You should read in the manual about groups - sections 3.3 and 8.1.1
> 
> I've read them -read both of them *now* again, actually- and still I see
> no mention of charge groups in the paragraph you linked -apart perhaps
> for this sentence in 3.3: " This is done separately for Lennard-Jones
> and Coulomb terms." -They are described in 3.4 ,and also in 4.6 for
> example, but I find no information on how they're differently described
> with respect to energy groups. It tells me they're described in the
> topology, but nothing more. If I look at my topology, I find no charge
> group definition.

Sorry, I was a bit incomplete last night. Charge groups are the 
fundamental unit for neighbour-searching (3.4.2) to construct lists of 
charge groups for nonbonded interactions, which determine lists of 
atom-atom interactions. However, the nonbonded interactions are 
evaluated as nested sums, first over energy groups. So for energy groups 
of Protein and SOL, the neighbour search finds all pairs of charge 
groups that are both Protein and inside the cutoffs, and lists them. 
Then all Protein-SOL similarly, then all SOL-SOL. This requires that the 
energy groups be a union of only complete charge groups (and I am not 
aware that this is spelled out anywhere in the manual!). So for energy 
groups of Calpha, Rest_of_Protein and SOL, it would be necessary to use 
an individual charge group for each Calpha. This would usually mean it 
has a net non-integral charge that is equal in magnitude of the charge 
of the group from which it is taken. It is well known that small charge 
groups of non-integral charge can then wander back and forth across 
cut-off boundaries and generate artefacts.

> The problem is that since I have a single molecule now, and the single
> molecule must be neutral, so it must be all a single charge group
> ("Therefore we have to keep groups of atoms with total charge 0
> together. These groups are called charge groups.", 4.6.2).
> 
>> Get the energy group stuff working before you worry about providing user
>> tables. Divide-and-conquer makes troubleshooting much easier.
> 
> Right, thanks...
> 
>>> Anyway, just for the sake of learning how to get gmx use a tabulated
>>> potential, I apply the tabulated potential to all the peptide chain,
>>> (and I rename the table table_protein_protein.xvg). At this point, it
>>> stops asking me for a "table.xvg" file. I realize that it probably wants
>>> to know how to define other non bonded interaction. Since there is only
>>> the protein, I rename the table_protein_protein.xvg as table.xvg and
>>> retry.
>> Group names are probably case sensitive, so your filename must match
>> your group names... the above doesn't look right.
> 
> The group name was in lower case, that is why I called it in lower case
> (for C-alpha I used upper case as in the group name). It recognized it
> somehow. But thanks for pointing it, I'll retry.
> 
>>> And now I cry because mdrun doesn't start, telling me:
>>>
>>>    Fatal error:
>>>    Out of range force value 2.44081e+41 in file 'table.xvg'
>>>
>>> The "out of range force value" error does not appear at all on google. I
>>> understand it's "kinda high" - but it's a very steep repulsion and so it
>>> is naturally like that at distances close to 0.
>> That looks like a bug with mdrun trying to read from a non-existent
>> file. It is always useful to provide your exact command lines, rather
>> than leave us to infer from text.
> 
> Sorry. This is from the script I use to launch the job on the cluster
> (it is not, however, a parallel job):
> 
> mdrun -v -s $STARTDIR/fullmd_xab10_01.tpr -o $STARTDIR/xab10_traj01.trr
> -x $STARTDIR/xab10_traj01.xtc -c $STARTDIR/xab10_final01.gro -e
> $STARTDIR/xab10_ener01.edr 2>&1 > $STARTDIR/output_20091127_xab01.txt
> 
> Now however, thanks to your suggestion, I checked my table and I found
> that I defined a derivative (the 2.44081e+41 value indeed) for r=0.
> 
> Adjusting that (putting it to 0, since r=0 makes little sense), it
> complains for the next, huge value, at r=0.002. The values it complains
> about are always real ones in the table, so it finds it and reads it. I
> can now regenerate it however such that it stays constant when energy is
> larger, say, than 50.000. This is probably the least problem, but some
> advice on it, just for the sake of completeness, could be nice.

If you're running in single precision, that precision cannot represent 
values as large as 10^41. Since in any simulation (but particularly 
coarse-grained one) non-bonded atoms aren't going to get this close, the 
values are next to irrelevant. Just choose 10^38 for anything larger 
than that.

>>> Now, my questions are:
>>> - What is the accepted range of values in tables?
>> I don't think this is the problem.
> 
> It is the least problem probably, given my confusion on energy-charge
> groups, but it seems it is too...
> 
>>> - How do I define a steep repulsion potential correctly?
>> It's terse, but manual 6.7 seems to have the necessary information.
> 
> 6.7 is one of the references I am obviously using, but it gives only
> general (even if essential!!) information, nothing speficic on "good" or
> "bad" potential shapes/values. But probably that's the least problem :)

Knowing a sensible shape is your problem, if you're choosing to 
unbalance the force field by changing one of the contributing potentials...

>>> - Is there a way to define a table for specific atoms in a chain,
>>> leaving the other atoms using the standard VdW potential? I guess it
>>> requires splitting the chain in two energy groups, but I don't
>>> understand if it a)makes sense/creates problems b)how to deal with the
>>> fact that the chain should be a single charge group.
>> Use energygrp-table like you've been trying. See 7.3.12
> 
> That is also something I keep under my pillow, but the charge/energy
> group confusion I must resolve, before going on...
> 
> Thanks a lot and sorry for my stubborness -learning MD almost by
> yourself is not the best of worlds, but that's how my laboratory seems
> to work...

Best of luck!

Mark