LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Temperature and Energy increases with langevin thermostat
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Temperature and Energy increases with langevin thermostat


From: Axel Kohlmeyer <akohlmey@...24...>
Date: Sun, 27 Aug 2017 16:00:59 -0400



On Sun, Aug 27, 2017 at 6:36 AM, surendra jain <jainsk.iitkgp@...24...> wrote:
Dear Lammps users,

  I am running MD simulations of mixture (N2 + CO2) in a confined
system. I encounter some problems. I have some questions :
(1) I used  "pair_style   lj/cut/coul/long/soft"  for CO2. and I have
N2 as rigid molecules (as it has a massless particle). I couldn't use
"pair_style   lj/cut/coul/long/soft" for N2.
(2) I used NVE/limit and langevin 0.0 0.0 100.0 to equilbrate the
system initially. I use different fix for rigid and non rigid
molecules. However, the temperature increases even if I use langevin
thermostat 0.0 0.0 100.0
(3) I have also found that ECOUL is +ve for using pair_style soft. Can
someone tell me why
 this is happening.


​most likely, your system is misbehaving because of a bad set of potential parameters.
many pairs of atom types have no LJ interaction, but only coulomb. so if you have atoms with opposite charge, they can get extremely close leading to very high negative energies and then divergence/crash or atoms being accelerated massively.

other issues:
- there doesn't seem to be a good reason to use a soft core potential here, at least not during equilibration, so using a hybrid potential makes no sense, and using hybrid/overlay even less
- it is a very bad idea to bypass errors over missing pair coefficients by initializing all lennard jones coefficients to zero. you don't turn off the coulomb part. if you explicitly want no interactions between certain atom types you should use pair style hybrid with either zero or none as a substyle for those pairs.
- to use a thermostat with a 0K target is not very meaningful. if you want to add friction use fix viscous for simulated annealing, otherwise use the minimize command.

that you complain about energy rising makes no sense, either, it is just your system ​relaxing into a lower energy state. if you do that during an MD run, potential energy will be converted to kinetic energy. that is how MD works.

​overall, you should probably first check your model and figure out how to correctly reproduce CO2 and N2 separately, and only then worry about how to do a combined simulation.

axel.​


 
I am attaching the input script :

dimension       3
units           real
boundary        p p p

kspace_style    pppm 1.0e-4
kspace_modify   order 3
atom_style      full
pair_style       hybrid/overlay lj/cut/coul/long 10.0 12.0
lj/cut/coul/long/soft 2.0 0.5 10.0 10.0 12.0
bond_style      harmonic
angle_style     harmonic
pair_modify     mix arithmetic
#special_bonds   lj/coul 0.0 0.0 0.5


read_data       lammps_data.dat

pair_coeff      * *   lj/cut/coul/long 0.0      0.0
pair_coeff      1 8   lj/cut/coul/long 6.307E-02  3.335
pair_coeff      2 8   lj/cut/coul/long 0.00       1.655
pair_coeff      3 8   lj/cut/coul/long 8.59E-02   3.53
pair_coeff      4 8   lj/cut/coul/long 0.00       1.655
pair_coeff      5 8   lj/cut/coul/long 0.12255    3.135
pair_coeff      6 8   lj/cut/coul/long 0.11029    3.155
pair_coeff      7 8   lj/cut/coul/long 4.629E-02  2.865

pair_coeff      1 10   lj/cut/coul/long/soft 0.0546   3.08    0.95
pair_coeff      2 10   lj/cut/coul/long/soft 0.00     1.4     0.95
pair_coeff      3 10   lj/cut/coul/long/soft 0.0744   3.275   0.95
pair_coeff      4 10   lj/cut/coul/long/soft 0.00     1.4     0.95
pair_coeff      5 10   lj/cut/coul/long/soft 0.106    2.88    0.95
pair_coeff      6 10   lj/cut/coul/long/soft 0.0955   2.9     0.95
pair_coeff      7 10   lj/cut/coul/long/soft 0.04     2.610   0.95

pair_coeff      1 11   lj/cut/coul/long/soft 0.0934   3.205   0.95
pair_coeff      2 11   lj/cut/coul/long/soft 0.00     1.525   0.95
pair_coeff      3 11   lj/cut/coul/long/soft 0.127    3.4     0.95
pair_coeff      4 11   lj/cut/coul/long/soft 0.00     1.525   0.95
pair_coeff      5 11   lj/cut/coul/long/soft 0.181    3.005   0.95
pair_coeff      6 11   lj/cut/coul/long/soft 0.163    3.025   0.95
pair_coeff      7 11   lj/cut/coul/long/soft 0.0685   2.735   0.95

pair_coeff      8 8    lj/cut/coul/long 7.152E-02  3.31
pair_coeff      10 10  lj/cut/coul/long/soft 0.053654  2.8    0.95
pair_coeff      10 11  lj/cut/coul/long/soft 0.0917    2.925  0.95
pair_coeff      11 11  lj/cut/coul/long/soft  0.1569    3.05  0.95
pair_coeff      8  10  lj/cut/coul/long 6.19e-02 3.055
pair_coeff      8  11   lj/cut/coul/long 0.105   3.18


group   matrix type 1 2 3 4 5 6 7
group   fluid type 8 9 10 11
group fluid1 type 8 9
group fluid2 type 10 11

#-----------------------------------------------------

 group clump type 8 9
 fix   99 clump rigid/small molecule

#----------------------------------------------------

neigh_modify    delay 0 every 1 check yes

fix             11 matrix setforce 0.0 0.0 0.0
fix             1 fluid2 nve/limit 0.1

fix              3 fluid2 langevin 0.0 0.0 10.0 498094
fix              31 clump langevin 0.0 0.0 10.0 498094
compute          mytemp1 fluid1 temp
compute          mytemp2 fluid2 temp

timestep        0.050

thermo_style    custom  step c_mytemp1 c_mytemp2 pe ke fmax evdwl
ecoul  epair ebond eangle elong etail
thermo          100
run             10000
unfix 3
unfix 1
unfix 31

fix             22 fluid2 nve

fix             24  fluid1 langevin 100.0 100.0 10.0 498094
fix             25  fluid2 langevin 100.0 100.0 10.0 498094

thermo      100
run         10000

unfix 24
unfix 25

fix             28  fluid1 langevin 300.0 300.0 10.0 498094
fix             29  fluid2 langevin 300.0 300.0 10.0 498094

thermo      100
run         10000

------------------------------
My output after some steps :

ERROR on proc 0: Bond atoms 7086 7088 missing on proc 0 at step 1146
(../ntopo_bond_all.cpp:63)  for N2andCO2 system

Step c_mytemp PotEng KinEng Fmax E_vdwl E_coul E_pair E_bond E_angle
E_long E_tail
       0            0   -41741.562            0          nan
-2208.7061   -12972.774   -41741.562 2.3415094e-05 6.0339129e-05
-26560.082            0
ERROR on proc 0: Out of range atoms - cannot compute PPPM (../pppm.cpp:1927)


Many thanks,
Surendra

------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...6297....sourceforge.net
https://lists.sourceforge.net/lists/listinfo/lammps-users



--
Dr. Axel Kohlmeyer  akohlmey@...92......  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.