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

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


From: surendra jain <jainsk.iitkgp@...24...>
Date: Thu, 31 Aug 2017 22:00:47 +0530

Hi Axel,

  Many thanks for your kind reply. I am modeling the Nitrogen as rigid
molecule (it  has a massless particle). My LJ coefficients are :

 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
 pair_coeff      1 8  6.307E-02  3.335
 pair_coeff      2 8  0.00       1.655
 pair_coeff      3 8  8.59E-02   3.53
 pair_coeff      4 8  0.00       1.655
 pair_coeff      5 8  0.12255    3.135
 pair_coeff      6 8  0.11029    3.155
 pair_coeff      7 8  4.629E-02  2.865

 pair_coeff      8 8   7.152E-02  3.31

atom type 1 to 7 are fixed with setforce.
 atom type (8 and 9) belongs to Nitrogen. atom type 9 is a massless
particle, but it has
  a charge.

The system is unstable when I run the simulation (KE, Fmax and
temperature increases).
  But if I turn off the coulombic interaction, the system is stable.
The force field has been taken from literature (used in GCMC simulations)

thanks for bearing with me,

Surendra


On 8/30/17, Axel Kohlmeyer <akohlmey@...24...> wrote:
> 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@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.
>> >
>>
>
>
>
> --
> 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.
>