
From:  Michael Winokur 
Subject:  [ESPResSousers] Re: Time step question 
Date:  Tue, 22 Mar 2011 11:40:39 0500 
Hello Ulf and Espresso uses, 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,
What kind of thermostat are you using?
Michael Winokur wrote:
> So after letting my Espresso project lay idle for awhile I've picked up
> the pieces of a connected GayBerne 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.
I'd suggest testing with a single particle whether the
> Does anyone have any suggestions as to how I should go about trouble
> shooting this?
fluctuationdissipation 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 (ICS2) Phone: +49 2461 611765
Forschungszentrum Jülich, Germany Fax: +49 2461 612850


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


2. Directly in tcl code:
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]"
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)]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.
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)]
}
}Here is my revised c code:
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
}#ifdef ROTATIONAL_INERTIA
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;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 */#endifenergy.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
[Prev in Thread]  Current Thread  [Next in Thread] 