|From:||Jacob Gissinger <Jacob.Gissinger@...780...>|
|Date:||Wed, 28 Mar 2018 12:04:02 -0600|
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.
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.
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.> 2. Both minimization procedures in case 1 and case 2 give the similar values
3) The "minimize" command is not compatible with constraints like "fix rigid". (Perhaps this is not an issue for you.)
> 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 http://lammps.sandia.gov/doc/D
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
>> 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
>> 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
>> 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 "rerun.data" by "write_data" command. Then
>> I start a new simulation taking the "rerun.data" 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?
>> Best wishes,
>> > -----Original Messages-----
>> > From: "Axel Kohlmeyer" <akohlmey@...24...>
>> > Sent Time: 2018-03-20 20:24:42 (Tuesday)
>> > To: chenwei@...7134...
>> > Cc: "LAMMPS Users Mailing List" <lammps-users@...42...
>> > 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...
-- 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://www.mpcs.cn/emms