structural relaxation for supercell with impurity

[fhssn1]

I'm a little unsure about the best strategy for structural relaxation of a supercell with a big impurity.

Let's say we start with a primitive cell of a pure crystal, with empirically determined lattice parameters. We carry out k point convergence, ecut convergence, then optcell 0, then optcell 2, and we get a very low total energy.

Next, we create a supercell out of this. As long as the material the same, I assume there is no need for further convergence studies (we can use the ngkpt and ecut values from the primitive calculation).

Next, we introduce an impurity, an atom which is somewhat larger than the existing atoms in the system. At this point, ignoring the issue of k and ecut convergence for the moment, it's clear that the overall size of the supercell is bound to increase. And yet the abinit documentation for optcell "STRONGLY suggest"s to first relax ionic positions without cell volume and shape using optcell 0, and only after that use optcell 2 to relax the whole system.

My question is, is this really a good strategy in this case? If we know for sure the system will expand, aren't we causing needless stress in the system by forcing the cell volume to not expand using optcell 0?

I'm in the middle of a similar calculation and I'm monitoring SCF energies in the first BFGS iteration, using 'ETOT' lines in the log file. Those ETOT error values are oscillating back and forth between 1e-1, 1e0, 1e-2, and 1e-3. Am I doing it wrong? Comments and suggestions highly appreciated. Thanks in advance.

PS: If I have to ignore the recommendation in optcell documentation, I think I would go straight to optcell 2 and ionmov 2, and I would keep the toldff and tolmxf very high (like 1e-4, and 1e-3) for two reasons:

- Supercell energy magnitudes are larger than primitive cell energies, so giving the calculation same absolute error tolerance is overkill. (I'm not sure why relative error tolerances aren't available/the-preferred-approach for all tol variables. I've only seen tolrff, but couldn't find relative alternative to, e.g., toldfe.)

