LDA pseudopotential for Carbon - good for molecules, bad for graphene?

Pseudopotentials and Basis Sets
1

Hello,

I have downloaded a LDA pseudopotential for Carbon in this link:
https://departments.icmab.es/leem/SIESTA_MATERIAL/Databases/Pseudopotentials/Pseudos_LDA_Abinit/C_html/C.html

I ran tests on molecules - pretty large ones - and it did well. I tried monolayer graphene and it also gave me good structural parameters. Now, for graphite or bilayer grahene, the energy keeps lowereing the smaller the interlayer distance (which from other codes I got 3.36 A as the optimal value.

My fdf file is:
SystemName bicamada
SystemLabel bicamada
NumberOfAtoms 4
NumberOfSpecies 1

%block ChemicalSpeciesLabel
1 6 C_LDA # Species index, atomic number, species label
%endblock ChemicalSpeciesLabel

LatticeConstant 1.0 Ang

%block LatticeVectors
1.22384691 -2.11932969 0.00000000
2.44723916 5.62061381E-04 0.00000000
0.00000000 0.00000000 30.0000000
%endblock LatticeVectors

K-sampling (alternative specification using kgrid_cutoff)

%block kgrid_Monkhorst_Pack
15 0 0 0.0
0 15 0 0.0
0 0 1 0.0
%endblock kgrid_Monkhorst_Pack

#Atomic coordinates
AtomicCoordinatesFormat Ang

%block AtomicCoordinatesAndAtomicSpecies
0.00000000 0.00000000 0.00000000 1
2.44739113 -1.41257062 0.00000000 1
2.44912297 -1.41157025 $a 1
1.22715986 -0.70437287 $a 1
%endblock AtomicCoordinatesAndAtomicSpecies

Basis set definition

PAO.EnergyShift 500 meV
PAO.SplitNorm 0.2
PAO.BasisSize DZP

#Real space grid
MeshCutoff 200.0 Ry

Convergence of SCF

MaxSCFIterations 0
DM.MixingWeight 0.4
DM.NumberPulay 4

Type of solution (diagon is the default for less than 100 atoms)

SolutionMethod diagon

#Geometrical optimization
MD.TypeOfRun CG
MD.NumCGsteps 0
#MD.MaxCGDispl 0.1 Bohr
#MD.MaxForceTol 0.04d0 eV/Ang
#continuation job
MD.UseSaveCG true
MD.UseSaveXV true
DM.UseSaveDM true
LongOutput true

The “$a” is because the script is changing these values to run several simulations in the search of the optimal structure. My energy curve is like that:
$a Etot (eV)
1.40 -631.562255
1.45 -629.984906
1.50 -628.463349
1.55 -627.018617
1.60 -625.663935
1.65 -624.406879
1.70 -623.251031
1.75 -622.196964
1.80 -621.243233
1.85 -620.386733
1.90 -619.623120
1.95 -618.946846
2.00 -618.351451
2.05 -617.830117
2.10 -617.375531
2.15 -616.980744
2.20 -616.638906
2.25 -616.343980
2.30 -616.090300
2.35 -615.872742
2.40 -615.686828
2.45 -615.528253
2.50 -615.393439
2.55 -615.278859
2.60 -615.181737
2.65 -615.099480
2.70 -615.029843
2.75 -614.971070
2.80 -614.921433
2.85 -614.879697
2.90 -614.844507
2.95 -614.814980
3.00 -614.790273
3.05 -614.769559
3.10 -614.752380
3.15 -614.738040
3.18 -614.730667
3.22 -614.722062
3.26 -614.714667
3.30 -614.708353
3.34 -614.702980
3.38 -614.698405
3.42 -614.694544
3.46 -614.691283
3.50 -614.688525
3.54 -614.686244
3.58 -614.684327
3.62 -614.682712
3.66 -614.681417
3.70 -614.680347

Have I made some mistake in my input or is it safe to assume that the pseudopotential is not good? In that last case, does anyone have a good LDA pseudopotential to share?

Thanks!
Cedric

2

Hi,

For graphite, you should use van der Waals functional (XC.Functional VDW), please check the manual.
You might want to try to use GGA-PBE pseudopotential + van der Waals functional.
You have also different implementations of van der Waals functional (XC.Authors), please check the manual.

Remember that you can optimized the geometry in SIESTA using “MD.TypeOfRun CG” and you can also set atom constraints (GeometryConstraints in SIESTA version 4.0, Geometry.Constraints in SIESTA version 4.1).

3

Thank you for the reply. I was checking the accuracy of the LDA pseudo that I downloaded, as one should always do. LDA is known for describing reasonably well the structural properties of systems with van der Waals interactions, even if that is usually attributed to error cancelation. GGA-PBE, on the other hand, is known for predicting interlayer distances grossly overestimated for VdW structures, which is attributed to the better description of correlation relative to LDA, whereas the exchange term remains equal, thus losing the effect of error cancelation in LDA. In any case, I would never expect a pseudopotential predicting such short vdW bonds as I obtained. Therefore I assumed the pseudo should be wrong. In that case, maybe it should be removed from the repository where I got it from. Otherwise there is something else that I have been overlooking so far. If anyone has comments on this it will be mostly appreciated. Best regards Cedric

4

I notice that you are using

The default value is 0.02 Ry (~272 meV) but you usually want to lower it.
Try to use a value below 50 meV, maybe 20 meV or so.

High values of energyshift generate a less extended basis set and can cause wrong results.
In order to set a proper energyshift one has to start from a reasonable value and modify it if necessary depending on the system and on the property of interest.

5

Thank you for point out the high Energy shift and the consequently short range of the basis set. Although this is definitely a problem, the reason for the abnormal behavior in my energy optimization for bilayer graphene was not allowing the code to perform the self consitent optimization of the electron density. The flag MaxSCFIterations was set to 0. I have obtained now an interlayer equilibriu distance of about 3.2 Ang for bilayer graphene using the pseudopotential from Siesta’s website. Good news is that the pseudo is good! Thanks again and if you think this post should be removed please go ahead and accept my apologies. Best regards. Cedric

6

You are right, I am glad you found and solved the problem.

For the future, should you have any other issues, you can also attach the fdf file and psf files to the post so it is easy to download the files and try to reproduce the issue.