Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2021 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_dfpt_dynmat
7 : USE m_types
8 : USE m_constants
9 :
10 : IMPLICIT NONE
11 :
12 : CONTAINS
13 0 : SUBROUTINE dfpt_dynmat_row(fi, stars, starsq, sphhar, xcpot, nococonv, hybdat, fmpi, qpts, iQ, iDtype_row, iDir_row, &
14 : eig_id, dfpt_eig_id, dfpt_eig_id2, enpara, results, results1, l_real, &
15 0 : rho, vTot, grRho3, grVext3, grVC3, denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im, dyn_row, &
16 0 : E2ndOrdII, sigma_ext, sigma_gext, q_eig_id)
17 : USE m_step_function
18 : USE m_convol
19 : USE m_dfpt_vgen
20 : USE m_vgen_coulomb
21 : USE m_dfpt_eii2
22 :
23 : TYPE(t_fleurinput), INTENT(IN) :: fi
24 : TYPE(t_stars), INTENT(IN) :: stars, starsq
25 : TYPE(t_sphhar), INTENT(IN) :: sphhar
26 : CLASS(t_xcpot), INTENT(IN) :: xcpot
27 : TYPE(t_nococonv), INTENT(IN) :: nococonv
28 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
29 : TYPE(t_mpi), INTENT(IN) :: fmpi
30 : TYPE(t_kpts), INTENT(IN) :: qpts
31 :
32 : TYPE(t_potden), INTENT(INOUT) :: rho, vTot
33 : TYPE(t_potden), INTENT(INOUT) :: grRho3(3), grVext3(3), grVC3(3)
34 : TYPE(t_potden), INTENT(INOUT) :: denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im
35 :
36 : TYPE(t_enpara), INTENT(INOUT) :: enpara
37 : TYPE(t_results), INTENT(INOUT) :: results, results1
38 :
39 : LOGICAL, INTENT(IN) :: l_real
40 :
41 : INTEGER, INTENT(IN) :: iQ, iDtype_row, iDir_row, eig_id, dfpt_eig_id, dfpt_eig_id2
42 :
43 : COMPLEX, INTENT(INOUT) :: dyn_row(:)
44 :
45 : COMPLEX, INTENT(INOUT) :: E2ndOrdII(:,:)
46 :
47 : COMPLEX, OPTIONAL, INTENT(IN) :: sigma_ext(2), sigma_gext(3,2)
48 :
49 : INTEGER, OPTIONAL, INTENT(IN) :: q_eig_id
50 :
51 0 : TYPE(t_fftgrid) :: fftgrid_dummy
52 0 : TYPE(t_potden) :: rho_dummy, rho1_dummy, vExt1, vExt1Im, grgrVCq, grgrVCqIm
53 0 : TYPE(t_hub1data) :: hub1data
54 :
55 : INTEGER :: col_index, row_index, iDtype_col, iDir_col, iType, iDir, iSpin
56 : COMPLEX :: tempval
57 : LOGICAL :: bare_mode
58 :
59 : REAL :: qvec(3)
60 0 : REAL :: e2_vm(fi%atoms%nat)
61 :
62 : complex :: sigma_loc(2)
63 :
64 0 : COMPLEX, ALLOCATABLE :: dyn_row_HF(:), dyn_row_eigen(:), dyn_row_int(:)
65 : COMPLEX, ALLOCATABLE :: theta1full(:, :, :), theta1full0(:, :, :)!, theta2(:, :, :)
66 0 : COMPLEX, ALLOCATABLE :: theta1_pw(:, :, :), theta1_pw0(:, :, :)!,theta2_pw(:, :, :)
67 0 : COMPLEX, ALLOCATABLE :: pww(:), pwwq(:), pww2(:), pwwq2(:)
68 0 : COMPLEX, ALLOCATABLE :: rho_pw(:), denIn1_pw(:), rho_vac(:,:,:), denIn1_vac(:,:,:)
69 0 : REAL, ALLOCATABLE :: rho_mt(:,:,:), grRho_mt(:,:,:), denIn1_mt(:,:,:), denIn1_mt_Im(:,:,:)
70 :
71 0 : bare_mode = .FALSE.
72 :
73 0 : ALLOCATE(dyn_row_HF(SIZE(dyn_row)), dyn_row_eigen(SIZE(dyn_row)), dyn_row_int(SIZE(dyn_row)))
74 0 : ALLOCATE(theta1full(0:27*starsq%mx1*starsq%mx2*starsq%mx3-1,fi%atoms%ntype,3))
75 0 : ALLOCATE(theta1full0(0:27*stars%mx1*stars%mx2*stars%mx3-1,fi%atoms%ntype,3))
76 0 : ALLOCATE(theta1_pw(starsq%ng3,fi%atoms%ntype,3))
77 0 : ALLOCATE(theta1_pw0(stars%ng3,fi%atoms%ntype,3))
78 0 : ALLOCATE(pww(stars%ng3),pwwq(starsq%ng3))
79 0 : ALLOCATE(pww2(stars%ng3),pwwq2(starsq%ng3))
80 :
81 0 : ALLOCATE(denIn1_mt(fi%atoms%jmtd,0:sphhar%nlhd,fi%atoms%ntype),denIn1_mt_Im(fi%atoms%jmtd,0:sphhar%nlhd,fi%atoms%ntype))
82 0 : ALLOCATE(denIn1_pw(starsq%ng3),rho_pw(stars%ng3))
83 0 : ALLOCATE(rho_mt(fi%atoms%jmtd,0:sphhar%nlhd,fi%atoms%ntype),grRho_mt(fi%atoms%jmtd,0:sphhar%nlhd,fi%atoms%ntype))
84 0 : IF (fi%input%film) THEN
85 0 : ALLOCATE(denIn1_vac(fi%vacuum%nmzd,starsq%ng2,2))
86 0 : ALLOCATE(rho_vac(fi%vacuum%nmzd,stars%ng2,2))
87 : END IF
88 :
89 0 : CALL rho_dummy%copyPotDen(rho)
90 0 : CALL rho1_dummy%copyPotDen(denIn1)
91 0 : CALL rho_dummy%resetPotDen()
92 0 : CALL rho1_dummy%resetPotDen()
93 :
94 0 : CALL grgrVCq%copyPotDen(vTot1)
95 0 : CALL grgrVCq%resetPotDen()
96 0 : CALL grgrVCqIm%copyPotDen(vTot1)
97 0 : CALL grgrVCqIm%resetPotDen()
98 :
99 0 : qvec = qpts%bk(:, iQ)
100 :
101 0 : theta1full = CMPLX(0.0,0.0)
102 0 : theta1full0 = CMPLX(0.0,0.0)
103 :
104 0 : CALL stepf_analytical(fi%sym, starsq, fi%atoms, fi%input, fi%cell, fmpi, fftgrid_dummy, qvec, iDtype_row, iDir_row, 1, theta1full)
105 0 : CALL stepf_analytical(fi%sym, stars, fi%atoms, fi%input, fi%cell, fmpi, fftgrid_dummy, [0.0,0.0,0.0], iDtype_row, iDir_row, 1, theta1full0)
106 : !CALL stepf_analytical(fi%sym, stars, fi%atoms, fi%input, fi%cell, fmpi, fftgrid_dummy, [0.0,0.0,0.0], iDtype_row, iDir_row, 2, theta2)
107 :
108 0 : DO iType = 1, fi%atoms%ntype
109 0 : DO iDir = 1, 3
110 0 : fftgrid_dummy%grid = theta1full(0:, iType, iDir)
111 0 : CALL fftgrid_dummy%takeFieldFromGrid(starsq, theta1_pw(:, iType, iDir))
112 0 : theta1_pw(:, iType, iDir) = theta1_pw(:, iType, iDir) * 3 * starsq%mx1 * 3 * starsq%mx2 * 3 * starsq%mx3
113 0 : CALL fftgrid_dummy%perform_fft(forward=.false.)
114 0 : theta1full(0:, iType, iDir) = fftgrid_dummy%grid
115 :
116 0 : fftgrid_dummy%grid = theta1full0(0:, iType, iDir)
117 0 : CALL fftgrid_dummy%takeFieldFromGrid(stars, theta1_pw0(:, iType, iDir))
118 0 : theta1_pw0(:, iType, iDir) = theta1_pw0(:, iType, iDir) * 3 * stars%mx1 * 3 * stars%mx2 * 3 * stars%mx3
119 0 : CALL fftgrid_dummy%perform_fft(forward=.false.)
120 0 : theta1full0(0:, iType, iDir) = fftgrid_dummy%grid
121 : END DO
122 : END DO
123 :
124 0 : row_index = 3 * (iDtype_row - 1) + iDir_row
125 :
126 0 : dyn_row = CMPLX(0.0,0.0)
127 0 : dyn_row_HF = CMPLX(0.0,0.0)
128 0 : dyn_row_eigen = CMPLX(0.0,0.0)
129 0 : dyn_row_int = CMPLX(0.0,0.0)
130 :
131 0 : denIn1_pw = (denIn1%pw(:,1)+denIn1%pw(:,fi%input%jspins))/(3.0-fi%input%jspins)
132 0 : denIn1_mt = (denIn1%mt(:,0:,:,1)+denIn1%mt(:,0:,:,fi%input%jspins))/(3.0-fi%input%jspins)
133 0 : IF (fi%input%film) denIn1_vac = (denIn1%vac(:,:,:,1)+denIn1%vac(:,:,:,fi%input%jspins))/(3.0-fi%input%jspins)
134 : ! Get "full" denIn1:
135 0 : denIn1_mt(:,0:,iDtype_row) = denIn1_mt(:,0:,iDtype_row) - (grRho3(iDir_row)%mt(:,0:,iDtype_row,1)+grRho3(iDir_row)%mt(:,0:,iDtype_row,fi%input%jspins))/(3.0-fi%input%jspins)
136 0 : denIn1_mt_Im = (denIn1Im%mt(:,0:,:,1)+denIn1Im%mt(:,0:,:,fi%input%jspins))/(3.0-fi%input%jspins)
137 :
138 0 : DO iDtype_col = 1, fi%atoms%ntype
139 0 : DO iDir_col = 1, 3
140 0 : IF (fmpi%irank==0) write(9989,*) "------------------"
141 0 : IF (fmpi%irank==0) write(9989,*) "Atom:", iDtype_col, "Direction", iDir_col
142 0 : tempval = CMPLX(0.0,0.0)
143 0 : col_index = 3 * (iDtype_col - 1) + iDir_col
144 :
145 : ! First calculate the HF contributions.
146 : ! \rho(1)V_{ext}(1) integral over whole unit cell
147 :
148 : ! Get V_{ext}(1) for \alpha, i with gradient cancellation
149 0 : CALL vExt1%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.TRUE.)
150 0 : CALL vExt1Im%init(starsq, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.FALSE.)
151 0 : sigma_loc = cmplx(0.0,0.0)
152 : !IF (iDir_col==3) sigma_loc = -sigma_ext
153 : CALL dfpt_vgen(hybdat,fi%field,fi%input,xcpot,fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
154 : fi%juphon,fi%cell,fmpi,fi%noco,nococonv,rho_dummy,vTot,&
155 0 : starsq,rho1_dummy,vExt1,.FALSE.,vExt1Im,rho1_dummy,iDtype_col,iDir_col,[0,0],sigma_loc)
156 :
157 : ! IR integral:
158 0 : pwwq2 = CMPLX(0.0,0.0)
159 0 : CALL dfpt_convol_big(1, starsq, stars, vExt1%pw(:,1), CMPLX(1.0,0.0)*stars%ufft, pwwq2)
160 0 : CALL dfpt_int_pw(starsq, fi%cell, denIn1_pw, pwwq2, tempval)
161 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
162 0 : IF (fmpi%irank==0) write(9989,*) "IR rho1 V1ext new ", tempval
163 0 : tempval = CMPLX(0.0,0.0)
164 :
165 0 : IF (fi%input%film) THEN
166 0 : CALL dfpt_int_vac(starsq,fi%vacuum,fi%cell,denIn1_vac,vExt1%vac(:,:,:,1),tempval)
167 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
168 0 : IF (fmpi%irank==0) write(9989,*) "VAC rho1 V1ext ", tempval
169 0 : tempval = CMPLX(0.0,0.0)
170 : END IF
171 :
172 : ! MT integral:
173 : ! If we use gradient cancellation, remove it from rho1
174 : IF (.NOT.bare_mode) denIn1_mt(:,0:,iDtype_row) = &
175 : denIn1_mt(:,0:,iDtype_row) + &
176 0 : (grRho3(iDir_row)%mt(:,0:,iDtype_row,1)+grRho3(iDir_row)%mt(:,0:,iDtype_row,fi%input%jspins))/(3.0-fi%input%jspins)
177 0 : DO iType = 1, fi%atoms%ntype
178 0 : IF (fmpi%irank==0) write(9989,*) "Loop atom:", iType
179 0 : CALL dfpt_int_mt(fi%atoms, sphhar, fi%sym, iType, denIn1_mt, denIn1_mt_Im, vExt1%mt(:,0:,:,1), vExt1Im%mt(:,0:,:,1), tempval)
180 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
181 0 : IF (fmpi%irank==0) write(9989,*) " MT rho1 V1ext ", tempval
182 0 : tempval = CMPLX(0.0,0.0)
183 : END DO
184 0 : IF (fmpi%irank==0) write(9989,*) "End atom loop"
185 : IF (.NOT.bare_mode) denIn1_mt(:,0:,iDtype_row) = &
186 : denIn1_mt(:,0:,iDtype_row) - &
187 0 : (grRho3(iDir_row)%mt(:,0:,iDtype_row,1)+grRho3(iDir_row)%mt(:,0:,iDtype_row,fi%input%jspins))/(3.0-fi%input%jspins)
188 :
189 0 : sigma_loc = cmplx(0.0,0.0)
190 : !IF (iDir_col==3) sigma_loc = sigma_gext(iDir_row,:)
191 : !IF (iDir_row==3) sigma_loc = sigma_gext(iDir_col,:)
192 : CALL vgen_coulomb(1, fmpi, fi%input, fi%field, fi%vacuum, fi%sym, fi%juphon, starsq, fi%cell, &
193 : & sphhar, fi%atoms, .TRUE., rho1_dummy, grgrVCq, sigma_loc, &
194 : & dfptdenimag=rho1_dummy, dfptvCoulimag=grgrVCqIm,dfptden0=rho1_dummy,stars2=stars,iDtype=iDtype_col,iDir=iDir_col,iDir2=iDir_row, &
195 0 : & sigma_disc2=MERGE(sigma_ext,[cmplx(0.0,0.0),cmplx(0.0,0.0)],iDir_col==3.AND.iDir_row==3.AND..FALSE.))
196 0 : IF (iDtype_col==iDtype_row) THEN
197 0 : e2_vm = 0.0
198 0 : CALL dfpt_e2_madelung(fi%atoms,fi%input%jspins,rho1_dummy%mt(:,0,:,:),grgrVCq%mt(:,0,:,1),e2_vm(:))
199 0 : E2ndOrdII(row_index,col_index) = E2ndOrdII(row_index,col_index) - e2_vm(iDtype_col)
200 0 : e2_vm = 0.0
201 0 : CALL dfpt_e2_madelung(fi%atoms,fi%input%jspins,rho1_dummy%mt(:,0,:,:),grgrVCqIm%mt(:,0,:,1),e2_vm(:))
202 0 : E2ndOrdII(row_index,col_index) = E2ndOrdII(row_index,col_index) - ImagUnit*e2_vm(iDtype_col)
203 : ELSE
204 0 : E2ndOrdII(row_index,col_index) = E2ndOrdII(row_index,col_index) + fi%atoms%zatom(iDtype_row)*(grgrVCq%mt(1,0,iDtype_row,1)+ImagUnit*grgrVCqIm%mt(1,0,iDtype_row,1))/sfp_const
205 : END IF
206 0 : CALL grgrVCq%resetPotDen()
207 0 : CALL grgrVCqIm%resetPotDen()
208 : ! Various V_ext integrals:
209 : ! IR:
210 0 : rho_pw = (rho%pw(:,1)+rho%pw(:,fi%input%jspins))/(3.0-fi%input%jspins)
211 0 : pwwq2 = CMPLX(0.0,0.0)
212 0 : CALL dfpt_convol_big(2, stars, starsq, rho_pw, theta1full(0:, iDtype_row, iDir_row), pwwq2)
213 0 : CALL dfpt_int_pw(starsq, fi%cell, pwwq2, vExt1%pw(:,1), tempval)
214 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
215 0 : IF (fmpi%irank==0) write(9989,*) "IR Theta1 rho V1ext new ", tempval
216 0 : tempval = CMPLX(0.0,0.0)
217 :
218 : ! MT:
219 : ! TODO: Apply baremode correction to MT correction and corresponding 2nd integral. They
220 : ! are the last two *big* terms.
221 : IF (.NOT.bare_mode) denIn1_mt(:,0:,iDtype_row) = &
222 : denIn1_mt(:,0:,iDtype_row) + &
223 0 : (grRho3(iDir_row)%mt(:,0:,iDtype_row,1)+grRho3(iDir_row)%mt(:,0:,iDtype_row,fi%input%jspins))/(3.0-fi%input%jspins)
224 0 : grRho_mt = grVC3(iDir_col)%mt(:,0:,:,1)
225 0 : CALL dfpt_int_mt(fi%atoms, sphhar, fi%sym, iDtype_col, denIn1_mt, denIn1_mt_Im, grRho_mt, 0*grRho_mt, tempval)
226 0 : dyn_row_int(col_index) = dyn_row_int(col_index) + tempval
227 0 : IF (fmpi%irank==0) write(9989,*) "MT rho1 grVC ", tempval
228 0 : tempval = CMPLX(0.0,0.0)
229 : IF (.NOT.bare_mode) denIn1_mt(:,0:,iDtype_row) = &
230 : denIn1_mt(:,0:,iDtype_row) - &
231 0 : (grRho3(iDir_row)%mt(:,0:,iDtype_row,1)+grRho3(iDir_row)%mt(:,0:,iDtype_row,fi%input%jspins))/(3.0-fi%input%jspins)
232 :
233 : IF (.NOT.bare_mode) vC1%mt(:,0:,iDtype_row,1) = &
234 : vC1%mt(:,0:,iDtype_row,1) + &
235 0 : grVC3(iDir_row)%mt(:,0:,iDtype_row,1)
236 0 : grRho_mt = -(grRho3(iDir_col)%mt(:,0:,:,1)+grRho3(iDir_col)%mt(:,0:,:,fi%input%jspins))/(3.0-fi%input%jspins)
237 0 : CALL dfpt_int_mt(fi%atoms, sphhar, fi%sym, iDtype_col, vC1%mt(:,0:,:,1), vC1Im%mt(:,0:,:,1), grRho_mt, 0*grRho_mt, tempval)
238 0 : dyn_row_int(col_index) = dyn_row_int(col_index) + tempval
239 0 : IF (fmpi%irank==0) write(9989,*) "MT grRho V1C ", tempval
240 0 : tempval = CMPLX(0.0,0.0)
241 : IF (.NOT.bare_mode) vC1%mt(:,0:,iDtype_row,1) = &
242 : vC1%mt(:,0:,iDtype_row,1) - &
243 0 : grVC3(iDir_row)%mt(:,0:,iDtype_row,1)
244 :
245 : IF (.NOT.bare_mode) THEN
246 0 : IF (.NOT.bare_mode) vExt1%mt(:,0:,iDtype_col,:) = vExt1%mt(:,0:,iDtype_col,:) + grVext3(iDir_col)%mt(:,0:,iDtype_col,:)
247 0 : grRho_mt = -(grRho3(iDir_row)%mt(:,0:,:,1)+grRho3(iDir_row)%mt(:,0:,:,fi%input%jspins))/(3.0-fi%input%jspins)
248 0 : CALL dfpt_int_mt(fi%atoms, sphhar, fi%sym, iDtype_row, grRho_mt, 0*grRho_mt, vExt1%mt(:,0:,:,1), vExt1Im%mt(:,0:,:,1), tempval)
249 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
250 0 : IF (fmpi%irank==0) write(9989,*) "MT correction grRho V1ext ", tempval
251 0 : tempval = CMPLX(0.0,0.0)
252 0 : IF (.NOT.bare_mode) vExt1%mt(:,0:,iDtype_col,:) = vExt1%mt(:,0:,iDtype_col,:) - grVext3(iDir_col)%mt(:,0:,iDtype_col,:)
253 : END IF
254 :
255 : ! SF:
256 0 : rho_mt = (rho%mt(:,0:,:,1)+rho%mt(:,0:,:,fi%input%jspins))/(3.0-fi%input%jspins)
257 0 : CALL dfpt_int_mt_sf(fi%atoms, sphhar, fi%sym, iDir_row, iDtype_row, rho_mt, vExt1%mt(:,0:,:,1), vExt1Im%mt(:,0:,:,1), tempval)
258 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
259 0 : IF (fmpi%irank==0) write(9989,*) "SF rho Vext1 ", tempval
260 0 : tempval = CMPLX(0.0,0.0)
261 :
262 : IF (.NOT.bare_mode) vC1%mt(:,0:,iDtype_row,1) = &
263 : vC1%mt(:,0:,iDtype_row,1) + &
264 0 : grVC3(iDir_row)%mt(:,0:,iDtype_row,1)
265 0 : CALL dfpt_int_mt_sf(fi%atoms, sphhar, fi%sym, iDir_col, iDtype_col, rho_mt, vC1%mt(:,0:,:,1), -vC1Im%mt(:,0:,:,1), tempval)
266 0 : dyn_row_int(col_index) = dyn_row_int(col_index) + tempval
267 0 : IF (fmpi%irank==0) write(9989,*) "SF rho VC1 ", tempval
268 0 : tempval = CMPLX(0.0,0.0)
269 : IF (.NOT.bare_mode) vC1%mt(:,0:,iDtype_row,1) = &
270 : vC1%mt(:,0:,iDtype_row,1) - &
271 0 : grVC3(iDir_row)%mt(:,0:,iDtype_row,1)
272 :
273 : ! Miscellaneous integrals:
274 0 : rho_pw = (rho%pw(:,1)+rho%pw(:,fi%input%jspins))/(3.0-fi%input%jspins)
275 0 : pwwq2 = CMPLX(0.0,0.0)
276 0 : CALL dfpt_convol_big(2, stars, starsq, rho_pw, theta1full(0:, iDtype_col, iDir_col), pwwq2)
277 0 : CALL dfpt_int_pw(starsq, fi%cell, vC1%pw(:,1), pwwq2, tempval)
278 0 : dyn_row_int(col_index) = dyn_row_int(col_index) + tempval
279 0 : IF (fmpi%irank==0) write(9989,*) "IR V1C rho Theta1 new ", tempval
280 0 : tempval = CMPLX(0.0,0.0)
281 :
282 0 : DO iSpin = 1, fi%input%jspins
283 0 : IF (fmpi%irank==0) write(9989,*) "Loop spin:", iSpin
284 : ! TODO: Ensure, that vTot/denIn1 is diagonal here, not 2x2.
285 0 : pwwq2 = CMPLX(0.0,0.0)
286 0 : CALL dfpt_convol_big(2, stars, starsq, vTot%pw(:, iSpin), theta1full(0:, iDtype_col, iDir_col), pwwq2)
287 0 : CALL dfpt_int_pw(starsq, fi%cell, denIn1%pw(:,iSpin), pwwq2, tempval)
288 0 : dyn_row_int(col_index) = dyn_row_int(col_index) + tempval
289 0 : IF (fmpi%irank==0) write(9989,*) " IR rho1 vTot Theta1 new ", tempval
290 0 : tempval = CMPLX(0.0,0.0)
291 : END DO
292 0 : IF (fmpi%irank==0) write(9989,*) "End spin loop"
293 :
294 0 : IF (iDtype_row==iDtype_col) THEN
295 0 : CALL vExt1%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.TRUE.)
296 0 : CALL vExt1Im%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT, l_dfpt=.FALSE.)
297 : ! Get V_{ext}(1) for \alpha, i, q=0 with gradient cancellation
298 0 : sigma_loc = cmplx(0.0,0.0)
299 : !IF (iDir_col==3) sigma_loc = -sigma_ext
300 : CALL dfpt_vgen(hybdat,fi%field,fi%input,xcpot,fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
301 : fi%juphon,fi%cell,fmpi,fi%noco,nococonv,rho_dummy,vTot,&
302 0 : stars,rho_dummy,vExt1,.FALSE.,vExt1Im,rho_dummy,iDtype_col,iDir_col,[0,0],sigma_loc)
303 :
304 : ! Integrals:
305 0 : rho_pw = (grRho3(iDir_row)%pw(:,1)+grRho3(iDir_row)%pw(:,fi%input%jspins))/(3.0-fi%input%jspins)
306 0 : pww2 = CMPLX(0.0,0.0)
307 0 : CALL dfpt_convol_big(1, stars, stars, vExt1%pw(:,1), CMPLX(1.0,0.0)*stars%ufft, pww2)
308 0 : CALL dfpt_int_pw(stars, fi%cell, rho_pw, pww2, tempval)
309 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
310 0 : IF (fmpi%irank==0) write(9989,*) "IR grRho V1ext0 new ", tempval
311 0 : tempval = CMPLX(0.0,0.0)
312 :
313 0 : IF (fi%input%film) THEN
314 0 : rho_vac = (grRho3(iDir_row)%vac(:,:,:,1)+grRho3(iDir_row)%vac(:,:,:,fi%input%jspins))/(3.0-fi%input%jspins)
315 0 : CALL dfpt_int_vac(stars,fi%vacuum,fi%cell,rho_vac,vExt1%vac(:,:,:,1),tempval)
316 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
317 0 : IF (fmpi%irank==0) write(9989,*) "VAC grRho V1ext0 ", tempval
318 0 : tempval = CMPLX(0.0,0.0)
319 : END IF
320 :
321 0 : rho_mt = (rho%mt(:,0:,:,1)+rho%mt(:,0:,:,fi%input%jspins))/(3.0-fi%input%jspins)
322 0 : rho_pw = -(rho%pw(:,1)+rho%pw(:,fi%input%jspins))/(3.0-fi%input%jspins)
323 0 : DO iType = 1, fi%atoms%ntype
324 0 : IF (fmpi%irank==0) write(9989,*) "Loop atom:", iType
325 0 : pww2 = CMPLX(0.0,0.0)
326 0 : CALL dfpt_convol_big(1, stars, stars, rho_pw, theta1full0(0:, iType, iDir_row), pww2)
327 0 : CALL dfpt_int_pw(stars, fi%cell, pww2, vExt1%pw(:,1), tempval)
328 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
329 0 : IF (fmpi%irank==0) write(9989,*) " IR Theta1 rho V1ext0 new ", tempval
330 0 : tempval = CMPLX(0.0,0.0)
331 :
332 0 : IF (.NOT.bare_mode) vExt1%mt(:,0:,iDtype_col,:) = vExt1%mt(:,0:,iDtype_col,:) + grVext3(iDir_col)%mt(:,0:,iDtype_col,:)
333 0 : grRho_mt = (grRho3(iDir_row)%mt(:,0:,:,1)+grRho3(iDir_row)%mt(:,0:,:,fi%input%jspins))/(3.0-fi%input%jspins)
334 0 : CALL dfpt_int_mt(fi%atoms, sphhar, fi%sym, iType, grRho_mt, 0*grRho_mt, vExt1%mt(:,0:,:,1), vExt1Im%mt(:,0:,:,1), tempval)
335 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
336 0 : IF (fmpi%irank==0) write(9989,*) " MT grRho V1ext0 ", tempval
337 0 : tempval = CMPLX(0.0,0.0)
338 0 : IF (.NOT.bare_mode) vExt1%mt(:,0:,iDtype_col,:) = vExt1%mt(:,0:,iDtype_col,:) - grVext3(iDir_col)%mt(:,0:,iDtype_col,:)
339 :
340 0 : CALL dfpt_int_mt_sf(fi%atoms, sphhar, fi%sym, iDir_row, iType, -rho_mt, vExt1%mt(:,0:,:,1), vExt1Im%mt(:,0:,:,1), tempval)
341 0 : dyn_row_HF(col_index) = dyn_row_HF(col_index) + tempval
342 0 : IF (fmpi%irank==0) write(9989,*) " SF rho V1ext0 ", tempval
343 0 : tempval = CMPLX(0.0,0.0)
344 : END DO
345 0 : IF (fmpi%irank==0) write(9989,*) "End atom loop"
346 :
347 0 : rho_pw = (rho%pw(:,1)+rho%pw(:,fi%input%jspins))/(3.0-fi%input%jspins)
348 0 : pww2 = CMPLX(0.0,0.0)
349 0 : CALL dfpt_convol_big(1, stars, stars, rho_pw, theta1full0(0:, iDtype_col, iDir_col), pww2)
350 0 : CALL dfpt_int_pw(stars, fi%cell, grVC3(iDir_row)%pw(:,1), pww2, tempval)
351 0 : dyn_row_int(col_index) = dyn_row_int(col_index) + tempval
352 0 : IF (fmpi%irank==0) write(9989,*) "IR grVC rho Theta1 new ", tempval
353 0 : tempval = CMPLX(0.0,0.0)
354 :
355 0 : DO iSpin = 1, fi%input%jspins
356 0 : IF (fmpi%irank==0) write(9989,*) "Loop spin:", iSpin
357 : ! TODO: Ensure, that vTot/gradrho is diagonal here, not 2x2 [NOCO].
358 0 : pww2 = CMPLX(0.0,0.0)
359 0 : CALL dfpt_convol_big(1, stars, stars, vTot%pw(:,iSpin), theta1full0(0:,iDtype_col,iDir_col), pww2)
360 0 : CALL dfpt_int_pw(stars, fi%cell, grRho3(iDir_row)%pw(:,iSpin), pww2, tempval)
361 0 : dyn_row_int(col_index) = dyn_row_int(col_index) + tempval
362 0 : IF (fmpi%irank==0) write(9989,*) " IR grRho vTot Theta1 new ", tempval
363 0 : tempval = CMPLX(0.0,0.0)
364 : END DO
365 0 : IF (fmpi%irank==0) write(9989,*) "End spin loop"
366 : END IF
367 :
368 0 : IF (fmpi%irank==0) write(9990,*) qvec, iDtype_row, iDir_row, iDtype_col, iDir_col
369 0 : IF (fmpi%irank==0) write(9990,*) "HF :", conjg(dyn_row_HF(col_index))
370 0 : IF (fmpi%irank==0) write(9990,*) "int :", conjg(dyn_row_int(col_index))
371 :
372 : ! Calculate the contributions to the dynamical matrix that stem
373 : ! from terms related to occupation numbers and the eigenenergies.
374 0 : IF (.NOT.PRESENT(q_eig_id)) THEN
375 : CALL dfpt_dynmat_eigen(fi, results, results1, fmpi, enpara, nococonv, &
376 : stars, starsq, sphhar, rho, hub1data, vTot, vTot, vTot1, vTot1Im, &
377 : eig_id, dfpt_eig_id, dfpt_eig_id2, iDir_col, iDtype_col, iDir_row, iDtype_row, &
378 : theta1_pw0(:,iDtype_col,iDir_col), theta1_pw(:,iDtype_col,iDir_col), &
379 0 : qvec, l_real, dyn_row_eigen(col_index),[1,1,1,1,1,1,1,1,1,1,1])
380 : ELSE
381 : CALL dfpt_dynmat_eigen(fi, results, results1, fmpi, enpara, nococonv, &
382 : stars, starsq, sphhar, rho, hub1data, vTot, vTot, vTot1, vTot1Im, &
383 : eig_id, dfpt_eig_id, dfpt_eig_id2, iDir_col, iDtype_col, iDir_row, iDtype_row, &
384 : theta1_pw0(:,iDtype_col,iDir_col), theta1_pw(:,iDtype_col,iDir_col), &
385 0 : qvec, l_real, dyn_row_eigen(col_index),[1,1,1,1,1,1,1,1,1,1,1],q_eig_id)
386 : END IF
387 0 : IF (fmpi%irank==0) write(9990,*) "eigen:", dyn_row_eigen(col_index)
388 0 : IF (fmpi%irank==0) write(9989,*) "Eii2:", conjg(E2ndOrdII(3*(iDtype_row-1)+iDir_row,3*(iDtype_col-1)+iDir_col))
389 0 : IF (fmpi%irank==0) write(9990,*) "Eii2:", conjg(E2ndOrdII(3*(iDtype_row-1)+iDir_row,3*(iDtype_col-1)+iDir_col))
390 : END DO
391 : END DO
392 :
393 0 : dyn_row = conjg(dyn_row_HF) + conjg(dyn_row_int) + dyn_row_eigen
394 :
395 0 : END SUBROUTINE dfpt_dynmat_row
396 :
397 0 : SUBROUTINE dfpt_int_pw(stars, cell, pw_conj, pw_pure, pw_int)
398 : TYPE(t_stars), INTENT(IN) :: stars
399 : TYPE(t_cell), INTENT(IN) :: cell
400 :
401 : COMPLEX, INTENT(IN) :: pw_conj(:), pw_pure(:)
402 :
403 : COMPLEX, INTENT(INOUT) :: pw_int
404 :
405 0 : pw_int = pw_int + cell%omtil * DOT_PRODUCT(pw_conj(:stars%ng3),pw_pure(:stars%ng3))
406 0 : END SUBROUTINE dfpt_int_pw
407 :
408 0 : SUBROUTINE dfpt_int_mt(atoms, sphhar, sym, nat, mt_conj, mt_conj_im, mt_pure, mt_pure_im, mt_int)
409 : USE m_intgr, ONLY: intgr3, intgr3LinIntp
410 :
411 : TYPE(t_atoms), INTENT(IN) :: atoms
412 : TYPE(t_sphhar), INTENT(IN) :: sphhar
413 : TYPE(t_sym), INTENT(IN) :: sym
414 :
415 : INTEGER, INTENT(IN) :: nat
416 :
417 : REAL, INTENT(IN) :: mt_conj(:,0:,:), mt_conj_im(:,0:,:), mt_pure(:,0:,:), mt_pure_im(:,0:,:)
418 :
419 : COMPLEX, INTENT(INOUT) :: mt_int
420 :
421 : REAL :: dpdot_re, dpdot_im
422 : COMPLEX :: tmt
423 : INTEGER :: j, lh
424 :
425 0 : REAL :: dpj_re(atoms%jmtd), dpj_im(atoms%jmtd)
426 :
427 0 : tmt = CMPLX(0.0,0.0)
428 0 : DO lh = 0, sphhar%nlh(sym%ntypsy(nat))
429 0 : DO j = 1, atoms%jri(nat)
430 0 : dpj_re(j) = mt_conj(j,lh,nat)*mt_pure(j,lh,nat)+mt_conj_im(j,lh,nat)*mt_pure_im(j,lh,nat)
431 0 : dpj_im(j) = mt_conj(j,lh,nat)*mt_pure_im(j,lh,nat)-mt_conj_im(j,lh,nat)*mt_pure(j,lh,nat)
432 : END DO
433 0 : CALL intgr3(dpj_re,atoms%rmsh(1,nat),atoms%dx(nat),atoms%jri(nat),dpdot_re)
434 0 : CALL intgr3(dpj_im,atoms%rmsh(1,nat),atoms%dx(nat),atoms%jri(nat),dpdot_im)
435 0 : tmt = tmt + CMPLX(dpdot_re,dpdot_im)*atoms%neq(nat)
436 : END DO
437 :
438 0 : mt_int = mt_int + tmt
439 :
440 0 : END SUBROUTINE dfpt_int_mt
441 :
442 0 : SUBROUTINE dfpt_int_mt_sf(atoms, sphhar, sym, iDir, nat, mt_conj, mt_pure, mt_pure_im, sf_int)
443 : USE m_gaunt, ONLY: gaunt1
444 :
445 : TYPE(t_atoms), INTENT(IN) :: atoms
446 : TYPE(t_sphhar), INTENT(IN) :: sphhar
447 : TYPE(t_sym), INTENT(IN) :: sym
448 :
449 : INTEGER, INTENT(IN) :: iDir, nat
450 :
451 : REAL, INTENT(IN) :: mt_conj(:,0:,:), mt_pure(:,0:,:), mt_pure_im(:,0:,:)
452 :
453 : COMPLEX, INTENT(INOUT) :: sf_int
454 :
455 : INTEGER :: lhPr, lh, lPr, l, iMemPr, iMem, mPr, m, mPrPr, nd
456 : REAL :: f_prod_re, f_prod_im, gaunt_coef
457 : COMPLEX :: f_prod, tmt
458 :
459 : COMPLEX :: trafo_mat(3,-1:1)
460 :
461 0 : nd = sym%ntypsy(nat)
462 :
463 0 : trafo_mat = 0.0
464 :
465 0 : trafo_mat(1,-1) = 1
466 0 : trafo_mat(1, 1) = -1
467 0 : trafo_mat(2,-1) = ImagUnit
468 0 : trafo_mat(2, 1) = ImagUnit
469 0 : trafo_mat(3, 0) = SQRT(2.0)
470 :
471 0 : trafo_mat = SQRT(2*pi_const/3) * trafo_mat
472 :
473 0 : tmt = CMPLX(0.0,0.0)
474 0 : DO lhPr = 0, sphhar%nlh(sym%ntypsy(nat))
475 0 : lPr = sphhar%llh(lhPr,nd)
476 0 : DO lh = 0, sphhar%nlh(sym%ntypsy(nat))
477 0 : l = sphhar%llh(lh,nd)
478 0 : f_prod_re = mt_pure(atoms%jri(nat),lh,nat) * mt_conj(atoms%jri(nat),lhPr,nat)
479 0 : f_prod_im = mt_pure_im(atoms%jri(nat),lh,nat) * mt_conj(atoms%jri(nat),lhPr,nat)
480 0 : f_prod = CMPLX(f_prod_re,f_prod_im)
481 0 : DO iMemPr = 1, sphhar%nmem(lhPr,nd)
482 0 : mPr = sphhar%mlh(iMemPr,lhPr,nd)
483 0 : DO iMem = 1, sphhar%nmem(lh,nd)
484 0 : m = sphhar%mlh(iMem,lh,nd)
485 0 : IF (lPr<ABS(l-1).OR.lPr>(l+1)) CYCLE ! |l-1| <= l' <= l+1
486 0 : IF (MOD(lPr+l,2)/=1) CYCLE ! l' + l + 1 even
487 0 : mPrPr = mPr - m ! m' = m'' + m
488 0 : IF (ABS(mPrPr)>1) CYCLE ! |m''| <= 1
489 0 : gaunt_coef = gaunt1(lPr,1,l,mPr,mPrPr,m,atoms%lmaxd)
490 : tmt = tmt + f_prod * CONJG(sphhar%clnu(iMemPr,lhPr,nd)) * sphhar%clnu(iMem,lh,nd) &
491 0 : * gaunt_coef * trafo_mat(iDir,mPrPr)
492 : END DO
493 : END DO
494 : END DO
495 : END DO
496 :
497 0 : sf_int = sf_int + tmt
498 :
499 0 : END SUBROUTINE dfpt_int_mt_sf
500 :
501 0 : SUBROUTINE dfpt_int_vac(stars,vacuum,cell,vac_conj,vac_pure,vac_int)
502 : USE m_types
503 : USE m_constants
504 : USE m_intgr, ONLY : intgz0
505 :
506 : IMPLICIT NONE
507 :
508 : TYPE(t_stars),INTENT(IN) :: stars
509 : TYPE(t_vacuum),INTENT(IN) :: vacuum
510 : TYPE(t_cell),INTENT(IN) :: cell
511 :
512 : COMPLEX, INTENT(IN) :: vac_conj(:,:,:), vac_pure(:,:,:)
513 :
514 : COMPLEX, INTENT(INOUT) :: vac_int
515 :
516 : REAL :: facv,tvacre,tvacim
517 : COMPLEX :: tvact
518 : INTEGER :: ip,ivac,k2,npz
519 : LOGICAL :: tail
520 : COMPLEX :: dpzc
521 :
522 0 : REAL :: dpzre(vacuum%nmzd), dpzim(vacuum%nmzd)
523 :
524 0 : npz = vacuum%nmz + 1
525 0 : tail = .TRUE.
526 0 : facv=2.0/vacuum%nvac
527 0 : tvacre = 0.
528 0 : tvacim = 0.
529 0 : tvact = CMPLX(0.0,0.0)
530 :
531 0 : dpzre=0.0
532 0 : dpzim=0.0
533 0 : DO ivac = 1,vacuum%nvac
534 0 : DO ip = 1,vacuum%nmz
535 0 : dpzre(npz-ip) = REAL(vac_conj(ip,1,ivac)) * REAL(vac_pure(ip,1,ivac)) + AIMAG(vac_conj(ip,1,ivac)) * AIMAG(vac_pure(ip,1,ivac))
536 0 : dpzim(npz-ip) = REAL(vac_conj(ip,1,ivac)) * AIMAG(vac_pure(ip,1,ivac)) - AIMAG(vac_conj(ip,1,ivac)) * REAL(vac_pure(ip,1,ivac))
537 : END DO
538 0 : DO k2 = 2,stars%ng2
539 0 : DO ip = 1,vacuum%nmzxy
540 0 : dpzc = stars%nstr2(k2) * CONJG(vac_conj(ip,k2,ivac)) * vac_pure(ip,k2,ivac)
541 0 : dpzre(npz-ip) = dpzre(npz-ip) + REAL(dpzc)
542 0 : dpzim(npz-ip) = dpzim(npz-ip) + AIMAG(dpzc)
543 : END DO
544 : END DO
545 0 : CALL intgz0(dpzre,vacuum%delz,vacuum%nmz,tvacre,tail)
546 0 : CALL intgz0(dpzim,vacuum%delz,vacuum%nmz,tvacim,tail)
547 0 : tvact = tvact + cell%area*tvacre*facv
548 0 : tvact = tvact + cell%area*tvacim*facv*ImagUnit
549 : END DO
550 :
551 0 : vac_int = vac_int + tvact
552 :
553 0 : END SUBROUTINE dfpt_int_vac
554 :
555 0 : SUBROUTINE dfpt_dynmat_eigen(fi, results, results1, fmpi, enpara, nococonv, &
556 : stars, starsq, sphhar, inden, hub1data, vx, v, v1real, v1imag, &
557 : eig_id, dfpt_eig_id, dfpt_eig_id2, iDir_col, iDtype_col, iDir_row, iDtype_row, &
558 0 : theta1_pw0, theta1_pw, bqpt, l_real, eigen_term, killcont, q_eig_id)
559 :
560 : USE m_types
561 : USE m_constants
562 : USE m_eigen_hssetup
563 : USE m_pot_io
564 : USE m_eigen_diag
565 : USE m_local_hamiltonian
566 : USE m_util
567 : USE m_eig66_io, ONLY : write_eig, read_eig
568 : USE m_xmlOutput
569 : USE m_types_mpimat
570 : USE m_dfpt_tlmplm
571 :
572 : ! TODO: One bright day, these things will also be relevant for DFPT.
573 : ! We cannot keep doing small systems on small CPUs forever.
574 : !#ifdef _OPENACC
575 : ! USE cublas
576 : !#define CPP_zgemv cublaszgemv
577 : !#else
578 : #define CPP_zgemv zgemv
579 : !#endif
580 :
581 : IMPLICIT NONE
582 :
583 : type(t_fleurinput), intent(in) :: fi
584 : TYPE(t_results),INTENT(INOUT):: results, results1
585 : TYPE(t_mpi),INTENT(IN) :: fmpi
586 : TYPE(t_enpara),INTENT(INOUT) :: enpara
587 : TYPE(t_nococonv),INTENT(IN) :: nococonv
588 : TYPE(t_stars),INTENT(IN) :: stars, starsq
589 : TYPE(t_sphhar),INTENT(IN) :: sphhar
590 : TYPE(t_potden),INTENT(IN) :: inden !
591 : TYPE(t_hub1data),INTENT(INOUT):: hub1data
592 : TYPE(t_potden), INTENT(IN) :: vx
593 : TYPE(t_potden),INTENT(IN) :: v, v1real, v1imag
594 :
595 : ! Scalar Arguments
596 : INTEGER, INTENT(IN) :: eig_id, dfpt_eig_id, iDir_col, iDtype_col, iDir_row, iDtype_row, dfpt_eig_id2
597 : COMPLEX, INTENT(IN) :: theta1_pw0(:), theta1_pw(:)
598 :
599 : REAL, INTENT(IN) :: bqpt(3)
600 : LOGICAL, INTENT(IN) :: l_real
601 :
602 : COMPLEX, INTENT(INOUT) :: eigen_term
603 :
604 : INTEGER, OPTIONAL, INTENT(IN) :: killcont(11), q_eig_id
605 :
606 : ! Local Scalars
607 : INTEGER jsp,nk
608 : INTEGER nk_i
609 : INTEGER err
610 : REAL :: q_loop(3)
611 : ! Local Arrays
612 : INTEGER :: ierr, nbands, nbands1, nbasfcn, nbasfcnq, noccbd, nbandsq, nbands2
613 :
614 0 : REAL, ALLOCATABLE :: bkpt(:)
615 0 : REAL, ALLOCATABLE :: eig(:), eig1(:), we(:), we1(:)
616 :
617 0 : TYPE(t_tlmplm) :: tdV1, tdmod, td
618 0 : TYPE(t_usdus) :: ud, uddummy
619 0 : TYPE(t_lapw) :: lapw, lapwq
620 0 : TYPE(t_hub1data) :: hub1datadummy
621 0 : TYPE (t_mat) :: zMat, zMat1, zMatq, zMat2
622 0 : CLASS(t_mat), ALLOCATABLE :: hmat1,smat1,hmat1q,smat1q,hmat2,smat2,vmat2
623 :
624 : ! Variables for HF or fi%hybinp functional calculation
625 : INTEGER :: comm(fi%kpts%nkpt), dealloc_stat
626 : character(len=300) :: errmsg
627 :
628 : INTEGER :: iEig
629 : INTEGER :: imlo, ilo, iklo, l, ikg, ikglo !!! Test!
630 : REAL :: gext(3)
631 : COMPLEX :: we_loop, we1_loop, eig_loop, eig1_loop, eigen_s1q, eigen_h1qs1q, eigen_bonus
632 : COMPLEX :: eigen_e1, eigen_w1e1, eigen_s2w1e1, eigen_h2s2w1e1, eigen_w1h2s2w1e1
633 :
634 0 : COMPLEX, ALLOCATABLE :: tempVec(:), tempVecq(:), z_loop(:), z1_loop(:), ztest_loop(:), zq_loop(:)
635 :
636 : COMPLEX zdotc
637 : EXTERNAL zdotc
638 :
639 0 : eigen_e1 = CMPLX(0.0,0.0)
640 0 : eigen_w1e1 = CMPLX(0.0,0.0)
641 0 : eigen_s2w1e1 = CMPLX(0.0,0.0)
642 0 : eigen_h2s2w1e1 = CMPLX(0.0,0.0)
643 0 : eigen_w1h2s2w1e1 = CMPLX(0.0,0.0)
644 0 : eigen_s1q = CMPLX(0.0,0.0)
645 0 : eigen_h1qs1q = CMPLX(0.0,0.0)
646 0 : eigen_bonus = CMPLX(0.0,0.0)
647 :
648 : ! Modify this from kpts only in DFPT case.
649 0 : ALLOCATE(bkpt(3))
650 :
651 0 : call ud%init(fi%atoms,fi%input%jspins)
652 0 : call uddummy%init(fi%atoms,fi%input%jspins)
653 0 : ALLOCATE(eig(fi%input%neig))
654 :
655 0 : CALL local_ham(sphhar,fi%atoms,fi%sym,fi%noco,nococonv,enpara,fmpi,v,vx,inden,fi%input,fi%hub1inp,hub1data,td,ud,0.0)
656 : ! Get matrix elements of perturbed potential and modified H/S in DFPT case.
657 0 : hub1datadummy = hub1data
658 :
659 0 : CALL dfpt_tlmplm(fi%atoms,fi%sym,sphhar,fi%input,fi%noco,enpara,fi%hub1inp,hub1data,v,fmpi,tdV1,v1real,v1imag,.FALSE.,iDtype_col)
660 :
661 0 : CALL local_ham(sphhar,fi%atoms,fi%sym,fi%noco,nococonv,enpara,fmpi,v,vx,inden,fi%input,fi%hub1inp,hub1datadummy,tdmod,uddummy,0.0,.true.)
662 :
663 0 : DO jsp = MERGE(1,1,fi%noco%l_noco), MERGE(1,fi%input%jspins,fi%noco%l_noco)
664 0 : k_loop:DO nk_i = 1,size(fmpi%k_list)
665 0 : nk = fmpi%k_list(nk_i)
666 0 : bkpt = fi%kpts%bk(:, nk)
667 0 : q_loop = bqpt
668 :
669 0 : CALL lapw%init(fi%input,fi%noco,nococonv,fi%kpts,fi%atoms,fi%sym,nk,fi%cell,fmpi)
670 0 : CALL lapwq%init(fi%input,fi%noco,nococonv,fi%kpts,fi%atoms,fi%sym,nk,fi%cell,fmpi,q_loop)
671 :
672 0 : we = results%w_iks(:,nk,jsp)
673 0 : we1 = results1%w_iks(:,nk,jsp)
674 0 : eig = results%eig(:,nk,jsp)
675 0 : eig1 = results1%eig(:,nk,jsp)
676 :
677 0 : noccbd = COUNT(we>1.e-8)
678 :
679 0 : nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*fi%atoms%nlotot,lapw%nv(1)+fi%atoms%nlotot,fi%noco%l_noco)
680 0 : nbasfcnq = MERGE(lapwq%nv(1)+lapwq%nv(2)+2*fi%atoms%nlotot,lapwq%nv(1)+fi%atoms%nlotot,fi%noco%l_noco)
681 :
682 0 : CALL zMat%init(l_real,nbasfcn,noccbd)
683 0 : CALL zMat1%init(.FALSE.,nbasfcnq,noccbd)
684 0 : CALL zMat2%init(.FALSE.,nbasfcnq,noccbd)
685 0 : IF (PRESENT(q_eig_id)) CALL zMatq%init(l_real,nbasfcnq,noccbd)
686 :
687 0 : ALLOCATE(tempVec(nbasfcn))
688 0 : IF (PRESENT(q_eig_id)) ALLOCATE(tempVecq(nbasfcnq))
689 0 : ALLOCATE(z_loop(nbasfcn),z1_loop(nbasfcnq))
690 0 : ALLOCATE(ztest_loop(nbasfcn))
691 0 : IF (PRESENT(q_eig_id)) ALLOCATE(zq_loop(nbasfcnq))
692 :
693 0 : CALL read_eig(eig_id, nk,jsp,neig=nbands,zmat=zMat)
694 0 : CALL read_eig(dfpt_eig_id, nk,jsp,neig=nbands1,zmat=zMat1)
695 0 : CALL read_eig(dfpt_eig_id2,nk,jsp,neig=nbands2,zmat=zMat2)
696 0 : IF (PRESENT(q_eig_id)) CALL read_eig(q_eig_id, nk,jsp,neig=nbandsq,zmat=zMatq)
697 :
698 : !!!! Test:
699 : !DO ikG = 1, lapw%nv(jsp)
700 : ! gExt = MATMUL(fi%cell%bmat,lapw%vk(:, ikG, jsp))
701 : ! IF (zMat%l_real) THEN
702 : ! zMat1%data_c(ikG,:) = -ImagUnit * gExt(idir_row) * zMat%data_r(ikG, :noccbd)
703 : ! ELSE
704 : ! zMat1%data_c(ikG,:) = -ImagUnit * gExt(idir_row) * zMat%data_c(ikG, :noccbd)
705 : ! END IF
706 : !END DO
707 :
708 : !DO ikG = lapw%nv(jsp) + 1, lapw%nv(jsp) + fi%atoms%nlo(iDtype_row)
709 : ! iLo = ikG-lapw%nv(jsp)
710 : ! l = fi%atoms%llo(iLo, iDtype_row)
711 : ! DO imLo = 1, 2*l+1
712 : ! ikLo = lapw%kvec(imLo,iLo,iDtype_row)
713 : ! ikGLo = lapw%nv(jsp) + lapw%index_lo(iLo,iDtype_row) + imLo
714 : ! gExt = MATMUL(fi%cell%bmat,lapw%vk(:,ikLo, jsp))
715 : ! IF (zMat%l_real) THEN
716 : ! zMat1%data_c(ikGLo,:) = -ImagUnit * gExt(idir_row) * zMat%data_r(ikGLo, :)
717 : ! ELSE
718 : ! zMat1%data_c(ikGLo,:) = -ImagUnit * gExt(idir_row) * zMat%data_c(ikGLo, :)
719 : ! END IF
720 : ! END DO
721 : ! gExt = MATMUL(fi%cell%bmat,lapw%vk(:, ikG, jsp))
722 : !END DO
723 :
724 0 : CALL timestart("Setup of H&S matrices")
725 0 : IF (.NOT.PRESENT(q_eig_id)) THEN
726 : CALL dfpt_dynmat_hssetup(jsp, fmpi, fi, enpara, nococonv, starsq, stars, &
727 : ud, tdmod, tdV1, lapw, lapwq, iDir_row, iDtype_row, iDir_col, iDtype_col, theta1_pw0, theta1_pw, &
728 0 : smat1, hmat1, smat1q, hmat1q, smat2, hmat2, nk, killcont)
729 : ELSE
730 : CALL dfpt_dynmat_hssetup(jsp, fmpi, fi, enpara, nococonv, starsq, stars, &
731 : ud, tdmod, tdV1, lapw, lapwq, iDir_row, iDtype_row, iDir_col, iDtype_col, theta1_pw0, theta1_pw, &
732 0 : smat1, hmat1, smat1q, hmat1q, smat2, hmat2, nk, killcont, vmat2)
733 : END IF
734 0 : CALL timestop("Setup of H&S matrices")
735 :
736 0 : DO iEig = 1, noccbd
737 0 : eig_loop = eig(iEig)
738 0 : eig1_loop = eig1(iEig)
739 0 : we_loop = (2.0/fi%input%jspins)*we(iEig)
740 0 : we1_loop = (2.0/fi%input%jspins)*we1(iEig)
741 0 : IF (l_real) THEN
742 0 : z_loop = CMPLX(1.0,0.0)*zMat%data_r(:,iEig)
743 0 : IF (PRESENT(q_eig_id)) zq_loop = CMPLX(1.0,0.0)*zMatq%data_r(:,iEig)
744 : ELSE
745 0 : z_loop = zMat%data_c(:,iEig)
746 0 : IF (PRESENT(q_eig_id)) zq_loop = zMatq%data_c(:,iEig)
747 : END IF
748 :
749 0 : z1_loop = zMat1%data_c(:,iEig)
750 :
751 0 : CALL CPP_zgemv('N',nbasfcn,nbasfcn,-we_loop*eig1_loop,smat1%data_c,nbasfcn,z_loop,1,CMPLX(0.0,0.0),tempVec,1)
752 0 : eigen_e1 = eigen_e1 + zdotc(nbasfcn,z_loop,1,tempVec,1)
753 0 : CALL CPP_zgemv('N',nbasfcn,nbasfcn,-we1_loop*eig_loop,smat1%data_c,nbasfcn,z_loop,1,CMPLX(1.0,0.0),tempVec,1)
754 0 : eigen_w1e1 = eigen_w1e1 + zdotc(nbasfcn,z_loop,1,tempVec,1)
755 0 : CALL CPP_zgemv('N',nbasfcn,nbasfcn,-we_loop*eig_loop,smat2%data_c,nbasfcn,z_loop,1,CMPLX(1.0,0.0),tempVec,1)
756 0 : eigen_s2w1e1 = eigen_s2w1e1 + zdotc(nbasfcn,z_loop,1,tempVec,1)
757 0 : CALL CPP_zgemv('N',nbasfcn,nbasfcn,we_loop,hmat2%data_c,nbasfcn,z_loop,1,CMPLX(1.0,0.0),tempVec,1)
758 0 : eigen_h2s2w1e1 = eigen_h2s2w1e1 + zdotc(nbasfcn,z_loop,1,tempVec,1)
759 0 : CALL CPP_zgemv('N',nbasfcn,nbasfcn,we1_loop,hmat1%data_c,nbasfcn,z_loop,1,CMPLX(1.0,0.0),tempVec,1)
760 0 : eigen_w1h2s2w1e1 = eigen_w1h2s2w1e1 + zdotc(nbasfcn,z_loop,1,tempVec,1)
761 :
762 0 : eigen_term = eigen_term + zdotc(nbasfcn,z_loop,1,tempVec,1)
763 :
764 0 : IF (PRESENT(q_eig_id)) THEN
765 0 : CALL CPP_zgemv('N',nbasfcnq,nbasfcn,we_loop,vmat2%data_c,nbasfcnq,z_loop,1,CMPLX(0.0,0.0),tempVecq,1)
766 0 : eigen_term = eigen_term + zdotc(nbasfcnq,zq_loop,1,tempVecq,1)
767 : END IF
768 :
769 0 : CALL CPP_zgemv('C',nbasfcnq,nbasfcn,-we_loop*eig_loop,smat1q%data_c,nbasfcnq,z1_loop,1,CMPLX(0.0,0.0),tempVec,1)
770 0 : eigen_s1q = eigen_s1q + 2.0*zdotc(nbasfcn,z_loop,1,tempVec,1)
771 0 : CALL CPP_zgemv('C',nbasfcnq,nbasfcn,we_loop,hmat1q%data_c,nbasfcnq,z1_loop,1,CMPLX(1.0,0.0),tempVec,1)
772 0 : eigen_h1qs1q = eigen_h1qs1q + 2.0*zdotc(nbasfcn,z_loop,1,tempVec,1)
773 :
774 0 : eigen_term = eigen_term + 2.0*zdotc(nbasfcn,z_loop,1,tempVec,1)
775 :
776 : ! Correction:
777 0 : CALL CPP_zgemv('C',nbasfcnq,nbasfcn,we_loop,smat1q%data_c,nbasfcnq,zMat2%data_c(:,iEig),1,CMPLX(0.0,0.0),tempVec,1)
778 0 : eigen_bonus = eigen_bonus + 2.0*zdotc(nbasfcn,z_loop,1,tempVec,1)
779 0 : eigen_term = eigen_term + 2.0 * zdotc(nbasfcn,z_loop,1,tempVec,1)
780 : END DO
781 :
782 0 : DEALLOCATE(tempVec)
783 0 : IF (PRESENT(q_eig_id)) DEALLOCATE(tempVecq)
784 0 : DEALLOCATE(z_loop,z1_loop)
785 0 : IF (PRESENT(q_eig_id)) DEALLOCATE(zq_loop)
786 0 : DEALLOCATE(ztest_loop)
787 0 : CALL smat1%free()
788 0 : CALL hmat1%free()
789 0 : CALL smat1q%free()
790 0 : CALL hmat1q%free()
791 0 : CALL smat2%free()
792 0 : CALL hmat2%free()
793 0 : IF (PRESENT(q_eig_id)) CALL vmat2%free()
794 0 : DEALLOCATE(hmat1,smat1,hmat1q,smat1q,hmat2,smat2, stat=dealloc_stat, errmsg=errmsg)
795 0 : IF (PRESENT(q_eig_id)) DEALLOCATE(vmat2)
796 0 : if(dealloc_stat /= 0) call juDFT_error("deallocate failed one of the matrices",&
797 0 : hint=errmsg, calledby="dfpt_dynmat.F90")
798 :
799 : ! Output results
800 0 : CALL timestart("EV output")
801 :
802 : #if defined(CPP_MPI)
803 : ! RMA synchronization
804 0 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
805 : #endif
806 0 : CALL timestop("EV output")
807 :
808 : !IF (allocated(zmat)) THEN
809 0 : CALL zMat%free()
810 0 : CALL zMat1%free()
811 0 : CALL zMat2%free()
812 : !deallocate(zMat)
813 : !ENDIF
814 : END DO k_loop
815 : END DO ! spin loop ends
816 : #ifdef CPP_MPI
817 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_term,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
818 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_e1,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
819 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_w1e1,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
820 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_s2w1e1,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
821 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_h2s2w1e1,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
822 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_w1h2s2w1e1,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
823 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_s1q,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
824 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_h1qs1q,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
825 0 : CALL MPI_ALLREDUCE(MPI_IN_PLACE,eigen_bonus,1,MPI_DOUBLE_COMPLEX,MPI_SUM,fmpi%mpi_comm,ierr)
826 : #endif
827 0 : IF (fmpi%irank==0) THEN
828 0 : write(9989,*) "eigen -w0e1s1 :", eigen_e1
829 0 : write(9989,*) "eigen -w1e0s1 :", eigen_w1e1-eigen_e1
830 0 : write(9989,*) "eigen -w0e0s2 :", eigen_s2w1e1-eigen_w1e1
831 0 : write(9989,*) "eigen w0h2 :", eigen_h2s2w1e1-eigen_s2w1e1
832 0 : write(9989,*) "eigen w1h1 :", eigen_w1h2s2w1e1-eigen_h2s2w1e1
833 0 : write(9989,*) "eigen -w0e0s1q:", eigen_s1q
834 0 : write(9989,*) "eigen w0h1q :", eigen_h1qs1q - eigen_s1q
835 0 : write(9989,*) "eigen corr :", eigen_bonus
836 : END IF
837 0 : END SUBROUTINE
838 :
839 0 : SUBROUTINE dfpt_dynmat_hssetup(isp, fmpi, fi, enpara, nococonv, starsq, stars, &
840 0 : ud, td, tdV1, lapw, lapwq, iDir_row, iDtype_row, iDir_col, iDtype_col, theta1_pw0, theta1_pw, &
841 : smat1_final, hmat1_final, smat1q_final, hmat1q_final, smat2_final, hmat2_final, nk, killcont, vmat2_final)
842 : USE m_types
843 : USE m_types_mpimat
844 : USE m_dfpt_hs_int
845 : USE m_dfpt_hsmt
846 : USE m_dfpt_eigen_redist_matrix
847 :
848 : IMPLICIT NONE
849 :
850 : INTEGER, INTENT(IN) :: isp
851 : TYPE(t_mpi), INTENT(IN) :: fmpi
852 : type(t_fleurinput), INTENT(IN) :: fi
853 : TYPE(t_stars), INTENT(IN) :: starsq, stars
854 : TYPE(t_enpara), INTENT(IN) :: enpara
855 : TYPE(t_nococonv), INTENT(IN) :: nococonv
856 : TYPE(t_usdus), INTENT(IN) :: ud
857 : TYPE(t_tlmplm), INTENT(IN) :: td, tdV1
858 : TYPE(t_lapw), INTENT(IN) :: lapw, lapwq
859 : INTEGER, INTENT(IN) :: iDir_row, iDtype_row, iDir_col, iDtype_col
860 : COMPLEX, INTENT(IN) :: theta1_pw0(:), theta1_pw(:)
861 : CLASS(t_mat), ALLOCATABLE, INTENT(INOUT) :: smat1_final, hmat1_final, smat1q_final, hmat1q_final, smat2_final, hmat2_final
862 : CLASS(t_mat), OPTIONAL, ALLOCATABLE, INTENT(INOUT) :: vmat2_final
863 : INTEGER, INTENT(IN) :: nk, killcont(11)
864 :
865 0 : CLASS(t_mat), ALLOCATABLE :: smat1(:, :), hmat1(:, :), smat1q(:, :), hmat1q(:, :), smat2(:, :), hmat2(:, :), vmat2(:,:)
866 :
867 : INTEGER :: i, j, nspins
868 :
869 0 : nspins = MERGE(2, 1, fi%noco%l_noco)
870 0 : IF (fmpi%n_size == 1) THEN
871 0 : ALLOCATE (t_mat::smat1(nspins, nspins), hmat1(nspins, nspins))
872 0 : ALLOCATE (t_mat::smat1q(nspins, nspins), hmat1q(nspins, nspins))
873 0 : ALLOCATE (t_mat::smat2(nspins, nspins), hmat2(nspins, nspins))
874 0 : IF (PRESENT(vmat2_final)) THEN
875 0 : ALLOCATE (t_mat::vmat2(nspins, nspins))
876 : END IF
877 : ELSE
878 0 : ALLOCATE (t_mpimat::smat1(nspins, nspins), hmat1(nspins, nspins))
879 0 : ALLOCATE (t_mpimat::smat1q(nspins, nspins), hmat1q(nspins, nspins))
880 0 : ALLOCATE (t_mpimat::smat2(nspins, nspins), hmat2(nspins, nspins))
881 0 : IF (PRESENT(vmat2_final)) THEN
882 0 : ALLOCATE (t_mpimat::vmat2(nspins, nspins))
883 : END IF
884 : END IF
885 :
886 0 : DO i = 1, nspins
887 0 : DO j = 1, nspins
888 0 : CALL smat1(i, j)%init(.FALSE., lapw%nv(i) + fi%atoms%nlotot, lapw%nv(j) + fi%atoms%nlotot, fmpi%sub_comm, .false.)
889 0 : CALL hmat1(i, j)%init(smat1(i, j))
890 0 : CALL smat1q(i, j)%init(.FALSE., lapwq%nv(i) + fi%atoms%nlotot, lapw%nv(j) + fi%atoms%nlotot, fmpi%sub_comm, .false.)
891 0 : CALL hmat1q(i, j)%init(smat1q(i, j))
892 0 : CALL smat2(i, j)%init(.FALSE., lapw%nv(i) + fi%atoms%nlotot, lapw%nv(j) + fi%atoms%nlotot, fmpi%sub_comm, .false.)
893 0 : CALL hmat2(i, j)%init(smat2(i, j))
894 0 : IF (PRESENT(vmat2_final)) THEN
895 0 : CALL vmat2(i, j)%init(.FALSE., lapwq%nv(i) + fi%atoms%nlotot, lapw%nv(j) + fi%atoms%nlotot, fmpi%sub_comm, .false.)
896 : END IF
897 : END DO
898 : END DO
899 :
900 0 : CALL timestart("Interstitial part")
901 : CALL dfpt_dynmat_hs_int(fi%noco, starsq, stars, lapwq, lapw, fmpi, fi%cell%bbmat, isp, theta1_pw0, theta1_pw, &
902 0 : smat1, hmat1, smat1q, hmat1q, killcont(1:4))
903 0 : CALL timestop("Interstitial part")
904 :
905 0 : CALL timestart("MT part")
906 0 : DO i = 1, nspins; DO j = 1, nspins
907 : !!$acc enter data copyin(hmat(i,j),smat(i,j))
908 : !!$acc enter data copyin(hmat(i,j)%data_r,smat(i,j)%data_r,hmat(i,j)%data_c,smat(i,j)%data_c)
909 : END DO; END DO
910 0 : IF (.NOT.PRESENT(vmat2_final)) THEN
911 : CALL dfpt_dynmat_hsmt(fi%atoms, fi%sym, enpara, isp, iDir_row, iDtype_row, iDir_col, iDtype_col, fi%input, fmpi, fi%noco, nococonv, fi%cell, &
912 0 : lapw, lapwq, ud, td, tdV1, hmat1, smat1, hmat1q, smat1q, hmat2, smat2, nk, killcont(5:11))
913 : ELSE
914 : CALL dfpt_dynmat_hsmt(fi%atoms, fi%sym, enpara, isp, iDir_row, iDtype_row, iDir_col, iDtype_col, fi%input, fmpi, fi%noco, nococonv, fi%cell, &
915 0 : lapw, lapwq, ud, td, tdV1, hmat1, smat1, hmat1q, smat1q, hmat2, smat2, nk, killcont(5:11), vmat2)
916 : END IF
917 0 : DO i = 1, nspins; DO j = 1, nspins; if (hmat1(1, 1)%l_real) THEN
918 : !!$acc exit data copyout(hmat(i,j)%data_r,smat(i,j)%data_r) delete(hmat(i,j)%data_c,smat(i,j)%data_c)
919 : !!$acc exist data delete(hmat(i,j),smat(i,j))
920 : ELSE
921 : !!$acc exit data copyout(hmat(i,j)%data_c,smat(i,j)%data_c) delete(hmat(i,j)%data_r,smat(i,j)%data_r)
922 : !!$acc exit data delete(hmat(i,j),smat(i,j))
923 : END IF; END DO; END DO
924 0 : CALL timestop("MT part")
925 :
926 : !Now copy the data into final matrix
927 : ! Collect the four fi%noco parts into a single matrix
928 : ! In collinear case only a copy is done
929 : ! In the parallel case also a redistribution happens
930 0 : ALLOCATE (smat1_final, mold=smat1(1, 1))
931 0 : ALLOCATE (hmat1_final, mold=smat1(1, 1))
932 0 : ALLOCATE (smat1q_final, mold=smat1q(1, 1))
933 0 : ALLOCATE (hmat1q_final, mold=smat1q(1, 1))
934 0 : ALLOCATE (smat2_final, mold=smat2(1, 1))
935 0 : ALLOCATE (hmat2_final, mold=smat2(1, 1))
936 0 : IF (PRESENT(vmat2_final)) ALLOCATE (vmat2_final, mold=vmat2(1, 1))
937 :
938 0 : CALL timestart("Matrix redistribution")
939 0 : CALL dfpt_eigen_redist_matrix(fmpi, lapw, lapw, fi%atoms, smat1, smat1_final)
940 0 : CALL dfpt_eigen_redist_matrix(fmpi, lapw, lapw, fi%atoms, hmat1, hmat1_final, smat1_final)
941 0 : CALL dfpt_eigen_redist_matrix(fmpi, lapwq, lapw, fi%atoms, smat1q, smat1q_final)
942 0 : CALL dfpt_eigen_redist_matrix(fmpi, lapwq, lapw, fi%atoms, hmat1q, hmat1q_final, smat1q_final)
943 0 : CALL dfpt_eigen_redist_matrix(fmpi, lapw, lapw, fi%atoms, smat2, smat2_final)
944 0 : CALL dfpt_eigen_redist_matrix(fmpi, lapw, lapw, fi%atoms, hmat2, hmat2_final, smat2_final)
945 0 : IF (PRESENT(vmat2_final)) CALL dfpt_eigen_redist_matrix(fmpi, lapw, lapw, fi%atoms, vmat2, vmat2_final)
946 0 : CALL timestop("Matrix redistribution")
947 0 : END SUBROUTINE
948 : END MODULE m_dfpt_dynmat
|
