espressomd-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [ESPResSo-users] Re: Time step question


From: Michael Winokur
Subject: Re: [ESPResSo-users] Re: Time step question
Date: Wed, 23 Mar 2011 23:27:40 -0500

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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]