angle harmonic with large angles

Dear all,

I'm simulating a simple spring-and-mass network under shear. I make use of longitudinal and angular springs among the atoms.

Sometimes the angle between two bonds is ~180 degrees, and in this case I am concerned about the way angle_harmonic.cpp evaluates the angle displacement:

dtheta = acos(c) - theta0[type];

Since acos(x)=[0..pi), an angle stretched from 179 to 181 degrees (or generally from 180-x to 180+x) will have dtheta=0, thus no force and no energy?

I cannot see any advice about this in the documentation, so I'm wondering if there's something I'm missing.

Thanks in advance to any reply,

Roberto

Am I missing some stupid detail? in that case please enlighten me!

Can linear molecules be simulated with harmonic angular springs in lammps? How?

Best,
Roberto

applying an angle potential that is hinged at the central atom to a near linear chain of atoms is a difficult issue. in case of a 180 degree angle, the direction of forces is undefined and close to it severely impacted by numerical noise and the projection of forces difficult. so angle style harmonic or anything similar is not a good choice for that case.

axel.

applying an angle potential that is hinged at the central atom to a near linear chain of atoms is a difficult issue. in case of a 180 degree angle, the direction of forces is undefined and close to it severely impacted by numerical noise and the projection of forces difficult. so angle style harmonic or anything similar is not a good choice for that case.

axel.

Dear Roberto

For what it’s worth, I have been using angle style harmonic to model long linear chains whose persistence length is up to 25 monomers. (equivalently, K = k/2 up to 25kBT). I have not noticed any issues with this yet. (Although this doesn’t mean there is not a problem with my method.)

I suggest that you run a simulation under the final conditions you intend to use (langevin temperature and damping parameters, etc). Then measure the final temperature of the system as well as the persistence length of your polymer. If they look okay, then the issue can be ignored. (There are so many issues that could in principal effect accuracy that I see no reason to choose one to focus on.)

(I will do the same.)

Thanks for the post, by the way. I did not think about this issue.

Andrew

Dear Axel and Andrew,

in my opinion the main problem is not the singularity at 180 degrees, but the fact that in the present implementation the potential has two minina, one at theta0, and a second at 2*pi-theta0. In some cases this can lead to unphysical dynamics. It should at least be documented.

Couldn’t a better definition be theta=atan2(r1 x r2, r1 * r2), where r1xr2 and r1*r2 are the cross product and dot product between the bond vectors?

That would at least remove the bistability.

I could provide a patch to bond_harmonic.cpp if you agree with this solution.

Best,
Roberto