LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] trouble using rigid bodies in simulation (lost atoms)
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] trouble using rigid bodies in simulation (lost atoms)


From: Axel Kohlmeyer <akohlmey@...24...>
Date: Sat, 29 Jul 2017 00:55:03 -0400



On Fri, Jul 28, 2017 at 11:53 PM, surendra jain <jainsk.iitkgp@...24...> wrote:
Dear Axel,

  Many thanks for your response. I ran the simulation again at time
step of 0.05fs and still got the error :

ERROR: Lost atoms: original 7577 current 7574 (../thermo.cpp:427)

When I looked at the log file, I found that the temperature and
kinetic energy (of the rigid bodies) was fluctuating wildly. I am
using langevin
thermostat.

​please define "wildly". please also recall, that the magnitude of temperature fluctuations ​varies with the number of particles.
 

I have 2 options in my mind
(1) should I use pair_style soft in the beginning of the simulation. Is it
   possible to change the pair_style later doing the simulation.

​switching pair styles is possible during a simulation. whether this is useful in this case, i am unable to say.​
 
(2) Is it correct to change the force to a smaller value (say to a
maximum of 100.0) for small distances. Will it affect the dynamics in
the long run.

​this is a change in the force field and thus affect the validity of your simulation.
​please recall, that i mentioned that the symptoms you are​ seeing can also be cause by incorrect force field input.
since a see, that you have both, coulomb and lj interactions, your "force spikes" could for example be caused by incorrect parameters, where atoms get too close due to the coulomb forces becoming very strong at short distance due to errors in the lennard jones resulting in it not repulsive enough. if this hypothesis is true, through capping the force, you would achieve the exact opposite of what you intend.
 

Below is the input  script :

​this input is a mess!! it is full of commented out stuff and commands that are not needed to debug and commands that make no sense.
overall, it is not something that inspires confidence in that you are doing a good job here. if i was your adviser (which i am *not*), i would encourage you to start over and build your system in small increments and, e.g., figure out the correct settings for a stable and correct simulation reproducing results of a pure CO2 system first. in fact, i would recommend you find yourself a competent local(!) tutor, as none of us on this list has the time to debug your input deck for you or teach you in detail how to do that.

axel.
 

----------------------------------------------------------

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
#special_bonds   lj/coul 0.0 0.0 0.5

#variable d index res_1 res_2 res_3 res_4 res_5 res_6 res_7 res_8
res_9 res_10 res_11 res_12 res_13 res_14 res_15 res_16

#shell cd $d

read_data       lammps_data.dat

pair_coeff      * *  0.0      0.0
pair_coeff      1 8  0.0546   3.08
pair_coeff      2 8  0.00     1.4
pair_coeff      3 8  0.0744   3.275
pair_coeff      4 8  0.00     1.4
pair_coeff      5 8  0.106    2.88
pair_coeff      6 8  0.0955   2.9
pair_coeff      7 8  0.04     2.610
pair_coeff      1 9  0.0934   3.205
pair_coeff      2 9  0.00     1.525
pair_coeff      3 9  0.127    3.4
pair_coeff      4 9  0.00     1.525
pair_coeff      5 9  0.181    3.005
pair_coeff      6 9  0.163    3.025
pair_coeff      7 9  0.0685   2.735

pair_coeff      8 8  0.053654  2.8
pair_coeff      8 9  0.0917    2.925
pair_coeff      9 9  0.1569    3.05


#log            log.$d.data

# 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

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

#neigh_modify    delay 4 every 1
neighbor   2.0  bin
neigh_modify    delay 0 every 1 check yes
neigh_modify   exclude molecule clump
#thermo_modify  lost warn
#-----------------------------------------------------------
velocity        clump create 100.0 97287

fix             31 clump viscous  1.0
fix             11 matrix setforce 0.0 0.0 0.0
#fix             12 fluid nvt temp 300.0 300.0 0.01
#fix             12 all nvt temp 300.0 300.0 100
#fix             1 fluid nve/limit 0.1
compute          mytemp clump temp

timestep        0.050
#--------------------------------------------
group subgroup  id  6303  6304  6305 2082  2065
dump            191 subgroup  custom 1 f.xyz id mol type xu yu zu ix iy iz
dump            192 subgroup  custom 1 f1.xyz id mol type vx vy vz fx fy fz

compute 1 clump rigid/local 99 mol x y z vx vy vz fx fy fz
dump 111 clump local 1 tmp.dump index c_1[1] c_1[2] c_1[3] c_1[4]
c_1[5] c_1[6] c_1[7] c_1[8] c_1[9] c_1[10]

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

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

#fix             22 fluid nve
fix             2  clump langevin 300.0 300.0 100.0 498094

#dump 4a clump custom 10 dump.myforce.* id type x y vx fx fy fz
#thermo_style    custom step  c_mytemp pe ke
#thermo_modify  lost warn

thermo      10
run         1000

# do only NVE from now on
#unfix 2

reset_timestep  0

compute     2 fluid vacf
fix        14 fluid ave/time 1 1 1 c_2[4] file tmp.vacf
fix         5 fluid vector 1 c_2[4]
variable    diff equal dt*trap(f_5)
#thermo_style  custom step v_diff

#compute          CCRDF fluid rdf 200 1 1
compute          CCRDF all rdf 200  8 8
fix              3 all ave/time 100 1 1000 c_CCRDF[*] file cc.rdf mode vector
compute          CMRDF all rdf 200  1 8
fix             32 all ave/time 100 1 1000 c_CMRDF[*] file cm.rdf mode vector

compute         13 fluid msd com yes
variable        twopoint equal c_13[4]/4/(step*dt+1.0e-6)
fix             9 all vector 10 c_13[4]
variable        fitslope equal slope(f_9)/4/(10*dt)


#thermo_style    custom step temp pe ke etotal fmax c_13[4] v_twopoint
v_fitslope  v_diff c_mytemp
thermo_style    custom step temp pe ke etotal fmax c_13[4] c_mytemp
evdwl ecoul  epair ebond eangle elong etail

restart         100000 temp.restart

dump            108 all xyz 10 tmp.xyz
dump            109 fluid custom 10 tmp5.xyz id mol type xu yu zu ix iy iz
dump            110 fluid custom 10 tmp6.xyz id mol type x y z
dump            311 fluid custom 10 tmp7.xyz id mol type xu yu zu

timestep        0.1

#timestep        0.500

thermo          10
run             1000

------------------------------------------------------------------------------------------------

Best Regards,
Surendra



On 7/28/17, Axel Kohlmeyer <akohlmey@...24...> wrote:
> On Fri, Jul 28, 2017 at 8:02 AM, surendra jain <jainsk.iitkgp@...24...>
> wrote:
>
>> Dear Lammps users,
>>
>>   I am studying the diffusion of CO2 molecules in a carbon model. In
>> my simulation I fix the carbon atoms (setforce 0.0 ..) and make the
>> CO2 molecules as a rigid body molecule (having 3 atoms). However I get
>> the following error : "lost atoms:". After debugging I found that one
>> of the atoms of CO2 molecule is very close to a carbon atom (rsq = 0.5
>> Angsgtrom). This causes huge force and hence huge velocity on the
>> atom. I have the following question :
>>
>> (1) Does the rigid body moves as a whole  or the coordinates of the
>> atoms are updated individually.
>>
>
> ​this question is answered in the documentation ​
> http://lammps.sandia.gov/doc/fix_rigid.html
>
>
>> (2) If one atom faces huge force then the whole rigid body also
>> experiene the huge force.
>>
>
> ​see (1)
> ​
>
>> (3) Since the molecule is rigid, is it possible for one of the atoms
>> to be detached from the molecule.
>>
>
> ​no.​
>
>
>>
>> I ran a simulation of timestep 0.5fs and got the lost atoms error.
>> Should I use a smaller timestep. Can you please suggest me a suitable
>> timestep for rigid body simulations.
>>
>
> ​please make tests to find what works for your simulation conditions and
> force field. there is no "the one right value(tm)".
> please also note, that the situation you describe may as well be cause by
> unsuitable or incorrectly entered force field parameters.​
>
> ​axel.​
>
>
>
>>
>> Best Regards,
>> 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.