espressomd-users
[Top][All Lists]
Advanced

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

[ESPResSo-users] Re: Time step question


From: Ulf Schiller
Subject: [ESPResSo-users] Re: Time step question
Date: Sun, 20 Mar 2011 21:28:00 -0500
User-agent: Mozilla/5.0 (Windows; U; Windows NT 6.1; en-US; rv:1.9.2.15) Gecko/20110303 Thunderbird/3.1.9

Michael,

I am away from my desk for the next couple of weeks and unable to look into the details right now. Meanwhile, please post your replies to the mailing list. It increases the likelihood of getting a swift response. Also, it will help other people who may be interested in this topic is well.

Thanks,
Ulf

On 3/18/2011 12:48 PM, Michael Winokur wrote:
Hello Ulf,

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:

     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, 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 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
   /* 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
}

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


On Fri, Mar 18, 2011 at 3:57 AM, Ulf Schiller <address@hidden
<mailto:address@hidden>> wrote:

    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
    <tel:%2B49%202461%2061-1765>
    Forschungszentrum J├╝lich, Germany           Fax: +49 2461 61-2850
    <tel:%2B49%202461%2061-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
    
------------------------------------------------------------------------------------------------
    
------------------------------------------------------------------------------------------------






reply via email to

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