[SOLVED] Is Coulomb truncation applicable in a 2D system?

[Robin]

Dear friends,

While applying the Coulomb truncation to a hexagonal surface in a GW calculation, Abinit 6.2.2 outputs the following error, which I can’t fix. Since this case (Coulomb truncation for non-orthorhombic surface) is not mentioned in the help file nor included in the test sets, I am wondering if this feature is available now.

Thanks for your help.

Sincerely,
Guangfu Luo

--------------------------------------- Relevant parameters------------------------------
xred 1.0271626370E-33 1.5124252423E-01 0.0000000000
-7.7037197775E-34 2.9871361846E-01 0.0000000000
0.0000000000E+00 4.2921413811E-01 0.0000000000
… …
angdeg 90 90 120
icutcoul 2
vcutgeo 1 1 0
rcut 10

---------------------------------------------Error message---------------------------------
optimal value for ng0sh = 1 1 0
m_coulombian.F90:383:ERROR
Surface must be in the x-y plane
-P-0000
-P-0000 leave_new : decision taken to exit ...

[Robin]

Dear friends,

According to the error message, I checked the m_coulombian.F90 file, which sits under abinit-6.2.2/src/67_common. The error occurs at line 383 because in this 2D system, the value of Vcp%pdir is not "1 1 0", though I do set vcutgeo to be "1 1 0" in the *.in file. When tracing back, one can find that the latest Vcp%pdir value comes from some codes between line 373 and 379, and they change the Vcp%pdir to be "1 1 1". This is further caused by the if-condition at line 375--the condition is always true because the value of test, which counts the nonzero number in vcutgeo, is always positive.

Analyses suggest that this if-condition may be "ABS(check)>zero", and consequently the other if-condition would be "check<zero".

After changing the two places and recompiling abinit, the calculation of Coulomb cut-off in 2D system finally works. But this surely needs further verification of the developers. (p.s. This has been comfirmed by Matthieu.)

Sincerely,
Guangfu Luo

------------------possible bug in m_coulombian.F90 (line 373-379)---------------
do ii=1,3
check=Vcp%vcutgeo(ii)
if (ABS(test)>zero) then ! Use Rozzi"s method with a finite surface along x-y
Vcp%pdir(ii)=1
if (test<zero) Vcp%alpha(ii)=normv(check*rprimd(:,ii),rmet,'R')
end if
end do

------------------possible revision of m_coulombian.F90 (line 373-379)----------
do ii=1,3
check=Vcp%vcutgeo(ii)
if (ABS(check)>zero) then
Vcp%pdir(ii)=1
if (check<zero) Vcp%alpha(ii)=normv(check*rprimd(:,ii),rmet,'R')
end if
end do