FW: a question in fix rigid

LAMMPS LAMMPS Mailing List Mirror
1

The cnt.dat file has been attached (a part of the studied structure to reduce the size of the dat file). The lammps dat file has been created through the writelammpsdata in VMD .

dimension 3
boundary p p f
units real
atom_style atomic
neighbor 0.4 bin
neigh_modify delay 5
lattice fcc 4.0

######Data file######################################
read_data cnt.dat

######Regions########################################
region upper_container block 3.477263 70.522736 5.777263 77.922737 80 150 side in units box

#####Define Gases####################################
create_atoms 2 random 2000 38790 upper_container #
create_atoms 3 random 2000 38792 upper_container #
group tube type 1
group CO2 type 2
group CO type 3
group gas union CO2 CO

############Potentials################################
pair_style lj/cut 20.0
pair_coeff * * 0.0 3.0 0.0
pair_coeff 1 1 0.086 3.4 0.0 ###if I change 0.0 to somehting else it gives Nan and Inf.
pair_coeff 2 2 0.468781 3.454 9.0
pair_coeff 3 3 0.065 3.920 9.0

pair_coeff 1 2 0.2007863 3.427 9.0
pair_coeff 1 3 0.0747663025 3.66 9.0
pair_coeff 2 3 0.174558772 3.687 9.0
###############INITIALIZATION########################
velocity CO2 create 300.0 4928459 dist gaussian
velocity CO create 300.0 4928459 dist gaussian

fix zUpperWall all wall/reflect zhi EDGE units box
fix zLowWall all wall/reflect zlo EDGE units box

cnt.dat (8.64 KB)

2

The cnt.dat file has been attached (a part of the studied structure to
reduce the size of the dat file). The lammps dat file has been created
through the writelammpsdata in VMD .

if i run this system it produces a pressure of
2.135556e+16 atm in step zero. don't you
think that this is more than a little bit suspicious?

axel.

3

Yes there are a lot of molecules after some timesteps it is reduced to a normal pressure. Obviously the two randoms have made some overlaps I guess. If I reduce the number of the molecules it doesn’t give such pressures.

4

Yes there are a lot of molecules after some timesteps it is reduced to a
normal pressure. Obviously the two randoms have made some overlaps I guess.

yes, and that needs to be dealt with. you have atoms "shooting"
through the wall at the beginning and lots of "dangerous builds"

also the fact that you run fix temp/rescale on a frequency
commensurate with the thermo output is hiding that your system
is not as well behaving as you claim it is. try setting its frequency
to 90 instead of 100 for example.

there are several more issues with your input.

you choice of neighborlist skin is odd for your
system of fast moving particles and temp/rescale
overall is a horribly bad way to adjust kinetic energy.

for efficiency reasons, you should ramp up the
timestep fairly quickly (otherwise your "equilibration"
is effectively doing nothing). e.g. like this:

timestep 0.0001
run 100
timestep 0.001
run 1000
thermo 1000
timestep 0.01
run 10000
timestep 0.1
run 100000

overall, i have no problem with NaNs or Infs. i can only suspect
that you get them by chance from those overlapping atoms
colliding with the CNT.

i also don't understand the point of computing the forces between
the CNT atoms. i would remove them with neigh_exclude, since
they don't contribute to the motions in any way.

axel.

5

Deat Axel ,

I put in my script your suggestions (or at least I tried). When I am running the simulation with a pretty small number of molecules everything works fine. As I increase the number, the overlappings create the same problem with tht overlappings. I created 4 different regions to solve at some point the problem but i just “reduced” the problem by two orders of magnitude (+18 to +16). So nothing. I know the number of molecules is big but I am trying to run a test case presented in a paper and I cannot even run it for half of the number of the molecules. In the lammps manual it is said that in case someone uses the create_random a lot of overlappings are created and it is suggested to use nve/limit but in my case it doesn’t solve my problem.

