[lammps-users] constant torque for rigid bodies

Dear lammps users

Is it possible to rotate rigid spheres with constant velocity (torque)? The problem is that I need to measure the (orientational) averaged interaction forces between 2 rough nanospheres at fixed distances in vacuum. Using fix rigid implies an NVE ensemble, i.e. in order to keep the total energy fixed, near a potential energy minimum the spheres rotate fast while close to the maxima they rotate slow which has no physical meaning. What can I do to receive a constant torque while measuring the binding energies?

Thank you for your help.

regards

Sabine

infile:

velocity all create 300 12348765

velocity all zero linear

neighbor 5.0 bin

neigh_modify delay 2

group sphere1 id <> 1 2808

group sphere2 id <> 2809 5616

velocity sphere2 zero linear

velocity sphere1 zero linear

fix 1 all rigid group 2 sphere1 sphere2 force * off off off

compute 1 sphere1 group/group sphere2

thermo 100

thermo_style custom step etotal ke c_1 c_1[1] c_1[2] c_1[3]

timestep 3

run 100000

Constant velocity and torque are not the same thing. I assume
you mean omega (rotational velocity) or angular momentum.
If you turn off torque with the fix rigid option, then the rigid
bodies will rotate with constant angular momentum. Is
that what you want?

From the fix rigid doc page:

By default, each rigid body is acted on by other atoms which induce an
external force and torque on its center of mass, causing it to
translate and rotate. Components of the external center-of-mass force
and torque can be turned off by the force and torque keywords. This
may be useful if you wish a body to rotate but not translate, or vice
versa, or if you wish it to rotate or translate continuously
unaffected by interactions with other particles

Steve

Dear Steve,

Thank you for the correction I confused the term torque with angular momentum, sorry. I have redone the simulations now fixing the rotational energies, but there are many things in the results I cannot understand.

  1. I use compute group/group to calculate the potential energies and pair forces between the two particles. I would expect that doing it for several constrained particle positions, the integral over the forces along their connecting axis should lead to the same free energy as given by the measured potential energies, but this is by far not the case.

  2. I think my calculations make no sense insofar as the NVE ensemble attempts to conserve the total energy, which causes problems when I fix the kinetic energies. Maybe a certain energy state of the system is then measured as it gives in average (sampled over several 100000 conf. at each constrained particle distance) completely different forces and pot. energies depending on whether I rotate the bodies with constant or changing velocities. This should not be the case if I calculated what I want: just to measure the interparticle forces and pot. energies while I rotate the particles. Is it possible at all to calculate the interaction energies for certain orientations, then to change the orientations (using run not to turn them by hand to give lammps a new data-file) and redo the calculation and so on till one has several 100000 measurement?

  3. The system (2 rigid bodies) is for sure to small to have a stable system. Would it be a way out to increase the number of particles, to constrain two of them and use the group/group command to calculate the forces and energies.

I hope this all was not to confusing but I am gracefull for any hint you can provide.

Best regards and thank you very much.

Sabine

Steve Plimpton [email protected] 8/25/2010 4:25 PM >>>
Constant velocity and torque are not the same thing. I assume
you mean omega (rotational velocity) or angular momentum.
If you turn off torque with the fix rigid option, then the rigid
bodies will rotate with constant angular momentum. Is
that what you want?

From the fix rigid doc page:

By default, each rigid body is acted on by other atoms which induce an
external force and torque on its center of mass, causing it to
translate and rotate. Components of the external center-of-mass force
and torque can be turned off by the force and torque keywords. This
may be useful if you wish a body to rotate but not translate, or vice
versa, or if you wish it to rotate or translate continuously
unaffected by interactions with other particles

Steve

Comments below.

Steve

Dear Steve,

Thank you for the correction I confused the term torque with angular
momentum, sorry. I have redone the simulations now fixing the rotational
energies, but there are many things in the results I cannot understand.

1) I use compute group/group to calculate the potential energies and pair
forces between the two particles. I would expect that doing it for several
constrained particle positions, the integral over the forces along their
connecting axis should lead to the same free energy as given by the
measured potential energies, but this is by far not the case.

If you're saying the interactions energy/force between 2 particles should
be the same for any relative orientation of the 2 particles, I don't
see why that would be the case.

2) I think my calculations make no sense insofar as the NVE ensemble
attempts to conserve the total energy, which causes problems when I fix the
kinetic energies. Maybe a certain energy state of the system is then
measured as it gives in average (sampled over several 100000 conf. at each
constrained particle distance) completely different forces and pot. energies
depending on whether I rotate the bodies with constant or changing
velocities. This should not be the case if I calculated what I want: just to
measure the interparticle forces and pot. energies while I rotate the
particles. Is it possible at all to calculate the interaction energies for
certain orientations, then to change the orientations (using run not to turn
them by hand to give lammps a new data-file) and redo the calculation and so
on till one has several 100000 measurement?

By zeroing the torque on the 2 particles (so they rotate freely) you are
certainly messing with energy conervation when the overall system is
doing NVE. It is possible to use a loop in a LAMMPS input script
and/or generate configs in different orientations using some external
program and run the set of configs thru LAMMPS. Look at the
"loop" command doc in LAMMPS or you could write a shell or Python
script to iterate over many short LAMMPS runs.

3) The system (2 rigid bodies) is for sure to small to have a stable system.
Would it be a way out to increase the number of particles, to constrain two
of them and use the group/group command to calculate the forces and
energies.

You can certainly constrain some subset of particles and you get
to define the 2 groups for group/group, so they could include
whatever atoms you want.