Probably a bug in "fix langevin"

Hi, Lammps users,

I implement the langevin method to control the temp of water flow. My system is something like water is confined by two plates, and the top plate can move up and down. So the volume is not fixed.

The problem is: I tested two versions of Lammps , one is 22Feb, the other one is 7June. Same .in file, same workstation, but the temp is different. See below:

This is the output from 7June:

Step Temperat PotEng
1002 282.25884 78034743
2002 588.39786 78042837
ERROR: Lost atoms: original 10400 current 9600 (…/thermo.cpp:389)

This is the output from 22Feb:

Step Temperat PotEng
1002 282.25884 78034743
2002 296.61572 78034247
Loop time of 17.4039 on 1 procs for 1000 steps with 10400 atoms.

Then this is the output from 7June, but by Nose-Hoover thermostat:

Step Temperat PotEng
1002 282.19582 78034743
2002 292.93137 78034475
Loop time of 17.1521 on 1 procs for 1000 steps with 10400 atoms.

Best wishes,

Wei

Hi,

I think it would be much better to provide small input file so that one could reproduce observed behaviour. There may be many reasons for this.

Also, take a look at the changes in fix langevin made since 22 Feb:
http://lammps.sandia.gov/bug.html
For example, there is a backward compatibility note (24 Apr):
BACKWARD COMPATIBILITY: The previous "angmom yes" syntax of the fix langevin command is no longer allowed. You must choose a numeric scale factor.

Regards,
Oleg.

20.06.2013, 12:41, "Chen Wei" <[email protected]>:

Hi,

Thanks for your reply.

These are the .in and .data file.

The system is SPC/Fw water confined by graphene sheets.

Best wishes,

Wei

spce_graphene.in (1.49 KB)

system.data (556 KB)

The problem is in fact reproduced with current version of lammps while does not take place with 14 May version (though I did not observe lost atoms)... Maybe it is related to some recent changes in fix langevin.

Anyway, I don't think your system is defined correctly. During first run with NVT there are pressures of 10^8 bar:
Step Temp E_pair E_mol TotEng Press
       2 206.29887 78096213 2882.9089 78105490 1.0320949e+08
    1002 197.3678 78031926 2796.4166 78040840 1.0318184e+08

You may have bad dynamics and lost atoms because of too dense system. Try to visualize your cell and see what atoms were lost and how.

Oleg.

20.06.2013, 13:18, "Chen Wei" <[email protected]>:

Maybe the problem is here:

http://lammps.sandia.gov/bug.html

7 Jun 2013

Aidan Thompson (Sandia) has added the new Gronbech-Jensen/Farago formulation of the Langevin thermostat to LAMMPS as a gjf option on the fix langevin command. It allows the use of much larger timesteps (e.g. 3 fs) for systems containing bonds with light atoms. See the fix langevin doc page for more info.

And thanks for your suggestion, I will check that.

Best wishes,

Wei

It is not used by default, but of course Dr. Thompson knows better.

Regards,
Oleg.

20.06.2013, 16:24, "Chen Wei" <[email protected]>:

hmmm... there are several problems with this input deck. it doesn't
quite act like you describe it and its neighborlist settings result in
*lots* of dangerous builds. that is not good. also you should exclude
the C-C interactions, they just make the output confusing and do a few
more tweaks.

diff -u spce_graphene.in spce_graphene_mod.in
--- spce_graphene.in 2013-06-20 13:18:30.000000000 +0200
+++ spce_graphene_mod.in 2013-06-20 15:25:09.050134908 +0200
@@ -3,7 +3,8 @@
# Created on Mar 11, 2013

#Basic parameters
-newton off
+newton on
+processors * * 1

variable T equal 298.0

@@ -52,13 +53,16 @@

group bot molecule <> 3201 4000

-neighbor 1.0 bin
+neighbor 2.0 bin
+neigh_modify every 2 delay 2 check yes exclude type 3 3

-neigh_modify every 5 delay 0 check yes

velocity spce create $T 4928459 mom yes rot yes dist gaussian

velocity top set 0 0 0 units box
+velocity bot set 0 0 0 units box
+compute wat spce temp
+thermo_modify temp wat

fix 1 bot setforce 0 0 0

ok. i think i found the bug. try this change to src/fix_langevin.cpp:

[[email protected] src]$ git diff fix_langevin.cpp
diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp
index d66efb9..6f11342 100644
--- a/src/fix_langevin.cpp
+++ b/src/fix_langevin.cpp
@@ -513,6 +513,8 @@ void FixLangevin::post_force_untemplated
     }
   }

+ if (Tp_BIAS) temperature->compute_scalar();

Thanks Axel. I just checked it in.

Hi, Axel,

Thank you very much for your reply.

Actually, I was running Lammps with "-suffix gpu", neither "newton on" nor "neigh_modify exclude" is allowed then.
Do you have a suggestion?

Thanks again,

Best wishes,

Wei

Hi, Axel,

Thank you very much for your reply.

Actually, I was running Lammps with "-suffix gpu", neither "newton on" nor "neigh_modify exclude" is allowed then.
Do you have a suggestion?

no.