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