Convergence issue with cation migration in bulk tio2

[ipcxj1]

Dear Prof. Henkelman,

I am new to using CINEB calculations. I am trying to study the pathway of Ti atom from 1 site to the adjacent interstitial site in Ti73O143. (The original structure is supercell composed of 3*3*2 anatase tio2 . I have created 1 O vacancy and relaxed the structure. IS was also relaxed with the same INCAR). In the past few weeks, I have read most of the posts about speeding up convergence （using conservative optimizer of iopt=3, cheap setup). But I failed and I don't know why.

I noticed that you mentioned the problem of image translation. I checked my images and nothing seems to be happening. Can you please confirm it for me? Also, I noticed your comment about freezing an atom for slab model, but not for the bulk materials. Therefore, I would like to ask you to confirm whether this method is suitable for the bulk materials. Although I have tried frozen an atom, it still failed to converge.

Another question, when constructing the bulk model with a vacancy, is there a difference between the O vacancy constructed in the center of the supercell and the area near the surface of the supercell?

Could you please help me understand why this is happening and what I must do to correct this ? I have attached the files here.

Thanks for your help and support.

Jing Xia

[graeme]

These calculations are generally looking ok. I have just a few suggestions:

1. relax the cell without a vacancy (ISIF=3) and then leave the cell fixed for your IS and FS optimizations (ISIF=2) so that the cell does not change over the NEB.

2. I see a small shift in the structures, which is not really a problem, but could be avoided by freezing an atom.

3. These are periodic calculations, so it doesn't matter where the vacancy is in the cell.

4. Your MIDDLE calculation looks reasonable. The max forces are below 0.1 eV/Ang, and with some additional iterations there will be tighter convergence. To speed things up, you can try IOPT=1 or IOPT=2. But what you are doing (IOPT=3) is a nice slow and steady optimizer.

5. I don't get what is going on in the HIGH calculation; it is not using our vtstcode. I guess you are trying to optimize the NEB with more kpoints, which is fine, but you will need to link in our vtstcode in your vasp binary as you did with vasp_gam.

[ipcxj1]

Dear Prof. Henkelman:
I didn't expect to receive your reply so quickly, thank you very much.
I haven't slept well for weeks because of this problem. I have tried almost all the methods mentioned in the forum.

In fact, in order to get most stable structure （IS and FS）, I have already tried to relax the cell without a vacancy (ISIF=3 or 2) and then leave the cell fixed for IS and FS optimizations (ISIF=2 or 3), but I got the same result failed.

" I don't get what is going on in the HIGH calculation; it is not using our vtstcode. I guess you are trying to optimize the NEB with more kpoints, which is fine, but you will need to link in our vtstcode in your vasp binary as you did with vasp_gam..". Yes ,I am trying to optimize the CINEB with more kpoints using the CONTCAR results from MIDDLE. I don't understand why you said I didn't link the two programs. How did you come to this conclusion. In fact, I have linked the vtstcode to our vasp5.4, no matter vasp_gam or vasp_std. Maybe my problem just comes from the linking bewteen the vasp_std (I renamed it as vasp) and vtstcode. But I don't know how to judge whether the linking is successful.

Thanks for your help and support.
Jing


Below is the script I used for cineb.

#!/bin/sh
##job1-LOW_cineb
cd LOW
pwd
rm cineb_time1
time mpirun -np 64 vasp_gam
T=$(grep "CPU" 01/OUTCAR | tail -1 | awk '{printf "%-12.6f \n", $6}')
T1=$(echo "scale=2; $T/3600" | bc)
echo cineb_time1=$T1 h >> cineb_time1
cd ..

##job2 MIDDLE_cineb
cd MIDDLE
pwd
rm cineb_time2
cp ../LOW/cineb_time1 cineb_time2
for i in {0..9};
do
rm -rf 0$i
cp -rf ../LOW/0$i 0$i
cd 0$i
mv CONTCAR POSCAR
cd ..
done
time mpirun -np 64 vasp_gam
T=$(grep "CPU" 01/OUTCAR | tail -1 | awk '{printf "%-12.6f \n", $6}')
T2=$(echo "scale=2; $T/3600" | bc)
echo cineb_time2=$T2 h >> cineb_time2
cd ..

##job3 HIGH cineb
cd HIGH
pwd
rm cineb_time3
cp ../MIDDLE/cineb_time2 cineb_time3
for i in {0..9};
do
rm -rf 0$i
cp -rf ../MIDDLE/0$i 0$i
cd 0$i
mv CONTCAR POSCAR
cd ..
done
time mpirun -np 64 vasp
T=$(grep "CPU" 01/OUTCAR | tail -1 | awk '{printf "%-12.6f \n", $6}')
T3=$(echo "scale=2; $T/3600" | bc)
echo cineb_time3=$T3 h >> cineb_time3
total_time=$(echo "scale=2; $T1+$T2+$T3" | bc)
echo total_time=$total_time h >> cineb_time3
cat cineb_time3

