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

[Axel Arnold]

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