GTH pseudopotentials have a special treatment, which unfortunately wa…

…s not

implemented in some parts of the code involving stress calculation and maybe
other cases as well. For the time being, I have added it to the cpu version
(but not in the GPU one, that stops with an error message). For the future:
I think that the more special cases we allow, the more we run into trouble.
In my opinion, GTH pseudopotentials should be treated on the same footing as
other pseudopotentials: computed on an interpolation table and interpolated.

Side remark: there is another special case, spline interpolation (as opposed to
the hand-made one). I think there should be just ONE type of interpolation.

upflib/gen_us_dy.f90

@@ -19,6 +19,7 @@ SUBROUTINE gen_us_dy_base &
   USE uspp,        ONLY: nkb, indv, nhtol, nhtolm
   USE uspp_data,   ONLY: nqx, tab, tab_d2y, dq, spline_ps
   USE uspp_param,  ONLY: upf, lmaxkb, nbetam, nh
+  USE m_gth,       ONLY : mk_ffnl_gth
   USE splinelib
   !
   IMPLICIT NONE
@@ -113,25 +114,29 @@ SUBROUTINE gen_us_dy_base &
   DO nt = 1, ntyp
      ! calculate beta in G-space using an interpolation table
      DO nb = 1, upf(nt)%nbeta
-        DO ig = 1, npw
-           IF ( spline_ps ) THEN
-             vkb0(ig,nb,nt) = splint( xdata, tab(:,nb,nt), &
-                                      tab_d2y(:,nb,nt), q(ig) )
-           ELSE
-             px = q(ig)/dq - INT(q(ig)/dq)
-             ux = 1.d0 - px
-             vx = 2.d0 - px
-             wx = 3.d0 - px
-             i0 = q(ig)/dq + 1
-             i1 = i0 + 1
-             i2 = i0 + 2
-             i3 = i0 + 3
-             vkb0(ig, nb, nt) = tab(i0, nb, nt) * ux * vx * wx / 6.d0 + &
-                                tab(i1, nb, nt) * px * vx * wx / 2.d0 - &
-                                tab(i2, nb, nt) * px * ux * wx / 2.d0 + &
-                                tab(i3, nb, nt) * px * ux * vx / 6.d0
-           ENDIF
-        ENDDO
+        IF ( upf(nt)%is_gth ) THEN
+           CALL mk_ffnl_gth( nt, nb, npw, omega, q, vkb0(1,nb,nt) )
+        ELSE
+           DO ig = 1, npw
+              IF ( spline_ps ) THEN
+                 vkb0(ig,nb,nt) = splint( xdata, tab(:,nb,nt), &
+                      tab_d2y(:,nb,nt), q(ig) )
+              ELSE
+                 px = q(ig)/dq - INT(q(ig)/dq)
+                 ux = 1.d0 - px
+                 vx = 2.d0 - px
+                 wx = 3.d0 - px
+                 i0 = q(ig)/dq + 1
+                 i1 = i0 + 1
+                 i2 = i0 + 2
+                 i3 = i0 + 3
+                 vkb0(ig, nb, nt) = tab(i0, nb, nt) * ux * vx * wx / 6.d0 + &
+                                    tab(i1, nb, nt) * px * vx * wx / 2.d0 - &
+                                    tab(i2, nb, nt) * px * ux * wx / 2.d0 + &
+                                    tab(i3, nb, nt) * px * ux * vx / 6.d0
+              ENDIF
+           ENDDO
+        ENDIF
      ENDDO
   ENDDO
   !

upflib/gen_us_dy_gpu.f90

@@ -20,7 +20,6 @@ SUBROUTINE gen_us_dy_gpu_ ( npw, npwx, igk_d, xk, nat, tau, ityp, ntyp, &
   USE upf_const,   ONLY: tpi
   USE uspp,        ONLY: nkb, indv_d, nhtol_d, nhtolm_d
   USE uspp_data,   ONLY: nqx, tab, tab_d2y, tab_d, dq, spline_ps
-  USE m_gth,       ONLY: mk_dffnl_gth, mk_dffnl_gth_gpu
   USE splinelib
   USE uspp_param,  ONLY: upf, lmaxkb, nbetam, nh, nhm
   USE device_fbuff_m,   ONLY: dev_buf
@@ -92,6 +91,8 @@ SUBROUTINE gen_us_dy_gpu_ ( npw, npwx, igk_d, xk, nat, tau, ityp, ntyp, &
   attributes(DEVICE) :: dvkb_d
 #endif
   !
+  IF ( ANY(upf(1:ntyp)%is_gth ) ) &
+       CALL upf_error( 'gen_us_dy_gpu',' GTH not implemented', 1)
   dvkb_d = (0._DP,0._DP)
   !
   IF (lmaxkb <= 0) RETURN

upflib/init_us_2_base_gpu.f90

@@ -1,5 +1,5 @@
 !
-! Copyright (C) 2001-2015 Quantum ESPRESSO group
+! Copyright (C) 2001-2021 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,