[ESPResSo-users] Fwd: Time step question

[Michael Winokur]

---------- 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