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

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


From: Axel Kohlmeyer <akohlmey@...24...>
Date: Thu, 31 Aug 2017 13:29:12 -0400



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

​you keep posting stuff without really making an effort to resolve your problems and systematically debugging what you are doing.​ i don't have time to remind you of the same things over and over again. if the simulation is still not behaving correctly, then the part that is at fault is still in there.

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.

with this ​you are confirming​ yourself, that your force field parameters are not correct. 
 
The force field has been taken from literature (used in GCMC simulations)

​but obviously not applied correctly. LAMMPS is like any other MD code and follows the G​I-GO (= garbage in, garbage out) principle. 
 

thanks for bearing with me,

but my patience has reached is maximum extent now. ​this is my last response on this subject. i am tired of giving advice that is not followed (or understood?) properly.
 
axel.



 

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

------------------------------------------------------------------------------
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@...24...  http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.