Problem with lattice optimisation

Elk Users
  • dtompsett

    dtompsett - 2014-04-01

    Dear All,

    I have been attempting to optimize the lattice parameter of cubic perovskite BaTiO3 using PBEsol. I am using the latest Elk 2.2.10. However, using the input file below I obtain a lattice parameter that is too small compared to both experiment and benchmarks with other codes (Wien2k and VASP).

    ! lattice vector optimisation while maintaining the unit cell volume

    ! speed the calculation up with Broyden mixing


    10 10 10

    7.576000000 0.000000000 0.000000000
    0.000000000 7.576000000 0.000000000
    0.000000000 0.000000000 7.576000000

    ! decrease the unit cell volume to increase the pressure


    3 : nspecies
    '' : spfname
    1 : natoms; atposl below
    0.00000000 0.00000000 0.00000000
    '' : spfname
    1 : natoms; atposl below
    0.50000000 0.50000000 0.50000000
    '' : spfname
    3 : natoms; atposl below
    0.50000000 0.50000000 0.00000000
    0.50000000 0.00000000 0.50000000
    0.00000000 0.50000000 0.50000000

    I have reduced the muffin tin radii to 2.0 bohr for both Ba and Ti, and for O to 1.6 bohr. The optimized lattice parameter is 7.403 bohr, while experiment and benchmarks on other codes are close to 7.57 bohr. The energy decreases nicely during optimisation as seen in TOTENERGY_OPT.OUT but the end lattice constant is not correct.

    I tried increasing rgkmax and gmaxvr, but that makes no difference. There is one warning during the SCF steps about "Warning(linengy): could not find 1 linearisation energies in s.c. loop".

    Any help much appreciated.


  • mfechner

    mfechner - 2014-04-04

    Dear david,

    I repeated you calculation and indeed you claims here are correct so with this method I obtained a=7.403 bohr indeed to small for PBEsol (my reference point is VASP pbesol a= 7.531 bohr). I add the following lines to check convergence.

    epsengy (I add this line already in the first run and you should to it too since the default is too small)






    using this setup I obtained a=7.452, so a little bit larger and going in the right direction. I then did equation of state fitting with this setup and further increase gmaxvr=18.

    so finally I end up with the following table

    code(method) a

    experimental 7.56
    vasp 7.53
    elk (via stress) 7.45
    elk (eos) 7.51
    elk (eos+gmaxvr18) 7.53

    you see the stress calculation give smaller lattice constant, whereas the the eos calculation gives reasonable results. I don't know the exact details of the stress calculation but I would guess that for BTO the stress converges slow with respect to gmaxvr/rgkmax. A similar behavior is also present in the stress calculation in vasp, but there with respect to the number of plane waves. I will test this but may take a while.

    In the mean time do you have some infos on the wien2k calculation did they do eos or stress relaxation ? Moreover what rgkmax/gmaxvr you use in wien2k ?


    two further remarks:

    -I increased the muffin tin radius of Ba to 2.2 bohr this will remove the linenergy warning. Actually, 2 is too small and result in a large core leakage (See INFO.OUT).

    -For the calculation with gmaxvr>17 you may have to add


    since otherwise the calculation will not converge due to the small Fourier components of Veff (See manual).

    Last edit: mfechner 2014-04-04
  • mfechner

    mfechner - 2014-04-07

    Okay I finished my convergence tests the result is that a further increase of RGKMAX/GMAXVR doesn't help. So I checked now the stress calculation itself. If you look into the manual you find that there is a parameter you can set for the stress calculation

    deltast -->size of the change in lattice vectors used for calculating the stress tensor

    I convinced myself that it is too large somehow, I in&decreased it for tests with the above mentioned settings and the decrease to 0.005 lead to the desired 7.51 lattice constant.


  • dtompsett

    dtompsett - 2014-04-07

    Hi Michael,

    Thank you very much for the help. I was considering looking at the same tweaks to deltast. Perhaps the default is not appropriate for this kind of system. The developers may wish to look at this since a cubic perovskite is one of easiest lattice optimisations possible. I suspect more testing of the stress calculation may have been done for metals. Metals are typically "softer" and this may allow a larger deltast in those cases.


  • John Kay Dewhurst

    Hi David and Michael,

    I've been running your BaTiO3 for a few weeks now and I can't get it to optimise in a stable way for PBESol. It works OK for PBE though.

    The reason is simple, but tricky to fix, namely that as the volume is varied the number of G-vectors and G+k-vectors change. This results in a small but discontinuous change in the total energy, which makes the finite-differences unstable. We'll get this fixed, although I'm not quite sure how yet. The idea would be to fix the number of vectors at the beginning and not change them during the optimisation.

    In the meantime, setting

     10 10 10

    ... gave pretty good results.


    Last edit: John Kay Dewhurst 2014-04-09
  • Ken Roberts

    Ken Roberts - 2014-04-09

    Dear Kay,
    Yes, the volume affecting the number of vectors is a concern. I encountered something similar with a manual exploration of as simple a structure as Au, using lattice constant as the parameter -- there is a curious sawtoothing effect on the graph of energy vs lattice as the volume change allows additional vectors to enter the computation.
    Ken R.

  • julian

    julian - 2016-04-07

    Hy everyone,

    concerning this and generally I have a question about the cutoff radii. It seems to me that the default cutoff radii are often too large. Can anyone tell me how they were determined?

    Then, if we take this example, the smallest Ti-O bond in the experimental structure should be 1.875 A. ELK default values (as far as I can extract them from the 'species' files) for Ti and O are 1.27 and 0.95 A, respectively, adding up to 2.22 A, thus, even exceeding the Ti-O bond length. Generally I would think the two radii should be chosen that they are like 10% smaller then this bond length, isn't it?

    Taking a O2 molecule as a second example, with a bond length of 1.208 A, makes the default O radius look even worth.

    Please let me know if I am missunderstanding/overlooking something here. Or am I right and the default cutoff radii are just not good and have to be adjusted for every system?


    Last edit: julian 2016-04-07
  • mfechner

    mfechner - 2016-04-09

    Dear Full Name,

    I guess with cutoff radii you mean the muffin tin radius of the atoms ? Indeed you are right for structural relaxation on should choose this radii to left some freedom for the atoms to move. This is also true for volume or other relaxations. In other situations however, where this degree of freedom is fixed its sometimes favorable to have bigger muffin tins, because they give you a better description of the system. So please note that for example spin orbit coupling and DFT+U potential is only applied within the muffin tin spheres.
    Finally, you are right one has to adjust this radii for each system regarding what you want to do.

    best regards

    PS: If you have been working before with pseudopotential codes this problem is indeed new for you since there the radii are fixed within the pseudopotentials itself.

  • julian

    julian - 2016-04-09

    Dear Michael,

    indeed, so far I was happy with pseudopotential codes and I am very new to ELK and its framework. Your remark about spin-orbit coupling and +U only acting within the MT sphere is indeed quite interesting (and a bit concerning). Since both can be implemented with plane-waves, may I ask what the reason for this is? Is it 'just not done yet' or is there specific reasoning behind it (I guess one could argue that both effects are less pronounced for the outer most valence bands/orbitals, but if there is no specific reason against it (probably the connection of MT and interstitial?) why not include it)?

    Considering this and the other aspects of your reply, I am not sure what to do with the MT spheres exactly. So I do understand now that the MT spheres are chosen as large as possible. Thus, the defaults are probably not applicable to all systems. Okay, so let's adjust them. But is there a good recipe you can share from your experience so far? For simplicty, let's stay with the O2 example. So there we need to decrease the MT sphere. Since I want to optimize it's structure, we should give it some freedom. I think you agreed so far and say: "for structural relaxation on should choose this radii to left some freedom". Is the value I gave (estimated bond length/2 -10%) resonable or does your experience show that this leads to too small MT spheres (i.e., can you specify 'some' further)?

    Next, after decreasing the MT sphere, I assume the basis and probably the separation into core/semi-core might also need readjustments (in the chosen example the oxygen MT radius would be required to be decreased by almost 1/2, from 0.95 A to ~0.54 A in order to prevent overlap and give an additional safety of 10%)? This makes it a bit tedious, but anyways, what is you experience here? What do we need to test look out for? (Probably the example is quite extreme, since one will mostly deal with solid systems, but in any case it would be nice to solve since being able to treat referece cases on the same footing as, e.g., metal oxides would be a necessity to have).

    Last, on the one hand, you are saying the MT spheres should be as large as possible, on the other hand, results should be independen of the chosen MT sphere (this is what is the advised test in this methotology as far as I understand). The first statement would suggest to chose a smaller MT sphere first for a optimization, rerunning the calculation with a maximal MT radius (up to touching MT spheres) on a fixed geometry for better results, once the geometry is optimized. The latter statement, however, contradicts this and suggests that results are only good if they are independent from a (small) change of the MT sphere and I really hope that this is the case, since the scheme becomes quite unreliable otherwhise. So again, is there a recipe or experience you can share?

    Overall, I am sorry if those questions are naive for people working in this field for long, but maybe you can help me out, turning me into a new fan of this methodology.



  • mfechner

    mfechner - 2016-04-10


    its seem I generated more confusion than clarfication. Okay first of all a good converged setting regardind the basis in LAPW should not lead to a variation of results if the MT is changed! So for your structural relaxation I would recommend you to go with the 10% reduction of the MT radius an go on with it. Should you encounter with your setting a situation where you have to decrease the MT again ( because the atoms prefere a even smaller bond length) you will have to lower the MT again. Always keep an eye on the ""core leakage that it does not become to large (>0.1) per atom" otherwise you will have to change core-valence configuration in the species file. By the way, if you encounter ''strange'' results or you want to verify your convergence with respect to basis a good check is always to modulate the MT and redo the calculation. Moreover, I prefer to use isgkmax=-2 so the plane wasve basis becomes independent on the MT.

    Regarding your question concerning the reduction of the MT of oxygen by 50 %, indeed this would requiere a readjustment of semicore and core states. In principle you could move states step by step from the core to the valence and rerun the calculation and check always the "core leakage" in INFO.OUT. If this becomes small (0.01) per each atom you are fine and can continue with this setting.

    I hope this clarified some points.


    PS: Regardin the ls coupling and DFT+U in the pseudopotential codes both are also only applied within the pseudopotential radius! Hence there we are facing the same situation as in LAPW, but actually even worse since ls coupling strengths have been generated with the pseudopotential and may not fit your situation and or depends on how the potential has been generated.


Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks