Make MP correction work with ibrav=0 and cubic cell

Doc/release-notes

@@ -2,6 +2,7 @@ Fixed in development version:
   * XSPECTRA gave incorrect results with k-point parallelization, since
     at least v. 6.6, due to missing broadcast of recomputed Fermi energy
     (found and fixed by Fanchen Meng, Brookhaven)
+  * Makov-Payne correction wasn't working with ibrav=0 and cubic cell (A. Fonari)
 
 New in 7.0 version:
   * GPU support for PWscf and CP significantly extended

PW/src/input.f90

@@ -339,10 +339,11 @@ SUBROUTINE iosys()
   CHARACTER(LEN=256), EXTERNAL :: trimcheck
   CHARACTER(LEN=256):: dft_
   !
-  INTEGER  :: ia, nt, tempunit, i, j
+  INTEGER  :: ia, nt, tempunit, i, j, ibrav_mp
   LOGICAL  :: exst, parallelfs, domag, stop_on_error
   REAL(DP) :: at_dum(3,3), theta, phi, ecutwfc_pp, ecutrho_pp, V
   CHARACTER(len=256) :: tempfile
+  INTEGER, EXTERNAL :: at2ibrav
   !
   ! MAIN CONTROL VARIABLES, MD AND RELAX
   !
@@ -1417,6 +1418,11 @@ SUBROUTINE iosys()
   !
   IF (.NOT. tqmmm) CALL qmmm_config( mode=-1 )
   !
+  ! CRYSTAL LATTICE (MP correction depends on at)
+  !
+  call cell_base_init ( ibrav, celldm, a, b, c, cosab, cosac, cosbc, &
+                        trd_ht, rd_ht, cell_units )
+  !
   ! BOUNDARY CONDITIONS, ESM
   !
   do_makov_payne  = .false.
@@ -1429,7 +1435,10 @@ SUBROUTINE iosys()
     CASE( 'makov-payne', 'm-p', 'mp' )
       !
       do_makov_payne = .true.
-      IF ( ibrav < 1 .OR. ibrav > 3 ) CALL errore(' iosys', &
+      !
+      ibrav_mp = ibrav
+      IF ( ibrav .EQ. 0 ) ibrav_mp = at2ibrav(at(:, 1), at(:, 2), at(:, 3))
+      IF ( ibrav_mp < 1 .OR. ibrav_mp > 3 ) CALL errore(' iosys', &
               'Makov-Payne correction defined only for cubic lattices', 1)
       !
     CASE( 'martyna-tuckerman', 'm-t', 'mt' )
@@ -1492,11 +1501,6 @@ SUBROUTINE iosys()
   IF ( tefield ) ALLOCATE( forcefield( 3, nat_ ) )
   IF ( gate ) ALLOCATE( forcegate( 3, nat_ ) ) 
   !
-  ! CRYSTAL LATTICE
-  !
-  call cell_base_init ( ibrav, celldm, a, b, c, cosab, cosac, cosbc, &
-                        trd_ht, rd_ht, cell_units )
-  !
   ! ... once input variables have been stored, read optional plugin input files
   !
   CALL plugin_read_input("PW")

PW/src/makov_payne.f90

@@ -84,7 +84,8 @@ SUBROUTINE write_dipole( etot, x0, dipole_el, quadrupole_el, qq )
   REAL(DP) :: dipole_ion(3), quadrupole_ion(3), dipole(3), quadrupole(3)
   REAL(DP) :: zvia, zvtot
   REAL(DP) :: corr1, corr2, aa, bb
-  INTEGER  :: ia, ip
+  INTEGER  :: ia, ip, ibrav_mp
+  INTEGER, EXTERNAL :: at2ibrav
   !
   ! Note that the definition of the Madelung constant used here:
   !   Lento, Mozos, Nieminen, J. Phys.: Condens. Matter 14 (2002), 2637-2645
@@ -152,7 +153,9 @@ SUBROUTINE write_dipole( etot, x0, dipole_el, quadrupole_el, qq )
   WRITE( stdout, '( 5X,"    Total quadrupole moment",F20.8," a.u. (Ha)")' ) &
       SUM(quadrupole(:))
   !
-  IF ( ibrav < 1 .OR. ibrav > 3 ) THEN
+  ibrav_mp = ibrav
+  IF ( ibrav .EQ. 0 ) ibrav_mp = at2ibrav(at(:, 1), at(:, 2), at(:, 3))
+  IF ( ibrav_mp < 1 .OR. ibrav_mp > 3 ) THEN
      call errore(' write_dipole', &
                   'Makov-Payne correction defined only for cubic lattices', 1)
      !
@@ -161,7 +164,7 @@ SUBROUTINE write_dipole( etot, x0, dipole_el, quadrupole_el, qq )
   ! ... Makov-Payne correction, PRB 51, 4014 (1995)
   ! ... Note that Eq. 15 has the wrong sign for the quadrupole term
   !
-  corr1 = - madelung(ibrav) / alat * qq**2 / 2.0D0 * e2
+  corr1 = - madelung(ibrav_mp) / alat * qq**2 / 2.0D0 * e2
   !
   aa = SUM(quadrupole(:))
   bb = dipole(1)**2 + dipole(2)**2 + dipole(3)**2
@@ -173,7 +176,7 @@ SUBROUTINE write_dipole( etot, x0, dipole_el, quadrupole_el, qq )
   WRITE( stdout, '(/,5X,"*********    MAKOV-PAYNE CORRECTION    *********")' )
   WRITE( stdout, &
          '(/5X,"Makov-Payne correction with Madelung constant = ",F8.4)' ) &
-      madelung(ibrav)
+      madelung(ibrav_mp)
   !
   WRITE( stdout,'(/5X,"Makov-Payne correction ",F14.8," Ry = ",F6.3, &
        &              " eV (1st order, 1/a0)")'   ) -corr1, -corr1*rytoev