[ESPResSo-users] Re: Time step question

[Michael Winokur]

Hi All,

I've done a bit more assessment and I believe I've identified the source of the problem with anomalous behavior between the time step and the kinetic energy.  The problem lies within the rotational energy portion of the code and any systematic drop in the rotational kinetic energy with decreasing time step is eventually reflected in an anomalous drop in the translational kinetic energy after equilibration.  

My local code modification was to extend the energy calculation to particles with an anisotropic moments of inertia and I pretty much assumed that the given form of the source code was correct.  However it is also possible that I introduced something that inadvertently broke the rotational motion pieces of Espresso.  My guess is that this isn't used by many Espresso users so they aren't overly interested or impacted by this issue.  

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.

Thanks,

Michael

On Sun, Mar 20, 2011 at 9:57 PM, Michael Winokur <address@hidden> wrote:

---------- Forwarded message ----------
From: Ulf Schiller <address@hidden>
Date: Fri, Mar 18, 2011 at 3:57 AM
Subject: Re: Time step question
To: "address@hidden" <address@hidden>
Cc: Michael Winokur <address@hidden>


Michael,

Michael Winokur wrote:
> So after letting my Espresso project lay idle for awhile I've picked up
> the pieces of a connected Gay-Berne particle interaction simulation.
> Briefly: I'm in the process of tuning the time stepping to prevent
> unphysical results.  I immediately noticed that after halving the time
> step the simulation's average kinetic energy per particle actually
> dropped even though I kept the temperature set point constant.  I
> haven't looked into the source code but this suggests something isn't
> quite right.

What kind of thermostat are you using?

> Does anyone have any suggestions as to how I should go about trouble
> shooting this?

I'd suggest testing with a single particle whether the
fluctuation-dissipation relation checks out. Then with two particles,
and so on. This will tell you whether it is an issue with the thermostat
or the particle interaction.

Cheers,
Ulf

--
Dr. Ulf D. Schiller                         Building 04.8, Room 231a
Institute of Complex Systems (ICS-2)        Phone:  +49 2461 61-1765
Forschungszentrum Jülich, Germany           Fax:    +49 2461 61-2850


------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
Forschungszentrum Juelich GmbH
52425 Juelich
Sitz der Gesellschaft: Juelich
Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
Vorsitzender des Aufsichtsrats: MinDirig Dr. Karl Eugen Huthmacher
Geschaeftsfuehrung: Prof. Dr. Achim Bachem (Vorsitzender),
Dr. Ulrich Krafft (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt,
Prof. Dr. Sebastian M. Schmidt
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------

Hello Ulf and Espresso uses,

Thanks for the suggestions and quick reply.  I'm using the Langevin thermostat.  I can reduce the number of particles to two but having a single particle breaks the tcl front end and that would be hard for me to fix. 

With two particles I obtain the kinetic energy two ways:

1. Through the call "set energies [analyze energy]"

2. Directly in tcl code:

    for {set k 0} { $k < $num_mol } {incr k} {
    for {set i $start($k)} { $i <= $stop($k) } {incr i 1} {
    foreach {vx vy vz} [part $i print v] break           
    set mass [part $i print mass]
    set energy [expr $energy+0.5*$mass*($vx*$vx+$

vy*$vy+$vz*$vz)]
    foreach {vx vy vz} [part $i print omega] break           
    foreach {ix iy iz} [part $i print rinertia] break           
    set renergy [expr $renergy+0.5*($ix*$vx*$vx+$iy*$vy*$vy+$iz*$vz*$vz)]
      }
     }

Both give the same answers and show, with 2 particles, what seems to be a linear response to the magnitude of the time step.  Note that I have significantly modified the code (I don't know if my revisions have made it into the regular distribution) to include an anisotropic moment of inertia for non spherical particles.  This effects both rotational and translation motions.  I note that there where some seemingly cryptic scaling of the mass and the time step in the source code to make things "scale invariant".  I do not believe I broke them when added the rotational pieces.

Here is the original code in energy.h:
*/
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   
}

Here is my revised c code:

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

#ifdef ROTATIONAL_INERTIA

  /* the rotational part is added to the total kinetic energy;

     Here we use the rotational inertia  */

  energy.data.e[0] += (   SQR(p1->m.omega[0])*p1->p.

rinertia[0]
              + SQR(p1->m.omega[1])*p1->p.rinertia[1]
              + SQR(p1->m.omega[2])*p1->p.rinertia[2])*time_step*time_step;
//printf("energy: %d %f %f \n",p1->p.identity,p1->m.omega[0],p1->p.rinertia[1]);
#else

  /* the rotational part is added to the total kinetic energy;

     at the moment, we assume unit inertia tensor I=(1,1,1) why isn't there the 0.5 factor  */

  energy.data.e[0] += (SQR(p1->m.omega[0]) + SQR(p1->m.omega[1]) + SQR(p1->m.omega[2]))*time_step*time_step;
#endif

#endif   
}

I'm really not certain why time_step isn't used for translation motion while it is used in rotational motion except to refer to the assumption of some sort of scale invariance in the case of translational motion.

I note that my _expression_ was seemingly necessary to yield the classical equipartition result of 1/2 k_b T per degree of freedom.   I should have caught the time step issue earlier....


Thanks again for any suggestions.

Michael