[lammps-users] minimization of graphitic layers


sorry the delay in following up with this: the computer
were i ran the calculations was down (Franklin, at nersc).

I have attached here three files:

first.log and second.log are the outputs for the run. The
difference between first and second is this:

in first.log, the minimize commands is

minimize 0.0 1.0e-5 100 100

whereas in second.log, the minimize command is

minimize 0.0 1.0e-5 100 1000

That is, the runs differ only in the value of maxeval (100
in first, 1000 in second).

As you can see, first.log converges in 9 steps. Second.log does

I also attached the file data.in, which contains the structure.


data.in (152 KB)

first.log (2.15 KB)

second.log (7.65 KB)

When you run with a limit of 100, you get this output

  Stopping criterion = max force evaluations
  Iterations, force evaluations = 9 104

meaning you hit the maxeval critterion.

When you run with a limit of 1000, it runs longer.
This is normal behavior. You may not like the
answer/structure you're getting as the minimization
proceeds, but the energy is dropping dramatically,
so the minimizer thinks it's doing its job.



Well, this is clear, but the answer lammps is given is still very
unclear to me. So to better understand where the source of the problem
could be, i did more tests.

I have taken the conventional cell of a graphite crystal (8 atoms per
cell; orthogonal cell) and expanded it 3 times in each direction, so as
to produce an orthogonal cell with 216 atoms.

I minimize this cell using the LJ potential

pair_style lj/cut 8.0
pair_coeff * * 0.002413 3.4
pair_modify shift yes

and the minimization when perfectly: It finishes well, very fast, and the final
structure is graphite.

Now, the same is not true if i put a vacuum (make z very large, i.e. 100)
and continue using periodic boundary conditions. In this case, the minimization
does finish, and converges, but the final structure is not graphite. Instead, it
looks like graphite got stretched in the z-direction. Why is this happening? i
don’t know.

I can see three reasons why this could be happening : 1) the LJ is
not working for a surface; 2) there’s a bug in the minimizer; 3) i’m not
using lammps properly.

I’m trying to run these tests with another MD code and the same LJ, just
to clarify question 1). For question 3) i still have not found a clue in
lammps’ manual. For question 2), unlikely for sure, i do not have
an answer, but i’d be more than willing to provide the files i’m
using to run these calculations.


I think the answer to your question is mostly (1). The LJ potential won’t hold the layers in a bulk structure because you have too few atoms in the system to form the necessary interactions. The first few layers at the surface of any material are often not the same structure as the bulk material. By making your system essentially all surface, you won’t see the expected bulk structure.

Therefore, if you need the bulk structure, you must impose periodic boundary conditions that mimic a bulk structure. Using a large vacuum layer mimics a surface because the PBC don’t help with the necessary interactions if the atoms are too far apart.


Cutting off the LJ potential will interfere with the true dispersion force as well. The cutoff may have been picked to work in a particular bulk environment, with a certain layer spacing. Layer spacings > 8 won’t have any contribution from the LJ potl given below.



Is there a way to impose a constraint minimization? In particular, can one minimize a certain group of atoms and
leave the other group frozen?

I can see this is possible with an MD, but it seems to not be possible with a minimization.


Did you try the fix setforce command (Full description at http://lammps.sandia.gov/doc/fix_setforce.html)??)

From the documentation:

Set each component of force on each atom in the group to the specified values fx,fy,fz. This erases all previously computed forces on the atom, though additional fixes could add new forces. This command can be used to freeze certain atoms in the simulation by zeroing their force, assuming their initial velocity zero.

The forces due to this fix are imposed during an energy minimization, invoked by the minimize command.


Read the doc pages on minimize and various fix commands that can
be invoked during minimization, including the setforce that Joanne mentioned.




I want to follow-up a bit on the minimization issue. The answer to the problem was
number 3), that is “i’m not using lammps properly”. Specifically, I was using a LJ alone
for the graphite, and that, of course, did not work. AIREBO, on the other hand, works
very nice.

But i’d like to ask another 3 small questions.

I’m interested on investigating a metallic nanoparticle on top of graphite. The metallic
nanoparticle (atom type 1) is described with an EAM potential; the graphite (C is atom
type 2) is described with AIREBO; the interaction between the nanoparticle and the
graphite is described with a LJ.

I’ve used this set of commands to describe the interactions. Are they correct? (I’m a
bit concerned i did not use the “pair_coeff * * airebo” properly.)

pair_style hybrid airebo 3.0 1 1 lj/cut 8.0 eam
pair_coeff * * airebo /scratch/scratchdirs/m96/eam/potentials/CH.airebo NULL C
pair_coeff 1 1 eam /scratch/scratchdirs/m96/eam/potentials/niu3
pair_coeff 1 2 lj/cut 0.023049 2.852
pair_modify shift yes

Another question is the following: In the pair_style lj/cut, is the cutoff the “(sigma scale factor)”,
as in the pair_style airebo.

And finally, i’m running a minimization with the system mentioned above, with the pair-
interactions described above as well. All is working, but not always. Specifically, I notice
that depending on the number of processors use, the code runs or crashes. I suppose this has
something do to the number of atoms are distributed among nodes.

In particular, i have 7169 atoms. I’m running tests in 4 processors. 8 and 16 processors
do not work. Is there any reason for this?


Hi Miguel,

Which version of LAMMPS are you using? There was a bug using AIREBO
with hybrid pair potentials in pre-20Aug08 versions of LAMMPS.


Hi Zhun-Yong,

thanks! I’m using this version: LAMMPS (21 May 2008)
and i assume the patch for airebo is not there.

Is this the reason why the code crashes when using different
set of processors?

This is the message i get when the code crashes:

Hi Miguel,

Yes. That's probably why your script crashed. Patching up your LAMMPS
code should fix things.


Your pair_coeff commands look correct.
The lj cutoff is what you set it to, 8 Angs.
The crash is probably due to the bug that
was subsequently patched.