Line data Source code
1 : !!****m* ABINIT/m_rf2_init
2 : !! NAME
3 : !! m_rf2_init
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2015-2026 ABINIT group (LB,MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! TODO
15 : !! Can be merged with m_rf2
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : module m_rf2_init
26 :
27 : use defs_basis
28 : use m_xmpi
29 : use m_errors
30 : use m_wfk
31 : use m_hamiltonian
32 : use m_cgtools
33 : use m_rf2
34 : use m_abicore
35 : use m_dtset
36 : use m_dtfil
37 :
38 : use m_time , only : timab
39 : use defs_abitypes, only : MPI_type
40 : use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_copy,pawcprj_get,pawcprj_free,pawcprj_output
41 : use m_cgprj, only : getcprj
42 :
43 : implicit none
44 :
45 : private
46 : !!***
47 :
48 : public :: rf2_init
49 : !!***
50 :
51 : contains
52 : !!***
53 :
54 : !!****f* ABINIT/rf2_init
55 : !!
56 : !! NAME
57 : !! rf2_init
58 : !!
59 : !! FUNCTION
60 : !! Compute terms needed for the 2nd order Sternheimer equation.
61 : !! All terms are stored in a rf2_t object.
62 : !!
63 : !! INPUTS
64 : !! cg(2,mpw*nspinor*mband*nsppol)=planewave coefficients of wavefunctions at k
65 : !! cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
66 : !! projected with non-local projectors: cprj=<p_i|Cnk>
67 : !! rf2 : the object we want to initialize (see m_rf2.F90 for more information)
68 : !! dtfil <type(datafiles_type)>=variables related to files
69 : !! dtset <type(dataset_type)>=all input variables for this dataset
70 : !! eig0_k(mband*nsppol)=GS eigenvalues at k (hartree)
71 : !! eig1_k(2*mband*mband*nsppol)=2nd-order eigenvalues at k,q (hartree)
72 : !! ffnl1=nonlocal form factors
73 : !! ffnl1_test=nonlocal form factors used for tests (i.e. when dtset%nonlinear_info>2)
74 : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
75 : !! ibg=shift to be applied on the location of data in the array cprj
76 : !! icg=shift to be applied on the location of data in the array cg
77 : !! idir=direction of the perturbation
78 : !! ikpt=number of the k-point
79 : !! ipert=type of the perturbation
80 : !! isppol=index of current spin component
81 : !! mkmem =number of k points trated by this node (GS data).
82 : !! mpi_enreg=information about MPI parallelization
83 : !! mpw=maximum dimensioned size of npw or wfs at k
84 : !! nband_k=number of bands at this k point for that spin polarization
85 : !! ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
86 : !! nspinor=number of spinorial components of the wavefunctions
87 : !! nsppol=1 for unpolarized, 2 for spin-polarized
88 : !! rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
89 : !! rf_hamk_dir2 <type(rf_hamiltonian_type)>= (used only when ipert=natom+11, so q=0)
90 : !! same as rf_hamkq, but the direction of the perturbation is different
91 : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k.
92 : !! rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
93 : !! ddk_f<wfk_t>=struct info for DDK file.
94 : !!
95 : !! OUTPUT
96 : !! rf2%RHS_Stern
97 : !!
98 : !! NOTES
99 : !!
100 : !! SOURCE
101 :
102 14616 : subroutine rf2_init(cg,cprj,rf2,dtset,dtfil,eig0_k,eig1_k,ffnl1,ffnl1_test,gs_hamkq,ibg,icg,idir,ikpt,ipert,isppol,mkmem,&
103 14616 : mpi_enreg,mpw,nband_k,nsppol,rf_hamkq,rf_hamk_dir2,occ_k,rocceig,ddk_f)
104 :
105 : ! *************************************************************************
106 : !Arguments -------------------------------
107 : !scalars
108 : integer,intent(in) :: ibg,icg,idir,ipert,isppol,ikpt
109 : integer,intent(in) :: mkmem,mpw,nband_k,nsppol
110 : type(datafiles_type),intent(in) :: dtfil
111 : type(dataset_type),intent(in) :: dtset
112 : type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
113 : type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq,rf_hamk_dir2
114 : type(MPI_type),intent(in) :: mpi_enreg
115 :
116 : !arrays
117 : real(dp),intent(in),target :: cg(2,mpw*gs_hamkq%nspinor*dtset%mband*mkmem*nsppol)
118 : real(dp),intent(in) :: eig0_k(dtset%mband)
119 : real(dp),intent(inout) :: eig1_k(2*dtset%mband**2) ! Here eig1_k contains 2nd order eigenvalues...
120 : real(dp),intent(in) :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
121 : real(dp),intent(in) :: occ_k(nband_k),rocceig(nband_k,nband_k)
122 : type(pawcprj_type),intent(in) :: cprj(gs_hamkq%natom,gs_hamkq%nspinor*dtset%mband*mkmem*nsppol*gs_hamkq%usecprj)
123 : type(rf2_t),intent(inout) :: rf2
124 : type(wfk_t),intent(inout) :: ddk_f(4)
125 : !
126 : !Local variables-------------------------------
127 : !scalars
128 : integer,parameter :: iorder_cprj=0
129 : integer :: choice_cprj,cpopt_cprj,iband,icpgr_loc,idir1,idir2,idir_cprj,ierr
130 : integer :: indb,ipert1,ipert2,iproc,jband,kdir1
131 : integer :: me,my_nband,natom,ncpgr_loc,nproc_band,debug_mode
132 : integer :: size_cprj,size_wf,shift_band1,shift_band2,shift_cprj_band1,shift_cprj_dir1,shift_proc
133 : integer :: shift_dir1_lambda,shift_dir2_lambda,shift_dir1,shift_dir1_loc,shift_dir2,shift_jband_lambda
134 : logical :: has_cprj_jband,has_dudkprj
135 : real(dp) :: doti,dotr,dot2i,dot2r,invocc,tol_final
136 : real(dp) :: factor,factor_in
137 : logical :: symmetric,antisymmetric,total
138 : character(len=500) :: msg
139 : !arrays
140 : integer :: file_index(2)
141 : real(dp) :: lambda_ij(2),tsec(2)
142 14616 : real(dp),allocatable :: cg_jband(:,:,:),ddk_read(:,:),dudkdk(:,:),dudk_dir2(:,:)
143 14616 : real(dp),allocatable :: eig1_read(:),gvnlx1(:,:),h_cwave(:,:),s_cwave(:,:),dsusdu_loc(:,:),dsusdu_gather(:,:)
144 14616 : real(dp),allocatable,target :: dsusdu(:,:),dudk(:,:),eig1_k_stored(:)
145 14616 : real(dp), ABI_CONTIGUOUS pointer :: cwave_dudk(:,:),cwave_i(:,:),cwave_j(:,:),eig1_k_jband(:)
146 14616 : real(dp),pointer :: rhs_j(:,:)
147 14616 : type(pawcprj_type),target :: cprj_empty(0,0)
148 14616 : type(pawcprj_type),allocatable,target :: cprj_jband(:,:),dudkprj(:,:)
149 14616 : type(pawcprj_type),pointer :: cprj_dudk(:,:),cprj_j(:,:)
150 : type(rf_hamiltonian_type),pointer :: rf_hamk_idir
151 :
152 : ! *********************************************************************
153 :
154 :
155 : DBG_ENTER("COLL")
156 :
157 14616 : call timab(514,1,tsec)
158 :
159 : !my mpi rank :
160 14616 : me=mpi_enreg%me_kpt
161 :
162 14616 : size_wf=gs_hamkq%npw_k*gs_hamkq%nspinor
163 14616 : size_cprj=gs_hamkq%nspinor
164 14616 : natom = gs_hamkq%natom
165 14616 : debug_mode = 0
166 14616 : if (dtset%nonlinear_info>2) debug_mode = 1 ! also active a lot of tests
167 :
168 : !Define some attributes of the rf2 object
169 14616 : rf2%nband_k = nband_k
170 14616 : rf2%size_wf = size_wf
171 14616 : rf2%size_cprj = size_cprj
172 :
173 14616 : if(ipert<natom+10.or.ipert>natom+11) then
174 0 : write(msg,'(a)') 'ipert must be equal to natom+10 or natom+11 for rf2 calculations.'
175 0 : ABI_BUG(msg)
176 : end if
177 :
178 : !Define perturbations and idirs
179 14616 : rf2%iperts(1) = natom+1
180 14616 : rf2%iperts(2) = natom+1
181 14616 : if (ipert==natom+11) rf2%iperts(2) = natom+2
182 :
183 14616 : if (ipert==natom+10.and.idir<=3) then ! One perturbation, one direction
184 4104 : rf2%ndir=1
185 4104 : rf2%idirs(1)=idir ; rf2%idirs(2)=idir
186 : else ! Two perturbations or/and two directions
187 10512 : rf2%ndir=2
188 10512 : call rf2_getidirs(idir,idir1,idir2)
189 10512 : rf2%idirs(1)=idir1
190 10512 : rf2%idirs(2)=idir2
191 : end if
192 :
193 : !Choose dkdk option
194 14616 : symmetric = .false.
195 14616 : antisymmetric = .false.
196 14616 : total = .false.
197 :
198 : ! even if rf2_dkdk is 5 (nonesense),
199 : ! as long as rf2_dkdk /= 0, symmetric = .true.
200 14616 : if (dtset%rf2_dkdk==2) then
201 : antisymmetric = .true.
202 14580 : else if (dtset%rf2_dkdk==3) then
203 : total = .true.
204 : else
205 5760 : symmetric = .true.
206 : endif
207 :
208 :
209 : ! **************************************************************************************************
210 : ! Get info from ddk files
211 : ! **************************************************************************************************
212 :
213 : !Allocate work spaces
214 43848 : ABI_MALLOC(eig1_read,(2*nband_k))
215 43848 : ABI_MALLOC(ddk_read,(2,size_wf))
216 146936 : eig1_read(:)=zero
217 5230020 : ddk_read(:,:)=zero
218 :
219 : ! "eig1_k_stored" contains dLambda_{nm}/dpert every bands n and m and ndir (=1 or 2) directions
220 : ! pert = k_dir (wavevector) or E_dir (electric field)
221 43848 : ABI_MALLOC(eig1_k_stored,(2*rf2%ndir*nband_k**2))
222 1263768 : eig1_k_stored=zero
223 :
224 : ! "dudk" contains du/dpert1 for every bands and ndir (=1 or 2) directions
225 43848 : ABI_MALLOC_OR_DIE(dudk,(2,rf2%ndir*nband_k*size_wf), ierr)
226 41021316 : dudk=zero
227 14616 : has_dudkprj=.false.
228 14616 : if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
229 52224 : ABI_MALLOC(dudkprj,(natom,rf2%ndir*nband_k*size_cprj))
230 960 : ncpgr_loc=1;if(ipert==natom+10.or.ipert==natom+11) ncpgr_loc=3
231 960 : call pawcprj_alloc(dudkprj,ncpgr_loc,gs_hamkq%dimcprj)
232 960 : choice_cprj=5 ; cpopt_cprj=0
233 960 : has_dudkprj=.true.
234 : else
235 27312 : ABI_MALLOC(dudkprj,(natom,0))
236 : end if
237 :
238 14616 : if (debug_mode/=0) then
239 240 : write(msg,'(4(a,i2))') 'RF2_INIT : ipert-natom = ',ipert-natom,' , idir = ',idir,&
240 480 : ' , ikpt = ',ikpt,' , isppol = ',isppol
241 240 : call wrtout(std_out,msg,'COLL')
242 : end if
243 :
244 14616 : file_index(1)=1 ! dir1
245 14616 : file_index(2)=2 ! dir2
246 14616 : if (ipert==natom+11) then ! see dfpt_looppert.F90
247 3456 : file_index(1)=3 ! dir1
248 : file_index(2)=2 ! dir2
249 : end if
250 :
251 39744 : do kdir1=1,rf2%ndir
252 25128 : idir1=rf2%idirs(kdir1)
253 25128 : ipert1=rf2%iperts(kdir1)
254 153936 : do iband=1,nband_k
255 114192 : call ddk_f(file_index(kdir1))%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read)
256 : ! Copy ddk_read in "dudk"
257 114192 : shift_band1=(iband-1)*size_wf
258 114192 : shift_dir1=(kdir1-1)*nband_k*size_wf
259 41120892 : dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=ddk_read(:,:)
260 : ! Copy eig1_read in "eig1_k_stored"
261 114192 : shift_band1=(iband-1)*2*nband_k
262 114192 : shift_dir1_lambda=2*(kdir1-1)*nband_k**2
263 1363344 : eig1_k_stored(1+shift_band1+shift_dir1_lambda:2*nband_k+shift_band1+shift_dir1_lambda)=eig1_read(:)
264 : ! Get this dudk projected on NL projectors
265 139320 : if (has_dudkprj.and.mpi_enreg%proc_distrb(ikpt,iband,isppol)==me) then
266 16128 : shift_cprj_band1=(iband-1)*size_cprj
267 16128 : shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj
268 : cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1: &
269 16128 : & size_cprj+shift_cprj_band1+shift_cprj_dir1)
270 16128 : idir_cprj=0;if (dudkprj(1,1)%ncpgr/=3) idir_cprj=idir1
271 : call getcprj(choice_cprj,cpopt_cprj,ddk_read,cprj_dudk,gs_hamkq%ffnl_k,idir_cprj,&
272 : & gs_hamkq%indlmn,gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kpg_k,gs_hamkq%kpt_k,&
273 : & gs_hamkq%lmnmax,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,&
274 : & gs_hamkq%nloalg,gs_hamkq%npw_k,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkxred,&
275 16128 : & gs_hamkq%ph1d,gs_hamkq%ph3d_k,gs_hamkq%ucvol,gs_hamkq%useylm)
276 : end if
277 : end do
278 : end do
279 :
280 14616 : ABI_MALLOC(dudkdk,(2,0))
281 14616 : ABI_MALLOC(dudk_dir2,(2,0))
282 :
283 : !Get dudkdk for ipert==natom+11
284 14616 : if(ipert==natom+11) then
285 3456 : ABI_FREE(dudkdk)
286 10368 : ABI_MALLOC(dudkdk,(2,nband_k*size_wf))
287 3456 : if (idir>3) then
288 2304 : ABI_FREE(dudk_dir2)
289 4608 : ABI_MALLOC(dudk_dir2,(2,nband_k*size_wf))
290 : end if
291 21120 : do iband=1,nband_k
292 17664 : call ddk_f(1)%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read)
293 17664 : shift_band1=(iband-1)*size_wf
294 3767424 : dudkdk(:,1+shift_band1:size_wf+shift_band1)=ddk_read(:,:)
295 : ! Check that < u^(0) | u^(2) > = - Re[< u^(1) | u^(1) >]
296 17664 : if (debug_mode/=0 .and. gs_hamkq%usepaw==0) then
297 : ! Compute < u^(0) | u^(2) >
298 0 : do jband=1,nband_k
299 0 : cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
300 0 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,ddk_read,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
301 0 : if (idir<=3 .and. iband==jband .and. abs(occ_k(iband))>tol8) then
302 : ! Compute < u^(1) | u^(1) > = Re[< u^(1) | u^(1) >]
303 0 : cwave_dudk => dudk(:,1+shift_band1:size_wf+shift_band1)
304 0 : call sqnorm_g(dot2r,gs_hamkq%istwf_k,size_wf,cwave_dudk,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
305 0 : dotr = dotr + dot2r
306 0 : dotr = sqrt(dotr**2+doti**2)
307 0 : if (dotr > tol7) then
308 0 : write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST dudkdk iband = ',iband,&
309 0 : ' : NOT PASSED. | < u^(0) | u^(2) > + Re[< u^(1) | u^(1) >] | = ',dotr
310 0 : call wrtout(std_out,msg)
311 0 : call wrtout(ab_out,msg)
312 : ! else
313 : ! write(msg,'(a,i2,a)') 'RF2 TEST dudkdk iband = ',iband,' : OK.'
314 : ! call wrtout(std_out,msg)
315 : end if
316 : end if ! idir<=3
317 : end do ! jband
318 : end if ! debug_mode
319 : ! Read ddk for idir2
320 21120 : if (idir>3) then
321 11776 : call ddk_f(4)%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read)
322 2511616 : dudk_dir2(:,1+shift_band1:size_wf+shift_band1)=ddk_read(:,:)
323 : end if
324 : end do !iband
325 : end if ! ipert=natom+11
326 :
327 14616 : ABI_FREE(ddk_read)
328 14616 : ABI_FREE(eig1_read)
329 :
330 : ! **************************************************************************************************
331 : ! COMPUTATION OF "dsusdu", A PART OF "A_mn" AND A PART OF "Lambda_mn" (see defs in m_rf2)
332 : ! **************************************************************************************************
333 :
334 : !Allocate work spaces for one band
335 29232 : ABI_MALLOC(h_cwave,(2,size_wf))
336 29232 : ABI_MALLOC(s_cwave,(2,size_wf))
337 29232 : ABI_MALLOC(gvnlx1,(2,size_wf))
338 5230020 : h_cwave(:,:) = zero
339 5230020 : s_cwave(:,:) = zero
340 5230020 : gvnlx1(:,:) = zero
341 :
342 : ! "dsusdu" contains dS/dpert_dir |u_band> + S|du_band/dpert1> for every bands and ndir (=1 or 2) directions
343 43848 : ABI_MALLOC_OR_DIE(dsusdu,(2,rf2%ndir*nband_k*size_wf), ierr)
344 41021316 : dsusdu=zero
345 :
346 43848 : ABI_MALLOC(rf2%amn,(2,nband_k**2))
347 1094616 : rf2%amn=zero
348 :
349 29232 : ABI_MALLOC(rf2%lambda_mn,(2,nband_k**2))
350 1094616 : rf2%lambda_mn(:,:)=zero
351 :
352 : !Allocate work spaces when debug_mode is activated
353 14616 : has_cprj_jband=.false.
354 14616 : if (debug_mode/=0) then ! Only for test purposes
355 960 : ABI_MALLOC(cg_jband,(2,size_wf*nband_k,2))
356 814560 : cg_jband(:,:,1) = cg(:,1+icg:size_wf*nband_k+icg)
357 240 : if (ipert==natom+11) then ! Note the multiplication by "i"
358 144 : if (idir<=3) then
359 54336 : cg_jband(1,:,2) = -dudk(2,1:size_wf*nband_k) ! for dir1
360 54336 : cg_jband(2,:,2) = dudk(1,1:size_wf*nband_k) ! for dir1
361 : else
362 108672 : cg_jband(1,:,2) = -dudk_dir2(2,1:size_wf*nband_k) ! for dir2
363 108672 : cg_jband(2,:,2) = dudk_dir2(1,1:size_wf*nband_k) ! for dir2
364 : end if
365 : end if
366 240 : if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
367 9600 : ABI_MALLOC(cprj_jband,(natom,size_cprj*nband_k))
368 240 : has_cprj_jband=.true.
369 : else
370 0 : ABI_MALLOC(cprj_jband,(natom,0))
371 : end if
372 : else
373 14376 : ABI_MALLOC(cg_jband,(2,0,2))
374 28752 : ABI_MALLOC(cprj_jband,(natom,0))
375 : end if
376 :
377 : ! define dkdk factor
378 14616 : factor = one
379 14616 : if(ipert==natom+10 .and. idir<=3 .and. symmetric) factor=two ! in order to not compute same terms twice
380 :
381 39744 : do kdir1=1,rf2%ndir
382 : ! First iteration (kdir1=1) :
383 : ! pert1 = rf2%iperts(1) along rf2%idirs(1)
384 : ! pert2 = rf2%iperts(2) along rf2%idirs(2)
385 : ! Second iteration (kdir1=2) :
386 : ! pert1 = rf2%iperts(2) along rf2%idirs(2)
387 : ! pert2 = rf2%iperts(1) along rf2%idirs(1)
388 25128 : idir1=rf2%idirs(kdir1)
389 25128 : ipert1=rf2%iperts(kdir1)
390 25128 : shift_dir1=(kdir1-1)*nband_k*size_wf
391 25128 : shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj
392 25128 : shift_dir1_lambda=(kdir1-1)*2*nband_k**2
393 25128 : if(ipert==natom+10 .and. idir<=3) then
394 4104 : shift_dir2=0
395 4104 : idir2 = idir1
396 4104 : ipert2 = ipert1
397 4104 : rf_hamk_idir => rf_hamkq
398 : else
399 21024 : shift_dir2=(2-kdir1)*nband_k*size_wf
400 21024 : idir2 = rf2%idirs(3-kdir1)
401 21024 : ipert2 = rf2%iperts(3-kdir1)
402 21024 : if (kdir1==1) rf_hamk_idir => rf_hamkq
403 21024 : if (kdir1==2) rf_hamk_idir => rf_hamk_dir2
404 : end if
405 :
406 : ! define zero factor instead of exiting the loop: we need to compute dsusdu
407 25128 : if (antisymmetric .and. rf2%ndir==1) factor = zero
408 25128 : if (total .and. rf2%ndir==2 .and. kdir1==2) factor = zero
409 :
410 : ! define factor for antisymmetric dkdk case
411 25128 : if (rf2%ndir==2 .and. kdir1==2 .and. antisymmetric) factor = -one
412 :
413 : ! define factor_in for accumulate_bands
414 25128 : factor_in = factor
415 25128 : if ((gs_hamkq%usepaw==1 .or. ipert/=natom+10) .and. rf2%ndir==1) factor_in = two
416 :
417 : ! Load projected WF according to ipert1 and idir1
418 25128 : cprj_j => cprj_empty ; cprj_dudk => cprj_empty
419 25128 : if (has_cprj_jband) then
420 432 : call pawcprj_free(cprj_jband)
421 432 : ncpgr_loc= 3;if(ipert1==natom+1.or.ipert1==natom+2) ncpgr_loc=1
422 432 : icpgr_loc=-1;if(ipert1==natom+1.or.ipert1==natom+2) icpgr_loc=idir1
423 432 : call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj)
424 : call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,&
425 : & isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,&
426 432 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc)
427 : end if
428 :
429 : ! LOOP OVER BANDS
430 153936 : do jband=1,nband_k ! = band n
431 :
432 : ! Skip bands not treated by current proc
433 114192 : if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
434 :
435 114192 : shift_band1=(jband-1)*size_wf
436 114192 : shift_cprj_band1=(jband-1)*size_cprj
437 114192 : shift_jband_lambda=(jband-1)*2*nband_k
438 :
439 139320 : if (abs(occ_k(jband))>tol8) then
440 :
441 : ! Extract first order wavefunction and eigenvalues for jband
442 102672 : eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir1_lambda:2*nband_k+shift_jband_lambda+shift_dir1_lambda)
443 102672 : cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
444 102672 : if (has_dudkprj) cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1:size_cprj+shift_cprj_band1+shift_cprj_dir1)
445 :
446 : ! Compute H^(0) | du/dpert1 > (in h_cwave) and S^(0) | du/dpert1 > (in s_cwave)
447 : call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_dudk,cprj_dudk,h_cwave,s_cwave,&
448 : & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnlx1,0,0,ikpt,isppol,mkmem,&
449 102672 : & mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf)
450 :
451 36022140 : if (gs_hamkq%usepaw==0) s_cwave(:,:)=cwave_dudk(:,:) ! Store | du/dpert1 > in s_cwave
452 :
453 : ! Copy infos in dsusdu
454 : dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=s_cwave(:,:)&
455 37976508 : +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
456 :
457 102672 : if (debug_mode/=0) then
458 1728 : write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 1 kdir1 = ',kdir1,' jband = ',jband
459 1728 : call wrtout(std_out,msg)
460 : end if
461 :
462 : ! For every occupied iband, we compute :
463 : ! < du/dpert2(iband) | H^(0) | du/dpert1(jband) > and add it to lambda_mn
464 : ! < du/dpert2(iband) | S^(0) | du/dpert1(jband) > and add it to amn
465 607440 : do iband=1,rf2%nband_k ! = band m
466 607440 : if (abs(occ_k(iband))>tol8) then
467 458688 : shift_band2=(iband-1)*size_wf
468 458688 : cwave_dudk => dudk(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2)
469 : call rf2_accumulate_bands(rf2,1,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
470 458688 : jband,debug_mode,cwave_dudk,h_cwave,s_cwave,factor_in)
471 : end if
472 : end do
473 :
474 : ! Extract GS wavefunction for jband
475 102672 : cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
476 102672 : if(has_cprj_jband) cprj_j => cprj_jband(:,1+shift_cprj_band1:size_cprj+shift_cprj_band1)
477 :
478 102672 : if (ipert1==natom+2) then
479 : ! Extract ddk and multiply by i :
480 13824 : if(idir<=3) then ! in this case : idir1=idir2
481 305216 : gvnlx1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1)
482 305216 : gvnlx1(2,:) = dudk(1,1+shift_band1:size_wf+shift_band1)
483 : else
484 610432 : gvnlx1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1)
485 610432 : gvnlx1(2,:) = dudk_dir2(1,1+shift_band1:size_wf+shift_band1)
486 : end if
487 : end if
488 :
489 : ! Compute dH/dpert1 | u^(0) > (in h_cwave) and dS/dpert1 | u^(0) > (in s_cwave)
490 : call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_j,cprj_j,h_cwave,s_cwave,&
491 : & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnlx1,idir1,ipert1,ikpt,isppol,&
492 102672 : & mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf)
493 :
494 : ! Copy infos in dsusdu
495 102672 : if (gs_hamkq%usepaw==1) then
496 : dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=s_cwave(:,:)&
497 1961280 : +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
498 : end if
499 :
500 102672 : if (debug_mode/=0) then
501 1728 : write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 2 kdir1 = ',kdir1,' jband = ',jband
502 1728 : call wrtout(std_out,msg)
503 : end if
504 :
505 : ! For every occupied iband, we compute :
506 : ! < du/dpert2(iband) | dH/dpert1 | u^(0)(jband) > and add it to lambda_mn
507 : ! < du/dpert2(iband) | dS/dpert1 | u^(0)(jband) > and add it to amn
508 607440 : do iband=1,rf2%nband_k ! = band m
509 607440 : if (abs(occ_k(iband))>tol8) then
510 458688 : shift_band2=(iband-1)*size_wf
511 458688 : cwave_dudk => dudk(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2)
512 : call rf2_accumulate_bands(rf2,2,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
513 458688 : jband,debug_mode,cwave_dudk,h_cwave,s_cwave,factor_in)
514 : end if
515 : end do
516 :
517 : end if ! empty band test
518 : end do ! jband
519 : end do ! idir1
520 :
521 : ! Allgather dsusdu
522 14616 : nproc_band = xmpi_comm_size(mpi_enreg%comm_band)
523 14616 : if (nproc_band>1) then
524 :
525 0 : my_nband = nband_k/nproc_band;if (mod(nband_k,nproc_band)/=0) my_nband=my_nband+1
526 0 : ABI_MALLOC(dsusdu_loc,(2,size_wf*my_nband*rf2%ndir))
527 0 : ABI_MALLOC(dsusdu_gather,(2,size_wf*my_nband*rf2%ndir*nproc_band))
528 0 : dsusdu_loc(:,:) = zero
529 0 : dsusdu_gather(:,:) = zero
530 :
531 0 : do kdir1=1,rf2%ndir
532 0 : indb = 1
533 0 : shift_dir1=(kdir1-1)*size_wf*nband_k
534 0 : shift_dir1_loc=(kdir1-1)*size_wf*my_nband
535 0 : do jband=1,nband_k
536 : ! Skip bands not treated by current proc
537 0 : if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
538 :
539 0 : shift_band1=(jband-1)*size_wf
540 : dsusdu_loc(:,indb+shift_dir1_loc:indb-1+size_wf+shift_dir1_loc) = &
541 0 : dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
542 0 : indb = indb + size_wf
543 : end do
544 : end do
545 :
546 0 : call xmpi_allgather(dsusdu_loc,2*size_wf*my_nband*rf2%ndir,dsusdu_gather,mpi_enreg%comm_band,ierr)
547 :
548 0 : do kdir1=1,rf2%ndir
549 0 : shift_dir1=(kdir1-1)*size_wf*nband_k
550 0 : shift_dir1_loc=(kdir1-1)*size_wf*my_nband
551 0 : do iproc=1,nproc_band
552 0 : shift_proc = (iproc-1)*size_wf*my_nband*rf2%ndir
553 0 : indb = 1
554 0 : do jband=1,my_nband
555 0 : iband = jband+(iproc-1)*my_nband
556 0 : if(iband<=nband_k) then
557 0 : shift_band1=(iband-1)*size_wf
558 : dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) = &
559 0 : dsusdu_gather(:,indb+shift_dir1_loc+shift_proc:indb-1+size_wf+shift_dir1_loc+shift_proc)
560 : end if
561 0 : indb = indb + size_wf
562 : end do
563 : end do
564 : end do
565 0 : ABI_FREE(dsusdu_loc)
566 0 : ABI_FREE(dsusdu_gather)
567 :
568 : end if
569 :
570 : ! **************************************************************************************************
571 : ! COMPUTATION OF "RHS_Stern", THE LAST PART OF "A_mn" AND A PART OF "Lambda_mn"
572 : ! **************************************************************************************************
573 :
574 43848 : ABI_MALLOC_OR_DIE(rf2%RHS_Stern,(2,nband_k*size_wf), ierr)
575 24118668 : rf2%RHS_Stern(:,:)=zero
576 :
577 : ! Define prefactor for the H^{(2)} term
578 14616 : factor = one
579 14616 : if (total) factor = half
580 :
581 : !Computation of terms containing H^(2)
582 14616 : if (ipert/=natom+11 .or. gs_hamkq%usepaw==1) then ! Otherwise H^(2) = 0
583 11736 : if (.not. antisymmetric) then ! check antisymmetric dkdk
584 :
585 : !Load projected WF according to ipert and idir
586 11700 : cprj_j => cprj_empty
587 11700 : if (has_cprj_jband) then
588 240 : call pawcprj_free(cprj_jband)
589 240 : ncpgr_loc= 3;if(ipert==natom+1.or.ipert==natom+2) ncpgr_loc=1
590 240 : icpgr_loc=-1;if(ipert==natom+1.or.ipert==natom+2) icpgr_loc=idir
591 240 : call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj)
592 : call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,&
593 : & isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,&
594 240 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc)
595 : end if
596 :
597 11700 : if (ipert==natom+10) then
598 : rf_hamk_idir => rf_hamkq ! all info are in rf_hamkq
599 576 : else if (ipert==natom+11) then
600 576 : rf_hamk_idir => rf_hamk_dir2 ! all info are in rf_hamk_dir2
601 : end if
602 :
603 65428 : do jband=1,nband_k
604 :
605 : ! Skip bands not treated by current proc
606 53728 : if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
607 :
608 65428 : if (abs(occ_k(jband))>tol8) then
609 48096 : shift_band1=(jband-1)*size_wf
610 48096 : shift_cprj_band1=(jband-1)*size_cprj
611 48096 : shift_jband_lambda=(jband-1)*2*nband_k
612 :
613 : ! Extract GS wavefunction
614 48096 : cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
615 48096 : if(has_cprj_jband) cprj_j => cprj_jband(:,1+shift_cprj_band1:size_cprj+shift_cprj_band1)
616 :
617 : ! Not used here, but a null pointer is not allowed... (called by rf2_apply_hamiltonian but not used)
618 48096 : eig1_k_jband => eig1_k_stored(1+shift_jband_lambda:2*nband_k+shift_jband_lambda)
619 :
620 48096 : if (ipert==natom+11) then
621 : ! Extract ddk and multiply by i :
622 2304 : if(idir<=3) then ! in this case : idir1=idir2
623 73152 : gvnlx1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1)
624 73152 : gvnlx1(2,:) = dudk(1,1+shift_band1:size_wf+shift_band1)
625 : else
626 146304 : gvnlx1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1)
627 146304 : gvnlx1(2,:) = dudk_dir2(1,1+shift_band1:size_wf+shift_band1)
628 : end if
629 : end if
630 :
631 : ! Compute : d^2H/(dpert1 dpert2)|u^(0)> (in h_cwave)
632 : ! and : d^2S/(dpert1 dpert2)|u^(0)> (in s_cwave)
633 : call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_j,cprj_j,h_cwave,s_cwave,&
634 : & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnlx1,idir,ipert,ikpt,isppol,&
635 : & mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf,&
636 48096 : & ffnl1=ffnl1,ffnl1_test=ffnl1_test)
637 :
638 48096 : if (debug_mode/=0) then
639 960 : write(msg,'(a,i2)') 'RF2 TEST before accumulate_bands choice = 3 jband = ',jband
640 960 : call wrtout(std_out,msg)
641 : end if
642 :
643 : ! For every occupied iband, we compute :
644 : ! < u^(0)(iband) | d^2H/(dpert1 dpert2) | u^(0)(jband) > and add it to lambda_mn
645 : ! < u^(0)(iband) | d^2S/(dpert1 dpert2) | u^(0)(jband) > and add it to amn
646 291808 : do iband=1,rf2%nband_k ! = band m
647 291808 : if (abs(occ_k(iband))>tol8) then
648 221184 : shift_band2=(iband-1)*size_wf
649 221184 : cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
650 221184 : if(ipert == natom+10) then
651 211968 : ipert1 = natom+1
652 211968 : ipert2 = natom+1
653 : else
654 9216 : ipert1 = natom+1
655 9216 : ipert2 = natom+2
656 : end if
657 221184 : call rf2_getidirs(idir,idir1,idir2)
658 : call rf2_accumulate_bands(rf2,3,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
659 221184 : jband,debug_mode,cwave_i,h_cwave,s_cwave,factor)
660 : end if
661 : end do
662 :
663 : ! Add d^2H/(dk_dir1 dk_dir2)|u^(0)> to RHS_Stern :
664 1133856 : if (gs_hamkq%usepaw==1) h_cwave(:,:)=h_cwave(:,:)-eig0_k(jband)*s_cwave(:,:) ! if PAW : we add H^(2)-eps^(0) S^(2)
665 48096 : rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
666 144288 : call cg_zaxpy(size_wf,(/factor,zero/),h_cwave,rhs_j)
667 :
668 : end if ! empty band test
669 : end do ! jband
670 : endif ! check antisymmetric dkdk
671 : end if ! H^(2) exists
672 :
673 : ! define dkdk factor
674 14616 : factor = one
675 14616 : if(ipert==natom+10 .and. idir<=3 .and. symmetric) factor=two ! in order to not compute same terms twice
676 :
677 : !Computation of terms containing H^(1)
678 39732 : do kdir1=1,rf2%ndir
679 : ! First iteration (kdir1=1) :
680 : ! pert1 = rf2%iperts(1) along rf2%idirs(1)
681 : ! pert2 = rf2%iperts(2) along rf2%idirs(2)
682 : ! Second iteration (kdir1=2) :
683 : ! pert1 = rf2%iperts(2) along rf2%idirs(2)
684 : ! pert2 = rf2%iperts(1) along rf2%idirs(1)
685 25128 : shift_dir1=(kdir1-1)*nband_k*size_wf
686 25128 : shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj
687 25128 : shift_dir1_lambda=(kdir1-1)*2*nband_k**2
688 25128 : idir1=rf2%idirs(kdir1)
689 25128 : ipert1=rf2%iperts(kdir1)
690 25128 : if(ipert==natom+10 .and. idir<=3) then
691 4104 : idir2=idir1
692 4104 : ipert2=ipert1
693 4104 : shift_dir2=0
694 4104 : shift_dir2_lambda=0
695 4104 : rf_hamk_idir => rf_hamkq
696 : else
697 21024 : idir2=rf2%idirs(2-kdir1+1)
698 21024 : ipert2=rf2%iperts(2-kdir1+1)
699 21024 : shift_dir2=(2-kdir1)*nband_k*size_wf
700 21024 : shift_dir2_lambda=(2-kdir1)*2*nband_k**2
701 21024 : if (kdir1==1) rf_hamk_idir => rf_hamk_dir2 ! dir2
702 21024 : if (kdir1==2) rf_hamk_idir => rf_hamkq ! dir1
703 : end if
704 :
705 : ! determine when we need to compute terms according to dkdk option
706 25128 : if (antisymmetric .and. rf2%ndir==1) exit
707 25116 : if (total .and. rf2%ndir==2 .and. kdir1==2) cycle
708 :
709 : ! define factor for antisymmetric dkdk case
710 19236 : if (rf2%ndir==2 .and. kdir1==2 .and. antisymmetric) factor = -one
711 :
712 : ! define factor_in for accumulate_bands
713 19236 : factor_in = factor
714 19236 : if ((gs_hamkq%usepaw==1 .or. ipert/=natom+10) .and. rf2%ndir==1) factor_in = two
715 :
716 : ! Load projected WF according to ipert2 and idir2
717 19236 : cprj_j => cprj_empty ; ; cprj_dudk => cprj_empty
718 19236 : if (has_cprj_jband) then
719 432 : call pawcprj_free(cprj_jband)
720 432 : ncpgr_loc= 3;if(ipert2==natom+1.or.ipert2==natom+2) ncpgr_loc=1
721 432 : icpgr_loc=-1;if(ipert2==natom+1.or.ipert2==natom+2) icpgr_loc=idir2
722 432 : call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj)
723 : call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,&
724 : & isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,&
725 432 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc)
726 : end if
727 :
728 123612 : do jband=1,nband_k
729 :
730 : ! Skip bands not treated by current proc
731 89760 : if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
732 :
733 114876 : if (abs(occ_k(jband))>tol8) then
734 78240 : shift_band1=(jband-1)*size_wf
735 78240 : shift_cprj_band1=(jband-1)*size_cprj
736 78240 : shift_jband_lambda=(jband-1)*2*nband_k
737 :
738 : ! Extract first order wavefunction | du/dpert1 > and eigenvalues
739 78240 : eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir2_lambda:2*nband_k+shift_jband_lambda+shift_dir2_lambda)
740 78240 : cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)
741 78240 : if (has_dudkprj) cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1:size_cprj+shift_cprj_band1+shift_cprj_dir1)
742 :
743 78240 : if (ipert2==natom+2) then
744 : ! Extract dkdk and multiply by i :
745 915648 : gvnlx1(1,:) = -dudkdk(2,1+shift_band1:size_wf+shift_band1)
746 915648 : gvnlx1(2,:) = dudkdk(1,1+shift_band1:size_wf+shift_band1)
747 : end if
748 :
749 : ! Compute dH/dpert2 | du/dpert1 > (in h_cwave) and dS/dpert2 | du/dpert1 > (in s_cwave)
750 : call rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave_dudk,cprj_dudk,h_cwave,s_cwave,&
751 : & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnlx1,idir2,ipert2,ikpt,isppol,&
752 78240 : & mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf)
753 :
754 78240 : if (debug_mode/=0) then
755 1728 : write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 4 kdir1 = ',kdir1,' jband = ',jband
756 1728 : call wrtout(std_out,msg)
757 : end if
758 :
759 : ! For every occupied iband, we compute :
760 : ! < u^(0)(iband) | dH/dpert2 | du/dpert1(jband) > and add it to lambda_mn
761 : ! < u^(0)(iband) | dS/dpert2 | du/dpert1(jband) > and add it to amn
762 466080 : do iband=1,rf2%nband_k ! = band m
763 466080 : if (abs(occ_k(iband))>tol8) then
764 341760 : shift_band2=(iband-1)*size_wf
765 341760 : cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
766 : call rf2_accumulate_bands(rf2,4,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
767 341760 : jband,debug_mode,cwave_i,h_cwave,s_cwave,factor_in)
768 : end if
769 : end do
770 :
771 : ! Add dH/dpert2 | du/dpert1 > to RHS_Stern :
772 2032608 : if (gs_hamkq%usepaw==1) h_cwave(:,:)=h_cwave(:,:)-eig0_k(jband)*s_cwave(:,:) ! if PAW : we add H^(1)-eps^(0) S^(1)
773 78240 : rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
774 234720 : call cg_zaxpy(size_wf,(/factor,zero/),h_cwave,rhs_j)
775 :
776 : ! Compute : -factor * sum_iband ( dLambda/dpert1_{iband,jband} * dsusdu_{iband} )
777 466080 : do iband=1,nband_k
778 466080 : if (abs(occ_k(iband))>tol8) then ! if empty band, nothing to do
779 :
780 : ! Extract lambda_ij(iband,jband) for dir1
781 341760 : lambda_ij(1)=eig1_k_stored(2*iband-1+shift_jband_lambda+shift_dir1_lambda)
782 341760 : lambda_ij(2)=eig1_k_stored(2*iband +shift_jband_lambda+shift_dir1_lambda)
783 :
784 : ! Extract dsusdu for iband and pert2 (in cwave_i)
785 341760 : shift_band2=(iband-1)*size_wf
786 341760 : cwave_i => dsusdu(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2)
787 :
788 : ! Compute Lambda_{iband,jband} * dsusdu_{iband} and add it to RHS_Stern
789 341760 : rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
790 : call cg_zaxpy(size_wf,(/-factor*lambda_ij(1), &
791 1025280 : & -factor*lambda_ij(2)/),cwave_i,rhs_j) !do not forget the minus sign!
792 :
793 : end if ! empty iband test
794 : end do ! iband
795 :
796 : end if ! empty jband test
797 : end do ! jband
798 :
799 : end do ! kdir1
800 :
801 14616 : ABI_FREE(gvnlx1)
802 14616 : ABI_FREE(h_cwave)
803 14616 : ABI_FREE(s_cwave)
804 14616 : ABI_FREE(cg_jband)
805 14616 : ABI_FREE(dudk)
806 14616 : ABI_FREE(dudkdk)
807 14616 : ABI_FREE(dudk_dir2)
808 14616 : ABI_FREE(dsusdu)
809 14616 : ABI_FREE(eig1_k_stored)
810 14616 : if (has_cprj_jband) call pawcprj_free(cprj_jband)
811 20376 : ABI_FREE(cprj_jband)
812 14616 : if (has_dudkprj) call pawcprj_free(dudkprj)
813 46872 : ABI_FREE(dudkprj)
814 :
815 : ! Compute the part of 2nd order wavefunction that belongs to the space of empty bands
816 80776 : do jband=1,nband_k
817 :
818 : ! Skip bands not treated by current proc
819 66160 : if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
820 :
821 66160 : shift_band1=(jband-1)*size_wf
822 80776 : if (abs(occ_k(jband))>tol8) then
823 59760 : invocc = one/occ_k(jband)
824 59760 : rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
825 353200 : do iband=1,nband_k
826 353200 : if (iband /= jband) then
827 233680 : if (debug_mode/=0) then
828 10560 : if (abs(occ_k(iband) - occ_k(jband)) > tol12 .and. occ_k(iband) > tol8) then
829 0 : write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband
830 0 : call wrtout(std_out,msg)
831 0 : call wrtout(ab_out,msg)
832 0 : write(msg,'(a)') 'ERROR : occ_k(iband) /= occ_k(jband) (and both are >0)'
833 0 : call wrtout(std_out,msg)
834 0 : call wrtout(ab_out,msg)
835 : end if
836 10560 : if (abs(eig0_k(iband) - eig0_k(jband)) < tol8 ) then
837 120 : write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband
838 120 : call wrtout(std_out,msg)
839 120 : write(msg,'(a,es22.13e3)') 'WARNING : DEGENERATE BANDS Eig(jband) = Eig(jband) = ',eig0_k(jband)
840 120 : call wrtout(std_out,msg)
841 : end if
842 10560 : if ( (eig0_k(iband) - eig0_k(jband) < -tol12) .and. (jband < iband) ) then
843 0 : write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband
844 0 : call wrtout(std_out,msg)
845 0 : call wrtout(ab_out,msg)
846 0 : write(msg,'(a)') 'ERROR : Eig(jband) < Eig(iband) with jband < iband'
847 0 : call wrtout(std_out,msg)
848 0 : call wrtout(ab_out,msg)
849 0 : write(msg,'(a,es22.13e3)') 'Eig(jband) = ',eig0_k(jband)
850 0 : call wrtout(std_out,msg)
851 0 : call wrtout(ab_out,msg)
852 0 : write(msg,'(a,es22.13e3)') 'Eig(iband) = ',eig0_k(iband)
853 0 : call wrtout(std_out,msg)
854 0 : call wrtout(ab_out,msg)
855 : end if
856 : end if ! end tests
857 233680 : if ( abs(occ_k(iband))<tol8 ) then ! for empty bands only
858 25600 : shift_band2=(iband-1)*size_wf
859 25600 : cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
860 25600 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_i,rhs_j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
861 : ! Store it in a_mn
862 : ! /!\ There is a factor "-2" to simplify the use of amn in the following.
863 : ! /!\ Occupied and empty bands will be treated in a same way.
864 : rf2%amn(:,iband+(jband-1)*nband_k)=-two*rocceig(iband,jband)*invocc*(/dotr,doti/)&
865 76800 : +rf2%amn(:,iband+(jband-1)*nband_k)
866 : end if ! empty band test
867 : end if ! iband \= jband
868 : end do ! iband
869 : end if ! empty band test
870 : end do ! jband
871 :
872 : ! **************************************************************************************************
873 : ! COMPUTATION OF "dcwavef" AND "Lambda_mn" FROM "A_mn"
874 : ! **************************************************************************************************
875 :
876 43848 : ABI_MALLOC_OR_DIE(rf2%dcwavef,(2,nband_k*size_wf), ierr)
877 24118668 : rf2%dcwavef=zero
878 :
879 80776 : do jband=1,nband_k
880 :
881 : ! Skip bands not treated by current proc
882 66160 : if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
883 :
884 66160 : shift_band1=(jband-1)*size_wf
885 80776 : if (abs(occ_k(jband))>tol8) then
886 353200 : do iband=1,nband_k
887 293440 : shift_band2=(iband-1)*size_wf
888 :
889 : ! Extract GS wavefunction for iband
890 293440 : cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg)
891 :
892 : call cg_zaxpy(size_wf,-half*rf2%amn(:,iband+(jband-1)*nband_k), &
893 880320 : & cwave_i,rf2%dcwavef(:,1+shift_band1))
894 :
895 353200 : if (abs(occ_k(iband))>tol8 .and. abs(occ_k(jband))>tol8) then
896 : rf2%lambda_mn(:,iband+(jband-1)*nband_k) = rf2%lambda_mn(:,iband+(jband-1)*nband_k) &
897 803520 : -half*(eig0_k(iband)+eig0_k(jband))*rf2%amn(:,iband+(jband-1)*nband_k)
898 :
899 267840 : eig1_k(2*iband-1+(jband-1)*2*nband_k) = rf2%lambda_mn(1,iband+(jband-1)*nband_k)
900 267840 : eig1_k(2*iband +(jband-1)*2*nband_k) = rf2%lambda_mn(2,iband+(jband-1)*nband_k)
901 :
902 : end if ! empty band test
903 : end do ! iband
904 : end if ! empty band test
905 : end do ! jband
906 :
907 : ! For the following, "rf2%lambda_mn" and "rf2%RHS_Stern" must be computed for every bands
908 14616 : call xmpi_barrier(mpi_enreg%comm_band)
909 :
910 : ! **************************************************************************************************
911 : ! FINAL TEST
912 : ! **************************************************************************************************
913 :
914 14616 : tol_final = tol6
915 14616 : if (debug_mode/=0) then
916 3120 : do jband=1,nband_k
917 :
918 : ! Skip bands not treated by current proc
919 2880 : if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle
920 :
921 3120 : if (abs(occ_k(jband))>tol8) then
922 : ! write(msg,'(3(a,i2))') 'RF2 TEST FINAL : ipert=',ipert-natom,' idir=',idir,' jband=',jband
923 : ! call wrtout(std_out,msg)
924 960 : shift_band1=(jband-1)*size_wf
925 960 : rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1)
926 960 : cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg)
927 960 : call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,rhs_j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
928 960 : dot2r = dotr - rf2%lambda_mn(1,jband+(jband-1)*nband_k)
929 960 : dot2i = doti - rf2%lambda_mn(2,jband+(jband-1)*nband_k)
930 960 : dot2r = sqrt(dot2r**2+dot2i**2)
931 960 : if (dot2r > tol_final) then
932 0 : write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST FINAL iband = ',jband,' : NOT PASSED dotr = ',dotr
933 0 : call wrtout(std_out,msg)
934 0 : call wrtout(ab_out,msg)
935 0 : write(msg,'(2(a,es22.13E3))') ' < cwave_j | rhs_j > =',dotr,',',doti
936 0 : call wrtout(std_out,msg)
937 0 : call wrtout(ab_out,msg)
938 0 : write(msg,'(2(a,es22.13E3))') ' lambda_jj =',&
939 0 : & rf2%lambda_mn(1,jband+(jband-1)*nband_k),',',rf2%lambda_mn(2,jband+(jband-1)*nband_k)
940 0 : call wrtout(std_out,msg)
941 0 : call wrtout(ab_out,msg)
942 : else
943 960 : if (dot2r<tol9) dot2r = zero ! in order to hide the numerical noise
944 : write(msg,'(a,i2,a,es22.13E3,a,es7.1E2)') &
945 960 : 'RF2 TEST FINAL iband = ',jband,' : OK. |test| = ',dot2r,' < ',tol_final
946 960 : call wrtout(std_out,msg)
947 : end if
948 : end if
949 : end do
950 : end if
951 :
952 : ! **************************************************************************************************
953 : ! JOB FINISHED
954 : ! **************************************************************************************************
955 :
956 : ! Deallocations of arrays
957 14616 : if (debug_mode==0) then
958 14376 : ABI_FREE(rf2%amn)
959 : end if
960 :
961 14616 : call timab(514,2,tsec)
962 :
963 : DBG_EXIT("COLL")
964 :
965 29232 : end subroutine rf2_init
966 : !!***
967 :
968 : end module m_rf2_init
969 : !!***
|
