## VMD-L Mailing List

**From:** Axel Kohlmeyer (*akohlmey_at_gmail.com*)

**Date:** Tue Jul 05 2016 - 08:37:15 CDT

**Next message:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Previous message:**Fernando Vallejos-Burgos: "Re: I want to change back ground color"**In reply to:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Next in thread:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Reply:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]

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.

**Next message:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Previous message:**Fernando Vallejos-Burgos: "Re: I want to change back ground color"**In reply to:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Next in thread:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Reply:**Peter Freddolino: "Re: On CHARMM angle bending potential"**Messages sorted by:**[ date ] [ thread ] [ subject ] [ author ] [ attachment ]