dpd/tstat behavior

I am using the dpd/tstat pair style with nve as a thermostat for an explicitly solvated system. I have two questions/issues that I would appreciate input on.

1. Using a thermostat gamma parameter on the order of ~1 kcal fsec /mol �^2 or less I have no problem running the initial simulation at the desired temperature. However, if I continue the simulation from a restart file the temperature drops rapidly close to zero. I don't understand why this is happening only after a restart. Perhaps someone can educate me and give me ideas on how to avoid this behavior?

2. I can successfully restart a run if I use a very small value for gamma (i.e. ~10^-8 kcal fsec /mol �^2). However in this case duplicate simulations from the same restart file yield exactly the same trajectories. I expected since the specific allocation of the random force is different after a restart that duplicate simulations should produce different trajectories? Why don't I get different trajectories? Is gamma too small in this case?

I am using lammps-30Aug12

Thank you,
Adam

I am using the dpd/tstat pair style with nve as a thermostat for an
explicitly solvated system. I have two questions/issues that I would
appreciate input on.

1. Using a thermostat gamma parameter on the order of ~1 kcal fsec /mol
Å^2 or less I have no problem running the initial simulation at the
desired temperature. However, if I continue the simulation from a
restart file the temperature drops rapidly close to zero. I don't
understand why this is happening only after a restart. Perhaps someone
can educate me and give me ideas on how to avoid this behavior?

hard to say. there is probably something wrong in the way how you do
the restarting. in these cases, it is usually best, if you could
provide a test system input that reproduces the behavior (the simpler
and the smaller, the better) so that people can have a look, try to
reproduce it and either point out your mistake or try to track down
the bug.

2. I can successfully restart a run if I use a very small value for
gamma (i.e. ~10^-8 kcal fsec /mol Å^2). However in this case duplicate
simulations from the same restart file yield exactly the same
trajectories. I expected since the specific allocation of the random
force is different after a restart that duplicate simulations should
produce different trajectories? Why don't I get different trajectories?
Is gamma too small in this case?

again, hard to say without a suitable test system to reproduce it.

axel.

Axel,
I am attaching a test system. It is just 256 flexible SPC water molecules.

SPCFw.in will run with the dpd/tstat gamma parameter set to either 1e-3 or 1e-8 (I'm using metal units). Although with the smaller value the system approaches the target temperature very slowly. (see log.1e-3 and log.1e-8)

However, when restarting (see SPCFw_restart.in) the temperature will drop if gamma=1e-3, but will stay about the same if gamma=1e-8 (see log.1e-3.restart and log.1e-8.restart). Also, as I mentioned in the last message, if I restart the job with gamma=1e-8 multiple times I get the same trajectory each time.

I'd appreciate any thoughts on how to correct this behavior, or how to correct mine as the case may be.
Thank you,
Adam

SPCFw_restart.in (477 Bytes)

SPCFw.in (487 Bytes)

log.1e-3 (40.3 KB)

log.1e-3.restart (40.2 KB)

log.1e-8 (40.3 KB)

log.1e-8.restart (40.2 KB)

SPCFw.data (104 KB)

hi adam,

thanks for the perfect test example. your input is correct, you have
run into a bug in dpd/tstat where the target temperature is not
properly initialized when reading from a restart. you can work around
it by moving the pair_style line *after* the read_restart command.

the real fix is this one line file to the source code:

diff --git a/src/pair_dpd_tstat.cpp b/src/pair_dpd_tstat.cpp
index f3f7ecc..9e9e3ee 100644
--- a/src/pair_dpd_tstat.cpp
+++ b/src/pair_dpd_tstat.cpp
@@ -229,6 +229,8 @@ void PairDPDTstat::read_restart_settings(FILE *fp)
   MPI_Bcast(&seed,1,MPI_INT,0,world);
   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);

+ temperature = t_start;

thanks for finding this - will be in the next patch …

Steve