Unintentional symmetry breaking in unit cell convergece

Hello everyone, I’m a newcomer to Abinit and just getting my feet wet with the code, though have a fair amount of experience with QE, VASP, and DFT in general. After converging the geometry of a few simple unit cells, I tried to test my understanding by converging a small yet “messy” unit cell with many degrees of freedom - monoclinic C2/m CuBr2. So far, whenever I allow unit cell geometry to change, there’s an error with the symmetry and I end up in space group P1. I am using the manual-suggested method of a two-step run that first optimizes internal degrees of freedom before allowing the cell to change. I suspect things might be due to some numerical problem with the choice of primitive unit cell vectors - I know that in other codes the software can be picky about exactly which ones are used. What’s strange is that I’m using the vectors suggested by running abistruct convert on a cif file, so I assumed that choice would be a good one.

I have tried:

  1. Increasing tolsym to 1d-4, but nothing much seemed to change
  2. Using a larger two formula unit cell with “nicer” lattice vectors and chkprim=0, but it seemed like Abinit wasn’t able to properly identify the space group
  3. When I ran the non-primitive cell without setting chkprim=0, it suggested a primitive unit cell with different lattice vectors than the ones given by abistruct. I tried those, and symmetry still broke in much the same way
  4. Using angdeg instead of setting rprim, but I got an unexpected error, “Number of bands nband= 18 > number of planewaves npw= 1”, which didn’t make sense to me since I had not changed ecut.

Things remain stable if I use optcell=1 and only allow volume changes to the cell, but that’s not a very satisfying solution. It’s also possible that the system is delicate and I might just need a higher ecut or kpoint mesh. I’ve included my .abi file below, which uses LDA and pseudopotentials from Pseudodojo (v4.1).

Thank you for any advice you are able to give! If you have any tips on how to generate the primitive cell from the .cif than abistruct (or what might be going wrong with angdeg) I am happy to hear them.

ndtset 2
natom 3
ntypat 2
typat 1 2 2
znucl 29 35
0.0000000000 0.0000000000 0.0000000000
0.5056360000 0.5056360000 0.7612900000
0.4943640000 0.4943640000 0.2387100000
acell 1.0 1.0 1.0
7.8496356240 0.0000000000 -3.3648092894
5.0741829977 5.9891106465 -3.3648092894
0.0000000000 0.0000000000 13.0609347069

pp_dirpath “.”
pseudos “Cu.psp8, Br.psp8”
ixc = 1

#Definition of the planewave basis set
ecut 20

#Definition of the k-point grid
kptopt 1
ngkpt 6 6 6
#nshiftk 1 # these are the defaults so we don’t need to include them
#shiftk 0.5 0.5 0.5

optcell1 0
optcell2 2
getxred2 -1
ionmov 2
ntime 10
dilatmx 1.05
ecutsm 0.5
tolmxf 1.0d-4
nstep 50

#Definition of the SCF procedure
toldfe 1.0d-6 #
diemac 12.0

Hi, I’m running the system you mentioned, with the cif file of Oeckler and Simon, and I am not seeing the problems you mention. However, I do have a number of comments on your input file that might be causing the issue you mentioned and others as well.

Most importantly, your ngkpt choice is almost certainly not consistent with the cell symmetry. This is definitely a case where you should be allowing abinit to do the work for you. Do a preparatory run with prtkpt 1 and kptrlen 70, and look at the various kpt grids proposed. They will be consistent with your cell symmetry, although ngkpt probably won’t be the input style you’ll end of using. prtkpt 1 gives choices using the kptrlatt formalism, which is the most general. When I did this I ended up using kptrlatt 5 -10 0 10 0 0 -5 5 5, which has kptrlen 65.6, which from experience is usually sufficient for ok stress convergence. You might have to go higher, for reasons discussed below. Notice that this kptrlatt is very unlike an ngkpt input, but that’s because the cell symmetry is far from cubic.

Next, you are converging with toldfe, that is, the total energy. But because you want to converge stresses (and forces) you should be converging stresses and forces, which is much more demanding than energies. Typically you’d use toldff 2.0D-6 and tolmxf 2.0D-5. This strategy fails, by the way, when the cell is so symmetric that no atoms ever have forces on them, such as diamond, then you have to use tolvrs and tolmxf.

Next, your system is metallic, it has an odd number of valence electrons. occopt 1, which is the default, tries to fill the bands 2-by-2, in this case it’s better to allow for metallic occupancy, by using for example occopt 7. The metallic occupancies also have a smearing variable which is very important for convergence, so don’t assume for “real” work that the default tsmear value is the best one. It should be checked.

In any case, with the various changes outlined above, the system converges ok and the “getxred -1” strategy works fine. Hope this helps–

Thanks, this ended up working. I’m still getting a handle on what the default behavior of the code is, but with respect to kpoints, I had expected the system to throw an error and stop if I used an incompatible grid, since I hadn’t turned off chksymbreak - am I misunderstanding the point of that flag?

And thank you for pointing out the mistake with my convergence criteria - it looks like I had confused the flag for force convergence tolerance within a single scf cycle and convergence within the structural relaxation loops. Also good to know that the default occupations are fixed - in VASP there is some default smearing, and I think I got used to it. I stress that I would converge these quantities for a publication - this was just a first stab at things for me.

No, I think you are understanding chksymbreak correctly, I’m not quite sure exactly why your run crashed, I did not get that behavior (but I didn’t run your kmesh either). As for the defaults, this is always a design problem. Many beginners (not you obviously, I mean the beginners who also don’t want to bother to learn about DFT) want everything to “just work” and so they want defaults for everything, which is impossible of course. Real experts usually want very few defaults. And, by setting certain defaults in certain ways, you can make your code look like it has superior performance, although it may actually not be well converged at all (we try not to do that in abinit).