dimension 3
boundary p p f
units real
atom_style atomic
pair_style lj/cut 20.0
lattice fcc 4.0
######Data file######################################
read_data cnt.dat
######Regions########################################
region upper_container block 3.477263 70.522736 5.777263 77.922737 80 150 side in units box
region co2_1 block 3.477263 70.522736 5.777263 77.922737 80 95 side in units box
region co_1 block 3.477263 70.522736 5.777263 77.922737 96 110 side in units box
region co2_2 block 3.477263 70.522736 5.777263 77.922737 111 125 side in units box
region co_2 block 3.477263 70.522736 5.777263 77.922737 126 150 side in units box
#####Define Gases####################################
create_atoms 2 random 500 38790 co2_1 #CO2 atoms
create_atoms 3 random 500 38792 co_1 #CO atoms
create_atoms 2 random 500 1790 co2_2 #CO2 atoms
create_atoms 3 random 500 38792 co_2 #CO atoms

group tube type 1
group CO2 type 2
group CO type 3
group gas union CO2 CO

neighbor 5.0 bin
neigh_modify delay 0 every 1 check yes
neigh_modify exclude group tube tube

############Potentials################################
pair_coeff * * 0.0 3.0 0.0
pair_coeff 1 1 0.086 3.4 1.0
pair_coeff 2 2 0.468781 3.454 9.0
pair_coeff 3 3 0.065 3.920 9.0
pair_coeff 1 2 0.2007863 3.427 9.0
pair_coeff 1 3 0.0747663025 3.66 9.0
pair_coeff 2 3 0.174558772 3.687 9.0
###############INITIALIZATION########################
velocity CO2 create 300.0 4928459 dist gaussian
velocity CO create 300.0 4928459 dist gaussian

fix zUpperWall all wall/reflect zhi EDGE units box
fix zLowWall all wall/reflect zlo EDGE units box

velocity tube set 0.0 0.0 0.0
velocity tube zero angular
velocity tube zero linear
fix 1a tube setforce 0.0 0.0 0.0
fix 1b tube momentum 1 linear 0 0 0 angular
fix 1c tube rigid single

#####################System’s Minimization####################
fix 2a gas nve/limit 20.0
fix 2b gas temp/rescale 50 300.0 300.0 20.0 0.8

thermo 1000
thermo_style custom step temp press pe evdwl epair etotal
dump 1 all atom 5000 mixt_co2_co.lammpstrj
timestep 0.0001
run 100
timestep 0.001
run 1000
thermo 1000
timestep 0.01
run 10000
timestep 0.1
run 100000

6

Deat Axel ,

I put in my script your suggestions (or at least I tried). When I am running
the simulation with a pretty small number of molecules everything works
fine. As I increase the number, the overlappings create the same problem
with tht overlappings. I created 4 different regions to solve at some point
the problem but i just "reduced" the problem by two orders of magnitude (+18
to +16). So nothing. I know the number of molecules is big but I am trying
to run a test case presented in a paper and I cannot even run it for half of
the number of the molecules. In the lammps manual it is said that in case
someone uses the create_random a lot of overlappings are created and it is
suggested to use nve/limit but in my case it doesn't solve my problem.

nve/limit is not the best solution to handle overlap.
using a soft potential, which has no cusp, in combination
with fix adapt to slowly crank up the prefector A is a better
solution for cases with many overlaps. you can switch to
your potential of choice after that without a problem.

axel.

7

So I will put it with a hybrid pair style and then when I see that the overlappings have been disolved I have to unfix the fix adapt. Concerning the change of the soft potential let’s say in 1-1 interaction is it a better practice to write a restart file or with pair modify or there is no difference?

Jim

8

So I will put it with a hybrid pair style and then when I see that the
overlappings have been disolved I have to unfix the fix adapt. Concerning
the change of the soft potential let's say in 1-1 interaction is it a better
practice to write a restart file or with pair modify or there is no
difference?

no hybrid pair style needed. just run for a bit with pair_style soft
and have the same parameters for everything. with A set to zero,
this would be like simulating non-interacting particles. you can then
use fix adapt to have push atoms away from where they overlap
in a smooth fashion. then you unfix fix adapt and change the
pair style(s) to what you want to have.

