From: Axel Kohlmeyer (akohlmey_at_gmail.com)
Date: Tue Jul 05 2016 - 08:37:15 CDT

On Mon, Jul 4, 2016 at 5:09 PM, Peter Freddolino <petefred_at_umich.edu> wrote:
> Why do you think the "force field" is proportional to 1/sin(a)? The magnitude of the force is proportional to (a-a0), as should be clear from the form of the potential.
> Best,
> Peter

peter,

it is not the "force field", but the fore projection on the atoms. for
the case where the two bonds forming an angle are (near) parallel, the
standard expression diverges. in fact, for exactly parallel bonds, one
could pick a random direction to bend the molecule.

to answer bernard's question. if you look into the source code, you'll
see that NAMD will use an approximation for sin(theta) < 1.0e-6 that
is suitable for theta equals 0 or pi (n.b. diff = theta - theta0 and
normal = 1 for CHARMM)

  // Calculate constant factor 2k(theta-theta0)/sin(theta)
  BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
  if (normal != 1) {
    diff *= (2.0* k);
  } else if ( sin_theta < 1.e-6 ) {
    // catch case where bonds are parallel
    // this gets the force approximately right for theta0 of 0 or pi
    // and at least avoids small division for other values
    if ( diff < 0. ) diff = 2.0 * k;
    else diff = -2.0 * k;
  } else {
    diff *= (-2.0* k) / sin_theta;
  }
  BigReal c1 = diff * d12inv;
  BigReal c2 = diff * d32inv;

  // Calculate the actual forces
  Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
  Force force2 = force1;
  Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
  force2 += force3; force2 *= -1;

another approach employed in other codes is to not allow sin_theta to
become smaller than a given value, say 1.0e-3 and thus curb the effect
of dividing by a few small number.

axel.

>
>
>
>> On Jul 3, 2016, at 11:56 PM, Mr Bernard Ramos <bgrquantum_at_REMOVE_yahoo.com> wrote:
>>
>>
>> Hi all!
>>
>> The potential for the angle-vibration potential in CHARMM is V = k(a - a0)^2. Its corresponding force field is proportional to 1/sin(a) which blows up when a approaches \pi or when the molecule is in linear arrangement. I am wondering how NAMD calculates the angle-bending potential of a linear molecule such as carbon dioxide when it the force vector is singular at a equal to \pi?
>>
>> Thank you.
>>
>> Bernard
>
>

-- 
Dr. Axel Kohlmeyer  akohlmey_at_gmail.com  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.