Regarding the convergence of structure optimization

[mukesh.phy]

Dear VASP users and admin,

I am exploring the different ways available in vasp for structure optimization as I have many structures ( Not so many that that I can call them high throughput), and convergence of each of them are taking time. A bit more concrete ways, I am performing calculation on single atom catalysts for Li-S battery. In these simulations, Li2Sx molecules, where x = 1, 2, 4, 6, 8 absorbs on Fe/Ni/Co/Mn doped in 2D materials. And I am looking forward the suggestions to perform fast structure optimization. And asking some questions?

Q1- At the vasp website:
https://vasp.at/wiki/Structure_optimization
the diagram says few atoms, and degree of freedom > 20. what is meaning of these two term in term of number of atoms.

In my opinion, if system has 7-8 atom or more, then it has > 20 degree of freedom. Is this correction?

Please tell few atoms means here?

Q 2. Again on the same page, it mention getting values of POTIM from CG algo, as this algo choose dynamically POTIM.
""
The quantity trialstep is the size of the current trialstep. Multiply the current value of POTIM with this number to obtain the optimal step size.
""
Let's say, I perform 50 iterations using CG algo. Now, which trialstep, should I consider? Minimum, maximum or average one? Any suggestion will be appreciated.

Q 3. Finally, regarding ibrion=3, which is damped molecular dynamics, SMASS and POTIM both are required. Suggestion on POTIM are specific (step for ibrion 3 = 2* step from CG). Further in a talk slide () by Prof. Kresse, he has mentioned that SMASS = 2* sqrt (Gamma_min/Gamma_max), where Gamma_min, Gammma_max are minimum and maximum value of Hessian. Do they are same as eigenvalue of mixing* dielectric matrix in OUTCAR, which are printed in OUTCAR as below:

eigenvalues of (default mixing * dielectric matrix)
average eigenvalue GAMMA= 1.5392
8.7355 5.6152 5.6152 3.2220 2.7562 2.7562 2.1360 2.1360 2.0139 1.6124
1.6124 1.3593 1.3593 1.1860 1.1860 0.2855 0.2855 1.0180 1.0180 0.9508
0.9508 0.8691 0.8691 0.7214 0.7214 0.7770 0.7770 0.7674 0.7674 0.7036
0.7036 0.6124 0.6124 0.6165 0.6165 0.5905 0.5905 0.5385 0.3662

Can minimum and maximum of these eigenvalues can be considered some approximation of minimum and maximum of eigenvalues values of Hessian matrix?

Thanks,
Mukesh Singh

[ahampel]

Dear Mukesh,

thank you for reaching out to us on the VASP forum. Let me address your questions one by one.

Q1:

The number of ionic degrees of freedom (DOF) for a standard relaxation (no cell relaxation) should be something like 3N - N_(constraints) where N is the number of atoms.
Therefore, technically 7 atoms → 21 DOF and 8 atoms → 24 DOF and you are correct that already 7–8 atoms correspond to more than 20 degrees of freedom.

However, this assumes that all atoms can move independently. If symmetry is enabled (ISYM ≠ 0), the number of independent degrees of freedom may be significantly reduced. For adsorption systems such as Li₂Sₓ on doped 2D materials, it is generally advisable to set ISYM = 0. Otherwise, symmetry can artificially constrain the relaxation. Just in case that is not clear. So using this choice of ISYM 8 atoms would be larger than 20 DOF.

Regarding the choice of algorithm: for systems of this type (tens to >100 atoms) I would still recommend. IBRION = 2 (conjugate gradient) if your starting structures are reasonably close to the energetic minimum. This algorithm is robust and performs very reliably for adsorption problems. The benefit is that it does not require careful parameter tuning. I would use it unless you have convergence problems. Then switch to IBRION=3, then to IBRION=2 and 1 once you are closer to the optimal structure.

Q2:

As described on the VASP Wiki, the quantity trialstep printed during IBRION=2 runs gives information about the optimal step length relative to the current POTIM. You should not use the minimum value nor the maximum value. Instead I would recommend the following:
1. Ignore the first 5–10 ionic steps (far from minimum, strongly non-linear regime).
2. Look at the last 10 steps before convergence.
3. Take a representative average of trialstep in that region.

Then set: POTIM= POTIM_(old) x avg(trialstep). Only in the near-quadratic regime (close to convergence) does this estimate meaningfully reflect the local curvature. That said, for most adsorption systems, a POTIM in the range 0.3–0.5 already works well with IBRION=2 and does not necessarily require further optimization.

Q3:

There is an important conceptual difference. The formula

$$ SMASS = 2 \sqrt{\Gamma{\min}/\Gamma{\max}} $$

refers to the eigenvalues of the ionic Hessian, i.e.,

$$ \frac{\partial2 E}{\partial R_i \partial R_j} $$

The eigenvalues printed in the OUTCAR under: "eigenvalues of (default mixing * dielectric matrix)" are not the ionic Hessian eigenvalues. They are related to electronic density mixing and SCF convergence. These quantities describe electronic response and are unrelated to the curvature of the ionic potential energy surface as far as I know.

Therefore, the printed GAMMA eigenvalues cannot be used as an approximation for the Hessian eigenvalues.

Regarding SMASS = 0:

For SMASS ≥ 0, a Nosé thermostat is used.
SMASS = 0 simply selects an automatic Nosé mass corresponding to a thermostat period of about 40 time steps. It does not estimate the phonon spectrum or the Hessian. It is a fixed heuristic choice. I think I would stick with my above recommendation that if IBRION = 3 is used, a reasonable empirical starting value is SMASS = 0.4 and adjust if strong oscillations (increase SMASS) or very slow relaxation (decrease SMASS) are observed.

I hope this clarifies the conceptual points and helps you choose an efficient and stable relaxation strategy.

Best regards,
Alexander