[gmx-users] Re: energy conservation: shift vs shifted user potential

[Mark Abraham]

On 12/06/2012 9:05 PM, Anja Kuhnhold wrote:
> Ok, that didn't help me much.
> Of course I changed both, potential and force to shifted ones.

That's not what you said you did with g(x) and h(x)... but whatever.

> And sure I used a little fortran code to generate the tables --
>
> I have found the problem:
> I thought GROMACS would ignore the table values above the cutoff distance,
> but it does not.
> I added zeros and now it works.

The tables have to extend past the cut-off because particles can diffuse 
past the cut-off between neighbour list updates. There's an .mdp setting 
controlling how far the extension goes.

Mark

>
> Thanks anyway
> Anja
>
>
>
>
>> ----- Original Message -----
>> ----- Original Message -----
>> From Mark Abraham<Mark.Abraham at anu.edu.au>
>> Date Sat, 09 Jun 2012 12:59:17 +1000
>> To Discussion list for GROMACS users<gmx-users at gromacs.org>
>> Subject Re: [gmx-users] Re: Re: energy conservation: shift vs
>> shifted	user potential
>> On 9/06/2012 2:05 AM, Anja Kuhnhold wrote:
>>> Ok, I try to explain it again :)
>>> I know the format of the user tables (x, f(x), -f'(x), g(x), -g'(x),
>> h(x), -h'(x);
>>> f is coulomb, g dispersion, h repulsion).
>>> For testing the simple cut-off version, my table looks like the
>>> example in the installation.
>>> So I run a simulation with vdwtype=cut-off
>>> and the same with vdwtype=user.
>>> The results are the same.
>>>
>>> The problem is, that I want to use shifted forces and potentials
>> within my tables.
>>> (Aim is better energy conservation)
>>> In the manual (ch. 4.1.5) is explained how the shift functions
>>> for forces and potentials have to be.
>>> So for -g(x) and h(x) I use the potential from the manual, Phi, with
>> alpha=6 and 12, respectivly.
>>> And for g'(x) and -h'(x) I use the shifted force from the manual,
>> Fs, with alpha=6 and 12, respectivly.
>>> But when I run the simulation with this table on the one hand
>>> and with vdwtype=Shift on the other hand,
>>> the results are not comparable, mostly the energy is not conserved
>> by using the table.
>>
>> The force is the negative derivative of the potential. If you change
>> only the force columns then you have made the interaction piecewise
>> discontinuous at every segment. So when GROMACS uses both the values
>> and
>> derivatives to convert your table into the internal format used
>> (section
>> 6.7.1), it comes up with junk. There should be a warning printed if
>> the
>> force column is not close enough to the numerical derivative of the
>> potential column, but maybe this change is close enough that you avoid
>>
>> that. Anyway, this makes energy basically impossible to conserve. You
>>
>> should generate all your table values from the actual function you
>> wish
>> to implement, e.g. with a spreadsheet or scripting language.
>>
>> Mark
>>
>>> Can you follow me?
>>>
>>> Anja
>>>
>>>
>>>
>>> ----- Ursprüngliche Nachricht -----
>>> Von: gmx-users-request at gromacs.org
>>> Datum: Freitag, Juni 8, 2012 15:55
>>> Betreff: gmx-users Digest, Vol 98, Issue 52
>>> An: gmx-users at gromacs.org
>>>
>>>> ----- Ursprüngliche Nachricht -----
>>>> Von Mark Abraham<Mark.Abraham at anu.edu.au>
>>>> Datum Fri, 08 Jun 2012 23:40:35 +1000
>>>> An Discussion list for GROMACS users<gmx-users at gromacs.org>
>>>> Betreff Re: [gmx-users] Re: energy conservation: shift vs shifted user
>>>> potential
>>>> On 8/06/2012 5:59 PM, Anja Kuhnhold wrote:
>>>>> Hello all,
>>>>>
>>>>> can someone give me a hint, please?
>>>>> Do you need more information?
>>>>> Has anyone had a similar problem.
>>>>>
>>>>> I really need to figure that out, because I will simulate other systems
>>>>> which can be run only with user tables.
>>>>>
>>>>> Anja
>>>>>
>>>>>
>>>>>
>>>>> -
>>>>>> ----- Original Message -----
>>>>>> ----- Original Message -----
>>>>>>   From Anja Kuhnhold<anja.kuhnhold at physik.uni-halle.de>
>>>>>> Date Wed, 06 Jun 2012 15:07:40 +0200
>>>>>> To gmx-users at gromacs.org
>>>>>> Subject [gmx-users] energy conservation: shift vs shifted user potential
>>>>>> Dear gmx-users,
>>>>>>
>>>>>> I have a problem concerning energy conservation when using user
>>>>>> potentials (tables).
>>>>>> I'm using gromacs 4.5.4
>>>>>> I simulate a simple Lennard-Jones(6-12) +Fene polymer melt (1600
>>>>>> chains a 10 beads in a 26x26x26 periodic box).
>>>>>>
>>>>>> I tried different vdwtypes (cutoff always 3.24):
>>>>>> The cut-off version does not conserve energy -- okay.
>>>>>> The shift and switch versions conserve energy -- fine.
>>>>>>
>>>>>> Now I wanted to do the same with user tables:
>>>>>> Simple Lennard-Jones table gives really the same results as the
>>>>>> cut-off version -- good.
>>>>>>
>>>>>> But if I use a table with shifted Lennard-Jones potential it is not
>>>>>> comparable to the shift version
>>>>>> and the energy is not conserved -- ?
>>>>>>
>>>>>> I use a shift function as written in the manual (chapter 4.1.5) --
>>>>>> there must be a factor alpha added in the constants A and B --
>>>>>> (r1=0).
>>>>>>
>>>>>> The parameters are the same for shift version and shifted user version.
>>>>>>
>>>>>> Has someone an idea why the shifted user potential doesn't work in
>>>>>> this way?
>>>> We've no real idea what you've done... manual 6.7.2 describes the
>>>> required format and there are (unshifted) examples in your
>>>> installation
>>>> under $GMXLIB/share/top/table*.xvg.
>>>>
>>>> Makr
>>>>
>>>>>> Here is the mdp:
>>>>>>
>>>>>>
>>>>>> integrator        	= md-vv
>>>>>> dt                	= 0.0035
>>>>>> nsteps                	= 1000
>>>>>> nstxout                	= 1
>>>>>> nstvout                	= 1
>>>>>> nstfout                	= 1
>>>>>> nstlog                	= 1
>>>>>> ns_type                	= grid
>>>>>> pbc                	= xyz
>>>>>> rvdw                	= 3.24
>>>>>> rlist                	= 3.6
>>>>>> tcoupl                	= no
>>>>>> tc-grps                	= System
>>>>>> tau_t                	= 2.0
>>>>>> ref_t                	= 127.2717
>>>>>> vdwtype                	= user;Shift
>>>>>> rcoulomb        	= 3.6;2.24;1.12
>>>>>> coulombtype        	= Cut-off
>>>>>> rvdw-switch        	= 0.0
>>>>>>
>>>>>> energygrps        	= bead
>>>>>> energygrp_table        	= bead bead
>>>>>>
>>>>>>
>>>>>> Thanks in advance
>>>>>> Anja
>>>>>>