1! 2! Copyright (C) 2001-2007 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!---------------------------------------------------------------------- 10subroutine dvloc_of_g (mesh, msh, rab, r, vloc_at, zp, tpiba2, ngl, gl, & 11 omega, dvloc) 12 !---------------------------------------------------------------------- 13 ! 14 ! dvloc = D Vloc (g^2) / D g^2 = (1/2g) * D Vloc(g) / D g 15 ! 16 USE kinds 17 USE constants , ONLY : pi, fpi, e2, eps8 18 USE Coul_cut_2D, ONLY: do_cutoff_2D 19 USE esm, ONLY : do_comp_esm, esm_bc 20 implicit none 21 ! 22 ! first the dummy variables 23 ! 24 integer, intent(in) :: ngl, mesh, msh 25 ! the number of shell of G vectors 26 ! max number of mesh points 27 ! number of mesh points for radial integration 28 29 real(DP), intent(in) :: zp, rab (mesh), r (mesh), vloc_at (mesh), & 30 tpiba2, omega, gl (ngl) 31 ! valence pseudocharge 32 ! the derivative of the radial grid 33 ! the radial grid 34 ! the pseudo on the radial grid 35 ! 2 pi / alat 36 ! the volume of the unit cell 37 ! the moduli of g vectors for each s 38 ! 39 real(DP), intent(out) :: dvloc (ngl) 40 ! the fourier transform dVloc/dG 41 ! 42 real(DP) :: vlcp, g2a, gx 43 real(DP), allocatable :: aux (:), aux1 (:) 44 real(DP), external :: qe_erf 45 46 integer :: i, igl, igl0 47 ! counter on erf functions or gaussians 48 ! counter on g shells vectors 49 ! first shell with g != 0 50 51 ! the G=0 component is not computed 52 if (gl (1) < eps8) then 53 dvloc (1) = 0.0d0 54 igl0 = 2 55 else 56 igl0 = 1 57 endif 58 59 ! Pseudopotentials in numerical form (Vloc contains the local part) 60 ! In order to perform the Fourier transform, a term erf(r)/r is 61 ! subtracted in real space and added again in G space 62 63 allocate (aux1( mesh)) 64 ! 65 ! This is the part of the integrand function 66 ! indipendent of |G| in real space 67 ! 68 do i = 1, msh 69 aux1 (i) = r (i) * vloc_at (i) + zp * e2 * qe_erf (r (i) ) 70 enddo 71 ! 72!$omp parallel private(aux, gx, vlcp, g2a) 73 ! 74 allocate (aux( mesh)) 75 ! 76!$omp do 77 do igl = igl0, ngl 78 gx = sqrt (gl (igl) * tpiba2) 79 ! 80 ! and here we perform the integral, after multiplying for the |G| 81 ! dependent part 82 ! 83 ! DV(g)/Dg = Integral of r (Dj_0(gr)/Dg) V(r) dr 84 do i = 1, msh 85 aux (i) = aux1 (i) * (r (i) * cos (gx * r (i) ) / gx - sin (gx & 86 * r (i) ) / gx**2) 87 enddo 88 call simpson (msh, aux, rab, vlcp) 89 ! DV(g^2)/Dg^2 = (DV(g)/Dg)/2g 90 vlcp = fpi / omega / 2.0d0 / gx * vlcp 91 92 ! for ESM stress 93 ! In ESM, vloc and dvloc have only short term. 94 IF ( ( .not. do_comp_esm ) .or. ( esm_bc .eq. 'pbc' ) ) THEN 95 ! subtract the long-range term 96 IF (.not.do_cutoff_2D) then ! 2D cutoff: do not re-add LR part here (re-added later in stres_loc) 97 g2a = gl (igl) * tpiba2 / 4.d0 98 vlcp = vlcp + fpi / omega * zp * e2 * exp ( - g2a) * (g2a + & 99 1.d0) / (gl (igl) * tpiba2) **2 100 ENDIF 101 END IF 102 dvloc (igl) = vlcp 103 enddo 104!$omp end do nowait 105 deallocate (aux) 106!$omp end parallel 107 ! 108 deallocate (aux1) 109 110 return 111end subroutine dvloc_of_g 112! 113!---------------------------------------------------------------------- 114subroutine dvloc_coul (zp, tpiba2, ngl, gl, omega, dvloc) 115 !---------------------------------------------------------------------- 116 ! 117 ! Fourier transform of the Coulomb potential - For all-electron 118 ! calculations, in specific cases only, for testing purposes 119 ! 120 USE kinds 121 USE constants , ONLY : fpi, e2, eps8 122 implicit none 123 ! 124 integer, intent(in) :: ngl 125 ! the number of shell of G vectors 126 real(DP), intent(in) :: zp, tpiba2, omega, gl (ngl) 127 ! valence pseudocharge 128 ! 2 pi / alat 129 ! the volume of the unit cell 130 ! the moduli of g vectors for each s 131 real(DP), intent(out) :: dvloc (ngl) 132 ! fourier transform: dvloc = D Vloc (g^2) / D g^2 = 4pi e^2/omegai /G^4 133 ! 134 integer :: igl0 135 ! first shell with g != 0 136 137 ! the G=0 component is 0 138 if (gl (1) < eps8) then 139 dvloc (1) = 0.0d0 140 igl0 = 2 141 else 142 igl0 = 1 143 endif 144 145 dvloc (igl0:ngl) = fpi * zp * e2 / omega / ( tpiba2 * gl (igl0:ngl) ) ** 2 146 147return 148end subroutine dvloc_coul 149 150