Energy minimization issue and the way LAMMPS works

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps user,

I am trying to simulate diffusion of simple 1-site CO2 molecules in a solid porous carbon structure. The structure is rigid and not supposed to move.

I ran MD for a couple of million steps in order to minimize potential energy of the system using NVT-Noose Hoover. However, I noticed that potential energy of the system doesn’t decrease that much (only 350 kcal/mol over 1 million steps) and becomes equilibrated at a very high positive energy level which is clearly not correct. In fact in this type of systems, at a reasonably equilibrated system, the potential energy should be negative, considering I have only employed simple LJ potential to model interatomic interactions. Moreover, initial configuration of system comes from GCMC simulation, as such, I expect energy of the system not to be very far from a local minimum on the PES surface.

Alternatively, I have used “minimize” command for equilibration phase, but it seems that this command doesn’t care about the fixed positions of the solid carbon atoms and relocates every single atoms in the system (regardless of being fluid or solid) which actually messes up the whole system.

So, my specific questions are as following:

1- How can I minimize energy of the fluid (CO2 molecules) while the solid structure is fixed and rigid (the solid structure is supposed to interact with the fluid molecules but itself doesn’t move or deform).

2- Is the way that I have defined/reported potential energy of the system is wrong that Lammps shows such a high potential energy?

I would be grateful if anybody could give me a few comments on these two issues.

Following I have attached a copy of my input file.

Best regards,

Amir

INPUT SCRIPT ========================================

#diffusion simulation

2

Dear Lammps user,

I am trying to simulate diffusion of simple 1-site CO2 molecules in a solid
porous carbon structure. The structure is rigid and not supposed to move.

I ran MD for a couple of million steps in order to minimize potential energy
of the system using NVT-Noose Hoover. However, I noticed that potential
energy of the system doesn’t decrease that much (only 350 kcal/mol over 1
million steps) and becomes equilibrated at a very high positive energy level
which is clearly not correct. In fact in this type of systems, at a
reasonably equilibrated system, the potential energy should be negative,
considering I have only employed simple LJ potential to model interatomic
interactions. Moreover, initial configuration of system comes from GCMC
simulation, as such, I expect energy of the system not to be very far from a
local minimum on the PES surface.

you missed one important item to check: your input deck declares
periodic boundary conditions. high potential energy can be result of
periodic images of atoms getting very close, if the box dimensions are
not specified correctly.

Alternatively, I have used “minimize” command for equilibration phase, but
it seems that this command doesn’t care about the fixed positions of the
solid carbon atoms and relocates every single atoms in the system
(regardless of being fluid or solid) which actually messes up the whole
system.

if you *do* have an equilibrated system, then there should be no need
to run a minimization. that being said, with using fix setforce, you
can make certain that forces on atoms you don't want to move are zero
since atoms without forces won't move. if they still do move, then
there is likely to be an error in your input or in your understanding
of how something is supposed to work when using LAMMPS.

So, my specific questions are as following:

1- How can I minimize energy of the fluid (CO2 molecules) while the
solid structure is fixed and rigid (the solid structure is supposed to
interact with the fluid molecules but itself doesn’t move or deform).

using the term "rigid" is probably not the best choice for what you
describe. "immobile" would be a more accurate choice and as noted
above you can immobilize atoms using fix setforce or by excluding them
from time integration

2- Is the way that I have defined/reported potential energy of the
system is wrong that Lammps shows such a high potential energy?

LAMMPS is a very flexible piece of software and - like all software -
will do what you tell it to do. it doesn't know what you think it
should do. so if those two differ you get unexpected results. as noted
above, high potential energy is due to atoms being too close. this can
be due to one or more of many reasons: error in the data file, typo in
the input deck, error in the choice of units, error in the conversion
of parameters for the functional form used in LAMMPS, a bad model and
so on.

I would be grateful if anybody could give me a few comments on these two
issues.

as with many other pieces of scientific software, it is almost always
a very good idea to not start with the full problem, but first confirm
your choice of realizing your simulation on a very small toy problem,
where it is much easier to spot individual problem and then also
"build" the overall input deck by adding commands step by step.

axel.