Best way to run dimer in Vasp?

[buber]

I'm trying to find if an oxygen di-interstitial in UO2 can move in a concerted way. I already tested an initial guess with NEB and found it decomposed into two single hops. So, I'm wanting to do some dimer searches.

I don't want to bias the searches in terms of guess a saddle point, though I can say that I know what atoms are most involved. What would be the best way to run a dimer search then? I don't necessarily want to do a MODECAR that points in some "intuitive" way. Would it be best to have a random vector that includes the atoms I think are important in the MODECAR?

Thanks,
Blas

[graeme]

If you're willing to invest a little time to try a new system, I think our akmc script is the best way to go. We could also use some feedback - if there are parts of it that don't make sense, or are confusing, we'll fix, change, or add documentation for them.

The way this script works is that you first specify your initial poscar, and the atoms that can move. This is done in a DISPLACECAR_sp file (the same as your dynmat code), indicating the atoms that are initially displaced and the standard deviation of the Gaussian displacements. Then, the script generates a set of initial configurations for dimer searches, and submits these jobs on a cluster. If you want to run this on our local cluster, the submission files are already generated, but it is also fairly easy to change them for a different queuing system.

You can also use the akmc script to generate the dimer jobs, and submit them yourself, manually. But this can get old rather quickly, so if you use the script, it will take some of the tediousness out of the calculation. Lijun has it set up to check if jobs are stopped, and restart them if they are not finished. It will also take saddle points that are found and start minimization jobs to see if they connect back to the initial state. Then, if a convergence criteria is reached for an initial state (enough saddle points are found) it will do a KMC step and simulate the state-to-state dynamics.

Finally, if you just want to try some random searches, you can do this in the same way that the akmc script does, using the diminit.pl script. This script will move your chosen atoms (in a DISPLACECAR_sp) file by an average amount given by the values in that file, and also set the MODCAR to point from the minimum to that initial configuration.

Good luck.

[buber]

Thanks for the advice Graeme. I may start with the simpler method you mentioned at the end, but if I find some time, I'll try out the akmc scripts.

[buber]

Quick question: I set up the DISPLACECAR_sp file to have two moving atoms. In the dimer displacements it makes, it seems that only one is displaced at a time. Is that what is meant to happen, or have I messed something up?

[xulijun2]

[quote="buber"]Quick question: I set up the DISPLACECAR_sp file to have two moving atoms. In the dimer displacements it makes, it seems that only one is displaced at a time. Is that what is meant to happen, or have I messed something up?[/quote]

Right now we are displacing atoms within a radius around a central atom. the central atom can be chosen via several algorithms.

USAGE: diminit.pl DIRECTORY OR NUMBER (OPTIONAL: DisplaceAlgo DisplaceRange Rcut MaxCord POSCAR DISPLACECAR_sp)

The first argument can be either a directory name or a number and either one or several dimers will be created , respectively.

The second argument is the displacing algorithm: 0 for choosing the most undercoordinated atom as the center atom; "1" for randomly choosing an atom whose coordination number is smaller than the specified maximum coordination number, as the central atom. "2" and "3" are just for surface events: it divides the surface atom into different islands and treat each island just it does with the option of "0" and "1". Using other numbers means a random atom (which is allowed to move in the displacecar) will be chosen as the central atom.

"Rcut" is the cutoff range in deciding the coordination number for each atom;

"Maxcord" is the specified maximum coordination number just mentioned.

Only those atoms set to be displaced in the DISPLACECAR_sp can be displaced.

example: diminit.pl 5 0 4.0 3.0 8 POSCAR DISPLACECAR_sp
This will create 5 dimers by randomly choosing one of the most under-coordinated atoms (cutoff is 4.0 Angstrom) from the to-be-displaced atoms in the DISPLACECAR_sp as the central atom and displacing atoms within the radius of 3.0 Angstrom.

