Branch data 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), 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 : : !!***
|
