[lammps-users] minimize with rigid water

Dear LAMMPS users,

The LAMMPS manual indicates that fix rigid is not supported by minimization.
Has there been any work on developing code for minimization with rigid
bodies? In theory, it is simple to minimize with variables of center of
mass and orientation, but I imagine it could be difficult to implement in
LAMMPS where the atomic representation is use. This is of great interest to
me, because right now LAMMPS can't do inherent structures of common water
models, SPC/E, for example.

Regards,
Harold

no work on this to my knowledge - if someone wants
to implement some code to do this - e.g. in fix rigid,
called by the minimizer - then that would be great.

Steve

Dear Steve (or others?),

Could you give me a hint as to how to call fix rigid within the minimizer I am writing? I have my own minimizer written already, so setting that up isn't the problem (it does enthalpy minimization).

That would be great if I could use the fix_rigid code already written to do two things:
1. return center of mass positions, euler angles, forces, torques of all rigid bodies
2. given new center of mass positions and euler angles, return forces, torques and total system pressure

I almost started coding my own version of fix_rigid within the minimizer, because with the atomic coordinates and forces, I can get all the rigid body info I need (except system pressure, ouch!). But this would take me some time, and my minimization code would only work for my specific system.

Thank you!
Harold

The existing minimizers can call a fix, like fix box/relax,
so there is an interface defined for that. The real question
is how to maintain a rigid body during a minimization. I believe
this is non-trivial to do correctly, since you really need to
do constrained minimization, else your energy and gradients
will be inconsistent. Fix rigid does not use Euler angles;
it uses quaternions. Euler angles are inaccurate for time
integration.

Steve

in addition to steve's remarks. i know several codes that
"allow" constraints during minimization, but they silently
replace them with very stiff harmonic potentials. since
you don't do dynamics, they are much less of a problem
and should give pretty much the same result.

if you keep in mind that constraints are somewhat of an
approximation in the first place, then this isn't such bad
thing at all.

cheers,
   axel.

Dear Steve and Axel,

Thank you very much for the tips.

I plan to avoid constrained minimization altogether. For something as small as water, it is simple to take the degrees of freedom of the water molecule as the center of mass position and euler angles (transformed from quaterions if necessary), as opposed to having atomic degrees of freedom + constraints. Then, the objective function is simply U(center of mass positions, orientations), and derivatives are the forces and torques of the molecules. Of course, if I wanted to constrain hydrogens on a polymer, for example, I'd have to start worrying about constrained minimization.

Best,
Harold

One problem you will have is getting the minimizer
itself to not update the positions of the atoms in the
rigid body. Since it controls the 3N degrees of freedom
it will want to do that. The DOF you want from fix rigid
are not really additional DOF, but replace some of the 3N
DOF.

Steve

I agree about the problem of stopping the minimizer from using the full 3N degrees of freedom, which doesn't allow me to use the current min styles as an interface. Instead, I made my own min style from min_cg.cpp, using different working vectors, etc. Unfortunately it is a messy design, and the current implementation probably won't be of much use to the LAMMPS community in general. =(