Re: [ESPResSo] Position-dependent Langevoin thermostat

[Lorenzo Isella]

Dear All,
I am now taking a deeper look at the espresso source code (not sure
whether this should be moved to the developer mailing list, please let
me know).
I really have two questions here, one dealing
 with understanding Espresso as is and another about a feature I would
 like to add.
(1) I am having a look mainly at the implementation of the Verlet
algorithm for the integrating the equations of motion.
I am sure it works fine, but I must be misunderstanding something in
the scaling applied by Espresso. Quoting from integrate.c

   /* Integration Steps: Step 1 and 2 of Velocity Verlet scheme:
       v(t+0.5*dt) = v(t) + 0.5*dt * f(t)
       p(t + dt)   = p(t) + dt * v(t+0.5*dt)
       NOTE 1: Prefactors do not occur in formulas since we use
               rescaled forces and velocities.
       NOTE 2: Depending on the integration method Step 1 and Step 2
               cannot be combined for the translation.
    */


and
    /* Integration Step: Step 3 of Velocity Verlet scheme:
       Calculate f(t+dt) as function of positions p(t+dt) ( and
velocities v(t+0.5*dt) ) */

and finally

    /* Integration Step: Step 4 of Velocity Verlet scheme:
       v(t+dt) = v(t+0.5*dt) + 0.5*dt * f(t+dt) */

This looks like a velocity Verlet scheme
http://www.fisica.uniud.it/~ercolessi/md/md/node21.html
but if I replace f in espresso with f_phys*dt*dt/2m and v with v_phys*dt
(where the pedix "phys" denotes the unscaled, "physical" quantities),
then I end up with an extra 1/2 factor with respect to the formulation
in the link
I provide above and also the dt's give me some troubles with the units.
What am I missing?
(2) At this point I believe I was a bit ambitious with my initial
plan, so I had better start from something easier to implement, but
useful for my simulations. Better to think of a cluster which does not
change its structure as time progresses (as those generated in my
simulations), so that the number of first neighbors of each monomer is
given once for all and I can calculate it at time t=0 without even
resorting to Espresso necessarily.
Now I am quoting from thermostat

/** overwrite the forces of a particle with
    the friction term, i.e. \f$ F_i= -\gamma v_i + \xi_i\f$.
*/
MDINLINE void friction_thermo_langevin(Particle *p)
{
  extern double langevin_pref1, langevin_pref2;

  int j;
#ifdef MASS
  double massf = sqrt(PMASS(*p));
#else
  double massf = 1;
#endif

  for ( j = 0 ; j < 3 ; j++) {
#ifdef EXTERNAL_FORCES
    if (!(p->l.ext_flag & COORD_FIXED(j)))
#endif
      p->f.f[j] = langevin_pref1*p->m.v[j]*PMASS(*p) +
langevin_pref2*(d_random()-0.5)*massf;
#ifdef EXTERNAL_FORCES
    else p->f.f[j] = 0;
#endif
  }

Here is how I would like to change the noise acting on any particle:

p->f.f[j] = langevin_pref1*eta1*p->m.v[j]*PMASS(*p) +
langevin_pref2*eta2*(d_random()-0.5)*massf;

where eta1 and eta2 are two time-independent factors which vary
according to the particle identity (and nothing else) which I have
already calculated elsewhere. For instance:

id =0  eta1=0.5 eta2=sqrt(0.6)
id =1  eta1=0.7 eta2=sqrt(0.4)

and so on.

I would like to create somewhere two text files to be read (only once)
by Espresso into two arrays eta1[] and eta[2] use them as explained above.
Maybe I should add that I aim at having e.g. two (or more) sets of 50
monomers each (each set would be an aggregate; let us call them  A and
B) and having eta1[] and eta2[]  as arrays of length 50. Then I would
use the same eta1(2) correction for the monomers in A and B since the
two aggregates are identical and their corresponding monomers have the
same number of first neighbors.
What is the easiest way of implementing this? My C is not that great
and Espresso is a big project.
Many thanks

Lorenzo



2008/9/17 Christoph Junghans <address@hidden>:
> Hi Lorenzo & Olaf,
>
> there is no problem in modifying the langevin thermostat as long as you keep
> the ration between noise and friction. In Espresso the noise and friction
> prefactors are saved in langevin_pref1 and langevin_pref2.
> In theory you just have to change them.
>
> From the implementation side this is more complicated, because in Espresso
> there is only one global friction constant for all particles.
> So inside of the function friction_thermo_langevin, you need to know the
> number of neighbors, the position of the aggregate or whatever.
>
> Cheers,
>
> Christoph
>
> Olaf Lenz schrieb:
>>
>> -----BEGIN PGP SIGNED MESSAGE-----
>> Hash: SHA1
>>