# Thermostat in Dissipative Particle Dynamics

I'm doing a computational experiment about dynamics of a solvated line polymer, with dissipative particle dynamics(DPD). However when I try to measure kinetic temperature of the whole system(monomers + solvent particles) it is always overestimated and convergent to values about 1.5 ~ 2 times higher than that of DPD theory. For example, if I choose $\sigma=\sqrt{2\gamma kT},\,kT=1.0$ then kinetic temperature of the experiment is measured as 1.5 ~ 2.0.

As far as I know, thermostat in DPD is governed by dissipatve and random forces, regardless of the form of conservative force or other variables such as $l$ (standard length of LJ potential), $k$ (constant in harmonic potential). Tactics included in this experiment are as follows.

• All particles(4000 solvents, 20 monomers) are in a cubic box, with periodic boundary condition(PBC).
• LJ and harmonic potential, and DPD conservative/random/dissipative forces exist between monomers. They are calculated before PBC is exerted.
• Solvent-solvent, solvent-monomer interactions are calculated after PBC is exerted.
• $F_{Rx,\,ij}=\sigma\,(1-\frac{r_{ij}}{r_c})\xi_{ij}\,\frac{(x_i-x_j)}{\sqrt{(\Delta x)^2 + (\Delta y)^2 + (\Delta z)^2}}$ where $F_{Rx,\,ij}$ is random force via $x$ direction exerted on $i$-th particle by $j$-th particle and $\xi$ is random variable with mean = 0, std = 1.
• $F_{Dx,\,ij}=-\frac{(\sigma)^2}{2kT}\,(1-\frac{r_{ij}}{r_c})^2\,\frac{(x_i-x_j)(\Delta x\Delta v_x+\Delta y\Delta v_y+\Delta z\Delta v_z)}{{(\Delta x)^2 + (\Delta y)^2 + (\Delta z)^2}}$ where $F_{Dx,\,ij}$ is dissipative force via $x$ direction exerted on $i$-th particle by $j$-th particle$(\frac{(\sigma)^2}{2kT}=\gamma)$.
• Velocity-Verlet integral scheme with $dt=0.01$.
• I use linked-cell list method, modified little bit. Forces are calculated between particles in same cell and nearest cells.

Even if, since I'm undergraduate, I don't have reliable advisers to check my C-based code I'm wondering if there are some checklists when I found the kinetic temperature(so, velocities) of the system is overestimated, no matter how I choose the pair $(\sigma,kT)$ (and so $\gamma$).

• Which DPD formulation are you using (can you link to a paper)? The original DPD scheme is non-energy conserving and so it is isothermal only. You cannot get temperature information back out. There are extensions for energy conservation, DPD-E methods, but if you aren't using those extensions I doubt the energy will come out right. Commented Jan 20, 2016 at 11:52
• The other thing -- I don't know what the numbers are in your system, but your time step seems large. I usually only work in dimensional systems, so maybe in a non-dimensional formulation it isn't too big. But check out my question on SciComp and make sure you are within your stability limits. Commented Jan 20, 2016 at 11:54
• @tpg2114 I read link1, link2 and paragraph-28 of link3. Article in link3 asserts that we can get desired temperature via DPD thermostat but it doesn't seems to be DPD-E scheme you are referring. So I can see convergence of kinetic temperature but cannot foresee what is correct temperature of the system regardless of the value of initial $kT$ and $\sigma$.... Am I right? Commented Jan 20, 2016 at 12:08
• From the paper in the 3rd link, they say -- "The fluctuation-dissipation relationship described in equations (6) and (7) ensures that the temperature calculated from the kinetic energy of the particles will remain constant." Does your temperature remain constant or does it change? It could just be that you are initializing it with 2x the kinetic energy you think you have and that's why you're off by a factor of 2, but otherwise everything remains constant. If you run you calculation to compute temperature on your initial conditions, do you get what you expect? Commented Jan 20, 2016 at 13:31
• @tpg2114 Oh, since you've reminded me the importance of conditions that seem to be tiny, but in fact not at all, I have solved my thermostat problem. It comes from random forces, because I did not consider the symmetry of them between $i$ and $j$-th particles. Rather I set them totally independently. Thanks. Commented Jan 25, 2016 at 4:24

## 1 Answer

Not sure exactly how you are carrying out the velocity Verlet integration but as stated in the classic paper by Groot and Warren, "Care must be taken with the random force". The random force must be divided by the square root of the time step. So if you don't have that, then you are unlikely to get the correct behavior. (Paper is J. Chem. Pays 107, 4423 (1997)). Again, maybe you have taken account of this but I can't tell based on the expression above. Your time step is the typical one for DPD units but without knowing how stiff your springs are for the polymer, we can't determine if that is the problem.

Does this help?

• The thermostat problem is not from $\sqrt{dt}$ but the symmetry of $\xi_{ij}$ in random forces, as mentioned in my comment to @tpg2114, but with your advice, I could find the errors in random forces. Well, according to other related articles, my spring constants seem to be reasonable. Thanks. Commented Jan 25, 2016 at 4:27