LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] How should I insert one molecule into the box?
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] How should I insert one molecule into the box?

From: Jacob Gissinger <Jacob.Gissinger@...780...>
Date: Wed, 28 Mar 2018 12:04:02 -0600

Hi Wei,

The force is a just a function of distances, angles, etc. between atoms (if interested, you can look up these formulas on the relevant bond_style, etc. lammps pages).

The error you received simply means you still have an unstable configuration. 

Tips for using nve/limit: use as large an 'xmax' value as possible while maintaining stability. You may also need to run the system in nve/limit for longer.

More tricky situations, which could arise e.g. if you have rings in your molecule, can only debugged by directly visualizing the misbehaving atoms.


On Wed, Mar 28, 2018 at 3:49 AM, Wei Chen <chenwei@...7134...> wrote:

Dear Jacob and Andrew,

I used "fix nve/limit"  instead of "minimize" to relax the system.

After several steps, the energy looks reasonable. But some atoms still have very large force, and when shift to normal dynamic, like "fix nvt", the LAMMPS  immediately crashes due to the large force.

Does "fix nve/limit" not update the force ? According to the manual, this command only resets the velocity, the force is kept unaltered. Am I right?

Or I have to reset the force by "fix setforce" following "fix nve/limit"?

I am looking forward to your reply.

Best wishes,


On 03/24/2018 07:22 AM, Jacob Gissinger wrote:

Wei: most importantly, you should figure out how to visualize the unstable atoms of your system, particularly every timestep near the simulation failure. this will likely clarify things quickly. you should look at Axel's 'readvarxyz' feature in VMD

Andrew: I would also add that 'nve/limit' seems particularly useful. this is quite stable, operates at finite-temperature (unlike minimize), and allows for groups


On Fri, Mar 23, 2018 at 4:59 PM, Andrew Jewett <jewett.aij@...24...> wrote:
On Wed, Mar 21, 2018 at 1:58 AM, Wei Chen <chenwei@...7134...> wrote:
> I don't simply want an initial starting structure, what I did is just a
> test. But insertion is an important point of my research, I would be
> grateful if someone would solve this issue.

My apologies.
What you are trying to do is difficult.

The "minimize" command, in its current form, is not ideal for what you are trying to do (unless you try modifying the forces in your system).

In my experience, the "minimize" command usually fails when there are extremely large forces in your system (for example when using 1/r^12 Lennard-Jones-style repulsion).  The "minimize" command will fail unless the distance between all pairs of atoms are at least 10%(??) of the sum of their their radii.  (It's not possible to estimate the exact distance when failure occurs, but it definitely happens at non-zero distances.)  This means if you insert a molecule randomly with respect to the other molecules, eventually you will crash LAMMPS.  Even if the system is dilute, eventually you will have a collision that causes a pair of atoms to be too close together.  I am not sure whether this can be alleviated by modifying the source code for the "minimize" command.  I suspect it is an inherent limitation of double precision floating point and 1/r^12.

To help get around these problems, for the moment I have been using pair_styles which do not necessarily diverge to ∞ as r -> 0.  (Such as pair_style soft, pair_style gauss, or pair_style lj/charmm/coul/charmm/inter, which you can currently download here.  See the example using "rsoftcore" here.  Let me know if that code no longer compiles.  It's kindof funky.)
Using soft potentials does not always solve the problem, but it seems to help.
In this case, it should be either possible to minimize the system using the "minimize" command, or run an ordinary MD simulation with a small timestep and high damping (using some combination of "fix langevin" and "fix nve").  If you attempt molecule insertion every 10000 timesteps, then you will have to write a loop which minimizes the system, right after the insertion every 10000 timesteps (perhaps changing the pair style each time).

---- other issues that make minimization difficult (after insertion) ----

2) The "minimize" command effects all of the atoms in your system.  (There is no "group" argument)  The entire simulation has to halt and wait for minimization to finish.  In your case, it could be more efficient to only minimize atoms which are located nearby where you inserted the molecule.   Doing this would require writing some kind of new fix, or (better) some function which other fixes can invoke after a sudden change in the system(such as fix bond/create, fix bond/break).  I'm not yet certain whether it would be worth the effort.  It depends on how frequently these insertions occur, and in how many different locations.

I run a lot of simulations here with bonds being created and destroyed.  I have to use small timesteps to avoid numeric explosions each time I add or delete a bond.  At some point, I may decide to attempt writing a minimizer which solves some of these issues, but I haven't had the time to focus on the issue yet.  I recently saw an interesting discussion regarding the new "fix bond/react" feature being contributed by Jacob Gissinger.  This fix a group of atoms which are effected every time a bond is created or destroyed.  (But I forget the details.)  Perhaps he eventually planned to implement a minimizer that would act on only this group?  In the distant future, I was considering creating a new fix based on my code or his code which would create or insert molecules at locations depending on the nearby environment.  Perhaps we should collaborate.

3) The "minimize" command is not compatible with constraints like "fix rigid".  (Perhaps this is not an issue for you.)

