Improper angle calculation

LAMMPS LAMMPS Mailing List Mirror
1

Hi,
I have been trying to determine the exact procedure used in improper angle calculation in LAMMPS.

I have picked up a part of the code from compute_improper_local.cpp and have tried to figure out what happens - step by step

vb1x = x[atom1][0] - x[atom2][0];

vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z);
~~ vb1 = coordinates(1) - coordinates(2)

vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z);
~~ vb2 = coordinates(3) - coordinates(2)

vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2];
domain->minimum_image(vb3x,vb3y,vb3z);
~~ vb3 = coordinates(4) - coordinates(3)
I am correct in assuming that the improper angle should be the angle between the planes of (vb1,vb2) and (vb2,vb3) ??

ss1 = 1.0 / (vb1xvb1x + vb1yvb1y + vb1zvb1z);
ss2 = 1.0 / (vb2xvb2x + vb2yvb2y + vb2zvb2z);
ss3 = 1.0 / (vb3xvb3x + vb3yvb3y + vb3z*vb3z);
~~ Inverse norm_squared

r1 = sqrt(ss1);
r2 = sqrt(ss2);
r3 = sqrt(ss3);
~~ Inverse norm of the angles

c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * r1 * r3;
~~ cosine of angle between vb1 and vb3? (say) alpha

c1 = (vb1x * vb2x + vb1y * vb2y + vb1z * vb2z) * r1 * r2;
~~ cosine of angle between vb2 and vb1? (say) beta

c2 = -(vb3x * vb2x + vb3y * vb2y + vb3z * vb2z) * r3 * r2;
~~ negative cosine of angle between vb2 and vb3? (say) gamma

s1 = 1.0 - c1*c1; ~~ sin^2 of vb2 and vb1?
if (s1 < SMALL) s1 = SMALL;
s1 = 1.0 / s1;

s2 = 1.0 - c2*c2; ~~ sin^2 of vb2 and vb3?
if (s2 < SMALL) s2 = SMALL;
s2 = 1.0 / s2;

s12 = sqrt(s1s2); ~~ 1/sin(angle vb2 and vb1) * 1/sin(angle vb2 and vb3)
c = (c1c2 + c0) * s12;

~~ IS c = [ 1-* cos(alpha)*cos(beta) + cos(gamma))/ sin(alpha)/sin(beta) ]??

if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
cbuf[n] = 180.0*acos©/MY_PI;

Question 1:
I do not understand how cos_inverse of ‘c’ should be the angle between the two planes? Is there a geometric explanation that I should look at?

Question 2:
Also, while specifying the improper angles in the input files as ( i j k l ) should give the same magnitude for the improper angle as ( l j k i ) . Is that right?

Question 3:
Finally, for improper style harmonic and cvff does lammps explicitly calculate the displacement of atom “i” from the plane j-k-l? My understanding is that it just calculates the angle between the ijk and jkl planes.

Thank you,

Regards,
Ambarish

2

Different improper styles define and do things differently. I think the
code in compute imp/local was taken from improper_harmonic. You can
double check that by looking at the 2 files.

Steve

3

Hi,
Thank you for the reply.

You are correct - the calculation for improper harmonic is very similar.
My basic question (from the original email) is how does c, where c = (c1*c2 + c0)*s12, give the cosine of the angle between two planes geometrically?

Thank you,

Regards,
Ambarish

4

My basic question (from the original email) is how does c, where c = (c1*c2
+ c0)*s12, give the cosine of the angle between two planes geometrically?

I don't know - been too long since I looked at this code in detail. I would
do it by finding the normals to the 2 planes (by crossing 2 vectors in
each plane),
then take the dot product of the 2 normals to get the cosine of the angle.

If you work that out, I'm guessing it will be the same final formula
as in the code,
but the code is written in some short-hand way for efficiency. And this
was likely a piece of code inherited from some optimized implementation
of impropers long ago ...

Steve