check out the micelle example in the examples folder for a demo.

axel.

9

Dear Axel,

I followed your suggestions and along with the micelle example I tried to solve my problem. I tried several options but not only the problem has not been solved but it got worse. My question is focused in why is this hapenning? It gives me Potential energy +73 in comparison with the nve/limit +16.

dimension 3
boundary p p f
units real
atom_style atomic
lattice fcc 4.0

read_data pillared24_co2_co.dat

10

Dear Axel,

I followed your suggestions and along with the micelle example I tried to
solve my problem. I tried several options but not only the problem has not
been solved but it got worse. My question is focused in why is this
hapenning? It gives me Potential energy +73 in comparison with the nve/limit
+16.

that is because you are not doing it right.

#==================================================
neighbor 5.0 bin
neigh_modify delay 0 every 1 check yes
neigh_modify exclude group tube tube
#================================================
# Potentials
#================================================
pair_style soft 10.0
pair_coeff * * 0.0 10.0

why a cutoff of 10.0?? that doesn't make any sense.
the soft (core) potential is a repulsive only potential,
so its cutoff should coincide with the minimum of the
final potential that you want to use. so for your case,
something in the neighborhood of 4. angstrom should
be a reasonable choice.

variable prefactor equal 1.0+elapsed*(20.0-1.0)/1000
#ramp(0,30)

why this change?. you *want* A to be zero at the first
step and a final value of 30.0 looks pretty good when
plotting the soft potential next to your LJ potentials.

[...]

#####################System's Minimization####################
fix 2a gas nve
fix 2b gas temp/rescale 80 250.0 350.0 20.0 0.8

i *hate* temp/rescale it is *really* bad overall and it is
a particularly bad choice for your system. try this.

fix 2b gas langevin 250.0 350.0 10.0 14252346

fix 2c all adapt 1 pair soft a * * v_prefactor

thermo 100
thermo_style custom step temp press pe evdwl epair etotal
dump 1 all atom 5000 mixt_co2_co.lammpstrj

timestep 0.0001
run 1000
timestep 0.001
run 10000
thermo 1000
timestep 0.01

no reason to do all this time step juggling
with pair_style soft.

you can easily afford a time step of 2.0 or more.

run 100000
#timestep 0.1
#run 100000
unfix 2a
unfix 2b
unfix 2c

#======================================
# Main Potential
#======================================
pair_style lj/cut 12
pair_coeff 1 1 0.086 3.4 1.0
pair_coeff 2 2 0.468781 3.454 9.0
pair_coeff 3 3 0.065 3.920 9.0
pair_coeff 1 2 0.2007863 3.427 9.0
pair_coeff 1 3 0.0747663025 3.66 9.0
pair_coeff 2 3 0.174558772 3.687 9.0
#======================================
# Minimization_style
#======================================
fix 3a gas nve
fix 3b gas temp/rescale 90 300.0 300.0 50.0 0.8

again, ditch time/rescale. and just loosen the
time constant on fix langevin to reduce the
dissipative part of it, e.g.

fix 3b gas langevin 300.0 300.0 100.0 63456

axel.

11

Yes I know that it was not implemented right. in fact a value of 5.0 seems to work fine. My question is that why if not implemented right it makes things so bad (+73 etc), while in case of nve/limit the “bad” things keep a steady trend +16,+18 …

12

Yes I know that it was not implemented right. in fact a value of 5.0 seems

it just requires a minimal amount of common
sense to figure it out and, e.g., use gnuplot
to plot and compare the potentials.

to work fine. My question is that why if not implemented right it makes
things so bad (+73 etc), while in case of nve/limit the "bad" things keep a
steady trend +16,+18 ...

it is pointless to discuss why doing something
"less wrong" isn't giving a "less wrong" result.
in your case, the solution is just straightforward
but requires a systematic and clear analysis
of what you are doing and why.

i will consider this issue fully resolved now.
most certainly this is not an issue of LAMMPS,
itself but how you use it.

thanks,
     axel.