Using dpd/tstat changes E_coul energy

Dear all,

I am using Lammps-1Feb2014 and I am having a problem with dpd/tstat: Including dpd/tstat with hybrid/overlay option changes the E_coul energy (typically to very large positive number for condensed phases). However all the other energy terms remain the same. Here I copied my input files for just 3 water molecules with and without the thermostat.

With dpd/tstat

----------------- Init Section -----------------

units real

atom_style full

pair_style hybrid/overlay lj/charmm/coul/long 9.0 10.0 10.0 dpd/tstat 300 300 9.0 34387

bond_style hybrid harmonic

angle_style hybrid harmonic

kspace_style pppm 0.0001

----------------- Atom Definition Section -----------------

LAMMPS Description

9 atoms

6 bonds

3 angles

0 dihedrals

0 impropers

2 atom types

1 bond types

1 angle types

0.0 41.000 xlo xhi

0.0 41.000 ylo yhi

0.0 101.000 zlo zhi

Masses

1 15.9994

2 1.008

Atoms

1 1 1 -0.8476 5.690 12.751 11.651

2 1 2 0.4238 4.760 12.681 11.281

3 1 2 0.4238 5.800 13.641 12.091

4 2 1 -0.8476 15.551 15.111 7.030

5 2 2 0.4238 14.981 14.951 7.840

6 2 2 0.4238 14.961 15.211 6.230

7 3 1 -0.8476 17.431 6.180 8.560

8 3 2 0.4238 17.761 7.120 8.560

9 3 2 0.4238 17.941 5.640 9.220

Bonds

1 1 1 2

2 1 1 3

3 1 4 5

4 1 4 6

5 1 7 8

6 1 7 9

Angles

1 1 2 1 3

2 1 5 4 6

3 1 8 7 9

----------------- Settings Section -----------------

bond_coeff 1 harmonic 1000.0 1.0

angle_coeff 1 harmonic 1000.0 109.47

pair_coeff * * dpd/tstat 2.4

pair_coeff 1 1 lj/charmm/coul/long 0.1553 3.166

pair_coeff 2 2 lj/charmm/coul/long 0.0 2.058

group spce type 1 2

fix fShakeSPCE spce shake 0.0001 10 100 b 1 a 1

----------------- Run Section -----------------

communicate single vel yes

thermo_style multi

thermo 15

unfix fShakeSPCE

minimize 1.0e-5 1.0e-7 1000 10000

The output for step 0 is like this:

---------------- Step 0 ----- CPU = 0.0000 (sec) ----------------

TotEng = -140.5061 KinEng = 0.0000 Temp = 0.0000

PotEng = -140.5061 E_bond = 0.0789 E_angle = 0.1960

E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = -0.0008

E_coul = -33.7868 E_long = -106.9934 Press = 2.3176

For the case without dpd/tstat, I only change the init and setting sections:

----------------- Init Section -----------------

units real

atom_style full

pair_style hybrid lj/charmm/coul/long 9.0 10.0 10.0

bond_style hybrid harmonic

angle_style hybrid harmonic

kspace_style pppm 0.0001

----------------- Settings Section -----------------

bond_coeff 1 harmonic 1000.0 1.0

angle_coeff 1 harmonic 1000.0 109.47

pair_coeff 1 1 lj/charmm/coul/long 0.1553 3.166

pair_coeff 2 2 lj/charmm/coul/long 0.0 2.058

group spce type 1 2

fix fShakeSPCE spce shake 0.0001 10 100 b 1 a 1

And I get these energy terms:

---------------- Step 0 ----- CPU = 0.0000 (sec) ----------------

TotEng = 0.2434 KinEng = 0.0000 Temp = 0.0000

PotEng = 0.2434 E_bond = 0.0789 E_angle = 0.1960

E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = -0.0008

E_coul = 106.9627 E_long = -106.9934 Press = 1.5867

I don’t know what I am doing wrong here. I appreciate any help on this problem.

Best regards,

Mohammad

Postdoctoral fellow

Department of Chemical and Petroleum Engineering

University of Wyoming

Phone: 307-761-0316

Dear all,

I am using Lammps-1Feb2014 and I am having a problem with dpd/tstat:

please first try with the very latest version of LAMMPS (9 October
2014). nobody has a desire to look into a possible bug, when it may
have already been fixed.

axel.

Dear Axel,

As you suggested I tried the very latest version (9 Oct 2014) but I am still experiencing the same problem.

In the newest version I had to use the command "comm_modify vel yes" instead of "communicate single vel yes", but everything else is the same.

Thanks,
Mohammad

Dear Axel,

As you suggested I tried the very latest version (9 Oct 2014) but I am still experiencing the same problem.

In the newest version I had to use the command "comm_modify vel yes" instead of "communicate single vel yes", but everything else is the same.

ok. upon closer inspection of the input, it turns out that the problem
is not a LAMMPS problem, but a problem of your input.
your specification of the pair coefficients depends on LAMMPS
performing mixing of the parameters and assigning them to the missing
pairs for you.
this works for the simple case of not specifying dpd/tstat, however,
if you *do* specify it, LAMMPS has no way of knowing whether you want
it to generated the parameters between atom types 1 and 2 or to just
keep it as it is. for hybrid overlay, this is quite normal. thus
LAMMPS only does mixing if the choice is obvious. for details, have
another look at the documentation of pair_style hybrid and
hybrid/overlay.

pair_coeff 1 2 lj/charmm/coul/long 0.0 2.553

axel

Thanks Axel, the problem is fixed.

Regards,
Mohammad