Re: [ESPResSo-users] Sistematically decreasing of Temperature with DPD thermostat

[Stefan Kesselheim]

Dear Salvador,

I'd not be sure Georgs or Ulfs points are valid, but let me give it a try.
a) It could be a bug. The line should read 
thermostat off 
and not thermostat set off
b) Thermalization of rotational degrees of freedom is entirely different in 
Langevin than what happens when combining rigid bodies with DPD (this so far 
probably nobody has done). In LD there is a friction proportional to the 
angular velocity and the the COM motion. When using DPD the "usual" DPD stuff 
happens to the particles, that is every particle experiences a friction with 
every other particle. Per se this includes also the other particles in the same 
dumbbell and this might eat up energy (In the non-transversal DPD this should 
however be exactly zero. This you probably have switched off).
When other people used rigid bodies with LB I was expecting trouble, because 
then friction should be correlated among all the particles (they are connected 
rigidly) and the noise however is uncorrelated. This however obviously worked 
correctly. For dimers (or trimers) this should however not reduce the 
temperature by almost a factor 100.

Before going into details with that I'd check the following things:
First: is possible, that you gamma is 0? Try using 1, please.

Second: Can you cut off the WF/TWF tail of your command? It could be a bug in 
some parser.

Third: I'd try is replacing the rigid bonds by simple harmonic bonds between 
all three particles of the dumbbell. This is something that should definitely 
yield the right temperature. If this still produces the wrong temperature there 
is another bug in you script. If this produces the correct temperature then we 
should think again.

Four: I'm not sure what the exclusions do to DPD. You should switch them off if 
possible.

Five: You could switch off all "conservative", i.e. LJ and bonded interactions. 
If this does not produce the right temperature there is more bugs.

I hope one of these things works.
Cheers
Stefan


On Jan 14, 2014, at 9:34 PM, Salvador H-V <address@hidden> wrote:

> Dear All,
> 
> I am performing simulations of a two dimensional system of hard dummbbells.
> 
> I am using a generalized LJ potential to simulate the hard interaction and 
> rigid bonds to constraint the two center of the spheres that composed one 
> dimer.
> 
>  I was using LD thermostat but as I am interested in dynamical properties I 
> decided to switch to DPD thermostat. Here it is my [code_info] compilation:
> 
>  Compilation status { FFTW } { CONSTRAINTS } { BOND_ANGLE_HARMONIC } { 
> GAUSSIAN } { HERTZIAN } { PARTIAL_PERIODIC } { LENNARD_JONES } { MODES } { 
> EXTERNAL_FORCES } { TABULATED } { SOFT_SPHERE } { LENNARD_JONES_GENERIC } { 
> MASS } { ELECTROSTATICS } { BOND_CONSTRAINT } { DPD } { TRANS_DPD } { 
> EXCLUSIONS } }"
> 
> 
> I am running the same program that I previously used with LD only with the 
> following modifications:
> 
> thermostat set off
> galilei_transform
> set dpd_cut [ expr { ( pow ( 2.0 , 1.0/6.0 ) )*$diam } ]
> thermostat dpd $temp_ave 0.0 $dpd_cut WF 0 1.0 $dpd_cut TWF 0
> (dpd 1.5 0.0 1.122462048309373 WF 1 1.0 1.122462048309373 TWF 0 )
> 
> When I switch form LD to DPD (standard or transversal) thermostat, I get a 
> sistematically decrease of the temperature: set temp [expr [analyze energy 
> kinetic]/($np)]
> un 100 at time=49.99999999997417 E = 1.464997510334239 & Temp = 
> 1.478570945756926^M
> run 200 at time=59.99999999995087 E = 1.445825470054056 & Temp = 
> 1.463115509638189^M
> .
> .
> .
> run 124200 at time=12460.00000207969 E = 0.01086010276683752 & Temp = 
> 0.02004077077749828^M
> run 124300 at time=12470.00000208172 E = 0.01324369725959244 & Temp = 
> 0.02005819667634172
> 
> I don't understand what is going on? Any idea or suggestion would be greatly 
> appreciated.
> 
> Thanks a  lot in advance,
> 
> Salvador
> 
> 
> 
> 
> 
>