[gmx-users] harmonic angle potential

[Grace Brannigan]

>In calculating the potential, gromacs takes the difference betwen the
>actual angle and the reference angle value in radians,  not in degrees.
>Maybe there is the difference.

>Marcos.

Yes, that occurred to me too (is it spelled out in the manual?) but that's
not the problem. I'm working in radians.

This is just sort of mystifying. Again, my program and g_angle agree on
the average angle (~166degrees).
If you plug that angle into
V = K/2 * (166 - 180)^2, where
K = 5 * (6 angles) * (pi rad/ 180 deg)^2

you get V = 0.9
but Gromacs tells me Vbond = 1.1

(reduced units)
Could it have anything to do with my constraints on the bond length?

Here are my .top and .gro files:

/*******************************************************************/
conf.gro:

Two Chains
   10
    1LIPG    PG    1   0.085   2.523   3.385 -2.6878  7.6361  9.8849
    1LIPG    SG    2   0.786   3.236   3.395  4.5329  8.0346 -1.3729
    1LIPG     G    3   1.351   4.046   3.258 -2.6881 -4.2935  7.4651
    1LIPG     G    4   2.181   4.593   3.152 13.5878  3.3816  3.9792
    1LIPG     G    5   3.006   5.130   2.963  4.6537 -3.6607  4.5302
    2LIPG    PG    6   3.619   6.709   3.556 -8.2967  4.0282 -0.0105
    2LIPG    SG    7   2.873   6.049   3.548  0.6960 -8.3603  2.5179
    2LIPG     G    8   2.155   5.375   3.707  6.6096  5.3158 -1.0018
    2LIPG     G    9   1.260   4.941   3.600 -2.3264  6.9400  1.5371
    2LIPG     G   10   0.472   4.361   3.387 -4.9276 -6.0048  0.1470
  20.00000  20.00000  20.00000

/********************************************************/
topol.top

;
;	File 'topol.top' was generated
;	By user: grace (508)
;	On host: centaur
;	At date: Fri Sep 17 13:47:11 2004
;
;	This is your topology file
;	Two Chains
;
; Include forcefield parameters
#include "ffgmx.itp"

[ moleculetype ]
; Name            nrexcl
Protein              5

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB
   chargeB      massB
     1         PG      1   LIPG     PG      1          0     	1   ; qtot 0
     2         SG      1   LIPG     SG      1          0     	1   ; qtot 0
     3          G      1   LIPG      G      1          0     	1   ; qtot 0
     4          G      1   LIPG      G      1          0        1   ; qtot 0
     5          G      1   LIPG      G      1          0        1   ; qtot 0
     6         PG      2   LIPG     PG      1          0     	1   ; qtot 0
     7         SG      2   LIPG     SG      1          0     	1   ; qtot 0
     8          G      2   LIPG      G      1          0     	1   ; qtot 0
     9          G      2   LIPG      G      1          0        1   ; qtot 0
    10          G      2   LIPG      G      1          0        1   ; qtot 0

[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
    1     2     1    1	100000
    2     3     1    1	100000
    3     4     1    1  100000
    4     5     1    1  100000
    6     7     1    1  100000
    7     8     1    1  100000
    8     9     1    1  100000
    9    10     1    1  100000
[ angles ]
;  ai    aj    ak funct            c0            c1            c2         
  c3
    1     2     3     1    180	5
    2     3     4     1    180  5
    3     4     5     1    180	5
    6     7     8     1    180  5
    7     8     9     1    180  5
    8     9    10     1    180  5

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include water topology
#include "spc.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include generic topology for ions
#include "ions.itp"

[ system ]
; Name
Monte Carlo Simulation of Lipids

[ molecules ]
; Compound        #mols
Protein             1

On Monday, September 20, 2004, at 09:23 AM, gmx-users-request at gromacs.org
wrote:

Message: 1
Date: Mon, 20 Sep 2004 08:38:23 -0300
From: Marcos Villarreal <arloa at dqb.fcq.unc.edu.ar>
Subject: Re: [gmx-users] harmonic angle potential
To: <gmx-users at gromacs.org>
Message-ID: <200409200838.23216.arloa at dqb.fcq.unc.edu.ar>
Content-Type: text/plain;  charset="iso-8859-1"

On Sunday 19 September 2004 17:56, Grace Brannigan wrote:
Hi,

I'm trying to use Gromacs for simulations of short Lennard-Jones type
chains. I already have home-made MC code. When I compare Gromacs'
assesment of the potential energy with that of my own code, for one
configuration of two 5mers, the values of the bond angle potential
disagree.

I was looking through the Gromacs source code to make sure I had
consistent potentials, and I noticed that in the function "angles", the
sum goes over the number of bonds. I guess this makes sense for the
forces, but it seems weird for the potential; it seems like it should go
over the number of angles (i.e., num_bonds - 1). Is it possible that one
angle is being counted twice? Maybe I'm misunderstanding how this function
works...

Has anybody else had this problem? I'm hoping there's a simple thing I'm
doing wrong but I haven't found it in my own code.

-Grace

In calculating the potential, gromacs takes the difference betwen the actual
angle and the reference angle value in radians,  not in degrees. Maybe there
is the difference.

Marcos.