Hi, Everyone,
I am running a simple energy minimization with lammps. I tried two different minimize style, quickmin and cg, the resultes are far different by 200 eV. I got confused, why they can be so different? And which one will be more trustable.
Thank you very much. I appreciate your help.
Minimization is a very tricky business if you have many degrees of freedom. Even more so, if your potential energy hypersurface is rugged. So how much you can trust a result depends on the latter and not so much on the minimizer.
Axel.
In addition there are several settings, like time step
for the dynamic minimizers (FIRE, QuickMin) that
affect the convergence. They are more designed
for use with the neb command.
Steve
Hi, everyone,
Thank you Axel and Steve,
Yap, I am actually running neb. What I am doing is I have found more than two energy minimum with ABC (Autonomous Basin Climb), and then I wanna use NEB to find the barriers with the help of lammps. However, my structure from ABC are minimized with cg. In my NEB simulations, by default in lammps, I need to do a QuickMin and then NEB. In the QuickMin in NEB (which is the first step in neb in lammps), there is a big jump in my energy about 200eV. So, this makes me confused, why the two minimizers have such a big difference? And because I wanna identify the potential landscape, so, I need to keep the energy continuous, the jumps make some gaps in my energy landscape, and also, the jump can cover the barriers I calculated in neb, because compare the energy just, the barriers I found with neb are kind of small. So, what should I do to reduce the jump. U have some suggestions, I appreciate your help.
Best
Your Q is not clear to me, but to use NEB effectively you typically want
the start and end point to be well-relaxed, e.g. with a CG method. Then
NEB can find a good barrier. The end-points of the relaxed NEB chain
of states should still be very close to your start and end points.
Steve
The original Q seems to have gone missing in this response thread thus excuse me if I am to restate someone else’s old point. What is exactly the system you’re trying to use NEB on? Lammps by default employs Cartesian coordinates and varies all the DOFs. If the three basins of your system are not located along narrow channels within the PES with the channels being enclosed by high energy walls, the choice of the spring constant within NEB could be crucial to your success. Unless you are applying NEB to a highly symmetric system (such as adatoms on periodic surfaces, etc), lateral motion of the elastic band should be expected. This lateral motion could result in identification of “newer” and trivial minima (like the rotation of a methyl group about the protein backbone in the alanine dipeptide case) or some “non-trivial” and undesired minima. In general, these problems do not arise when using collective variable formulations or working with smooth energy landscapes.
Anyways, my words were aimed to keep you aware of other potential complications you could encounter.
Carlos
Hi, everyone,
Sorry, I was out for several days. Now, I am back. Thank everyone for the response.
Steve, In my simulation, I did the relax for all the energy minima with cg in advance and I am sure they are well relaxed. My problem is in NEB, the start and end points are not fixed, they still have a fluctuation. So, in further TST (Transition state theory), the variation of start and end points will influence the calculation of time. So, How to set the start and end points fixed. Thank you very much.
Carlos, what do you mean by “located along narrow channel”? How could I know that. Forgive me if my question is really stupid. I am new to this area. I tried spring constant as 0.5 and 1, no big difference. 0.5 case converged a little bit faster.
I appreciate your help.
Carlos, what do you mean by "located along narrow channel"? How could I
know that. Forgive me if my question is really stupid. I am new to this
area. I tried spring constant as 0.5 and 1, no big difference. 0.5 case
converged a little bit faster.
Xin,
Your Q is not stupid. In general one doesn't know a-priori what the
morphology of your reaction path is. By narrow channel I meant that the
energy of the system increases substantially when the band deviates from
the minimum energy path. This is usually the case in the majority of the 2D
toy-surfaces that people use to theoretically "benchmark" path sampling
techniques. The use of collective variables also helps to smooth out the
shape of the PES (as one is implicitly integrating out many DOFs) which
makes NEB equilibration less sensitive to the initial choice of the spring
constant. The problem arises when you try to apply NEB to ,for example, a
molecular reaction or molecular conformational changes involving many
DOFs. Anyways, from your brief comment sounds like that is not your case.
Again, I only intended to make you aware of the limitations of the method
(which one only gets to fully learn after implementing it and trying it on
a variety of systems).
Carlos
Hi, Carlos,
Thank you again for your reply. The research work I am doing now is to find the potential landscape. I used ABC (Autonomous Basin Climb) method to find the minima (I already finish this part of work), now I am working on using neb to find all the barriers between every two minima. That the reason I need to make the potential of all the minima fixed in NEB parts. Otherwise my potential surface will not be continuous. As you said it is really different when reading the tutorial and doing the simulation myself. I do obtain some benefits from discussion with you guys. I appreciate your help very much.
In NEB as implemented in LAMMPS, there are no
spring forces acting on atoms in the first and last replicas.
Hence if they are already at an energy min, those
atoms should not move. Only that atoms in
replicas between the first and last experience
spring forces with adjacent replicas.
Steve
Hi, Steve,
Thank you very much for your answer. It is very helpful. So, for NEB in lammps, the energy of the two ends my change in the first part which is energy minimization. But in the second climbing part, the energy of the two ends won’t change right? Am I right?
However, in my simulation. I already did the cg energy minimize (EnMin) before NEB. So, because I wanna make the two ends fixed. So, I set the iteration of quickmin EnMin in NEB is 1. But, what I observed is the energy of two ends changed in climbing part.
I have another question of NEB. The Nevery in neb command only affect printing the results or it also affect the calculation. I think it won’t affect the calculation. But I found for the same simulation, nothing changed in the structure files, input file, only the value. The simulation result changed. If I set Nevery larger or equal to 6, the iterations in EnMin and climbing are both 6 and then they both converged. If I set Nevery less the 6, interestingly, it run more iterations up to thousands. I feel confused about it. I guess that reveals I don’t understand how the background cods works. I just curious about it. I suddenly found it during my simulations. Thank you very much for your help.
Xin Yan
Hi, Steve,
I search online and find your reply to another user a couple of weeks ago. Here is the information of your reply.
If both end points are already minimized, they should
not move during either stage of NEB. Of course they
could move slightly if the system minimum-energy
state is not perfect (non-zero forces). You can verify
this by having one dump file per replica and looking at
the values at the beginning versus subsequent steps
for the first and last replica.
Changing nevery
should not change the simulation, only the output.
I suggest you do these tests with the example in
examples/neb/in.neb to verify what I am saying.
Then try it for your problem.
Steve
Thank you very much. I already tried it with the example in lammps. Everything works as you said. the ends do fluatuate slightly but it is acceptable. Based on my understanding, in that example, only one atom move from the start to the final minima. However, in my problem, almost all the atoms (1020 atoms) are moving, is that because of accumulation of small errors to make the total difference which is not acceptable? Do u have some neb examples with more atoms moving?
Thank you very much for your help. I benefit a lot.
The issue you originally asked about should be independent of how many
atoms are free to move during NEB. Again, if the first and last
states are relaxed, then no atoms in those states should
move (much) during NEB. In the inbetween states any atoms
that are free to move, will move to adust to the inter-replica
and intra-replica forces they feel.
Steve
Hi, Steve,
Thank you for your reply and patience. I really appreciate your help in my research. Sorry for keeping bothering you recently. I just very curious and wanna figure out this problem. Because it related more than one projects in my PhD. I reread the papers referred in lammps manual (Henkelman’s papers). I noticed it seems that the total force add on the images include two parts, one is the parallel component of spring force and the perpendicular component of the force due to potential gradient.
In my problem, the NEB calculation is conducted under a constant tensile loading. I also found that with and without the consideration of the tensile, the NEB calculation gives quite different results. Then I want to know which one is correct, with or without tensile. I also want to know is it proper to apply external force in NEB.
Sorry for the long Email again. If my question is not clear, please let me know. I will be very glad to provide further explanation.
Thank you very much for your help all the time.
Best.
Xin Yan
I don’t know what the tensile loading aspect means. I assume at
some point during the dynamics (some amount of load), you want to
perform NEB between 2 states (with the same load). At that point
NEB knows nothing about load. It just has 2 end points, creates
the inbetween states, and does it’s calculation. So long as there are
no external forces in your script driving the end states to different
positions, there should be no problem.
Steve
Hi, Steve,
Thank you for your reply. Yes, you are correct, the two ends are under the same loads. Here us my NEB input codes:
#1. Initialization
clear
units metal
boundary p p p
atom_style atomic
atom_modify map array sort 0 0.0
#2. Atom definition
read_data Initial_frame86.dat
#3,Setting
pair_style meam
pair_coeff * * library.meam Li Si LiSi.meam Li Si
neighbor 2.0 bin
neigh_modify delay 0 every 20 check no
region 1 block -0.330193 24.2724 -0.428579 30.6005 -25 50 units box
region lower block INF INF INF INF INF 2.2 units box
region upper block INF INF INF INF 19 INF units box
group lower region lower
group upper region upper
group boundary union lower upper
group mobile subtract all boundary
fix 91 upper addforce 0 0 0.5
fix 92 lower addforce 0 0 -0.5
#4. NEB
fix 2 all neb 1
thermo 1000
min_style quickmin
timestep 0.0001
neb 1e-10 0.001 3000 10000 5 Final_frame87.dat
The two structure are energy minimized with the same load (fix 91 and fix 92) with quickmin in advanced. So, In this way, I tell NEB about my load. Is it proper from you point of view? Based on my own understanding about NEB, the energy landscapes are different for the cases with and withoud load. I also tried both cases, the results are far different. If I don’t have loads in energy minimization and NEB part, the two ends are fixed perfectly as you describe, but if I include the loads, they start to fluctuate in a unacceptable range. So, I get stuck by this problem for long time and have no idea about how to figure it out. So, do you have some suggestion if I wanna include the load in NEB? Thank you very much for your help all the time.
Xin Yan
I’ve never tried to add forces when using QuickMin. You should
read the minimize section of the fix addforce doc page. You certainly
want to specify fix_modify energy as it discusses. What I would do
is dump the per-atom forces for the first and last replica,
every timestep of minimization. If those 2 replicas are already in
minimal states as you say, then all the forces should be negligible and
hence the atoms should not move. In other words, if the
atoms in those replicas are moving, figuure out why.
Steve
Hi, Steve,
Thank you very much for the previous discussion. They are very helpful.
I have another question would like to ask. It is about fix_modify energy. This command is only affect the output of potential energy which include the contribution of fixes or not? Or this command will affect both output, I mean the printed potential value in the out put file, as well as the algorithm of the calculation.
Thank you very much for your help.
Best.
Xin Yan