[QE-developers] Coulomb Cutoff 2D

[Andreussi, Oliviero]

Dear Developers,

I am working on the code to make the calculation of the Hartree potential a bit more clean and code independent (PW vs CP). It may not be very useful per se, but it may help to reduce the cyclic dependencies of Environ.

Having said this, while looking at the current version of the code I realized that there is a new feature related to the electrostatic interactions and the pbc corrections, the Coulomb Cutoff 2D. I double checked the formulas, and it seems to me that they are the same of the Martyna-Tuckermann approach in 2D (Mináry, P.; Tuckerman, M. E.; Pihakari, K. a.; Martyna, G. J. A New Reciprocal Space Based Treatment of Long Range Interactions on Surfaces. J. Chem. Phys. 2002, 116 (13), 5351.), when an appropriate smoothing is added to avoid the cusp at the boundary, similar to what is done in the 0D Martyna-Tuckermann implementation in PW (the additional factor exp(-tpiba2*gg(:)*beta/4._dp)**2 in wg_corr )

Apart from the stress subroutines, which are totally new and very useful, the other subroutines in the Coul_cut_2D.f90 module can be merged inside the martyna-tuckermann.f90 module with very little changes. As a matter of fact, I did it already a few years ago, but did not submit it because I struggled -and failed- to make the 1D case work. My impression is that one just needs to change the way the wg_corr(ig) coefficients are computed, which for the 2D case would be something like

     axis = at(3,3) * alat
     DO ig =1, ngm
       q2 = gg(ig)*tpiba2
       IF (q2>1.e-6_dp) THEN
         gxy = SQRT(g(1,ig)**2 + g(2,ig)**2)
         wg_corr(ig) = - COS( 0.5D0 * g(3,ig) * tpiba * axis ) * &
                     & EXP( -0.5D0 * gxy * tpiba * axis ) * &
                     & fpi / q2
       ELSE
         wg_corr(ig) = 0.D0
       END IF
    END DO
    wg_corr(:) =  wg_corr(:) * exp(-tpiba2*gg(:)*beta/4._dp)**2
    if (gamma_only) wg_corr(gstart:ngm) = 2.d0 * wg_corr(gstart:ngm)

I think having the two correction schemes with a more uniform implementation would be desirable. In particular, I find the implementation of the MT approach to be very clean and compact, I would actually suggest to move it to Modules, the only dependencies that forbid this are in a subroutines to print out debug stuff, which could be easily commented out.

If the developers of this new feature want to fix this by themselves, I would not mind at all. Otherwise, if there are strong opinions on why the Coulomb Cutoff 2D should be implemented in a totally different way from the MT approach, despite being both G-space-based corrections, please let me know. At least we should change the capitalization of the module’s name, it is really hArd tO looK At.

Best,

Oliviero Andreussi
--
Assistant Professor
Department of Physics
University of North Texas
Email: oliviero.andreussi at unt.edu<mailto:oliviero.andreussi at unt.edu>
Skype: olivieroandreussi
Web: https://sites.google.com/site/olivieroandreussi

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/developers/attachments/20180510/1a7ce477/attachment.html>