Line data Source code
1 : MODULE m_dfpt_gradient
2 : USE m_juDFT
3 : USE m_constants
4 : USE m_types
5 :
6 : IMPLICIT NONE
7 : CONTAINS
8 0 : subroutine mt_gradient_new(atoms, sphhar, sym, r2FlhMt, GrFshMt)
9 :
10 : use m_gaunt, only : gaunt1
11 :
12 : type(t_atoms), intent(in) :: atoms
13 : type(t_sphhar), intent(in) :: sphhar
14 : type(t_sym), intent(in) :: sym
15 :
16 : real, intent(in) :: r2FlhMt(:, 0:, :)
17 : complex, intent(out) :: GrFshMt(:, :, :, :)
18 :
19 : real :: pfac
20 : real :: tGaunt
21 : integer :: itype
22 : integer :: imesh
23 : integer :: mqn_m
24 : integer :: oqn_l
25 : integer :: mqn_mpp
26 : integer :: lm
27 : integer :: symType
28 : integer :: ilh
29 : integer :: imem
30 :
31 0 : real, allocatable :: rDerFlhMt(:)
32 0 : complex, allocatable :: r2GrFshMtNat(:, :, :, :)
33 :
34 0 : allocate( r2GrFshMtNat(atoms%jmtd, ( atoms%lmaxd + 1)**2, atoms%nat, 3) )
35 0 : allocate( rDerFlhMt(atoms%jmtd) )
36 0 : GrFshMt = cmplx(0., 0.)
37 0 : r2GrFshMtNat = cmplx(0., 0.)
38 0 : rDerFlhMt = 0.
39 :
40 : pfac = sqrt( fpi_const / 3. )
41 0 : do mqn_mpp = -1, 1
42 0 : do itype = 1, atoms%ntype
43 0 : symType = sym%ntypsy(itype)
44 0 : do ilh = 0, sphhar%nlh(symType)
45 0 : oqn_l = sphhar%llh(ilh, symType)
46 0 : do imem = 1, sphhar%nmem(ilh,symType)
47 0 : mqn_m = sphhar%mlh(imem,ilh,symType)
48 :
49 : ! l + 1 block
50 : ! oqn_l - 1 to l, so oqn_l should be < lmax not <= lmax
51 0 : if ( ( abs(mqn_m - mqn_mpp) <= oqn_l + 1 ) .and. ( abs(mqn_m) <= oqn_l ) .and. (oqn_l < atoms%lmax(itype)) ) then
52 0 : lm = ( oqn_l + 1 ) * ( oqn_l + 2 ) + 1 + mqn_m - mqn_mpp
53 0 : call derivative_loc( r2FlhMt(:, ilh, itype), itype, atoms, rDerFlhMt )
54 0 : tGaunt = Gaunt1( oqn_l + 1, oqn_l, 1, mqn_m - mqn_mpp, mqn_m, -mqn_mpp, atoms%lmaxd )
55 0 : do imesh = 1, atoms%jri(itype)
56 : r2GrFshMtNat(imesh, lm, itype, mqn_mpp + 2) = r2GrFshMtNat(imesh, lm, itype, mqn_mpp + 2) + pfac * (-1)**mqn_mpp &
57 : &* tGaunt * (rDerFlhMt(imesh) * sphhar%clnu(imem,ilh,symType) &
58 0 : &- ((oqn_l + 2) * r2FlhMt(imesh, ilh, itype) * sphhar%clnu(imem,ilh,symType) / atoms%rmsh(imesh, itype)))
59 : end do ! imesh
60 : end if ! ( abs(mqn_m - mqn_mpp) <= oqn_l + 1 ) .and. ( abs(mqn_m) <= oqn_l )
61 :
62 : ! l - 1 block
63 0 : if ( ( abs(mqn_m - mqn_mpp) <= oqn_l - 1 ) .and. ( abs(mqn_m) <= oqn_l ) ) then
64 : if ( oqn_l - 1 == -1 ) then
65 : write (*, *) 'oqn_l too low'
66 : end if
67 0 : lm = (oqn_l - 1) * oqn_l + 1 + mqn_m - mqn_mpp
68 : ! This is also a trade of between storage and performance, because derivative is called redundantly, maybe store it?
69 0 : call derivative_loc( r2FlhMt(:, ilh, itype), itype, atoms, rDerFlhMt )
70 0 : tGaunt = Gaunt1( oqn_l - 1, oqn_l, 1, mqn_m - mqn_mpp, mqn_m, -mqn_mpp, atoms%lmaxd )
71 0 : do imesh = 1, atoms%jri(itype)
72 : r2GrFshMtNat(imesh, lm, itype, mqn_mpp + 2) = r2GrFshMtNat(imesh, lm, itype, mqn_mpp + 2) + pfac * (-1)**mqn_mpp &
73 : & * tGaunt * (rDerFlhMt(imesh) * sphhar%clnu(imem,ilh,symType) &
74 0 : & + ((oqn_l - 1) * r2FlhMt(imesh, ilh, itype) * sphhar%clnu(imem,ilh,symType) / atoms%rmsh(imesh, itype)))
75 : end do ! imesh
76 : end if ! ( abs(mqn_m - mqn_mpp) <= oqn_l - 1 ) .and. ( abs(mqn_m) <= oqn_l )
77 : end do ! imem
78 : end do ! ilh
79 : end do ! itype
80 : end do ! mqn_mpp
81 :
82 : ! Conversion from natural to cartesian coordinates
83 0 : do itype = 1, atoms%ntype
84 0 : do oqn_l = 0, atoms%lmax(itype)
85 0 : do mqn_m = -oqn_l, oqn_l
86 0 : lm = oqn_l * (oqn_l + 1) + 1 + mqn_m
87 0 : do imesh = 1, atoms%jri(itype)
88 0 : grFshMt(imesh, lm, itype, 1:3) = matmul( Tmatrix0(1:3, 1:3), r2GrFshMtNat(imesh, lm, itype, 1:3) ) / atoms%rmsh(imesh, itype)**2
89 : end do
90 : end do ! mqn_m
91 : end do ! oqn_l
92 : end do ! itype
93 :
94 0 : end subroutine mt_gradient_new
95 :
96 0 : subroutine derivative_loc(f, itype, atoms, df)
97 :
98 : integer, intent(in) :: itype
99 : type(t_atoms), intent(in) :: atoms
100 : real, intent(in) :: f(atoms%jri(itype))
101 : real, intent(out) :: df(atoms%jri(itype))
102 : real :: h, r, d21, d32, d43, d31, d42, d41, df1, df2, s
103 : real :: y0, y1, y2
104 : integer :: i, n
105 :
106 0 : n = atoms%jri(itype)
107 0 : h = atoms%dx(itype)
108 0 : r = atoms%rmsh(1, itype)
109 :
110 : ! use Lagrange interpolation of 3rd order (and averaging) for points 3 to n
111 0 : d21 = r * (exp(h)-1) ; d32 = d21 * exp(h) ; d43 = d32 * exp(h)
112 0 : d31 = d21 + d32 ; d42 = d32 + d43
113 0 : d41 = d31 + d43
114 : df(1) = d31*d41 / (d21*d32*d42) * f(2) + ( -1d0/d21 - 1d0/d31 - 1d0/d41) * f(1)&
115 0 : & - d21*d41 / (d31*d32*d43) * f(3) + d21*d31 / (d41*d42*d43) * f(4)
116 : df(2) = - d32*d42 / (d21*d31*d41) * f(1) + ( 1d0/d21 - 1d0/d32 - 1d0/d42) * f(2)&
117 0 : & + d21*d42 / (d31*d32*d43) * f(3) - d21*d32 / (d41*d42*d43) * f(4)
118 : df1 = d32*d43 / (d21*d31*d41) * f(1) - d31*d43 / (d21*d32*d42) * f(2) +&
119 0 : & ( 1d0/d31 + 1d0/d32 - 1d0/d43 ) * f(3) + d31*d32 / (d41*d42*d43) * f(4)
120 0 : do i = 3, n - 2
121 0 : d21 = d32 ; d32 = d43 ; d43 = d43 * exp(h)
122 0 : d31 = d42 ; d42 = d42 * exp(h)
123 0 : d41 = d41 * exp(h)
124 : df2 = - d32*d42 / (d21*d31*d41) * f(i-1) + ( 1d0/d21 - 1d0/d32 - 1d0/d42) * f(i) + &
125 0 : & d21*d42 / (d31*d32*d43) * f(i+1) - d21*d32 / (d41*d42*d43) * f(i+2)
126 0 : df(i) = ( df1 + df2 ) / 2
127 : df1 = d32*d43 / (d21*d31*d41) * f(i-1) - d31*d43 / (d21*d32*d42) * f(i) +&
128 0 : & ( 1d0/d31 + 1d0/d32 - 1d0/d43 ) * f(i+1) + d31*d32 / (d41*d42*d43) * f(i+2)
129 : enddo
130 0 : df(n-1) = df1
131 : df(n) = - d42*d43 / (d21*d31*d41) * f(n-3) + d41*d43 / (d21*d32*d42) * f(n-2) -&
132 0 : & d41*d42 / (d31*d32*d43) * f(n-1) + ( 1d0/d41 + 1d0/d42 + 1d0/d43 ) * f(n)
133 : ! for first two points use Lagrange interpolation of second order for log(f(i))
134 : ! or, as a fall-back, Lagrange interpolation with the conditions f(1), f(2), f(3), f'(3).
135 0 : s = sign(1d0,f(1))
136 0 : if(sign(1d0,f(2)) /= s .or. sign(1d0,f(3)) /= s .or. any(abs(f(:3)) < 1e0)) then
137 0 : d21 = r * (exp(h)-1)
138 0 : d32 = d21 * exp(h)
139 0 : d31 = d21 + d32
140 0 : s = df(3) / (d31*d32) - f(1) / (d21*d31**2) + f(2) / (d21*d32**2) - f(3) / (d31**2*d32) - f(3) / (d31*d32**2)
141 0 : df(1) = - (d21+d31) / (d21*d31) * f(1) + d31 / (d21*d32) * f(2) - d21 / (d31*d32) * f(3) + d21*d31 * s
142 :
143 0 : df(2) = - (d21-d32) / (d21*d32) * f(2) - d32 / (d21*d31) * f(1) + d21 / (d31*d32) * f(3) - d21*d32 * s
144 : else
145 0 : y0 = log(abs(f(1)))
146 0 : y1 = log(abs(f(2)))
147 0 : y2 = log(abs(f(3)))
148 0 : df(1) = ( - 3*y0/2 + 2*y1 - y2/2 ) * f(1) / (h*r)
149 0 : df(2) = (y2-y0)/2 * f(2) / (h*r*exp(h))
150 : endif
151 0 : end subroutine derivative_loc
152 :
153 0 : SUBROUTINE sh_to_lh(sym, atoms, sphhar, jspins, radfact, rhosh, rholhreal, rholhimag)
154 :
155 : ! WARNING: This routine will not fold back correctly for activated sym-
156 : ! metry and gradients (rho in l=0 and lattice harmonics do not
157 : ! allow l=1 --> gradient in l=1 is lost)
158 :
159 : TYPE(t_sym), INTENT(IN) :: sym
160 : TYPE(t_atoms), INTENT(IN) :: atoms
161 : TYPE(t_sphhar), INTENT(IN) :: sphhar
162 : INTEGER, INTENT(IN) :: jspins, radfact
163 : COMPLEX, INTENT(IN) :: rhosh(:, :, :, :)
164 : REAL, INTENT(OUT) :: rholhreal(:, 0:, :, :), rholhimag(:, 0:, :, :)
165 :
166 : INTEGER :: iSpin, iType, iAtom, ilh, iMem, ilm, iR
167 : INTEGER :: ptsym, l, m
168 : REAL :: factor
169 :
170 0 : rholhreal = 0.0
171 0 : rholhimag = 0.0
172 :
173 0 : DO iSpin = 1, jspins
174 0 : DO iType = 1, atoms%ntype
175 0 : iAtom = atoms%firstAtom(iType)
176 0 : ptsym = sym%ntypsy(iAtom)
177 0 : DO ilh = 0, sphhar%nlh(ptsym)
178 0 : l = sphhar%llh(iLH, ptsym)
179 0 : DO iMem = 1, sphhar%nmem(ilh, ptsym)
180 0 : m = sphhar%mlh(iMem, ilh, ptsym)
181 0 : ilm = l * (l+1) + m + 1
182 0 : DO iR = 1, atoms%jri(iType)
183 0 : IF ((radfact.EQ.0).AND.(l.EQ.0)) THEN
184 0 : factor = atoms%rmsh(iR, iType) / sfp_const
185 0 : ELSE IF (radfact.EQ.2) THEN
186 0 : factor = atoms%rmsh(iR, iType)**2
187 : ELSE
188 : factor = 1.0
189 : END IF
190 : rholhreal(iR, ilh, iType, iSpin) = &
191 : & rholhreal(iR, ilh, iType, iSpin) + &
192 0 : & real(rhosh(iR, ilm, iatom, iSpin) * conjg(sphhar%clnu(iMem, ilh, ptsym))) * factor
193 : rholhimag(iR, ilh, iType, iSpin) = &
194 : & rholhimag(iR, ilh, iType, iSpin) + &
195 0 : & aimag(rhosh(iR, ilm, iatom, iSpin) * conjg(sphhar%clnu(iMem, ilh, ptsym))) * factor
196 : END DO
197 : END DO
198 : END DO
199 : END DO
200 : END DO
201 :
202 0 : END SUBROUTINE sh_to_lh
203 : END MODULE m_dfpt_gradient
|