> 2. Both minimization procedures in case 1 and case 2 give the similar values
> of those energies. The program crashes during the dynamic procedure.
> In case 2, the minimization after insertion gives E_pair = -12414.855 at
> 10707 step, but in the dynamic, the E_pair suddenly changes to 1450.7209 at
> 10707 step. In case 1, E_pair doesn't change after minimization.

I'm not sure why.  Again, weird behavior like what you are reporting is common in LAMMPS when using the "minimize" command with overlapping particles.
It used to bother me, but now I have gotten used to it.

I would guess that E_pair = -12414.855 is true at the beginning of the timestep (Epair(x(t)) and E_pair=1450.7209 at the end of that timestep (Epair(x(t+Δt)).

Displaying the energy (printing thermo to terminal or log.lammps) is not done until after the timestep has completed, which includes code which causes the particle to move.  (See page 14 of

To test whether this is the reason, try setting the velocity of all of the atoms to 0 before using the "run" command", and try using the simplest possible integrator such as "fix nve" (not "fix nvt" or "fix npt").  Also try reducing the timestep to 0 (or a very small value if LAMMPS won't allow timestep=0).


I don't know if this email reply helped you with your issue.  Hopefully it is comforting to know that other people suffer too.  Again, what you are trying to do is not easy.


> On 03/21/2018 02:26 PM, Andrew Jewett wrote:
> On Tue, Mar 20, 2018, 8:57 PM <chenwei@...7134...> wrote:
>> Dear Axel,
>> Thanks very much for your reply. I have tested LAMMPS in several ways. I
>> think it is probably a bug, as the following experiments show:
>> There are 3000 SPC/E water molecules and one surfactant in the simulation
>> box in its initial state. And the atoms in the initial state are highly
>> overlapped.
>> Case 1. I insert another surfactant after the "read_data" command, and
>> minimize the energy of the system, then run the dynamics. Everything is OK.
>> Case 2. The dynamic has run for 10000 steps from the initial state, then I
>> insert another surfactant and minimize the energy of the system. Fail!!!
>> And I find that the Minimization gives E_pair = -12414.855, however when
>> the program shifts to the dynamic procedure, the E_pair is 1450.7209.
>> Step Temp E_pair E_mol TotEng Press
>>    10530    289.40022 3.9184693e+10    1610.0117 3.9184697e+10
>> 8.5978858e+10
>>    10600    289.40022   -11863.682    927.10885   -8245.9764   -213.04207
>>    10700    289.40022   -12423.715    924.48632   -8808.6327    39.802331
>>    10707    289.40022   -12414.855    911.33202   -8812.9266    144.03736
>> Loop time of 3.70799 on 20 procs for 177 steps with 3120 atoms
>> 99.0% CPU use with 20 MPI tasks x no OpenMP threads
>> Minimization stats:
>>   Stopping criterion = energy tolerance
>>   Energy initial, next-to-last, final =
>>          39184694708.2     -11503.4662055     -11503.5228676
>>   Force two-norm initial, final = 1.55345e+12 100.973
>>   Force max component initial, final = 8.83876e+11 43.1717
>>   Final line search alpha, max atom move = 0.00269996 0.116562
>>   Iterations, force evaluations = 177 319
>> Step Temp E_pair E_mol TotEng Press
>>    10707    289.40022    1450.7209    911.33202    5052.6492    30604.512
>> ERROR on proc 1: Bond atoms 3069 3095 missing on proc 1 at step 10708
>> (../ntopo_bond_partial.cpp:64)
>> E_pair is not the same at 10707 step!
>> Case 3. In Case 2, I delete the overlapped atoms. Even if I delete all
>> water molecules, only keep two surfactants. Still fail!  Have the same
>> problem on E_pair.
>> Case 4. 3000 SPC/E water molecules and one surfactant run for 10000 steps,
>> and then write a data file named "" by "write_data" command. Then
>> I start a new simulation taking the "" as the initial state,
>> insert another surfactant after the "read_data" command like the case 1.
>> Everything is OK!  Actually I combine Case 1 and Case 2 by splitting the
>> procedure into two simulations.
> Just a quick email for now.
> The fact that it does not behave the same way it did in Case 2 might be
> because writing a data file and reading it again does not garuantee it will
> be in the same state.  (Rounding error.  The velocities are probably
> forgotten.)  Try using write_restart and read_restart.  But this will not
> resolve the main issue.  (See below...)
>> In summary, it is OK to insert one molecule before the dynamic starts, but
>> it fails to do that when the dynamic has run for a while.
> The conclusion I draw from this is that the current minimizers are not
> robust enough to reliably cope inserting a molecule suddenly into a dense
> system.  At least that's been my experience with highly overlapping
> polymers.  Minimization of overlapping systems is difficult at best.  I
> sometimes minimize using pair_style soft or pair_style gauss.
> If your goal is simply to creat an initial starting structure, why don't you
> use one of the molecule builders to control exactly where you want each
> molecule to go?
> Andrew
>> Best wishes,
>> Wei
>> > -----Original Messages-----
>> > From: "Axel Kohlmeyer" <akohlmey@...24...>
>> > Sent Time: 2018-03-20 20:24:42 (Tuesday)
>> > To: chenwei@...7134...
>> > Cc: "LAMMPS Users Mailing List" <>
>> > Subject: Re: [lammps-users] How should I insert one molecule into the
>> > box?
>> >
>> > On Tue, Mar 20, 2018 at 6:04 AM,  <chenwei@...7134...> wrote:
>> > > Dear LAMMPS users,
>> > >
>> > >
>> > > My system is quite simple, there are 3000 SPC/E water molecules and
>> > > one
>> > > surfactant in the simulation box in its initial state. After the
>> > > simulation
>> > > has been performed for a while, I want to insert another surfactant
>> > > into the
>> > > box.
>> > >
>> > >
>> > > I used the "create_atoms" to do this job. But I always got the error
>> > > "Bond
>> > > atoms # # missing on proc 1 at step #".
>> > >
>> > >
>> > > I test my script. If I insert another surfactant before the "run"
>> > > command,
>> > > there is no problem. But if I do this after the "run" command, I get
>> > > the
>> > > error.
>> >
>> > create_atoms does not check for overlaps. so you may have to delete
>> > atoms/molecules that overlap with the newly inserted molecule.
>> >
>> > axel.
>> >
>> >
>> > >
>> > > Please find all the input files from the attachments.
>> > >
>> > >
>> > > Thank you!
>> > >
>> > >
>> > > Best wishes,
>> > >
>> > >
>> > > Wei
> --
> State Key Laboratory of Multi-Phase Complex Systems
> Institute of Process Engineering (IPE)
> Chinese Academy of Sciences (CAS)
> P.O. Box 353, Beijing 100190, China
> Tel: +86-10-8254 4839
> Fax: +86-10-6255 8065
> Email: chenwei@...7134...
> Http://


State Key Laboratory of Multi-Phase Complex Systems
Institute of Process Engineering (IPE)
Chinese Academy of Sciences (CAS)
P.O. Box 353, Beijing 100190, China
Tel: +86-10-8254 4839
Fax: +86-10-6255 8065
Email: chenwei@...7134...