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