[Pw_forum] Multiple issues with graphene supercell calculations

[DHRUV SINGH]

Dear PWSCF users,

I use a hexagonal unit cell with 2 basis atoms and lattice parameter (C-C bond length) of 2.6868 a.u. - this value was arrived at by energy minimization over a range of lattice parameters near this value. With a 16x16x1 kpoint grid it is seen that the energy per atom converges to -11.395192 Ry and the forces on the 2 atoms are zero (atleast until a good precision). The ultrasoft pseudopotential with the RRKJ modification is used for C and the results obtained are in excellent agreement with Marzari and Mounet PRB 71, 205214, 2005. 
The essentials of the input file and the forces from the 2 atom unit cell are shown below
--------------------------------------------------------------------------------
 &system
     ibrav=0,
     celldm(1)=2.6868,
     nat=2,
     ntyp=1,
     ecutwfc=40.0,
     ecutrho=480.0,
     degauss=0.03,
     nbnd=8
 /
 &electrons
     mixing_beta = 0.5
     conv_thr =  1.0d-8
 /
ATOMIC_SPECIES
C 12.0 C.pbe-rrkjus.UPF

CELL_PARAMETERS {hexagonal}
 1.5 0.866025404 0.0
 1.5 -0.866025404 0.0
 0.0   0.0   20
ATOMIC_POSITIONS {crystal}
 C  0.0000000 0.0000000 0.0000000
 C  1/3.0 1/3.0 0.0000000

atom   1 type  1   force =    -0.00000000   -0.00000000    0.00000000
atom   2 type  1   force =     0.00000000    0.00000000    0.00000000
--------------------------------------------------------------------------------

Using this lattice parameter I then try to do supercell calculations with rectangular and hexagonal supercells and I am facing the following problems in the supercell scf calculations :

i) With a 4 atom rectangular supercell and the same lattice parameter of 2.6868 au, atomic coordinates and input shown below, I find that the energy/atom coverges to -11.395182 Ry(very close to the 2 atom unit cell) for a 8x8x1 k point grid but there seem to be non-negligible net forces on all 4 atoms as shown below, 
--------------------------------------------------------------------------------

&system
     ibrav=0,
     celldm(1)=2.6868,
     nat=4,
     ntyp=1,
     ecutwfc=40.0,
     ecutrho=480.0,
     degauss=0.03,
     nbnd=16
 /
 &electrons
     mixing_beta = 0.5
     conv_thr =  1.0d-8
 /

CELL_PARAMETERS {cubic}
 3.0 0.0 0.0
 0.0 1.732050808 0.0
 0.0 0.0 14
ATOMIC_POSITIONS {alat}
C       0       0       0
C       1       0       0
C       1.5     0.866025404     0
C       2.5     0.866025404     0
--------------------------------------------------------------------------------

 atom   1 type  1   force =     0.00310476    0.00000000    0.00000000
 atom   2 type  1   force =    -0.00310476    0.00000000    0.00000000
 atom   3 type  1   force =     0.00310476    0.00000000    0.00000000
 atom   4 type  1   force =    -0.00310476    0.00000000    0.00000000

ii) With a 48 atom rectangular supercell and 2x2x1 kpoint grid I find energy/atom of -11.395152 Ry (again very close to the 2 atom unit cell) but non negligible forces: Forces on atom 1 are shown below
atom   1 type  1   force =     0.00274639    0.00000000    0.00000000
The story is the same as I keep increasing the supercell size ... for 60 atoms the forces on atom 1 are
atom   1 type  1   force =    -0.00395744    0.00000000    0.00000000
The forces do not change even if I increase ecutwfc to 80.0 and ecutrho to 800.0. They also do not change if I increase the kpoint grid resolution

iii) Hexagonal supercells: I get slightly better answers with a hexagonal supercell in terms of forces. For example a 32 atom hexagonal supercell of graphene gives me an energy/atom of  11.395194 Ry for a 4x4x1 kpoint grid. The force on 1 of the atoms is 
atom   2 type  1   force =     0.00001419    0.00000887    0.00000000
and the magnitude is similar for most other atoms. Once again I find that increasing ecutrho or ecutwfc doesnt matter. 

However a more serious problem is seen with the hexagonal supercell. I find that scf calculations do not converge for a 72 atom supercell, with the following trend: 

estimated scf accuracy    <       7.42866516 Ry
     estimated scf accuracy    <       2.47702692 Ry
     estimated scf accuracy    <       0.05551173 Ry
     estimated scf accuracy    <       0.02263262 Ry
     estimated scf accuracy    <       0.01534456 Ry
     estimated scf accuracy    <       0.00799064 Ry
     estimated scf accuracy    <       0.00499060 Ry
     estimated scf accuracy    <       0.00272207 Ry
     estimated scf accuracy    <       0.00292201 Ry
     estimated scf accuracy    <       0.00255637 Ry
     estimated scf accuracy    <       0.00240574 Ry
     estimated scf accuracy    <       0.00239154 Ry
     estimated scf accuracy    <       0.00227859 Ry
     estimated scf accuracy    <       0.00214441 Ry
     estimated scf accuracy    <       0.00221359 Ry

I also tried changing the parameters ecutrho to 600.0; nband to 400 (from 72*4); decreasing degauss to 0.03 and trying mixing_mode='local-TF' as suggested but that does not help convergence one bit.

Please let me know if there is something else I could try which would help with these issues - decreasing the forces on each atom and helping convergence in the hexagonal supercell calculations.

Thank you very much
Dhruv Singh
Graduate Student, Purdue University ME