Problems with cell optimization of graphene

[Riccardo Pezzetta]

Dear CP2K community,

I am trying to optimize the lattice constant of a graphene cell with RPA at the Gamma point. I performed energy calculations with different lattice constants, but I get a very strange behavior close to the minimum. As shown in the image attached to this mail, the RPA energy shows an unphysical drop close to the minimum. With PBE the curve is regular, while the HF@PBE energy (i.e., RPA without the correlation term) shows the same irregular profile as RPA, so the issue seems to originate from the Hartree-Fock exchange term. In my input, that you find attached to the mail, the only parameter I set in the HF section is the cutoff of the truncated Coulomb operator, I set it equal to 4 Angstrom so that it does not exceed half the cell size for all values of the lattice constant considered. Have you ever encountered a similar problem? What can be the cause of this unphysical behavior? Thank you in advance for your answer.
Kind regards,

Riccardo Pezzetta

[Augustin Bussy]

Dear Riccardo,

your HFX input is not ideal, because a 4 Angstrom cutoff is probably too short. There are 2 things to consider when selecting a value for the cutoff:

1) is should not exceed half the cell size (OK in your case)

2) it should be long enough. In most cases, a value of ~6 Ang is deemed sufficient, but technically, the ideal cutoff depends on the system band gap (for smaller gap systems, a longer cutoff might be required)

I suggest you take a supercell, and extend the cutoff.

Note that the value of EPS_PGF_ORB you use is way too small (1.0E-16), and it will dramatically impact performance. Please read section 2.2 of https://pubs.acs.org/jpcbfk/article/130/4/1237/5088250/The-CP2K-Program-Package-Made-Simple for an explanation.

If using a supercell makes the RPA part of your calculation too expensive (RPA scales as N^4 with system size), consider using the low-scaling implementation; it is particularly efficient for 2D systems.

Best,
Augustin

[Riccardo Pezzetta]

Thank you very much for the answer.

Kind regards,

Riccardo Pezzetta 

[Frederick Stein]

Dear Riccardo,

You are using smearing. Can you turn it off? Are all orbitals fully occupied? How does the band-gap change with cell volume? Beware that our RPA-implementation does NOT consider partial occupation. Even worse, partially occupied orbitals are considered as fully occupied by our implementation such that the energies may even differ with the available orbitals.

In addition, I am not entirely sure how well the 2D-periodicity is covered in your executable of CP2K and I would rather use PERIODIC XYZ instead.

Best,

Frederick

[Riccardo Pezzetta]

Dear Professor Stein,

Thank you very much for your reply.

I checked the convergence with respect to the smearing temperature. I found that the energy is very similar for electronic temperatures of 100, 200, and 300 K, while a noticeable change starts to appear at 500 K. However, if I completely remove the smearing, the PBE SCF calculation does not converge.

I also checked the convergence with respect to the size of the cell along the z direction.

Following your suggestion, I added the following section to print the orbital occupations:

&PRINT
&MO
OCCUPATION_NUMBERS
&END MO
&END PRINT

I found that the occupation numbers are either 2 or 0 for the different lattice constants considered, with no partially occupied orbitals.

Finally, I tried to perform a pure HF calculation to investigate whether the problem originates entirely from the Hartree-Fock exchange term. However, the HF calculation does not converge, even when using a relatively high electronic temperature or when mixing HF exchange with PBE exchange or when using PERIODIC XYZ instead of PERIODIC XY.

Do you have any suggestions on what I could try next?

Thank you again for your help.

Best regards,

Riccardo Pezzetta

[Frederick Stein]

Dear Riccardo,

Usually, I restart the calculation from a preconverged one, usually PBE. and point CP2K to the reference wfn-file using WFN_RESTART_FILE_NAME <name of the restart file>. If that one does not converge, I perform a PBE0-calculation starting from PBE and restart from PBE0.

A few more notes on your input file:

- EPS_PGF_ORB is too tight (just comment it out to use the default)

- make sure to converge EPS_DEFAULT

- make sure to use a suitable pseudopotential (PBE or GGA for PBE, PBE0 or HYB for PBE0 and HF for HF)

- Why are you not applying an exchange correction in ADMM? Just use the default (omit EXCH_CORRECTION_FUNC)

- As already mentioned by Augustin, your cutoff for HF might be too low. this introduces artifacts. Can you extend your system in the x- and y-directions? You need at least 5-6 A for converged results.

- start with a smaller ADMM basis set (AUX_FIT) of dz or dzp-quality. The one you are using will not introduce a siginificant acceleration.

- the RI-basis set (RI_AUX) was not optimized for your primary basis set and also not tested with it. You may remove them and set AUTO_BASIS RI_AUX [SMALL|MEDIUM|LARGE|HUGE] (in the DFT-section) instead.

- Are you sure that you need such a high cutoff in the WF_CORRELATION section ? I commonly apply 200 Ry. I have never seen value larger than 300 Ry in this method.

Best,

Frederick