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