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: surendra jain <jainsk.iitkgp@...24...>
Date: Wed, 30 Aug 2017 10:42:16 +0530

Hi Axel,

  Many thanks for your kind reply. I ran a Nitrogen only diffusion
simulation in a carbon model. The carbon matrix atoms are fixed. The
nitrogen molecules are rigid (as it has a massless particle). I still
get error. I am attaching the input script and output

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      lj/cut/coul/long 10.0 12.0
bond_style      harmonic
angle_style     harmonic
pair_modify     mix arithmetic

read_data       lammps_data.dat

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

pair_coeff      8 8   7.152E-02  3.31  1.0

# Make the group of HRMC matrix (carbon + hydrogen + oxygen) as one group
# these group of atoms are fixed and do not move during MD simulaitons

group   matrix type 1 2 3 4 5 6 7
group   fluid type 8 9

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

neighbor   2.0  bin
neigh_modify    delay 0 every 1 check yes

fix              31 clump langevin 0.0 0.0 100.0 498094
fix             11 matrix setforce 0.0 0.0 0.0

compute          mytemp clump temp

timestep        0.050

thermo_style    custom  step c_mytemp pe ke fmax evdwl ecoul  epair
ebond eangle elong etail
thermo          100
run             1000
unfix 31

timestep        0.5
fix             2  clump langevin 100.0 100.0 100.0 498094


thermo      100
run         1000
----------------------------------

The output is :

Step c_mytemp PotEng KinEng Fmax E_vdwl E_coul E_pair E_bond E_angle
E_long E_tail
       0    97.025098   -106313.98    44.056839    516.92036
-285.34782   -87624.447   -106313.98 1.7179598e-06 2.8477085e-05
-18404.183            0
     100    89.482797   -106314.59    40.632056    516.57792
-285.49175   -87624.901   -106314.59 1.7179598e-06 2.8477075e-05
-18404.195            0
     200    85.722973   -106316.85    38.924807    516.13706
-286.42579   -87626.218   -106316.85 1.7179598e-06 2.8477083e-05
-18404.203            0
     300    85.178394   -106320.43    38.677527    515.59328
-288.04518   -87628.175   -106320.43 1.7179598e-06 2.8477076e-05
-18404.207            0
     400    87.130394   -106325.11    39.563884    514.95716
-290.19447   -87630.712   -106325.11 1.7179598e-06 2.8477072e-05
-18404.208            0
     500    90.861841   -106330.95    41.258248    514.24475
-292.67597    -87634.07   -106330.95 1.7179598e-06 2.8477074e-05
-18404.206            0
     600    95.744801   -106337.31    43.475486    513.46945
-295.29649    -87637.81   -106337.31 1.7179598e-06 2.8477079e-05
-18404.201            0
     700    101.31344   -106344.07    46.004073    513.74009
-297.89599   -87641.981   -106344.07 1.7179598e-06 2.8477079e-05
-18404.194            0
     800     107.3097   -106351.79     48.72684    513.76597
-300.31339   -87647.292   -106351.79 1.7179598e-06 2.8477072e-05
-18404.185            0
     900    113.70898   -106359.51    51.632603    513.34306
-302.46022   -87652.879   -106359.51 1.7179598e-06 2.8477074e-05
-18404.176            0
    1000    120.75433   -106367.99    54.831729    512.49498
-304.2267   -87659.596   -106367.99 1.7179598e-06 2.8477083e-05
-18404.166            0
    1100    129.04406    -106377.3    58.595903    511.26018
-305.53933   -87667.606    -106377.3 1.7179598e-06 2.8477077e-05
-18404.159            0
    1200    139.62714   -106387.93    63.401435      509.687
-306.27819   -87677.501   -106387.93 1.7179598e-06 2.8477074e-05
-18404.153            0
    1300    154.30551   -106401.29    70.066543    509.16382
-306.20051   -87690.934   -106401.29 1.7179598e-06 2.8477073e-05
-18404.151            0
    1400    176.39472    -106418.8     80.09674    509.08376
-304.76943    -87709.88    -106418.8 1.7179598e-06 2.8477079e-05
-18404.151            0
    1500    213.38352   -106444.42    96.892493    508.90096
-300.43412   -87739.836   -106444.42 1.7179598e-06 2.8477073e-05
-18404.155            0
    1600      291.777    -106491.2    132.48915    508.63867
-287.71178   -87799.326    -106491.2 1.7179598e-06 2.8477079e-05
-18404.159            0
    1700    799.98828   -106741.41    363.25606    1891.8344
-226.18088    -88111.07   -106741.41 1.7179598e-06 2.8477082e-05
-18404.158            0
    1800    6966443.5   -106258.06    3163299.9    3454.9425
-129.72925   -87724.374   -106258.06 1.7179598e-06 2.8477076e-05
-18403.958            0
ERROR: Lost atoms: original 11052 current 11049 (../thermo.cpp:433)
Last command: run             10000
----------------------------------------------------------

Can you please tell me what I am missing here.

Many thanks,
Surendra



On 8/28/17, Axel Kohlmeyer <akohlmey@...24...> wrote:
> 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@lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/lammps-users
>>
>
>
>
> --
> Dr. Axel Kohlmeyer  akohlmey@...24...  http://goo.gl/1wk0
> College of Science & Technology, Temple University, Philadelphia PA, USA
> International Centre for Theoretical Physics, Trieste. Italy.
>