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: Wed, 30 Aug 2017 13:38:11 -0400



On Wed, Aug 30, 2017 at 1:12 AM, surendra jain <jainsk.iitkgp@...24...> wrote:
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

​practically everything i commented on your previous setup still applies to this​ one. your force field parameters look highly bogus and the behavior of your system is consistent with having disabled LJ interactions between almost all atom types (and the ones that have feature a ridiculously small cutoff of 1 angstrom). this *cannot* be correct.
the behavior you see is the typical "coulomb collapse" that i would expect from such a setup.

so you need to repeat and review more thoroughly what i recommend previously.

axel.

 

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@...396...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.
>



--
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.