- When a supercell is far from a relaxed configuration (as in, when you introduce a large impurity), aiming for super low errors sounds like asking for SCF convergence failures. The system is probably not near the "convex" minimum, and so is likely to oscillate back and forth, without actually converging (probably what's happening in my case).

Once this technique (full relaxation with high tol values) works and we're able to get some kind of "crude" BFGS convergence, we can say that the system is probably closer to the relaxed configuration than it was when it started out. And now we can tighten toldff, and toldmxf, to get a fully relaxed configuration.

Again, kindly advise on the best strategy for this scenario. And, going back to the k and ecut convergence, what would be the best time to introduce such convergence calculation for the supercell-with-impurity: before or after structural relaxation? (if before, we again have the issue of the system far from a relaxed config, so maybe we should use a high value of toldfe?)

[ebousquet]

Dear fhssn1,
At this point I would say that you have to make your own experience while dealing with systems that are a bit away from classical ones (e.g. with defects/impurity ), but let me give some comments.
The advise to first relax the internal coordinates and then relax the cell is because in many systems doing so will help the algo to converge the relaxation by reducing the number of degrees of freedom. Starting the cell relaxation after will be less far from the equilibrium than from scratch, though you can start by relaxation everything together.
When you add impurities/defects in your system, you break some/many symmetries and it can also add new electronic states than can mess up both the electronic SCF and structure relaxation convergence.
Regarding the SCF, it often comes from the fact that the impurity/defect adds extra electronic states in the gap (if gaped system) or in the conduction states (if metallic system) that can be strongly "oscillating" or/and with many metastable states close in energy. You can play with diemix (reducing it), with diemac, or other preconditionning stuff (please look at https://www.abinit.org/sites/default/files/last/input_variables/generated_files/varset_gstate.html#iprcel) to help in converging.
For the relaxation, the strategy you proposed could be good, please use restartxf -1 to restart a structure relaxation by keeping the history of previous relaxations. I've used once tolrff but the result was not necessarily better, sometimes even worst than regular toldff...
In any case, often impurities/defects are messing up the convergence and you need more patience than bulk cases.

For the convergence criteria, as I said in another post, adding a defect could be like running a new crystal. For example, taking an extreme case, if you have an insulator and your defect/impurity makes your system metallic, the later will require much more k-points than the insulator reference pure case. So you can run by keeping the bulk convergence parameters and check if increasing k-points/ecut makes any change in the electronic structure and structural relaxation. Again, you have to make your own experience for your system.

Hope this can help,
Cheers,
Eric

[fhssn1]

Thanks Eric,

I have a quick question. Can I re-use WFK from previous run with lower number of k-points and lower ecut, lower pawecutdg, in a new run with higher number of k points and higher ecut, higher pawecutdg?

I just started a run and it didn't complain. But at the same time I'm wondering if it silently ignored the WFK or used it but generated wave-function components for the extra k points internally.

edit: I just realized the log file generated a ton of warnings, like k-points from input file will have to be translated, PAW ecut is different, rprim is different. Probably it had to ignore the WFK file and generate a new one?

[ebousquet]

Hi fhssn1,
I think if Abinit cannot re-use a WFK file it should complain and stops (if you used irdwfk flag), so this means he managed to re-start from the WF file you gave to him.
All the best,
Eric

[fhssn1]

The advise to first relax the internal coordinates and then relax the cell is because in many systems doing so will help the algo to converge the relaxation by reducing the number of degrees of freedom. Starting the cell relaxation after will be less far from the equilibrium than from scratch, though you can start by relaxation everything together.

I hope you wouldn't mind clarifying this point further. When we use optcell 0 (i.e., relax the internal coordinates), and we try to make it converge, using for example a value of tolmxf, there is a high likelihood that it will never converge! unless we pick a very high value of tolmxf, i.e., 1.0e-2, or something. I say that because if we're not relaxing the cell at the same time, we're keeping the system in a high stress state. It will reduce the stress due to re-positioning of atomic sites, I agree, but even in the least stressed state, max of forces could still be high enough that it'll never go below a given tolmxf (because the only way to release that stress is to let the cell itself expand or contract).

So when I do a optcell 0 run and the Broyd iterations don't converge, I should think about increasing the tolmxf value, (not just increasing ntime and doing more and more runs with restartxf). But if that's really the case, are we really doing convergence?

Please let me know if this is the correct line of thinking (I'm still learning and I might be off the mark here. I think I'm having problem with tolmxf as a convergence criteria for Broyd optimization. If it were 'difference of max forces' instead of 'max of forces' I think it would be a much better criteria. Do we have something like toldmxf?).

[ebousquet]

Please let me know if this is the correct line of thinking (I'm still learning and I might be off the mark here. I think I'm having problem with tolmxf as a convergence criteria for Broyd optimization. If it were 'difference of max forces' instead of 'max of forces' I think it would be a much better criteria. Do we have something like toldmxf?).

When you use optcell to 0, the stress does not enter into account in tolmxf, such as tolmxf will correspond to a "toldff" as if on the forces only.
If it does not converge, means your system is tricky to converge.
ionmov=22 could be worse to try too in case it is better for your case.
Best wishes,
Eric

[ebousquet]

PS: is your SCF well converged for each relaxation step?

[fhssn1]

PS: is your SCF well converged for each relaxation step?

Ok I guess my understanding of optcell 0 was incorrect (especially when the log file calls it 'max grad (forces/stress)' not just 'max (forces/stress)'. I guess if it's a grad (gradient) then it makes more sense to me.

As for your question, that's another issue I'm having, how to pick the best toldfe (or toldff) for scf convergence during Broyd iterations? In abinit help files, I think there is a mention that toldff is better than toldfe (for relaxation) and should be at least an order of magnitude smaller than tolmxf.

I'm doing a 32 atom GaN supercell with one impurity and if I pick tolmxf as 1e-4 and toldff as 1e-5 then it takes more than 24 hours worth of computation (on a 16 or 24 core setup). And for toldff of 1e-5 I keep track of what happens to toldfe (I guess I should call it dfe, difference of energies in 'ETOT' lines in the log file). It ends up going very very small, like 1e-10, or 1e-11 before toldff (dff) is able to go under 1e-5. But then toldmxf of 1e-4 is not that low.

If I pick toldmxf 1e-5 and toldff of 1e-6, I had to run it for multiple days, and with multiple restarts (with restartxf -2), and dfe had to be 1e-12 or 1e-13 or 1e-14 for convergence, and in the end it actually blew up! The whole experiment failed because the last scf run resulted in ETOT lines that kept increasing instead of decreasing: 1e+8, 1e+10, 1e+12, etc, etc.

I was told it may have something to do with a not well-tested PAW dataset of the impurity (because I have difficulties only when I add the impurity). My impurity is one of rare earth elements, and I recently came across VLab PAW datasets (and the paper by Topsakal) of rare earth elements, and I'm going to try that now.

[ebousquet]

As for your question, that's another issue I'm having, how to pick the best toldfe (or toldff) for scf convergence during Broyd iterations? In abinit help files, I think there is a mention that toldff is better than toldfe (for relaxation) and should be at least an order of magnitude smaller than tolmxf.

If you do a relaxation, it is better to use toldff for the SCF, because to stop the relaxation you need to reach the tolmxf precision, which is related to the forces. You can imagine converging with another tol like toldfe, but suppose that at this toldfe the residual on the forces is larger than tolmxf, this means you will never be able to reach tolmxf in your relaxation since the noise on the forces of your SCF is larger than tolmxf.
This is why it is recommended to use toldff when performing relaxation and to have a toldff smaller than tolmxf to have noise on the force that will be smaller than the stopping criterion tolmxf.

I'm doing a 32 atom GaN supercell with one impurity and if I pick tolmxf as 1e-4 and toldff as 1e-5 then it takes more than 24 hours worth of computation (on a 16 or 24 core setup). And for toldff of 1e-5 I keep track of what happens to toldfe (I guess I should call it dfe, difference of energies in 'ETOT' lines in the log file). It ends up going very very small, like 1e-10, or 1e-11 before toldff (dff) is able to go under 1e-5. But then toldmxf of 1e-4 is not that low.

If I pick toldmxf 1e-5 and toldff of 1e-6, I had to run it for multiple days, and with multiple restarts (with restartxf -2), and dfe had to be 1e-12 or 1e-13 or 1e-14 for convergence, and in the end it actually blew up! The whole experiment failed because the last scf run resulted in ETOT lines that kept increasing instead of decreasing: 1e+8, 1e+10, 1e+12, etc, etc.

This sounds to be more related to the fact that your system is troublesome rather than tolmxf/toldff problem.

I was told it may have something to do with a not well-tested PAW dataset of the impurity (because I have difficulties only when I add the impurity). My impurity is one of rare earth elements, and I recently came across VLab PAW datasets (and the paper by Topsakal) of rare earth elements, and I'm going to try that now.

Haha, I see, your impurity is a rare earth. Check if you really need to treat the f-electrons, if yes, at the present stage of the JTH table, I would recommand to use the VLab pseudos, they are better optimized. This might also explain your troubles in converging, thought not sure.
If you don't need the f-electrons (means you just want to have the steric effect of the rare earth atom), a pseudo with f-electrons treated as core electrons (means absent from the valence in the calculation) could do the job, but you really need to know what you are doing.

All the best,
Eric

[fhssn1]

I was told it may have something to do with a not well-tested PAW dataset of the impurity (because I have difficulties only when I add the impurity). My impurity is one of rare earth elements, and I recently came across VLab PAW datasets (and the paper by Topsakal) of rare earth elements, and I'm going to try that now.

Haha, I see, your impurity is a rare earth. Check if you really need to treat the f-electrons, if yes, at the present stage of the JTH table, I would recommand to use the VLab pseudos, they are better optimized. This might also explain your troubles in converging, thought not sure.
If you don't need the f-electrons (means you just want to have the steric effect of the rare earth atom), a pseudo with f-electrons treated as core electrons (means absent from the valence in the calculation) could do the job, but you really need to know what you are doing.

I would need to further look into whether I need f-electrons or not. However, I looked at both the JTH-PAW as well as VLab datasets for the RE (rare earth) elements and they are including a large number of valence electrons, more or less 30 from a single RE atom alone (3p6, 4s2, 4d10, the f electrons, plus at least one unoccupied band).

I don't remember coming across the term 'steric effects' before but I looked it up at Wikipedia, and it looks like I would be more interested in electronic effects, not just the steric effects. I am studying electronic properties of my system, and in a crystal structure form (my guess is steric-effect-only investigation is more common in chemistry/molecular side of DFT application, whereas I'm on the solid-state/crystal side of it). I could be wrong though.

I guess I have two points regarding the f electron inclusion/exclusion advice:

- The advice is a bit surprising in that it means we're moving the "outermost" electrons into core (4f) while those in the inner compared to it (3p, 4s, 4d) are kept as valence. But I guess looking at the bottom of JTH-PAW xml datasets, where the core and valence orbitals are listed as c, c, v, c, etc, it appears that sometimes lower energy orbitals are in valence, while higher in core (e.g., in Ga, it's listed as c, c, c, v, c, c, v, v. Meaning 3p6 is in core while 3s2 is in valence, which doesn't make sense to me: 3s is not only at a lower energy compared to 3p, it also comes earlier in the s-p-d-f ranking). So I guess my understanding about the theory part about how to decide what is in core vs valence is a bit lacking. I would appreciate if you can help me with that (or point out some reading resource).

- Even if I decide after careful examination that I can put f electrons into core, that still leaves me with the need to compute a new PAW dataset for the RE element(s), and me being completely inexperienced in PAW generation, would most likely generate something that would be way less stable/convergent than the JTH-PAW one, much less the VLab one.