1! 2! Copyright (C) 2001-2018 Quantum ESPRESSO Foundation 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!----------------------------------------------------------------------- 9SUBROUTINE dqvan2( ih, jh, np, ipol, ngy, g, tpiba, qmod, ylmk0, dylmk0,& 10 dqg ) 11 !----------------------------------------------------------------------- 12 !! This routine computes the derivatives of the Fourier transform of 13 !! the Q function needed in stress assuming that the radial fourier 14 !! transform is already computed and stored in table qrad. 15 !! The formula implemented here is: 16 17 !! dq(g,i,j) = sum_lm (-i)^l ap(lm,i,j) * 18 !! ( yr_lm(g^) dqrad(g,l,i,j) + dyr_lm(g^) qrad(g,l,i,j)) 19 ! 20 USE kinds, ONLY: DP 21 USE us, ONLY: dq, qrad 22 USE uspp_param, ONLY: lmaxq, nbetam 23 USE uspp, ONLY: nlx, lpl, lpx, ap, indv, nhtol, nhtolm 24 ! 25 IMPLICIT NONE 26 ! 27 INTEGER, INTENT(IN) :: ngy 28 !! the number of G vectors to compute 29 INTEGER, INTENT(IN) :: ih 30 !! first index of Q(ih,jh) 31 INTEGER, INTENT(IN) :: jh 32 !! second index of Q(ih,jh) 33 INTEGER, INTENT(IN) :: np 34 !! pseudopotential index 35 INTEGER, INTENT(IN) :: ipol 36 !! the polarization of the derivative 37 REAL(DP), INTENT(IN) :: g(3,ngy) 38 !! G vectors 39 REAL(DP), INTENT(IN) :: tpiba 40 !! 2pi/a factor, multiplies G vectors 41 REAL(DP), INTENT(IN) :: qmod(ngy) 42 !! moduli of q+G vectors 43 REAL(DP), INTENT(IN) :: ylmk0(ngy,lmaxq*lmaxq) 44 !! spherical harmonics 45 REAL(DP), INTENT(IN) :: dylmk0(ngy,lmaxq*lmaxq) 46 !! derivetives of spherical harmonics 47 COMPLEX(DP), INTENT(OUT) :: dqg(ngy) 48 !! the fourier transform of interest 49 ! 50 ! ... local variables 51 ! 52 COMPLEX(DP) :: sig 53 ! (-i)^L 54 INTEGER :: nb, mb, ijv, ivl, jvl, ig, lp, l, lm, i0, i1, i2, i3 55 ! the atomic index corresponding to ih 56 ! the atomic index corresponding to jh 57 ! combined index (nb,mb) 58 ! the lm corresponding to ih 59 ! the lm corresponding to jh 60 ! counter on g vectors 61 ! the actual LM 62 ! the angular momentum L 63 ! the possible LM's compatible with ih,j 64 ! counters for interpolation table 65 ! 66 REAL(DP) :: sixth, dqi, qm, px, ux, vx, wx, uvx, pwx, work, work1, qm1 67 ! 1 divided by six 68 ! 1 divided dq 69 ! qmod/dq 70 ! measures for interpolation table 71 ! auxiliary variables for intepolation 72 ! auxiliary variable 73 ! auxiliary variable 74 ! 75 ! compute the indices which correspond to ih,jh 76 ! 77 sixth = 1.d0 / 6.d0 78 dqi = 1 / dq 79 nb = indv(ih, np) 80 mb = indv(jh, np) 81 IF (nb >= mb) THEN 82 ijv = nb * (nb - 1) / 2 + mb 83 ELSE 84 ijv = mb * (mb - 1) / 2 + nb 85 ENDIF 86 ivl = nhtolm(ih, np) 87 jvl = nhtolm(jh, np) 88 ! 89 IF (nb > nbetam .OR. mb > nbetam) & 90 CALL errore (' dqvan2 ', ' wrong dimensions (1)', MAX(nb,mb)) 91 IF (ivl > nlx .OR. jvl > nlx) & 92 CALL errore (' dqvan2 ', ' wrong dimensions (2)', MAX(ivl,jvl)) 93 ! 94 dqg(:) = (0.d0,0.d0) 95 ! 96 ! and make the sum over the non zero LM 97 ! 98 DO lm = 1, lpx(ivl, jvl) 99 lp = lpl(ivl, jvl, lm) 100 ! 101 ! extraction of angular momentum l from lp: 102 ! 103 IF (lp==1) THEN 104 l = 1 105 ELSEIF ( (lp>=2) .AND. (lp<=4) ) THEN 106 l = 2 107 ELSEIF ( (lp>=5) .AND. (lp<=9) ) THEN 108 l = 3 109 ELSEIF ( (lp>=10) .AND. (lp<=16) ) THEN 110 l = 4 111 ELSEIF ( (lp>=17) .AND. (lp<=25) ) THEN 112 l = 5 113 ELSEIF ( (lp>=26) .AND. (lp<=36) ) THEN 114 l = 6 115 ELSEIF ( (lp>=37) .AND. (lp<=49) ) THEN 116 l = 7 117 ELSE 118 CALL errore (' dqvan2 ', ' lp.gt.49 ', lp) 119 ENDIF 120 ! 121 sig = (0.d0, -1.d0)**(l - 1) 122 sig = sig * ap(lp, ivl, jvl) 123 ! 124 qm1 = -1.0_dp ! any number smaller than qmod(1) 125 ! 126!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(qm,px,ux,vx,wx,i0,i1,i2,i3,uvx,pwx,work,work1) 127 DO ig = 1, ngy 128 ! 129 ! calculate quantites depending on the module of G only when needed 130 ! 131#if !defined(_OPENMP) 132 IF ( ABS( qmod(ig) - qm1 ) > 1.0D-6 ) THEN 133#endif 134 qm = qmod (ig) * dqi 135 px = qm - INT(qm) 136 ux = 1.d0 - px 137 vx = 2.d0 - px 138 wx = 3.d0 - px 139 i0 = qm + 1 140 i1 = qm + 2 141 i2 = qm + 3 142 i3 = qm + 4 143 uvx = ux * vx * sixth 144 145 pwx = px * wx * 0.5d0 146 147 work = qrad(i0, ijv, l, np) * uvx * wx + & 148 qrad(i1, ijv, l, np) * pwx * vx - & 149 qrad(i2, ijv, l, np) * pwx * ux + & 150 qrad(i3, ijv, l, np) * px * uvx 151 work1 = - qrad(i0, ijv, l, np) * (ux*vx + vx*wx + ux*wx) * sixth & 152 + qrad(i1, ijv, l, np) * (wx*vx - px*wx - px*vx) * 0.5d0 & 153 - qrad(i2, ijv, l, np) * (wx*ux - px*wx - px*ux) * 0.5d0 & 154 + qrad(i3, ijv, l, np) * (ux*vx - px*ux - px*vx) * sixth 155 156 work1 = work1 * dqi 157 158#if !defined(_OPENMP) 159 qm1 = qmod(ig) 160 ENDIF 161#endif 162 163 dqg(ig) = dqg(ig) + sig * dylmk0(ig, lp) * work / tpiba 164 IF (qmod(ig) > 1.d-9) dqg(ig) = dqg(ig) + & 165 sig * ylmk0(ig, lp) * work1 * tpiba * g(ipol, ig) / qmod(ig) 166 ENDDO 167!$OMP END PARALLEL DO 168 ! 169 ENDDO 170 ! 171 RETURN 172 ! 173END SUBROUTINE dqvan2 174 175