# [gmx-users] Tabulated potential, derivative errors and interpolation

David van der Spoel spoel at xray.bmc.uu.se
Thu Jan 7 20:10:33 CET 2010

```On 1/7/10 8:02 PM, ms wrote:
> David van der Spoel ha scritto:
>> ms wrote:
>>> David van der Spoel ha scritto:
>>>> ms wrote:
>>>>> Hi,
>>>>>
>>>>> I am trying to use a tabulated potential but I have a hard time getting
>>>>> the derivatives right (apparently), which puzzles me.
>>>>>
>>>>> I tried at first to use the analytical expression for the derivative,
>>>>> and I got warnings ranging from:
>>>>>
>>>>> "WARNING: For the 2168 non-zero entries for table 1 in
>>>>> table_10_70_0.4_5_C-alpha_C-alpha.xvg the forces deviate on average
>>>>> 133%
>>>>> from minu
>>>>> s the numerical derivative of the potential
>>>>>
>>>>> WARNING: For the 257 non-zero entries for table 2 in
>>>>> table_10_70_0.4_5_C-alpha_C-alpha.xvg the forces deviate on average
>>>>> 300%
>>>>> from minus
>>>>>   the numerical derivative of the potential"
>>>>>
>>>>> to
>>>>>
>>>>> "WARNING: For the 268 non-zero entries for table 2 in
>>>>> table_10_70_0.4_1_C-alpha_C-alpha.xvg the forces deviate on average
>>>>> 120138% from minus the numerical derivative of the potential"
>>>>>
>>>>>
>>>>> At first I thought it was because I am artificially "cutting" the
>>>>> highest values of the repulsive term, putting everything above 10^17 as
>>>>> 10^17 with a derivative of 0, but playing with that didn't help; also I
>>>>> found that there are errors too in the repulsive part, which hasn't
>>>>> out-of-range values.
>>>>>
>>>>> So I thought it could have been some mysterious error in the math, and
>>>>> tried estimating it numerically, but again with a numerically-estimated
>>>>> table I have errors, and even worse ones:
>>>>>
>>>>> "WARNING: For the 257 non-zero entries for table 2 in
>>>>> table_10_70_0.4_5.xvg the forces deviate on average 120773% from minus
>>>>> the numerical derivative of the potential"
>>>>>
>>>>> I used mdrun with -debug 1 to see what tables the guy is actually
>>>>> using,
>>>>> and I am also concerned, because there is something weird which could
>>>>> explain the problem: why are *positive* values routinely interspersed
>>>>> with the right *negative* values in the derivative?
>>>>>
>>>>> For example, this is an excerpt of rtab.xvg:
>>>>>
>>>>> 3.0992913246e-01  2.7511889648e+02  7.0945126953e+03
>>>>> 3.1012925506e-01  2.7884115601e+02  -1.2228130859e+04
>>>>> 3.1032931805e-01  2.7019714355e+02  -6.7522789062e+04
>>>>> 3.1052938104e-01  2.5448614502e+02  -8.2872140625e+04
>>>>> 3.1072947383e-01  2.3970016479e+02  -5.8276152344e+04
>>>>> 3.1092956662e-01  2.3383134460e+02  6.2780151367e+03
>>>>> 3.1112965941e-01  2.3702958679e+02  -1.0651265625e+04
>>>>> 3.1132972240e-01  2.2957331848e+02  -5.8167523438e+04
>>>>> 3.1152981520e-01  2.1604493713e+02  -7.1314695312e+04
>>>>> 3.1172990799e-01  2.0332643127e+02  -5.0073011719e+04
>>>>> 3.1192997098e-01  1.9829878235e+02  5.5352783203e+03
>>>>> 3.1213006377e-01  2.0104492188e+02  -9.2590859375e+03
>>>>> 3.1233015656e-01  1.9461547852e+02  -5.0073472656e+04
>>>>> 3.1253021955e-01  1.8298132324e+02  -6.1317703125e+04
>>>>> 3.1273031235e-01  1.7205114746e+02  -4.2998656250e+04
>>>>> 3.1293040514e-01  1.6774539185e+02  4.8916777344e+03
>>>>> 3.1313049793e-01  1.7010066223e+02  -8.0468198242e+03
>>>>> 3.1333056092e-01  1.6456481934e+02  -4.3058406250e+04
>>>>> 3.1353065372e-01  1.5456475830e+02  -5.2672746094e+04
>>>>> 3.1373074651e-01  1.4518064880e+02  -3.6885503906e+04
>>>>> 3.1393080950e-01  1.4149731445e+02  4.2984921875e+03
>>>>> 3.1413090229e-01  1.4351583862e+02  -6.9780966797e+03
>>>>> 3.1433099508e-01  1.3875308228e+02  -3.6990878906e+04
>>>>> 3.1453105807e-01  1.3016751099e+02  -4.5200566406e+04
>>>>> 3.1473115087e-01  1.2211988831e+02  -3.1613373047e+04
>>>>> 3.1493124366e-01  1.1897129822e+02  3.7776948242e+03
>>>>>
>>>>> Thank you for help!
>>>>> Massimo
>>>> Which gromacs version?
>>>>
>>>
>>> Sorry! 4.0.5
>>>
>> I can not reproduce this. Can you give exact command line and details
>> about the system? vdwtype etc.?
>>
>>
>
> The system is a single peptide-like polymer in vacuum. I am trying to
> build a CG model.
> That said, the command line is (as launched from a python script I use):
>
> /home/ms872/software/bin/mdrun -v -s
> \$STARTDIR/'''+top_name+'_fullmd_'+jobname+'''.tpr -o
> \$STARTDIR/'''+top_name+'_'+jobname+'''.trr -x
> \$STARTDIR/'''+top_name+'_'+jobname+'''.xtc -c
> \$STARTDIR/'''+top_name+'_final_'+jobname+'''.gro -e
> \$STARTDIR/'''+top_name+'_minbox_'+jobname+'''.edr -table '''+table+'''
> -tablep '''+tablep+''' 2>&1>  \$STARTDIR/output_'''+jobname+'''.txt -debug 1
>
> The .mdp is below:
> (note: the email continues after the .mdp)
> ------------------------------
> title                    = MD simulation of xaa10
>
> ;Preprocessor
>
> cpp                      = /lib/cpp
>
> ;Directories to include in the topology format
>
> include                  = -I../top
>
> ;Run control: stochastic dynamics
>
> integrator               = sd
>
> ;Total simulation time: 10,000 ps
>
> :time step in picoseconds
>
> dt                       = 0.002
>
> ;number of steps
>
> nsteps                   = 5000000
>
> ;nsteps                  = 15000000
>
> ;frequency to write coordinates to output trajectory file
>
> nstxout                  = 1000
>
> ;frequency to write velocities to output trajectory file
>
> nstvout                  = 1000
>
> ;frequency to write energies to log file
>
> nstlog                   = 200
>
> ;frequency to write energies to energy file
>
> nstenergy                = 200
>
> ;frequency to write coordinates to xtc trajectory
>
> nstxtcout                = 1000
>
> ;group(s) to write to xtc trajectory
>
> xtc_grps                 = protein
>
> ;group(s) to write to energy file
>
> energygrps               = C-alpha
>
> ;Frequency to update the neighbor list (and the long-range forces,
>
> ;when using twin-range cut-off's).
>
> nstlist                  = 10
>
> ;Make a grid in the box and only check atoms in neighboring grid cells
>
> ;when constructing a new neighbor list every nstlist steps.
>
> ns_type                  = grid
>
> cut-off distance for the short-range neighbor list
>
> rlist                    = 0.8
>
>
> ;treatment of electrostatic interactions
> coulombtype              = cut-off
> rcoulomb                 = 1.4
>
> ;treatment of van der waals interactions
> vdwtype                  = user
> rvdw                     = 1.4
>
> ;table stuff
> energygrp                = C-alpha
> energygrp_table          = C-alpha C-alpha
>
> ; Periodic boudary conditions in all the directions
> pbc                      = xyz
> ;Temperature coupling
> tcoupl                   = v-rescale ;(unused with sd integrator)
> tc-grps                  = protein
> tau_t                    = 2    ;see manual
> ref_t                    = 300
> ;Pressure coupling - ignore without solvent
> Pcoupl                   = no
> ;Pcoupltype               = isotropic
> ;tau_p                   = 1.0
> ;compressibility         = 4.5e-5
> ;ref_p                   = 1.0
> ;Velocity generation
> gen_vel                  = yes ;not meaningful with sd integrator
> gen_temp                 = 300
> gen_seed                 = 173529
> ;Constrain all bonds
> constraints              = all-bonds
>
> ;comm_mode               = angular
>
>
> define                   = -DMIXMORSE
> ---------------------
>
>
> I sent another mail where I attached a table input file I use and the
> corresponding full rtab.xvg -it is currently waiting mod approval
> because of the size :)

Let's try to see what's going on without it...

>
> I noticed that in the rtab14 (which should be the default LJ potential)
> the values are, instead, apparently all negative, as expected.
This is what I saw. Could it be that you specify the table with the
wrong sign for the derivative?

>
> If you need more information, feel free to ask.
>
> Thank you very much,
> M.
>
>

--
David.
________________________________________________________________________
David van der Spoel, PhD, Professor of Biology
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

```