[graeme]

Oh goodness, you shouldn't lose sleep over the calculation of vacancy diffusion! This is not a difficult calculation - it is just a matter of making a few adjustments.

First, you can check if the vtstcode is linked in by greping for VTST in any OUTCAR file and look for a line such as: "VTST: version 3.2, (02/03/18)". You can find it in your MIDDLE calculation, using the vasp_gam binary, but not the HIGH calculation, using the vasp binary.

I also can't help but make a few suggestions which everyone in my group has to repeatedly endure. I strongly recommend starting with cheap computational settings and then reconverging to test your parameters. In your case, starting from a gamma point calculation and then going to higher k-point sampling is exactly this idea. But there are other such system parameters including the cell size, choice of pseudopotential and energy cutoff. In your case, I would recommend starting with a regular Ti potential (not _sv), as well as O_s (not O) and a cutoff of about 274 eV. I would also cut down your cell in the z-direction or find another way to find a smaller roughly cubic cell. This changes would save you more than an order of magnitude in computational time, and make it so much faster to do your saddle calculations. Once you have a nicely converged NEB calculation, it is trivial to try a different pseudopotential or increase the cutoff or add kpoints or whatever parameters you are worried about, and also then learn if those parameters are important to the quantity that you care about, such as the diffusion barrier. Starting with the most expensive parameters that you can afford, on the other hand, makes all the calculation slow so that you might not be able to converge things, and you don't have any sense for which settings are important for the precision/accuracy that you have. Putting this simply: start cheap and get quick convergence and then test your settings in terms of the quantity that you care about.

Also, in terms of the science, it appears to me that you are moving a vacancy to a site that is ~2eV less stable in the FS as compared to the IS and there is just a small reverse barrier. So thinking about what is important here, it may be primarily the energy difference between the IS and FS.

[ipcxj1]

Dear Professor Graeme,

Thanks a lot for your reply and help. It did work!

This morning, I checked VTST in the OUTCAR of HIGH. You are right, there is no link bewteen the vasp and the vtst code.

So I try to ran CINEB in HIGH alone and get a new OUTCAR. Interestingly, I can find ' "VTST: version 3.2.." in the new OUTCAR. I tried several times and I can find the VTST characters in the new OUTCAR. I don't know why this happens, but I know the source of the problem and I think I can fix it.

Once again, I would like to express my sincerest thanks for your help.

Jing

[ipcxj1]

Dear Prof. Henkelman,
Following your suggestion, I rebuilt the supercells（tio2,4*4*1 and 3*3*1）. In the supercell 4*4*1, a Ti atom is frozen and no atom is fixed in supercell 3*3*1. I check the vtst, it is linked to vasp. Unfortunately， both cases can not converge.

Do you think the reason for not converging is because there are only shoulders? For the case where only the shoulder exists, and no transition state, the calculations are meaningless?

Could you please help me ? I have attached the files her.

Thanks very much for your help.

Jing

[graeme]

I think there may be a misunderstanding about convergence. Take, for example, your 331/CINEB_NO_FIX/LOW , you can see a single barrier around 4.9 eV. It is true that this is not converged in the sense that the forces are around 0.4 eV/Ang and you have requested maximum forces of 0.05 eV/Ang. But that is because you are using a conservative optimizer (IOPT=3) and allowed only 200 iterations. So you just need to continue the calculation. You copy the CONTCAR files to POSCAR and restart vasp. You might use one of our scripts: "vfin.pl dir" , which will move the calculation to dir and have the main directory set for a new calculation.

Then, once your LOW calculation is converged, then you can go to MIDDLE. Since LOW will be closer to the minimum energy path, MIDDLE should converge faster. But again, if it does not in 200 iterations, start it again until it does. Then you can move on to HIGH. Of course, you can also increase NSW, but the point is that you can't demand that the calculation converge in a fixed number of iterations; instead you need to check that the calculation is reasonable and converging and then continue them until they do converge.

[ipcxj1]

Dear Professor Graeme,

Thank you very much for your reply. In my experience, I think 200 is enough for cineb. If not , the mode, path or setup may be not good. However, I will continue the
calculation and tell you the results.
Thanks again.

Jing

[graeme]

No, it really depends upon the path, and in particularly how linear it is or how close your initial path is. I have seen calculations of hundreds of paths and some range of simple one which converge in a few iterations with one image, and ones that require 10-00 images and many more iterations with a careful optimizer to converge. A lot depends upon the ratio of soft to stiff modes, i.e. the stiffness of the problem. But anyway, the bottom line is that you have to check for convergence - you can just run a fixed number of iterations and conclude that things are not working if you don't get to a force that you are looking for.

[ipcxj1]

Dear Professor Graeme,

I got it. Thank you for your reply.

Have a good day。

Jing