example: diminit.pl 5 1 4.0 3.0 8 POSCAR DISPLACECAR_sp
This will create 5 dimers by randomly choosing an atom whose coordination number is smaller than 8 (cutoff is 4.0 Angstrom) from the to-be-displaced atoms in the DISPLACECAR_sp as the central atom and displacing atoms within the radius of 3.0 Angstrom.

example: diminit.pl 5 4 4.0 3.0 8 POSCAR DISPLACECAR_sp
This will create 5 dimers by randomly choosing an atom from the to-be-displaced atoms in the DISPLACECAR_sp as the central atom and displacing atoms within the radius of 3.0 Angstrom.

If you just want to simultaneously displace all the atoms marked as to be displaced in the DISPLACECAR, you can set the cutoff and displacement radius to a high value so that all of them can be displaced no matter which atom is chosen as the central atom.

You can try again and please let me know if it works. Thanks.
i know this script can be certainly improved; your feedback is also appreciated. -- Lijun

[buber]

Thanks. I think I got something to work. I appreciate the help.

So, I started a couple of runs and got about 45 ionic iterations before running out of time. During the runs, atoms moved at most 0.3 angstroms. Forces are still pretty high. Starting with random displacements, how many ionic iterations might I expect to have to do to reach convergence? Just wanted an idea.

Thanks!

[graeme]

Can you post the DIMCAR - this is a short file that will give a good idea of what is going on. There is some overhead to get the dimer oriented at the start, but I'm a little surprised that it moved so little after 45 force calls. Are you using conjugate gradients (IBRION=3, POTIM=0, IOPT=2)? To converge onto a saddle (convergence of 0.002 eV/Ang), it typically takes 400-500 force evaluations.

[buber]

Here is the DIMCAR from one of my runs:

Step Force Torque Energy Curvature Angle
1 6.74769 17.61536 -1002.94132 14.73217 44.27777
1 6.74769 7.41450 -1002.94132 3.57505 42.55853
1 6.74769 4.97267 -1002.94132 2.32063 33.29783
1 6.74769 6.96414 -1002.94132 3.73984 17.54976
2 6.54296 5.95716 -1002.75720 1.27683 35.16600
2 6.54296 7.20395 -1002.75720 2.66888 22.88683
2 6.54296 7.00910 -1002.75720 2.29922 19.87855
2 6.54296 5.45499 -1002.75720 1.12817 17.22603
3 6.39192 5.80065 -1002.40962 -2.10497 25.84620
3 6.39192 7.21336 -1002.40962 0.47974 14.44385
3 6.39192 7.38443 -1002.40962 -1.21334 21.49611
3 6.39192 9.43741 -1002.40962 0.11785 25.42363
4 6.37884 7.81631 -1001.98600 -5.72580 20.65006
4 6.37884 9.63237 -1001.98600 -1.99931 27.49938
4 6.37884 8.71196 -1001.98600 -2.51490 22.79489
4 6.37884 6.82406 -1001.98600 -4.32456 16.19636
5 3.76125 9.79255 -1002.49825 -1.81659 35.59228
5 3.76125 8.61818 -1002.49825 1.13851 16.66160
5 3.76125 7.89576 -1002.49825 -1.16182 21.04477
5 3.76125 9.02551 -1002.49825 -0.42740 23.85460

and the INCAR from that run:

IBRION=3
NSW=100
PREC=med
LCHARG=F
LWAVE=F
EDIFFG=-0.01
ISYM=0

ICHAIN=2
POTIM=0
IOPT=2

So, I am using the option you mentioned. This run did 50 force evaluations.

[graeme]

It looks to me as if the forces are not accurate enough to properly evaluate the curvature and rotational force. When things are working, the Torque should drop during the 4 rotation iterations that you are allowing per step, and the Curvature should also drop systematically. In steps 3-5, this is clearly not happening.

Setting ediff=1e-6 or 1-e7 will probably do the trick; the default value will not work. Also, since this looks like a large system, you might be able to go a lot quicker setting LREAL=.AUTO.

[buber]

Thanks Graeme. I'll try those. I guess I should have known about the EDIFF issue.