Re: [ESPResSo] Position-dependent Langevoin thermostat

[Lorenzo Isella]

2008/11/10 Axel Arnold <address@hidden>:
> Hi Lorenzo!
>
>> diff  /media/disk/espresso_modified/espresso-desperation/thermostat.c
>> /media/disk/espresso_modified/espresso-2.1.0/thermostat.c
>>
>> < /* Add the eta array */
>> <
>> < double eta[2] ;
>> <
>> <
>> 270,271c265,266
>> <   langevin_pref1 = -eta[0]*langevin_gamma/time_step;
>> <   langevin_pref2 =
>> sqrt(eta[1]*24.0*temperature*langevin_gamma/time_step);
>> ---
>>>
>>>  langevin_pref1 = -langevin_gamma/time_step;
>>>  langevin_pref2 = sqrt(24.0*temperature*langevin_gamma/time_step);
>>
>> I did not touch thermostat.h, while I did some work in particle data.c
>> (one routine to print and one parse eta[]  and one to set it without
>> resorting to  MPI, as you suggested)
>>
>> diff  /media/disk/espresso_modified/espresso-desperation/particle_data.c
>> /media/disk/espresso_modified/espresso-2.1.0/particle_data.c
>> 92,97d91
>> <   part->p.eta[0]   = 1.0;
>> <   part->p.eta[1]   = 1.0;
>
> That cannot work - you put the values of eta into the particle structure,
> but calculate the langevin prefactors from a local variable eta in
> thermostat.c. This variable is uninitialized and therefore 0, which explains
> why you don't see movement - the thermostat is off. What you need to do is
> in fact to not touch the langevin prefactors as they are global for all
> particles, but rather
> friction_thermo_langevin(Particle *p) in thermostat.h. Here, you introduce
> your new prefactors directly, as
>        p->f.f[j] = p->eta1*langevin_pref1*p->m.v[j]*PMASS(*p) +
> sqrt(p->eta2)*langevin_pref2*(d_random()-0.5)*massf;
>
> If you want to save some computation time, just save sqrt(eta2) with the
> particle, by changing the set/get routines in particle_data.c accordingly.
>
> Cheers,
> Axel

Hello Axel,
Yes, you are right and it all makes sense. In fact now I see that the
thermostat is active.
I hope I will not have to take it back, but it looks like I am now
done with the Espresso scripting for a while.
Many thanks

Lorenzo