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
7 : USE m_juDFT
8 : USE m_constants
9 : USE m_types
10 :
11 : IMPLICIT NONE
12 :
13 : CONTAINS
14 0 : SUBROUTINE dfpt(fi, sphhar, stars, nococonv, qpts, fmpi, results, enpara, &
15 : & rho, vTot, vxc, eig_id, xcpot, hybdat, mpdata, forcetheo)
16 : USE m_dfpt_check
17 : !USE m_dfpt_test
18 : USE m_dfpt_sternheimer
19 : USE m_dfpt_dynmat
20 : USE m_dfpt_eii2, only : CalcIIEnerg2, genPertPotDensGvecs, dfpt_e2_madelung
21 : USE m_dfpt_dynmat_eig
22 : USE m_juDFT_stop, only : juDFT_error
23 : USE m_vgen_coulomb
24 : USE m_dfpt_vgen
25 : USE m_fleur_init
26 : USE m_eigen
27 : USE m_fermie
28 : USE m_grdchlh
29 : USE m_dfpt_dynmat_sym
30 : USE m_types_eigdos
31 : USE m_make_dos
32 : USE m_dfpt_gradient
33 : USE m_npy
34 :
35 : TYPE(t_mpi), INTENT(IN) :: fmpi
36 : TYPE(t_fleurinput), INTENT(IN) :: fi
37 : TYPE(t_sphhar), INTENT(IN) :: sphhar
38 : TYPE(t_stars), INTENT(IN) :: stars
39 : TYPE(t_nococonv), INTENT(IN) :: nococonv
40 : TYPE(t_enpara), INTENT(INOUT) :: enpara
41 : TYPE(t_results), INTENT(INOUT) :: results
42 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
43 : TYPE(t_mpdata), INTENT(INOUT) :: mpdata
44 :
45 : CLASS(t_xcpot), INTENT(IN) :: xcpot
46 : CLASS(t_forcetheo), INTENT(INOUT) :: forcetheo
47 :
48 : TYPE(t_kpts), INTENT(IN) :: qpts !Possibly replace this by fi_nosym%kpts [read correctly!]
49 :
50 : TYPE(t_potden), INTENT(INOUT) :: rho
51 : TYPE(t_potden), INTENT(IN) :: vTot, vxc
52 : INTEGER, INTENT(IN) :: eig_id
53 :
54 0 : TYPE(t_hub1data) :: hub1data
55 0 : TYPE(t_potden) :: grvextdummy, imagrhodummy, rho_nosym, vTot_nosym, vext_dummy, vC_dummy
56 0 : TYPE(t_potden) :: grRho3(3), grVtot3(3), grVC3(3), grVext3(3)
57 0 : TYPE(t_potden) :: grgrVC3x3(3,3), grgrvextnum(3,3)
58 0 : TYPE(t_potden) :: denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im, vTot1m, vTot1mIm ! q-quantities
59 0 : TYPE(t_results) :: q_results, results1, qm_results, results1m
60 0 : TYPE(t_kpts) :: qpts_loc
61 0 : TYPE(t_kpts) :: kqpts, kqmpts ! basically kpts, but with q added onto each one.
62 0 : TYPE(t_dos),TARGET :: dos
63 :
64 0 : TYPE(t_eigdos_list),allocatable :: eigdos(:)
65 :
66 : ! Desymmetrized type variables:
67 0 : TYPE(t_mpi) :: fmpi_nosym
68 0 : TYPE(t_fleurinput) :: fi_nosym
69 0 : TYPE(t_sphhar) :: sphhar_nosym
70 0 : TYPE(t_stars) :: stars_nosym, starsq, starsmq
71 0 : TYPE(t_nococonv) :: nococonv_nosym
72 0 : TYPE(t_enpara) :: enpara_nosym
73 0 : TYPE(t_results) :: results_nosym
74 0 : TYPE(t_wann) :: wann_nosym
75 0 : TYPE(t_hybdat) :: hybdat_nosym
76 0 : TYPE(t_mpdata) :: mpdata_nosym
77 :
78 0 : CLASS(t_xcpot), ALLOCATABLE :: xcpot_nosym
79 0 : CLASS(t_forcetheo), ALLOCATABLE :: forcetheo_nosym
80 :
81 : ! Full symmetrized type variables:
82 0 : TYPE(t_mpi) :: fmpi_fullsym
83 0 : TYPE(t_fleurinput) :: fi_fullsym
84 0 : TYPE(t_sphhar) :: sphhar_fullsym
85 0 : TYPE(t_stars) :: stars_fullsym
86 0 : TYPE(t_nococonv) :: nococonv_fullsym
87 0 : TYPE(t_enpara) :: enpara_fullsym
88 0 : TYPE(t_results) :: results_fullsym
89 0 : TYPE(t_wann) :: wann_fullsym
90 0 : TYPE(t_hybdat) :: hybdat_fullsym
91 0 : TYPE(t_mpdata) :: mpdata_fullsym
92 :
93 0 : CLASS(t_xcpot), ALLOCATABLE :: xcpot_fullsym
94 0 : CLASS(t_forcetheo), ALLOCATABLE :: forcetheo_fullsym
95 :
96 0 : INTEGER, ALLOCATABLE :: recG(:, :)
97 : INTEGER :: ngdp2km
98 0 : complex, allocatable :: E2ndOrdII(:, :)
99 0 : complex, allocatable :: eigenFreqs(:)
100 0 : real, allocatable :: eigenVals(:), eigenValsFull(:,:,:)
101 0 : complex, allocatable :: eigenVecs(:, :)
102 :
103 0 : COMPLEX, ALLOCATABLE :: grrhodummy(:, :, :, :, :)
104 :
105 0 : COMPLEX, ALLOCATABLE :: dyn_mat(:,:,:), dyn_mat_r(:,:,:), dyn_mat_q_full(:,:,:), dyn_mat_pathq(:,:), sym_dynvec(:,:,:), sym_dyn_mat(:,:,:)
106 0 : REAL, ALLOCATABLE :: e2_vm(:,:,:)
107 :
108 : INTEGER :: ngdp, iSpin, iQ, iDir, iDtype, nspins, zlim, iVac, lh, iDir2, sym_count
109 : INTEGER :: iStar, xInd, yInd, zInd, q_eig_id, ikpt, ierr, qm_eig_id, iArray
110 : INTEGER :: dfpt_eig_id, dfpt_eig_id2, dfpt_eigm_id, dfpt_eigm_id2
111 : LOGICAL :: l_real, l_minusq, l_dfpt_scf, l_cheated
112 :
113 : LOGICAL :: l_dfpt_band, l_dfpt_dos, l_dfpt_full
114 :
115 : CHARACTER(len=20) :: dfpt_tag
116 : CHARACTER(len=4) :: dynfiletag
117 : CHARACTER(len=100) :: inp_pref, trash
118 :
119 0 : INTEGER, ALLOCATABLE :: q_list(:)
120 0 : INTEGER, ALLOCATABLE :: sym_list(:) ! For each q: Collect, which symmetries leave q unchanged.
121 :
122 : ! Desym-tests:
123 : INTEGER :: grid(3), iread
124 0 : REAL :: dr_re(fi%vacuum%nmzd), dr_im(fi%vacuum%nmzd), drr_dummy(fi%vacuum%nmzd), numbers(3*fi%atoms%nat,6*fi%atoms%nat)
125 : complex :: sigma_loc(2), sigma_ext(2), sigma_coul(2), sigma_gext(3,2)
126 :
127 0 : ALLOCATE(e2_vm(fi%atoms%nat,3,3))
128 :
129 0 : l_dfpt_scf = fi%juPhon%l_scf
130 :
131 0 : l_dfpt_band = fi%juPhon%l_band
132 0 : l_dfpt_full = fi%juPhon%l_intp
133 0 : l_dfpt_dos = fi%juPhon%l_dos
134 :
135 0 : l_cheated = .FALSE.
136 :
137 0 : l_real = fi%sym%invs.AND.(.NOT.fi%noco%l_soc).AND.(.NOT.fi%noco%l_noco).AND.fi%atoms%n_hia==0
138 :
139 : ! l_minusq is a hard false at the moment. It can be used to ignore +-q symmetries and
140 : ! run the Sternheimer loop etc etc for both +q and -q.
141 : ! This was only used for testing but may become relevant again for SOC systems with broken
142 : ! inversion symmetry
143 0 : l_minusq = .FALSE.
144 :
145 : nspins = MERGE(2, 1, fi%noco%l_noco)
146 :
147 0 : IF (fi%juPhon%l_jpCheck) THEN
148 : ! This function will be used to check the validity of juPhon's
149 : ! input. I.e. check, whether all prohibited switches are off and,
150 : ! once there is more expertise considering this topic, check whether
151 : ! the cutoffs are chosen appropriately.
152 0 : CALL dfpt_check(fi_nosym, xcpot_nosym)
153 : END IF
154 :
155 0 : IF (fi%sym%nop>1) THEN
156 0 : WRITE(*,*) "Desymmetrization needed. Going ahead!"
157 : ! Grid size for desym quality test:
158 0 : grid = 21
159 :
160 0 : inp_pref = ADJUSTL("desym_")
161 0 : fmpi_nosym%l_mpi_multithreaded = fmpi%l_mpi_multithreaded
162 0 : fmpi_nosym%mpi_comm = fmpi%mpi_comm
163 :
164 : CALL dfpt_desym(fmpi_nosym,fi_nosym,sphhar_nosym,stars_nosym,nococonv_nosym,enpara_nosym,results_nosym,wann_nosym,hybdat_nosym,mpdata_nosym,xcpot_nosym,forcetheo_nosym,rho_nosym,vTot_nosym,grid,inp_pref,&
165 0 : fi,sphhar,stars,nococonv,enpara,results,rho,vTot)
166 : ELSE
167 0 : fmpi_nosym = fmpi
168 0 : fi_nosym = fi
169 0 : sphhar_nosym = sphhar
170 0 : stars_nosym = stars
171 0 : nococonv_nosym = nococonv
172 0 : forcetheo_nosym = forcetheo
173 0 : enpara_nosym = enpara
174 0 : xcpot_nosym = xcpot
175 0 : results_nosym = results
176 0 : hybdat_nosym = hybdat
177 0 : mpdata_nosym = mpdata
178 0 : rho_nosym = rho
179 0 : vTot_nosym = vTot
180 : END IF
181 :
182 : ! TODO: Maybe rather replace this by switches filtered into dfpt_routines.
183 : ! IF (fi%juPhon%l_jpTest) THEN
184 : ! This function will be used to run (parts of) the test suite for
185 : ! OG juPhon, as provided by CRG.
186 : !CALL dfpt_test(fi, sphhar, stars, fmpi, rho, grRho, rho0, grRho0, xcpot, ngdp, recG, grVxcIRKern, ylm, dKernMTGPts, gausWts, hybdat)
187 : !END IF
188 :
189 : ! q_results saves the eigen-info for k+q, results1 for the perturbed quantities
190 0 : CALL q_results%init(fi%input, fi%atoms, fi%kpts, fi%noco)
191 0 : CALL results1%init(fi_nosym%input, fi_nosym%atoms, fi_nosym%kpts, fi_nosym%noco)
192 :
193 : IF (l_minusq) THEN
194 : CALL qm_results%init(fi%input, fi%atoms, fi%kpts, fi%noco)
195 : CALL results1m%init(fi_nosym%input, fi_nosym%atoms, fi_nosym%kpts, fi_nosym%noco)
196 : END IF
197 :
198 0 : IF (.NOT.fi%juPhon%qmode==0) THEN
199 : ! Read qpoints from fullsym_inp.xml and fullsym_kpts.xml
200 0 : inp_pref = ADJUSTL("fullsym_")
201 0 : fmpi_fullsym%l_mpi_multithreaded = fmpi%l_mpi_multithreaded
202 0 : fmpi_fullsym%mpi_comm = fmpi%mpi_comm
203 : CALL fleur_init(fmpi_fullsym, fi_fullsym, sphhar_fullsym, stars_fullsym, nococonv_fullsym, forcetheo_fullsym, &
204 : enpara_fullsym, xcpot_fullsym, results_fullsym, wann_fullsym, hybdat_fullsym, mpdata_fullsym, &
205 0 : inp_pref)
206 0 : qpts_loc = fi_fullsym%kpts
207 :
208 0 : ALLOCATE(q_list(SIZE(qpts_loc%bk,2)))
209 0 : q_list = (/(iArray, iArray=1,SIZE(qpts_loc%bk,2), 1)/)
210 :
211 0 : ALLOCATE(sym_list(fi_fullsym%sym%nop))
212 0 : sym_list = 0
213 : ELSE
214 : ! Read qpoints from the juPhon qlist in inp.xml
215 0 : qpts_loc = qpts
216 0 : qpts_loc%bk(:, :SIZE(fi%juPhon%qvec,2)) = fi%juPhon%qvec
217 :
218 0 : ALLOCATE(q_list(SIZE(fi%juPhon%qvec,2)))
219 0 : q_list = (/(iArray, iArray=1,SIZE(fi%juPhon%qvec,2), 1)/)
220 : END IF
221 :
222 : ! Generate the gradients of the density and the various potentials, that will be used at different points in the programm.
223 : ! The density gradient is calculated by numerical differentiation, while the potential gradients are constructed (from the
224 : ! density gradient) by a Weinert construction, just like the potentials are from the density.
225 : ! This is done to ensure good continuity.
226 0 : CALL timestart("Gradient generation")
227 0 : ALLOCATE(grrhodummy(fi_nosym%atoms%jmtd, (fi_nosym%atoms%lmaxd+1)**2, fi_nosym%atoms%nat, SIZE(rho_nosym%mt,4), 3))
228 :
229 0 : CALL imagrhodummy%copyPotDen(rho_nosym)
230 0 : CALL imagrhodummy%resetPotDen()
231 0 : CALL grvextdummy%copyPotDen(rho_nosym)
232 :
233 0 : CALL vext_dummy%copyPotDen(vTot_nosym)
234 0 : CALL vext_dummy%resetPotDen()
235 0 : CALL vC_dummy%copyPotDen(vTot_nosym)
236 0 : CALL vC_dummy%resetPotDen()
237 0 : sigma_loc = cmplx(0.0,0.0)
238 0 : sigma_ext = cmplx(0.0,0.0)
239 0 : sigma_coul = cmplx(0.0,0.0)
240 0 : sigma_gext = cmplx(0.0,0.0)
241 : CALL vgen_coulomb(1, fmpi_nosym, fi_nosym%input, fi_nosym%field, fi_nosym%vacuum, fi_nosym%sym, fi%juphon, stars_nosym, fi_nosym%cell, &
242 0 : & sphhar_nosym, fi_nosym%atoms, .FALSE., imagrhodummy, vext_dummy, sigma_ext)
243 : CALL vgen_coulomb(1, fmpi_nosym, fi_nosym%input, fi_nosym%field, fi_nosym%vacuum, fi_nosym%sym, fi%juphon, stars_nosym, fi_nosym%cell, &
244 0 : & sphhar_nosym, fi_nosym%atoms, .FALSE., rho_nosym, vC_dummy, sigma_coul)
245 0 : DO iDir = 1, 3
246 0 : CALL grRho3(iDir)%copyPotDen(rho_nosym)
247 0 : CALL grRho3(iDir)%resetPotDen()
248 0 : DO iDir2 = 1, 3
249 0 : CALL grgrvextnum(iDir2,iDir)%copyPotDen(vTot_nosym)
250 0 : CALL grgrvextnum(iDir2,iDir)%resetPotDen()
251 0 : CALL grgrVC3x3(iDir2,iDir)%copyPotDen(vTot_nosym)
252 0 : CALL grgrVC3x3(iDir2,iDir)%resetPotDen()
253 : END DO
254 0 : CALL grVext3(iDir)%copyPotDen(vTot_nosym)
255 0 : CALL grVext3(iDir)%resetPotDen()
256 0 : CALL grVtot3(iDir)%copyPotDen(vTot_nosym)
257 0 : CALL grVtot3(iDir)%resetPotDen()
258 0 : CALL grVC3(iDir)%copyPotDen(vTot_nosym)
259 0 : CALL grVC3(iDir)%resetPotDen()
260 : ! Generate the external potential gradient.
261 0 : write(oUnit, *) "grVext", iDir
262 0 : sigma_loc = cmplx(0.0,0.0)
263 : !IF (iDir==3) sigma_loc = sigma_ext
264 : CALL vgen_coulomb(1, fmpi_nosym, fi_nosym%input, fi_nosym%field, fi_nosym%vacuum, fi_nosym%sym, fi%juphon, stars_nosym, fi_nosym%cell, &
265 : & sphhar_nosym, fi_nosym%atoms, .FALSE., imagrhodummy, grVext3(iDir), sigma_loc, &
266 0 : & dfptdenimag=imagrhodummy, dfptvCoulimag=grvextdummy,dfptden0=imagrhodummy,stars2=stars_nosym,iDtype=0,iDir=iDir)
267 0 : IF (iDir==3) sigma_gext(iDir,:) = sigma_loc
268 : END DO
269 : !CALL vext_dummy%copyPotDen(vTot_nosym)
270 : !CALL vext_dummy%resetPotDen()
271 : ! Density gradient
272 0 : DO iSpin = 1, SIZE(rho_nosym%mt,4)
273 0 : CALL mt_gradient_new(fi_nosym%atoms, sphhar_nosym, fi_nosym%sym, rho_nosym%mt(:, :, :, iSpin), grrhodummy(:, :, :, iSpin, :))
274 : END DO
275 0 : DO zInd = -stars_nosym%mx3, stars_nosym%mx3
276 0 : DO yInd = -stars_nosym%mx2, stars_nosym%mx2
277 0 : DO xInd = -stars_nosym%mx1, stars_nosym%mx1
278 0 : iStar = stars_nosym%ig(xInd, yInd, zInd)
279 0 : IF (iStar.EQ.0) CYCLE
280 0 : grRho3(1)%pw(iStar,:) = rho_nosym%pw(iStar,:) * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
281 0 : grRho3(2)%pw(iStar,:) = rho_nosym%pw(iStar,:) * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
282 0 : grRho3(3)%pw(iStar,:) = rho_nosym%pw(iStar,:) * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
283 : grgrvextnum(1,1)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
284 0 : * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
285 : grgrvextnum(1,2)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
286 0 : * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
287 : grgrvextnum(1,3)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
288 0 : * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
289 : grgrvextnum(2,1)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
290 0 : * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
291 : grgrvextnum(2,2)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
292 0 : * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
293 : grgrvextnum(2,3)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
294 0 : * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
295 : grgrvextnum(3,1)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
296 0 : * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
297 : grgrvextnum(3,2)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
298 0 : * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
299 : grgrvextnum(3,3)%pw(iStar,:) = vext_dummy%pw(iStar,:) * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat))) &
300 0 : * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
301 : END DO
302 : END DO
303 : END DO
304 :
305 0 : IF (fi_nosym%input%film) THEN
306 0 : DO yInd = -stars_nosym%mx2, stars_nosym%mx2
307 0 : DO xInd = -stars_nosym%mx1, stars_nosym%mx1
308 0 : iStar = stars_nosym%ig(xInd, yInd, 0)
309 0 : IF (iStar.EQ.0) CYCLE
310 0 : iStar = stars_nosym%ig2(iStar)
311 0 : grRho3(1)%vac(:,iStar,:,:) = rho_nosym%vac(:,iStar,:,:) * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,0]),fi_nosym%cell%bmat)))
312 0 : grRho3(2)%vac(:,iStar,:,:) = rho_nosym%vac(:,iStar,:,:) * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,0]),fi_nosym%cell%bmat)))
313 0 : DO iVac = 1, fi_nosym%vacuum%nvac
314 0 : DO iSpin = 1, SIZE(rho_nosym%vac,4)
315 0 : zlim = MERGE(fi_nosym%vacuum%nmz,fi_nosym%vacuum%nmzxy,iStar==1)
316 0 : CALL grdchlh(fi_nosym%vacuum%delz, REAL(rho_nosym%vac(:zlim,iStar,iVac,iSpin)),dr_re(:zlim),drr_dummy(:zlim))
317 0 : CALL grdchlh(fi_nosym%vacuum%delz,AIMAG(rho_nosym%vac(:zlim,iStar,iVac,iSpin)),dr_im(:zlim),drr_dummy(:zlim))
318 0 : grRho3(3)%vac(:,iStar,iVac,iSpin) = (3-2*iVac)*(dr_re + ImagUnit * dr_im)
319 : END DO
320 : END DO
321 : END DO
322 : END DO
323 : END IF
324 :
325 : ! Coulomb/Effective potential gradients
326 0 : DO iDir = 1, 3
327 0 : CALL sh_to_lh(fi_nosym%sym, fi_nosym%atoms, sphhar_nosym, SIZE(rho_nosym%mt,4), 2, grrhodummy(:, :, :, :, iDir), grRho3(iDir)%mt, imagrhodummy%mt)
328 0 : CALL imagrhodummy%resetPotDen()
329 0 : write(oUnit, *) "grVeff", iDir
330 0 : sigma_loc = cmplx(0.0,0.0)
331 0 : IF (iDir==3) sigma_loc = sigma_coul
332 : CALL dfpt_vgen(hybdat_nosym, fi_nosym%field, fi_nosym%input, xcpot_nosym, fi_nosym%atoms, sphhar_nosym, stars_nosym, fi_nosym%vacuum, fi_nosym%sym, &
333 : fi%juphon, fi_nosym%cell, fmpi_nosym, fi_nosym%noco, nococonv_nosym, rho_nosym, vTot_nosym, &
334 0 : stars_nosym, imagrhodummy, grVtot3(iDir), .TRUE., grvextdummy, grRho3(iDir), 0, iDir, [0,0], sigma_loc)
335 0 : write(oUnit, *) "grVC", iDir
336 0 : sigma_loc = cmplx(0.0,0.0)
337 0 : IF (iDir==3) sigma_loc = sigma_coul
338 : CALL dfpt_vgen(hybdat_nosym, fi_nosym%field, fi_nosym%input, xcpot_nosym, fi_nosym%atoms, sphhar_nosym, stars_nosym, fi_nosym%vacuum, fi_nosym%sym, &
339 : fi%juphon, fi_nosym%cell, fmpi_nosym, fi_nosym%noco, nococonv_nosym, rho_nosym, vTot_nosym, &
340 0 : stars_nosym, imagrhodummy, grVC3(iDir), .FALSE., grvextdummy, grRho3(iDir), 0, iDir, [0,0], sigma_loc)
341 : END DO
342 :
343 0 : DO iDir2 = 1, 3
344 0 : DO iDir = 1, 3
345 0 : CALL imagrhodummy%resetPotDen()
346 0 : sigma_loc = cmplx(0.0,0.0)
347 :
348 : !IF (iDir2==3) sigma_loc = sigma_gext(iDir,:)
349 : !IF (iDir==3) sigma_loc = sigma_gext(iDir2,:)
350 : CALL vgen_coulomb(1, fmpi_nosym, fi_nosym%input, fi_nosym%field, fi_nosym%vacuum, fi_nosym%sym, fi%juphon, stars_nosym, fi_nosym%cell, &
351 : & sphhar_nosym, fi_nosym%atoms, .TRUE., imagrhodummy, grgrVC3x3(iDir2,iDir), sigma_loc, &
352 : & dfptdenimag=imagrhodummy, dfptvCoulimag=grvextdummy,dfptden0=imagrhodummy,stars2=stars_nosym,iDtype=0,iDir=iDir,iDir2=iDir2, &
353 0 : & sigma_disc2=MERGE(sigma_ext,[cmplx(0.0,0.0),cmplx(0.0,0.0)],iDir2==3.AND.iDir==3.AND..FALSE.))
354 0 : CALL dfpt_e2_madelung(fi_nosym%atoms,fi_nosym%input%jspins,imagrhodummy%mt(:,0,:,:),grgrVC3x3(iDir2,iDir)%mt(:,0,:,1),e2_vm(:,iDir2,iDir))
355 : END DO
356 : END DO
357 0 : CALL save_npy("radii.npy",fi_nosym%atoms%rmsh(:,1))
358 0 : DO iDir2 = 1, 3
359 0 : DO iDir = 1, 3
360 0 : CALL save_npy("grgrVC_"//int2str(idir2)//int2str(idir)//"_pw.npy",grgrVC3x3(iDir2,iDir)%pw(:,1))
361 0 : CALL save_npy("grgrVCnum_"//int2str(idir2)//int2str(idir)//"_pw.npy",grgrvextnum(iDir2,iDir)%pw(:,1))
362 : END DO
363 : END DO
364 :
365 0 : CALL grRho3(1)%distribute(fmpi%mpi_comm)
366 0 : CALL grRho3(2)%distribute(fmpi%mpi_comm)
367 0 : CALL grRho3(3)%distribute(fmpi%mpi_comm)
368 0 : CALL grVext3(1)%distribute(fmpi%mpi_comm)
369 0 : CALL grVext3(2)%distribute(fmpi%mpi_comm)
370 0 : CALL grVext3(3)%distribute(fmpi%mpi_comm)
371 0 : CALL timestop("Gradient generation")
372 :
373 0 : CALL test_vac_stuff(fi_nosym,stars_nosym,sphhar_nosym,rho_nosym,vTot_nosym,grRho3,grVtot3,grVC3,grVext3,grrhodummy,grid)
374 :
375 : ! Old CRG-jp Routine to get the vectors G+q for Eii2
376 0 : CALL genPertPotDensGvecs( stars_nosym, fi_nosym%cell, fi_nosym%input, ngdp, ngdp2km, [0.0,0.0,0.0], recG )
377 :
378 : ! The eig_ids, where the stuff of k+q, the perturbed stuff and some extra dynmat stuff will be saved.
379 : q_eig_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
380 0 : .NOT.fi%INPUT%eig66(1), fi%input%l_real, fi%noco%l_soc, fi%input%eig66(1), .FALSE., fmpi%n_size)
381 : dfpt_eig_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
382 0 : .NOT.fi%INPUT%eig66(1), .FALSE., fi%noco%l_soc, fi%INPUT%eig66(1), .FALSE., fmpi%n_size)
383 : dfpt_eig_id2 = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
384 0 : .NOT.fi%INPUT%eig66(1), .FALSE., fi%noco%l_soc, fi%INPUT%eig66(1), .FALSE., fmpi%n_size)
385 :
386 : IF (l_minusq) THEN
387 : qm_eig_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
388 : .NOT.fi%INPUT%eig66(1), fi%input%l_real, fi%noco%l_soc, fi%input%eig66(1), .FALSE., fmpi%n_size)
389 : dfpt_eigm_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
390 : .NOT.fi%INPUT%eig66(1), .FALSE., fi%noco%l_soc, fi%INPUT%eig66(1), .FALSE., fmpi%n_size)
391 : dfpt_eigm_id2 = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
392 : .NOT.fi%INPUT%eig66(1), .FALSE., fi%noco%l_soc, fi%INPUT%eig66(1), .FALSE., fmpi%n_size)
393 : END IF
394 :
395 0 : ALLOCATE(dyn_mat(SIZE(q_list),3*fi_nosym%atoms%ntype,3*fi_nosym%atoms%ntype))
396 0 : dyn_mat = cmplx(0.0,0.0)
397 0 : ALLOCATE(sym_dyn_mat(SIZE(q_list),3*fi_nosym%atoms%ntype,3*fi_nosym%atoms%ntype))
398 0 : sym_dyn_mat = cmplx(0.0,0.0)
399 0 : IF (l_dfpt_scf) THEN
400 : ! Do the self-consistency calculations for each specified q, for all atoms and for
401 : ! all three cartesian directions.
402 : ! TODO: The effort here should be greatly reducible by symmetry considerations.
403 0 : write(*,*) fi%juPhon%startq/=0, fi%juPhon%stopq, size(q_list)
404 :
405 0 : DO iQ = fi%juPhon%startq, MERGE(fi%juPhon%stopq,SIZE(q_list),fi%juPhon%stopq/=0)
406 0 : CALL timestart("q-point")
407 0 : IF (.NOT.fi%juPhon%qmode==0) THEN
408 0 : CALL make_sym_list(fi_fullsym%sym, qpts_loc%bk(:,q_list(iQ)),sym_count,sym_list)
409 0 : ALLOCATE(sym_dynvec(3*fi_nosym%atoms%ntype,3*fi_nosym%atoms%ntype-1,sym_count))
410 : END IF
411 0 : kqpts = fi%kpts
412 : ! Modify this from kpts only in DFPT case.
413 0 : DO ikpt = 1, fi%kpts%nkpt
414 0 : kqpts%bk(:, ikpt) = kqpts%bk(:, ikpt) + qpts_loc%bk(:,q_list(iQ))
415 : END DO
416 :
417 : IF (l_minusq) THEN
418 : kqmpts = fi%kpts
419 : ! Modify this from kpts only in DFPT case.
420 : DO ikpt = 1, fi%kpts%nkpt
421 : kqmpts%bk(:, ikpt) = kqmpts%bk(:, ikpt) - qpts_loc%bk(:,q_list(iQ))
422 : END DO
423 : END IF
424 :
425 0 : CALL timestart("Eii2")
426 0 : CALL CalcIIEnerg2(fi_nosym%atoms, fi_nosym%cell, qpts_loc, stars_nosym, fi_nosym%input, q_list(iQ), ngdp, recG, E2ndOrdII)
427 0 : CALL timestop("Eii2")
428 :
429 0 : write(9991,*) "Eii2 old:", E2ndOrdII
430 0 : E2ndOrdII = CMPLX(0.0,0.0)
431 0 : DO iDtype = 1, fi_nosym%atoms%ntype
432 0 : DO iDir2 = 1, 3
433 0 : DO iDir = 1, 3
434 0 : E2ndOrdII(3*(iDtype-1)+iDir2,3*(iDtype-1)+iDir) = e2_vm(iDtype,iDir2,iDir)
435 : END DO
436 : END DO
437 : END DO
438 0 : CALL timestart("Eigenstuff at k+q")
439 :
440 : ! This was an additional eig_id to test a specific shift from k to k+q in the dynmat setup. We leave it in
441 : ! for now, as it might be tested again in the future.
442 : !q_eig_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
443 : ! .NOT.fi%INPUT%eig66(1), fi%input%l_real, fi%noco%l_soc, fi%input%eig66(1), .FALSE., fmpi%n_size)
444 :
445 : ! Get the eigenstuff at k+q
446 0 : CALL q_results%reset_results(fi%input)
447 :
448 : CALL eigen(fi, fmpi, stars, sphhar, xcpot, forcetheo, enpara, nococonv, mpdata, &
449 : hybdat, 1, q_eig_id, q_results, rho, vTot, vxc, hub1data, &
450 0 : qpts_loc%bk(:,q_list(iQ)))
451 :
452 : ! Fermi level and occupancies
453 0 : CALL timestart("determination of fermi energy")
454 0 : CALL fermie(q_eig_id, fmpi, kqpts, fi%input, fi%noco, enpara%epara_min, fi%cell, q_results)
455 0 : CALL timestop("determination of fermi energy")
456 :
457 : #ifdef CPP_MPI
458 0 : CALL MPI_BCAST(q_results%ef, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
459 0 : CALL MPI_BCAST(q_results%w_iks, SIZE(q_results%w_iks), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
460 : #endif
461 :
462 0 : CALL timestop("Eigenstuff at k+q")
463 :
464 : IF (l_minusq) THEN
465 : CALL timestart("Eigenstuff at k-q")
466 : !qm_eig_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, fi%input%jspins, fi%noco%l_noco, &
467 : ! .NOT.fi%INPUT%eig66(1), fi%input%l_real, fi%noco%l_soc, fi%input%eig66(1), .FALSE., fmpi%n_size)
468 :
469 : CALL qm_results%reset_results(fi%input)
470 :
471 : CALL eigen(fi, fmpi, stars, sphhar, xcpot, forcetheo, enpara, nococonv, mpdata, &
472 : hybdat, 1, qm_eig_id, qm_results, rho, vTot, vxc, hub1data, &
473 : -qpts_loc%bk(:,q_list(iQ)))
474 :
475 : ! Fermi level and occupancies
476 : CALL timestart("determination of fermi energy")
477 : CALL fermie(qm_eig_id, fmpi, kqmpts, fi%input, fi%noco, enpara%epara_min, fi%cell, qm_results)
478 : CALL timestop("determination of fermi energy")
479 :
480 : #ifdef CPP_MPI
481 : CALL MPI_BCAST(qm_results%ef, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
482 : CALL MPI_BCAST(qm_results%w_iks, SIZE(qm_results%w_iks), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
483 : #endif
484 :
485 : CALL timestop("Eigenstuff at k-q")
486 : END IF
487 :
488 0 : DO iDtype = 1, fi_nosym%atoms%ntype
489 0 : CALL timestart("Typeloop")
490 0 : DO iDir = 1, 3
491 0 : CALL timestart("Dirloop")
492 0 : IF (.NOT.fi%juPhon%qmode==0.AND.fmpi%irank==0) THEN
493 0 : IF (iDtype==1.AND.iDir==2) sym_dyn_mat(iQ, 1, :) = dyn_mat(iQ, 1, :)
494 0 : IF (3 *(iDtype-1)+iDir>1) THEN
495 0 : CALL cheat_dynmat(fi_fullsym%atoms, fi_fullsym%sym, fi_fullsym%cell%amat, qpts_loc%bk(:,q_list(iQ)), iDtype, iDir, sym_count, sym_list(:sym_count), sym_dynvec, dyn_mat(iQ,:,:), sym_dyn_mat(iQ,:,:), l_cheated)
496 : END IF
497 0 : IF (l_cheated) WRITE(*,*) "Following row was cheated!"
498 0 : IF (l_cheated) write(*,*) sym_dyn_mat(iQ,3 *(iDtype-1)+iDir,:)
499 : END IF
500 0 : dfpt_tag = ''
501 0 : WRITE(dfpt_tag,'(a1,i0,a2,i0,a2,i0)') 'q', q_list(iQ), '_b', iDtype, '_j', iDir
502 :
503 0 : IF (fmpi%irank==0) THEN
504 0 : WRITE(*,*) 'Starting calculation for:'
505 0 : WRITE(*,*) ' q = ', qpts_loc%bk(:,q_list(iQ))
506 0 : WRITE(*,*) ' atom = ', iDtype
507 0 : WRITE(*,*) ' direction = ', iDir
508 : END IF
509 :
510 : !IF (fmpi_nosym%irank==0) THEN
511 0 : CALL starsq%reset_stars()
512 : IF (l_minusq) CALL starsmq%reset_stars()
513 0 : CALL denIn1%reset_dfpt()
514 0 : CALL denIn1Im%reset_dfpt()
515 0 : CALL vTot1%reset_dfpt()
516 0 : CALL vTot1Im%reset_dfpt()
517 : IF (l_minusq) CALL vTot1m%reset_dfpt()
518 : IF (l_minusq) CALL vTot1mIm%reset_dfpt()
519 0 : CALL vC1%reset_dfpt()
520 0 : CALL vC1Im%reset_dfpt()
521 0 : CALL results1%reset_results(fi_nosym%input)
522 : !END IF
523 :
524 0 : IF (fmpi%irank==0) WRITE(*,*) '-------------------------'
525 : ! This is where the magic happens. The Sternheimer equation is solved
526 : ! iteratively, providing the scf part of dfpt calculations.
527 : IF (l_minusq) THEN
528 : CALL timestart("Sternheimer with -q")
529 : CALL dfpt_sternheimer(fi_nosym, xcpot_nosym, sphhar_nosym, stars_nosym, starsq, nococonv_nosym, qpts_loc, fmpi_nosym, results_nosym, q_results, enpara_nosym, hybdat_nosym, &
530 : rho_nosym, vTot_nosym, grRho3(iDir), grVtot3(iDir), grVext3(iDir), q_list(iQ), iDtype, iDir, &
531 : dfpt_tag, eig_id, l_real, results1, dfpt_eig_id, dfpt_eig_id2, q_eig_id, &
532 : denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im, MERGE(sigma_ext,[cmplx(0.0,0.0),cmplx(0.0,0.0)],iDir==3), &
533 : MERGE(sigma_coul,[cmplx(0.0,0.0),cmplx(0.0,0.0)],iDir==3),&
534 : starsmq, qm_results, dfpt_eigm_id, dfpt_eigm_id2, qm_eig_id, results1m, vTot1m, vTot1mIm)
535 : CALL timestop("Sternheimer with -q")
536 : ELSE
537 0 : CALL timestart("Sternheimer")
538 : CALL dfpt_sternheimer(fi_nosym, xcpot_nosym, sphhar_nosym, stars_nosym, starsq, nococonv_nosym, qpts_loc, fmpi_nosym, results_nosym, q_results, enpara_nosym, hybdat_nosym, &
539 : rho_nosym, vTot_nosym, grRho3(iDir), grVtot3(iDir), grVext3(iDir), q_list(iQ), iDtype, iDir, &
540 : dfpt_tag, eig_id, l_real, results1, dfpt_eig_id, dfpt_eig_id2, q_eig_id, &
541 : denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im, MERGE(sigma_ext,[cmplx(0.0,0.0),cmplx(0.0,0.0)],iDir==3), &
542 0 : MERGE(sigma_coul,[cmplx(0.0,0.0),cmplx(0.0,0.0)],iDir==3))
543 0 : CALL timestop("Sternheimer")
544 : END IF
545 :
546 0 : IF (fmpi%irank==0) WRITE(*,*) '-------------------------'
547 0 : CALL timestart("Dynmat")
548 : ! Once the first order quantities are converged, we can construct all
549 : ! additional necessary quantities and from that the dynamical matrix.
550 : IF (.TRUE.) THEN
551 : CALL dfpt_dynmat_row(fi_nosym, stars_nosym, starsq, sphhar_nosym, xcpot_nosym, nococonv_nosym, hybdat_nosym, fmpi_nosym, qpts_loc, q_list(iQ), iDtype, iDir, &
552 : eig_id, dfpt_eig_id, dfpt_eig_id2, enpara_nosym, results_nosym, results1, l_real,&
553 : rho_nosym, vTot_nosym, grRho3, grVext3, grVC3, &
554 0 : denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im, dyn_mat(iQ,3 *(iDtype-1)+iDir,:), E2ndOrdII, sigma_ext, sigma_gext)
555 : ELSE
556 : CALL dfpt_dynmat_row(fi_nosym, stars_nosym, starsq, sphhar_nosym, xcpot_nosym, nococonv_nosym, hybdat_nosym, fmpi_nosym, qpts_loc, q_list(iQ), iDtype, iDir, &
557 : eig_id, dfpt_eig_id, dfpt_eig_id2, enpara_nosym, results_nosym, results1, l_real,&
558 : rho_nosym, vTot_nosym, grRho3, grVext3, grVC3, &
559 : denIn1, vTot1, denIn1Im, vTot1Im, vC1, vC1Im, dyn_mat(iQ,3 *(iDtype-1)+iDir,:), E2ndOrdII, sigma_ext, sigma_gext, q_eig_id)
560 : END IF
561 0 : CALL timestop("Dynmat")
562 0 : dyn_mat(iQ,3 *(iDtype-1)+iDir,:) = dyn_mat(iQ,3 *(iDtype-1)+iDir,:) + conjg(E2ndOrdII(3 *(iDtype-1)+iDir,:))
563 0 : IF (.NOT.fi%juPhon%qmode==0) THEN
564 0 : CALL make_sym_dynvec(fi_fullsym%atoms, fi_fullsym%sym, fi_fullsym%cell%amat, qpts_loc%bk(:,q_list(iQ)), iDtype, iDir, sym_count, sym_list(:sym_count), dyn_mat(iQ,3 *(iDtype-1)+iDir,:), sym_dynvec)
565 : END IF
566 :
567 0 : IF (fmpi%irank==0) write(*,*) "dynmat row for ", dfpt_tag
568 0 : IF (fmpi%irank==0) write(*,*) dyn_mat(iQ,3 *(iDtype-1)+iDir,:)
569 0 : IF (fmpi%irank==0.AND.l_cheated) write(*,*) "The cheat:"
570 0 : IF (fmpi%irank==0.AND.l_cheated) write(*,*) sym_dyn_mat(iQ,3 *(iDtype-1)+iDir,:)
571 0 : l_cheated = .FALSE.
572 0 : IF (fmpi%irank==0) WRITE(9339,*) dyn_mat(iQ,3 *(iDtype-1)+iDir,:)
573 0 : IF (fmpi_nosym%irank == 0 .AND. fi_nosym%juphon%l_rm_qhdf) call system("rm "//TRIM(dfpt_tag)//".hdf")
574 0 : CALL timestop("Dirloop")
575 : END DO
576 0 : CALL timestop("Typeloop")
577 :
578 : #if defined(CPP_MPI)
579 0 : CALL MPI_BARRIER(fmpi%MPI_COMM,ierr)
580 : #endif
581 : END DO
582 :
583 0 : IF (fmpi%irank==0) THEN
584 0 : WRITE(*,*) '-------------------------'
585 0 : CALL timestart("Dynmat diagonalization")
586 0 : CALL DiagonalizeDynMat(fi_nosym%atoms, qpts_loc%bk(:,q_list(iQ)), fi%juPhon%calcEigenVec, dyn_mat(iQ,:,:), eigenVals, eigenVecs, q_list(iQ),.TRUE.,"raw",fi_nosym%juphon%l_sumrule)
587 0 : CALL timestop("Dynmat diagonalization")
588 :
589 0 : CALL timestart("Frequency calculation")
590 0 : CALL CalculateFrequencies(fi_nosym%atoms, q_list(iQ), eigenVals, eigenFreqs,"raw")
591 0 : CALL timestop("Frequency calculation")
592 0 : write(9991,*) "Eii2 new:", E2ndOrdII
593 0 : DEALLOCATE(eigenVals, eigenVecs, eigenFreqs, E2ndOrdII)
594 : END IF
595 : !CALL close_eig(q_eig_id)
596 : !IF (l_minusq) CALL close_eig(qm_eig_id)
597 0 : IF (.NOT.fi%juPhon%qmode==0) THEN
598 0 : DEALLOCATE(sym_dynvec)
599 : END IF
600 0 : CALL timestop("q-point")
601 :
602 : END DO
603 : END IF
604 :
605 : ! If the Dynmats-Files were already created, we can read them in and do postprocessing.
606 : ! a) Transform the q-Mesh onto real space.
607 : ! b) Transform it back onto a dense q-path.
608 : ! c) Transform it back to a denser grid
609 : ! d) Perform a DOS calculation for the denser grid.
610 0 : IF (fmpi%irank==0) THEN ! Band/Dos stuff
611 : ! 0) Read
612 0 : DO iQ = 1, fi_fullsym%kpts%nkpt ! Loop over dynmat files to read
613 0 : IF (iQ<=9) THEN
614 0 : OPEN( 3001, file="dynMatq=000"//int2str(iQ), status="old")
615 : ELSE
616 0 : OPEN( 3001, file="dynMatq=00"//int2str(iQ), status="old")
617 : END IF
618 0 : DO iread = 1, 3 + 3*fi%atoms%nat ! Loop over dynmat rows
619 0 : IF (iread<4) THEN
620 0 : READ( 3001,*) trash
621 0 : write(*,*) iread, trash
622 : ELSE
623 0 : READ( 3001,*) numbers(iread-3,:)
624 0 : write(*,*) iread, numbers(iread-3,:)
625 0 : dyn_mat(iQ,iread-3,:) = CMPLX(numbers(iread-3,::2),numbers(iread-3,2::2))
626 : END IF
627 : END DO ! iread
628 0 : CLOSE(3001)
629 : END DO ! iQ
630 :
631 : ! a) Real space transformation
632 0 : ALLOCATE(dyn_mat_r(fi_fullsym%kpts%nkptf,3*fi%atoms%nat,3*fi%atoms%nat))
633 0 : CALL ft_dyn(fi_fullsym%atoms, fi_fullsym%kpts, fi_fullsym%sym, fi_fullsym%cell%amat, dyn_mat, dyn_mat_r, dyn_mat_q_full)
634 :
635 : ! b/c) reciprocal space transformation for bands/dense grid
636 0 : IF (l_dfpt_band.OR.l_dfpt_full) THEN
637 0 : IF (l_dfpt_band) THEN
638 0 : dynfiletag = "band"
639 0 : ELSE IF (l_dfpt_full) THEN
640 0 : dynfiletag = "full"
641 : ELSE
642 0 : dynfiletag = "intp"
643 : END IF
644 0 : IF (l_dfpt_dos) ALLOCATE(eigenValsFull(3*fi%atoms%nat,fi_nosym%kpts%nkpt,fi_nosym%input%jspins))
645 0 : DO iQ = 1, fi_nosym%kpts%nkpt
646 0 : CALL ift_dyn(fi_fullsym%atoms,fi_fullsym%kpts,fi_fullsym%sym,fi_fullsym%cell%amat,fi_nosym%kpts%bk(:,iQ),dyn_mat_r,dyn_mat_pathq)
647 0 : WRITE(*,*) '-------------------------'
648 0 : CALL timestart("Dynmat diagonalization")
649 0 : CALL DiagonalizeDynMat(fi_nosym%atoms, fi_nosym%kpts%bk(:,iQ), fi%juPhon%calcEigenVec, dyn_mat_pathq, eigenVals, eigenVecs, iQ,.FALSE.,TRIM(dynfiletag),fi_nosym%juphon%l_sumrule)
650 0 : CALL timestop("Dynmat diagonalization")
651 :
652 0 : CALL timestart("Frequency calculation")
653 0 : CALL CalculateFrequencies(fi_nosym%atoms, iQ, eigenVals, eigenFreqs,TRIM(dynfiletag))
654 0 : CALL timestop("Frequency calculation")
655 :
656 0 : IF (l_dfpt_dos) eigenValsFull(:,iQ,1) = eigenFreqs(:)
657 :
658 0 : DEALLOCATE(eigenVals, eigenVecs, eigenFreqs, dyn_mat_pathq)
659 : END DO ! iQ
660 : END IF ! bands/interpolation
661 0 : IF (l_dfpt_dos) THEN
662 0 : fi_nosym%banddos%dos = .TRUE.
663 0 : CALL dos%init(fi_nosym%input,fi_nosym%atoms,fi_nosym%kpts,fi_nosym%banddos,eigenValsFull)
664 0 : allocate(eigdos(1))
665 0 : eigdos(1)%p=>dos
666 0 : CALL make_dos(fi_nosym%kpts,fi_nosym%atoms,fi_nosym%vacuum,fi_nosym%input,fi_nosym%banddos,fi_nosym%sliceplot,fi_nosym%noco,fi_nosym%sym,fi_nosym%cell,results,eigdos,fi%juPhon )
667 : END IF ! dos
668 : END IF ! Band/Dos stuff
669 :
670 0 : CALL close_eig(q_eig_id)
671 0 : CALL close_eig(dfpt_eig_id)
672 0 : CALL close_eig(dfpt_eig_id2)
673 : IF (l_minusq) THEN
674 : CALL close_eig(qm_eig_id)
675 : CALL close_eig(dfpt_eigm_id)
676 : CALL close_eig(dfpt_eigm_id2)
677 : END IF
678 :
679 0 : DEALLOCATE(recG)
680 :
681 0 : WRITE (oUnit,*) '------------------------------------------------------'
682 :
683 0 : CALL juDFT_end("Phonon calculation finished.",fmpi%irank)
684 :
685 0 : END SUBROUTINE dfpt
686 :
687 0 : SUBROUTINE dfpt_desym(fmpi_nosym,fi_nosym,sphhar_nosym,stars_nosym,nococonv_nosym,enpara_nosym,results_nosym,wann_nosym,hybdat_nosym,mpdata_nosym,xcpot_nosym,forcetheo_nosym,rho_nosym,vTot_nosym,grid,inp_pref,&
688 : fi,sphhar,stars,nococonv,enpara,results,rho,vTot)
689 : USE m_desymmetrizer
690 : USE m_outcdn
691 : USE m_plot
692 : USE m_fleur_init
693 :
694 : TYPE(t_mpi), INTENT(INOUT) :: fmpi_nosym
695 : TYPE(t_fleurinput), INTENT(INOUT) :: fi_nosym
696 : TYPE(t_sphhar), INTENT(INOUT) :: sphhar_nosym
697 : TYPE(t_stars), INTENT(INOUT) :: stars_nosym
698 : TYPE(t_nococonv), INTENT(INOUT) :: nococonv_nosym
699 : TYPE(t_enpara), INTENT(INOUT) :: enpara_nosym
700 : TYPE(t_results), INTENT(INOUT) :: results_nosym
701 : TYPE(t_wann), INTENT(INOUT) :: wann_nosym
702 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat_nosym
703 : TYPE(t_mpdata), INTENT(INOUT) :: mpdata_nosym
704 :
705 : CLASS(t_xcpot), ALLOCATABLE, INTENT(INOUT) :: xcpot_nosym
706 : CLASS(t_forcetheo), ALLOCATABLE, INTENT(INOUT) :: forcetheo_nosym
707 :
708 : TYPE(t_potden), INTENT(INOUT):: rho_nosym, vTot_nosym
709 : TYPE(t_fleurinput), INTENT(IN) :: fi
710 : TYPE(t_sphhar), INTENT(IN) :: sphhar
711 : TYPE(t_stars), INTENT(IN) :: stars
712 : TYPE(t_nococonv), INTENT(IN) :: nococonv
713 : TYPE(t_enpara), INTENT(IN) :: enpara
714 : TYPE(t_results), INTENT(IN) :: results
715 : TYPE(t_potden), INTENT(IN) :: rho, vTot
716 :
717 : INTEGER, INTENT(IN) :: grid(3)
718 :
719 : CHARACTER(len=100), INTENT(IN) :: inp_pref
720 :
721 : INTEGER :: ix, iy, iz, iv_old, iflag_old, iv_new, iflag_new
722 : INTEGER :: iType_old, iAtom_old, iType_new, iAtom_new!, inversionOp
723 : REAL :: old_point(3), new_point(3), pt_old(3), pt_new(3), xdnout_old, xdnout_new!, atom_shift(3)
724 : LOGICAL :: test_desym
725 : CALL fleur_init(fmpi_nosym, fi_nosym, sphhar_nosym, stars_nosym, nococonv_nosym, forcetheo_nosym, &
726 : enpara_nosym, xcpot_nosym, results_nosym, wann_nosym, hybdat_nosym, mpdata_nosym, &
727 0 : inp_pref)
728 :
729 0 : CALL rho_nosym%init(stars_nosym,fi_nosym%atoms,sphhar_nosym,fi_nosym%vacuum,fi_nosym%noco,fi%input%jspins,POTDEN_TYPE_DEN)
730 0 : CALL vTot_nosym%init(stars_nosym,fi_nosym%atoms,sphhar_nosym,fi_nosym%vacuum,fi_nosym%noco,fi%input%jspins,POTDEN_TYPE_POTTOT)
731 :
732 : ! TODO: Correctly account for such a shift in the desymmetrization.
733 : ! For now: Just build input, that does not necessitate a shift.
734 : ! inversionOp = -1
735 : ! symOpLoop: DO iSym = 1, sym%nop
736 : ! IF (ALL(sym%mrot(:,:,iSym)==invs_matrix)) THEN
737 : ! inversionOp = iSym
738 : ! EXIT symOpLoop
739 : ! END IF
740 : ! END DO symOpLoop
741 :
742 : ! atom_shift = 0.0
743 : ! IF (inversionOp.GT.0) THEN
744 : ! IF(ANY(ABS(sym%tau(:,inversionOp)).GT.eps7).and..not.(film.and.ABS(sym%tau(3,inversionOp))>eps7)) THEN
745 : ! atom_shift = 0.5*sym%tau(:,inversionOp)
746 : ! END IF
747 : ! END IF
748 :
749 0 : ALLOCATE(vTot_nosym%pw_w, mold=vTot_nosym%pw)
750 0 : vTot_nosym%pw_w = CMPLX(0.0,0.0)
751 :
752 0 : CALL desymmetrize_pw(fi%sym, stars, stars_nosym, rho%pw, rho_nosym%pw)
753 0 : CALL desymmetrize_pw(fi%sym, stars, stars_nosym, vTot%pw, vTot_nosym%pw, vTot%pw_w, vTot_nosym%pw_w)
754 0 : CALL desymmetrize_mt(fi%sym, fi_nosym%sym, fi%cell, fi%atoms, fi_nosym%atoms, sphhar, sphhar_nosym, rho%mt, rho_nosym%mt)
755 0 : CALL desymmetrize_mt(fi%sym, fi_nosym%sym, fi%cell, fi%atoms, fi_nosym%atoms, sphhar, sphhar_nosym, vTot%mt, vTot_nosym%mt)
756 :
757 : CALL desymmetrize_types(fi%input, fi_nosym%input, fi%atoms, fi_nosym%atoms, fi%noco, &
758 0 : nococonv, nococonv_nosym, enpara, enpara_nosym, results, results_nosym)
759 :
760 0 : test_desym = .FALSE.
761 : IF (test_desym) THEN
762 : DO iz = 0, grid(3)-1
763 : DO iy = 0, grid(2)-1
764 : DO ix = 0, grid(1)-1
765 : old_point = fi%cell%amat(:,1)*REAL(ix)/(grid(1)-1) + &
766 : fi%cell%amat(:,2)*REAL(iy)/(grid(2)-1) + &
767 : fi%cell%amat(:,3)*REAL(iz)/(grid(3)-1)
768 :
769 : new_point = fi%cell%amat(:,1)*REAL(ix)/(grid(1)-1) + &
770 : fi%cell%amat(:,2)*REAL(iy)/(grid(2)-1) + &
771 : fi%cell%amat(:,3)*REAL(iz)/(grid(3)-1)! - &
772 : !atom_shift
773 :
774 : ! Set region specific parameters for point
775 : ! Get MT sphere for point if point is in MT sphere
776 : CALL getMTSphere(fi%input,fi%cell,fi%atoms,old_point,iType_old,iAtom_old,pt_old)
777 : CALL getMTSphere(fi_nosym%input,fi_nosym%cell,fi_nosym%atoms,new_point,iType_new,iAtom_new,pt_new)
778 : IF (iAtom_old.NE.0) THEN
779 : iv_old = 0
780 : iflag_old = 1
781 : ELSE
782 : iv_old = 0
783 : iflag_old = 2
784 : pt_old(:) = old_point(:)
785 : END IF
786 :
787 : IF (iAtom_new.NE.0) THEN
788 : iv_new = 0
789 : iflag_new = 1
790 : ELSE
791 : iv_new = 0
792 : iflag_new = 2
793 : pt_new(:) = new_point(:)
794 : END IF
795 :
796 : ! Old point:
797 : CALL outcdn(pt_old,iType_old,iAtom_old,iv_old,iflag_old,1,.FALSE.,stars,&
798 : fi%vacuum,sphhar,fi%atoms,fi%sym,fi%cell ,&
799 : rho,xdnout_old)
800 : ! New point:
801 : CALL outcdn(pt_new,iType_new,iAtom_new,iv_old,iflag_old,1,.FALSE.,stars_nosym,&
802 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
803 : rho_nosym,xdnout_new)
804 :
805 : WRITE(9004,*) xdnout_new-xdnout_old
806 :
807 : ! Old point:
808 : CALL outcdn(pt_old,iType_old,iAtom_old,iv_old,iflag_old,1,.TRUE.,stars,&
809 : fi%vacuum,sphhar,fi%atoms,fi%sym,fi%cell ,&
810 : vTot,xdnout_old)
811 : ! New point:
812 : CALL outcdn(pt_new,iType_new,iAtom_new,iv_old,iflag_old,1,.TRUE.,stars_nosym,&
813 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
814 : vTot_nosym,xdnout_new)
815 : WRITE(9005,*) xdnout_new-xdnout_old
816 : END DO !x-loop
817 : END DO !y-loop
818 : END DO !z-loop
819 :
820 : CALL save_npy("sym_on_rhopw.npy",rho%pw)
821 : CALL save_npy("sym_off_rhopw.npy",rho_nosym%pw)
822 : CALL save_npy("sym_on_rhomt.npy",rho%mt)
823 : CALL save_npy("sym_off_rhomt.npy",rho_nosym%mt)
824 : CALL save_npy("sym_on_vpw.npy",vTot%pw)
825 : CALL save_npy("sym_off_vpw.npy",vTot_nosym%pw)
826 : CALL save_npy("sym_on_vmt.npy",vTot%mt)
827 : CALL save_npy("sym_off_vmt.npy",vTot_nosym%mt)
828 : END IF
829 0 : END SUBROUTINE
830 :
831 0 : SUBROUTINE test_vac_stuff(fi_nosym,stars_nosym,sphhar_nosym,rho_nosym,vTot_nosym,grRho3,grVtot3,grVC3,grVext3,grrhodummy,grid)
832 : USE m_npy
833 : USE m_outcdn
834 : USE m_grdchlh
835 : USE m_dfpt_gradient
836 :
837 : TYPE(t_fleurinput), INTENT(IN) :: fi_nosym
838 : TYPE(t_stars), INTENT(IN) :: stars_nosym
839 : TYPE(t_sphhar), INTENT(IN) :: sphhar_nosym
840 : TYPE(t_potden), INTENT(IN) :: rho_nosym, vTot_nosym, grRho3(3), grVtot3(3), grVC3(3), grVext3(3)
841 : INTEGER, INTENT(IN) :: grid(3)
842 : COMPLEX, INTENT(INOUT) :: grrhodummy(:, :, :, :, :)
843 :
844 : INTEGER :: ix, iy, iVac, iStar, iSpin, zlim, xInd, yInd, zInd
845 : REAL :: xdnout_grrho_up_pw, xdnout_grrho_up_vac, xdnout_grrho_down_pw, xdnout_grrho_down_vac
846 : REAL :: xdnout_grvc_up_pw, xdnout_grvc_up_vac, xdnout_grvc_down_pw, xdnout_grvc_down_vac
847 : REAL :: point_plus(3), point_minus(3)
848 0 : REAL :: dr_re(fi_nosym%vacuum%nmzd), dr_im(fi_nosym%vacuum%nmzd), drr_dummy(fi_nosym%vacuum%nmzd)
849 :
850 0 : COMPLEX, ALLOCATABLE :: grVtotvac(:,:,:,:), grVtotpw(:,:)
851 :
852 0 : ALLOCATE(grVtotpw(stars_nosym%ng3,3))
853 0 : ALLOCATE(grVtotvac(fi_nosym%vacuum%nmz,stars_nosym%ng2,2,3))
854 0 : DO iSpin = 1, SIZE(rho_nosym%mt,4)
855 0 : CALL mt_gradient_new(fi_nosym%atoms, sphhar_nosym, fi_nosym%sym, vTot_nosym%mt(:, :, :, iSpin), grrhodummy(:, :, :, iSpin, :))
856 : END DO
857 :
858 0 : DO zInd = -stars_nosym%mx3, stars_nosym%mx3
859 0 : DO yInd = -stars_nosym%mx2, stars_nosym%mx2
860 0 : DO xInd = -stars_nosym%mx1, stars_nosym%mx1
861 0 : iStar = stars_nosym%ig(xInd, yInd, zInd)
862 0 : IF (iStar.EQ.0) CYCLE
863 0 : grVtotpw(iStar,1) = vTot_nosym%pw(iStar,1) * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
864 0 : grVtotpw(iStar,2) = vTot_nosym%pw(iStar,1) * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
865 0 : grVtotpw(iStar,3) = vTot_nosym%pw(iStar,1) * cmplx(0.0,dot_product([0.0,0.0,1.0],matmul(real([xInd,yInd,zInd]),fi_nosym%cell%bmat)))
866 : END DO
867 : END DO
868 : END DO
869 :
870 0 : IF (fi_nosym%input%film) THEN
871 0 : DO yInd = -stars_nosym%mx2, stars_nosym%mx2
872 0 : DO xInd = -stars_nosym%mx1, stars_nosym%mx1
873 0 : iStar = stars_nosym%ig(xInd, yInd, 0)
874 0 : IF (iStar.EQ.0) CYCLE
875 0 : iStar = stars_nosym%ig2(iStar)
876 0 : grVtotvac(:,iStar,:,1) = vTot_nosym%vac(:,iStar,:,1) * cmplx(0.0,dot_product([1.0,0.0,0.0],matmul(real([xInd,yInd,0]),fi_nosym%cell%bmat)))
877 0 : grVtotvac(:,iStar,:,2) = vTot_nosym%vac(:,iStar,:,1) * cmplx(0.0,dot_product([0.0,1.0,0.0],matmul(real([xInd,yInd,0]),fi_nosym%cell%bmat)))
878 0 : DO iVac = 1, fi_nosym%vacuum%nvac
879 0 : DO iSpin = 1, SIZE(rho_nosym%vac,4)
880 0 : zlim = MERGE(fi_nosym%vacuum%nmz,fi_nosym%vacuum%nmzxy,iStar==1)
881 0 : CALL grdchlh(fi_nosym%vacuum%delz, REAL(vTot_nosym%vac(:zlim,iStar,iVac,1)),dr_re(:zlim),drr_dummy(:zlim))
882 0 : CALL grdchlh(fi_nosym%vacuum%delz,AIMAG(vTot_nosym%vac(:zlim,iStar,iVac,1)),dr_im(:zlim),drr_dummy(:zlim))
883 0 : grVtotvac(:,iStar,iVac,3) = (3-2*iVac)*(dr_re + ImagUnit * dr_im)
884 : END DO
885 : END DO
886 : END DO
887 : END DO
888 : END IF
889 :
890 : IF (.FALSE.) THEN!!!!! Test grRho/grVC on real space
891 : DO iy = 0, grid(2)-1
892 : DO ix = 0, grid(1)-1
893 : point_plus = fi_nosym%cell%amat(:,1)*REAL(ix)/(grid(1)-1) + &
894 : fi_nosym%cell%amat(:,2)*REAL(iy)/(grid(2)-1) + &
895 : [0.0,0.0,fi_nosym%cell%z1]
896 :
897 : point_minus = fi_nosym%cell%amat(:,1)*REAL(ix)/(grid(1)-1) + &
898 : fi_nosym%cell%amat(:,2)*REAL(iy)/(grid(2)-1) - &
899 : [0.0,0.0,fi_nosym%cell%z1]! - &
900 : !atom_shift
901 :
902 : ! IR rho:
903 : CALL outcdn(point_plus,1,0,0,2,1,.FALSE.,stars_nosym,&
904 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
905 : rho_nosym,xdnout_grrho_up_pw)
906 :
907 : ! Vac rho:
908 : CALL outcdn(point_plus,1,0,1,0,1,.FALSE.,stars_nosym,&
909 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
910 : rho_nosym,xdnout_grrho_up_vac)
911 :
912 : ! IR rho:
913 : CALL outcdn(point_minus,1,0,0,2,1,.FALSE.,stars_nosym,&
914 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
915 : rho_nosym,xdnout_grrho_down_pw)
916 :
917 : ! Vac rho:
918 : CALL outcdn(point_minus,1,0,2,0,1,.FALSE.,stars_nosym,&
919 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
920 : rho_nosym,xdnout_grrho_down_vac)
921 :
922 : write(5395,*) "Gridx/y:", ix, iy
923 : write(5395,*) "Upper rho:", xdnout_grrho_up_vac, xdnout_grrho_up_pw
924 : write(5395,*) "Lower rho:", xdnout_grrho_down_vac,xdnout_grrho_down_pw
925 :
926 : ! IR grrho:
927 : CALL outcdn(point_plus,1,0,0,2,1,.FALSE.,stars_nosym,&
928 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
929 : grRho3(3),xdnout_grrho_up_pw)
930 : ! IR grvc:
931 : CALL outcdn(point_plus,1,0,0,2,1,.FALSE.,stars_nosym,&
932 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
933 : grVC3(3),xdnout_grvc_up_pw)
934 :
935 : ! IR grrho:
936 : CALL outcdn(point_minus,1,0,0,2,1,.FALSE.,stars_nosym,&
937 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
938 : grRho3(3),xdnout_grrho_down_pw)
939 : ! IR grvc:
940 : CALL outcdn(point_minus,1,0,0,2,1,.FALSE.,stars_nosym,&
941 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
942 : grVC3(3),xdnout_grvc_down_pw)
943 :
944 : ! Vac grrho:
945 : CALL outcdn(point_plus,1,0,1,0,1,.FALSE.,stars_nosym,&
946 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
947 : grRho3(3),xdnout_grrho_up_vac)
948 : ! Vac grvc:
949 : CALL outcdn(point_plus,1,0,1,0,1,.FALSE.,stars_nosym,&
950 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
951 : grVC3(3),xdnout_grvc_up_vac)
952 :
953 : ! Vac grrho:
954 : CALL outcdn(point_minus,1,0,2,0,1,.FALSE.,stars_nosym,&
955 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
956 : grRho3(3),xdnout_grrho_down_vac)
957 : ! Vac grvc:
958 : CALL outcdn(point_minus,1,0,2,0,1,.FALSE.,stars_nosym,&
959 : fi_nosym%vacuum,sphhar_nosym,fi_nosym%atoms,fi_nosym%sym,fi_nosym%cell ,&
960 : grVC3(3),xdnout_grvc_down_vac)
961 :
962 : write(5395,*) "Upper grrho:", xdnout_grrho_up_vac,xdnout_grrho_up_pw
963 : write(5395,*) "Lower grrho:", xdnout_grrho_down_vac,xdnout_grrho_down_pw
964 : write(5395,*) "Upper grvc: ", xdnout_grvc_up_vac,xdnout_grvc_up_pw
965 : write(5395,*) "Lower grvc: ", xdnout_grvc_down_vac,xdnout_grvc_down_pw
966 : END DO !x-loop
967 : END DO !y-loop
968 : END IF!!!!!
969 :
970 0 : IF (fi_nosym%input%film)CALL save_npy("rhovac.npy",rho_nosym%vac(:,:,:,1))
971 0 : IF (fi_nosym%input%film)CALL save_npy("rhogr1vac.npy",grRho3(1)%vac(:,:,:,1))
972 0 : IF (fi_nosym%input%film)CALL save_npy("rhogr2vac.npy",grRho3(2)%vac(:,:,:,1))
973 0 : IF (fi_nosym%input%film)CALL save_npy("rhogr3vac.npy",grRho3(3)%vac(:,:,:,1))
974 0 : CALL save_npy("rhogr3pw.npy",grRho3(3)%pw(:,1))
975 0 : IF (fi_nosym%input%film)CALL save_npy("vcgr1vac.npy",grVC3(1)%vac(:,:,:,1))
976 0 : IF (fi_nosym%input%film)CALL save_npy("vcgr2vac.npy",grVC3(2)%vac(:,:,:,1))
977 0 : IF (fi_nosym%input%film)CALL save_npy("vcgr3vac.npy",grVC3(3)%vac(:,:,:,1))
978 0 : IF (fi_nosym%input%film)CALL save_npy("vextgr1vac.npy",grVext3(1)%vac(:,:,:,1))
979 0 : IF (fi_nosym%input%film)CALL save_npy("vextgr2vac.npy",grVext3(2)%vac(:,:,:,1))
980 0 : IF (fi_nosym%input%film)CALL save_npy("vextgr3vac.npy",grVext3(3)%vac(:,:,:,1))
981 0 : CALL save_npy("vextgr1pw.npy",grVext3(1)%pw(:,1))
982 0 : CALL save_npy("vextgr2pw.npy",grVext3(2)%pw(:,1))
983 0 : CALL save_npy("vextgr3pw.npy",grVext3(3)%pw(:,1))
984 0 : IF (fi_nosym%input%film)CALL save_npy("vtotgr1vac.npy",grVtot3(1)%vac(:,:,:,1))
985 0 : IF (fi_nosym%input%film)CALL save_npy("vtotgr2vac.npy",grVtot3(2)%vac(:,:,:,1))
986 0 : IF (fi_nosym%input%film)CALL save_npy("vtotgr3vac.npy",grVtot3(3)%vac(:,:,:,1))
987 0 : IF (fi_nosym%input%film)CALL save_npy("vtotgr1vacnum.npy",grVtotvac(:,:,:,1))
988 0 : IF (fi_nosym%input%film)CALL save_npy("vtotgr2vacnum.npy",grVtotvac(:,:,:,2))
989 0 : IF (fi_nosym%input%film)CALL save_npy("vtotgr3vacnum.npy",grVtotvac(:,:,:,3))
990 0 : CALL save_npy("vtotgr1pw.npy",grVtot3(1)%pw(:,1))
991 0 : CALL save_npy("vtotgr2pw.npy",grVtot3(2)%pw(:,1))
992 0 : CALL save_npy("vtotgr3pw.npy",grVtot3(3)%pw(:,1))
993 0 : CALL save_npy("vtotgr1pwnum.npy",grVtotpw(:,1))
994 0 : CALL save_npy("vtotgr2pwnum.npy",grVtotpw(:,2))
995 0 : CALL save_npy("vtotgr3pwnum.npy",grVtotpw(:,3))
996 0 : END SUBROUTINE
997 :
998 : END MODULE m_dfpt
|
