Re: [ESPResSo-users] Time step question

[Michael Winokur]

Hi Axel,

It has taken far longer than I would have wished but I've tracked down and fixed a number of errors in my modified versions of thermostat.c, thermostat.h and, additionally, rotation.h
I have now fixed the Langenvin thermostat so that it will give the same kinetic energy for both translational motion and rotational motion when explicit masses and moments of inertia are included.  The main issue had to do with the lack of an explicit mass term in the langevin thermostat with rotational motion.  I now provide for this if the code is compiled with rotational inertia (via my myconfig.h).  However I did not add this in the case that only rotation is used and, here, it just retains the original form of the thermostat.  Thanks for pointing out the test tcl code area.

I think these are now good enough to add to the base package assuming one of the core Espresso developers vets them.

In regards to my simulation of chemically linked Gay-Berne particles there are some peculiar quirks that will require some further customization.  I note that the molecule I'm modeling forces two intramolecular GB particles to be in very close proximity (but they are not nearest neighbors).  Thus the simulation has a combination of very strong hard core  forces between these pairs and the much weaker Van der Waal intermolecular interactions.  This seems to a problem because when  I choose damping and time steps that are appropriate for the intermolecular dynamics then the intramolecular motion generates unphysical outcomes which, at higher temperatures, create large energies which must be then dissipated.  

I need to somehow bias the system to reduce the likelihood of the this problematic intramolecular interaction.  My first thought would be to artificially introduce a larger damping coefficient in the case of this intramolecular interaction.  I don't have enough experience with MD to really know what one could/should do to deal with this artifact.

Thanks again,

Michael


On Wed, Mar 23, 2011 at 11:27 PM, Michael Winokur <address@hidden> wrote:

Hi Axel,

You were correct.  My modification of the thermostat was at fault.  I'm not sure what the exact form of the rotational thermostat should be but I've looked at one paper and I've made some tentative corrections and retested the code.   The behavior is much better but I still think scaling isn't fully correct.  I'll forward the link to the paper tomorrow and additional details to see if I can get the form to be 100% correct. 

Many thanks,

Michael

On Wed, Mar 23, 2011 at 10:28 AM, Axel Arnold <address@hidden> wrote:

On Tuesday 22 March 2011 23:06:14 Ulf Schiller wrote:
> Hi Michael,
>
> On 3/22/2011 11:40 AM, Michael Winokur wrote:
> > One thing to note in the original code of energy.h is that there is an
> > overt difference in the _expression_ so that the variable time_step only
> > appears in the rotational energy term:
> > */
> > MDINLINE void add_kinetic_energy(Particle *p1)
> > {
> >
> >    /* kinetic energy */
> >    energy.data.e[0] += (SQR(p1->m.v[0]) + SQR(p1->m.v[1]) +
> >
> > SQR(p1->m.v[2]))*PMASS(*p1);
> >
> > #ifdef ROTATION
> >
> >    /* the rotational part is added to the total kinetic energy;
> >
> >       at the moment, we assume unit inertia tensor I=(1,1,1)  */
> >
> >    energy.data.e[0] += (SQR(p1->m.omega[0]) + SQR(p1->m.omega[1]) +
> >
> > SQR(p1->m.omega[2]))*time_step*time_step;
> > #endif
> > }
> >
> > If anyone can comment in regards to this I would be happy to hear it.
>
> This depends on the units the velocities and energies are stored. You
> have to make sure that the units for the translational and rotational
> parts are consistent.

Hi!

For historic reasons, the translational velocities are stored internally with
a rescaling by the time step, while the rotational degrees are not. Therefore,
these factors are correct.

There is a testcase also for the rotational thermostat in
testsuite/thermostat.tcl, which shows that the classical langevin thermostat
works also correct for inertia 1,1,1. So chances are good that something with
your additional code breaks the thermostat.

Axel

--
JP Dr. Axel Arnold
ICP, Universität Stuttgart
Pfaffenwaldring 27
70569 Stuttgart, Germany
Email: address@hidden
Tel: +49 711 685 67609