UPDATE: SOLVED (See my final comment at the bottom of this page)
Apologies for the long post, but I feel that I need to give some context/details about the problem at hand before getting to the problem. Thank you for your time in advance!
I am minimizing the electronic structure and atomic positions (ISIF = 2) of a surface slab of CdSe with pseudo-hydrogen passivation in specific places (not important for the discussion below), and the desired surface is in the plane perpendicular to the z-axis. Whenever I reach the desired force cutoff for an ISIF = 2 calculation like this one, I always perform a couple of follow-up calculations where I read from the WAVECAR and keep the relaxation tags turned on (usually IBRION = 2, ISIF = 2, and NSW = 100). If the structure does not relax any further and if the force drift remains approximately the same, then I feel confident that my structure is at the desired position.
For this relaxation, I require the forces to be no greater than 1.8 meV/Angstrom, which I know is a strict convergence criterion for a large surface slab (this one is 140 atoms). However, with parameters such as PREC and ENAUG set in appropriate ways, I found this force convergence could be achieved with a force drift that adheres to VASP manual’s recommendation --- that is, the force drift is less than the desired force cutoff.
After relaxing the structure, I performed two further “relaxations” as mentioned in the second paragraph above. The output files corresponding to those jobs are ‘OUTCAR.1208664’ and ‘OUTCAR.1211575’. Both jobs give reasonable values in the drift, namely:
total drift: 0.00000045 0.00159302 0.00053071
and
total drift: 0.00000011 0.00080124 0.00101411
respectively. There are some differences between these, but nothing out of the ordinary for such a large slab in my experience. Plus, these drifts are acceptable with respect to my desired force cutoff. Both of these runs were performed on my university’s supercomputer on 6 nodes (24 cores-per-node), and I set NCORE = 24 to match this architecture.
Now, here is the issue: If I turn off the relaxation tags (IBRION, ISIF, and NSW) and change the calculation to run on a high-memory node (48 cores-per-node) with no NCORE set, I get the following drift (see ‘OUTCAR.1232305’):
total drift: -0.00000005 0.00045769 -0.01691489
The x- and y- components are still fine, but the z-component is much larger than I would like. It’s still “small” relative to most other calculations of this scale, but it bothers me that the drift is now about an order of magnitude larger than my desired EDIFFG value. I’m not sure if it’s due to the change to the high-memory node, or turning NCORE off, or both, or neither. But it seems very strange. I mention the high-memory/NCORE requirements because I perform the VASP2WANNIER conversion, otherwise I wouldn't mention those details.
Any help on this would be greatly appreciated. I have included the aforementioned OUTCAR* files, as well as the corresponding INCAR* files (appended with the same numbers as the OUTCAR* files), KPOINTS file, POSCAR file, and POTCAR file. Also, I have slightly modified the OUTCAR* files to give me more digits of precision in certain places, so there may be some areas that seem a little more cluttered with numbers than usual.
I know this issue is not that big of a deal in the grand scheme of things, but to see the drift jump from < 1.8 meV/Angstrom to 10 times this amount between different runs is troubling to me.
UPDATE: I performed a follow-up calculation reading from the CHGCAR instead, this time back on 6 nodes with NCORE turned on. The drifts this time were excellent, on the order of 0.1 meV/Angstrom. So it seems as if the atomic positions in my slab are in good places. Also, the total energy across all calculations is nearly the same. Why is it, then, that the drift is so different for that one calculation in my post? It seems strange that the forces would be "less accurate" in 1/4 calculations for no reason, other than I changed how many cores I was using and I turned NCORE off...
I always use NPAR instead of NCORE, and I never see such a problem. On the other hand, something related to your case was observed when the same calculation was performed on different computers with different compilations. One crucial point is the cut-off energy, k-points, and sigma values, and if those values were well optimized, then you won't observe such differences. Also, if you start from a pre converged CHGCAR you will reach faster the minimum energy.
Sinhué López, thanks for the reply. I've not noticed any difference between NPAR and NCORE on my end if they are set equivalently. For example, if I run on six 24-core nodes, I can either set NCORE=24 or NPAR=6, and they each yield comparable results.
About your mention of the CHGCAR - I actually restart from the WAVECAR in these calculations, which should be just as reliable (and faster) than starting from the CHGCAR. However, I could try restarting from the CHGCAR more often and see if that helps. It's tough to do that, though, on my high-memory VASP2WANNIER conversion jobs because those already run pretty close to the 7-day max walltime when I restart from the WAVECAR.
Next, as far as encut/k-points/sigma are concerned, I have those set as accurately as they can be set without my calculations breaching the max walltime. If those are the issue, there's not much I can do about that. But, I've generally converged my settings to
UPDATE: I have an update for those who are reading this thread in the future, and I am now confident this problem is SOLVED.
My issue was that I kept ALGO = Normal (IALGO = 38) and re-minimized the electronic structure in the same step as the vasp2wannier conversion step. Instead of this, what you should do is set IALGO = 2 or IALGO = 3, and also set NELMIN = 0 in your INCAR when you want to invoke the vasp2wannier conversion process. This should go without saying, but make sure you converge everything beforehand and save that WAVECAR.
I recommend setting IALGO = 2 because this will fix the eigenvalues and keep the orbitals the same, which are both important for Wannier. I have noticed some inconsistencies with IALGO = 3 and the final *hr.dat file generated by Wannier90. These inconsistencies "average out" in the right circumstances (e.g. if you average equivalent orbitals' hopping matrix elements, they match the individual hopping matrix elements of an IALGO = 2 calculation), but this is really inconvenient to deal with. As such, if you care about Hamiltonian matrix elements, I would just use IALGO = 2.
Note to the above: You may notice total energy (TOTEN) differences for IALGO = 2 because I have noticed the Hartree energy, DENC, can be calculated incorrectly for some reason. I wish I knew why this was, but I do not. However, as long as your orbitals are preserved, the Wannier conversion process will be OK. The Hartree energy does not play any role whatsoever in Wannier calculations.
One final comment: I did compile a new version of VASP on the high-memory nodes I use. Therefore, I have 2 different compilations of VASP, and I will use one or the other depending on which type of node I'm on. This will eliminate any minor inconsistencies between the architectures of the two types of nodes, and it is something I recommend for anyone else with similar computing resources.
- Badges
- Science topic
- Similar topics
- Geoscience
- Atmospheric Sciences
- Atmosphere
- Air