Rigid molecule with multiple correlated constraints

[CHANG YUN SON]

Is it possible to use constraints to fix a molecule completely rigid and planar?
I'm testing a polarizable nitrate force field (NO3-) for certain ionic liquids.
For non-polarizable ones I can get away with putting some improper dihedrals, but the polarizable FF was pulling the dihedral always too strongly which eventually leads into polarization catastrophe, so I wanna just fix the geometry of the ion.

I initially tried to modify the forcefield.py in the python wrapper to put bond and angle constraints for all the bonds and angles in it, and exclude the impropers of this molecule. I think my system class is well-defined, but when I set the velocities for the system using setVelocitiesToTemperature, it does not assign a value for nitrate ions. This results in nan for both the velocities and kinetic energy of those particles, so the system crashes. I tried to dig into more details in the source code but it was a bit hard to figure out how the coupled constraints are handled with ccma.

Thanks,
ChangYun

[Peter Eastman]

It's possible, but you have to be a bit careful about how you do it. Constraints work by applying a force along the constrained direction, then solving for the values of all constraint forces needed to keep the distances from changing. That doesn't work for planar molecules. All the constraints lie in a plane, so none of them can apply a force out of the plane. The workaround is to add an extra "dummy particle" out of the plane. Set its charge and epsilon to 0 so it doesn't interact with anything. You can then add constraints between that particle and all the real particles, which allows constraint forces to be applied out of the plane.

The next problem is that you need to be careful not to add too many constraints. When making a molecule completely rigid, its easy to add more constraints than you need, which causes the constraint matrix to be singular and everything blows up. For example, consider a methane molecule. It has five atoms, so 15 degrees of freedom. To make it fully rigid, you need to constrain 9 of them so you're left with 6 degrees of freedom, the number that a rigid body has (3 translations and 3 rotations). But if you naively add a constraint between every atom and every other, that's 10 constraints. The constraints are overdetermined and the matrix is singular. You need to omit one of them to get the correct number of constraints.

Peter

[CHANG YUN SON]

Thank you for the suggestion of adding dummy particle out of plane. I agree that it's more robust way of constraining planar molecule. Although I still don't understand why the applying constraints to velocities would not assign a value even during initialization. In practice, constrained distances often deviate slightly from the equilibrium distance, and there should be some out of plane force as well. Here I'm setting six constraints for every atom pair distances, and this is the right number of constraints to keep the molecule rigid (4*3-6=6). This set of constraints work okay with SHAKE or even LINCS in GROMACS, so I'm not sure why the cmma module would not initiate the velocities at all. Do you have any idea on this?

ChangYun

[Peter Eastman]

When you call setVelocitiesToTemperature(), it first initializes all the velocities to random values, and then applies constraints to them to remove any component of velocity along a constrained direction. If the matrix is singular, that can produce nans.

SHAKE only deals with one constraint at a time, so if the matrix is overdetermined it just doesn't notice. (But if the constraints are inconsistent, it won't converge.) CCMA starts by explicitly constructing and inverting the full matrix, so it's more sensitive to problems of this sort.

Peter

[CHANG YUN SON]

Another issue with adding a dummy particle for the out of plane constraint is that one cannot assign constraints for massless particles. So if I do the constraints in this way, I have to set a small mass for the dummy particle which is a bit strange and deviates from the actual mass distribution of the molecule.

And do you think that the singularity of the ccma matrix comes from the nature that the molecule is planar? I haven't gone through the exact math of the ccma algorithm, but I thought given that the number of constraints are consistent with the molecular degree of freedom, it should be fine. Is there any mathematical reason why planar molecule would make singularity in the matrix even with the right number of constraints?

[Peter Eastman]

Think about it this way. You're solving a system of equations to find the values of all constraint forces that give the correct distances. But if the constraint forces all lie in a plane, and an out-of-plane force is needed to get the correct distances, then that's impossible. There does not exist any set of values you could pick that would solve the equations. Mathematically, that means the matrix is singular.

Here's another way to think about it. The whole calculation of degrees of freedom is different in 2D than in 3D. For a planar molecule moving in a plane, you have 2N degrees of freedom, and you need to remove all but 3. For a 3D molecule moving in space, you have 3N degrees of freedom and you need to remove all but 6. So it's a different problem.

Peter

[CHANG YUN SON]

Thanks for the explanation.
Yes I thought it might not satisfy the dof in 2D, and it makes sense now.

I'll probably play around with different strengths of impropers and/or some bond constraints without dummy particles (ex, constraining O-O bond lengths and putting non-constrained bonds for the internal N-O distance) since I don't want to change the inertia of the molecule.

Thanks,
ChangYun