[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [ESPResSo] Position-dependent Langevoin thermostat
From: |
Lorenzo Isella |
Subject: |
Re: [ESPResSo] Position-dependent Langevoin thermostat |
Date: |
Mon, 10 Nov 2008 18:03:30 +0100 |
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