[gmx-users] Tabulated potential, derivative errors and interpolation
ms
devicerandom at gmail.com
Thu Jan 7 20:02:58 CET 2010
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 :)
I noticed that in the rtab14 (which should be the default LJ potential)
the values are, instead, apparently all negative, as expected.
If you need more information, feel free to ask.
Thank you very much,
M.
More information about the gromacs.org_gmx-users
mailing list