TABLE OF CONTENTS
ABINIT/m_rf2_init [ Modules ]
NAME
m_rf2_init
FUNCTION
COPYRIGHT
Copyright (C) 2015-2024 ABINIT group (LB,MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
TODO
Can be merged with m_rf2
SOURCE
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
ABINIT/rf2_init [ Functions ]
NAME
rf2_init
FUNCTION
Compute terms needed for the 2nd order Sternheimer equation. All terms are stored in a rf2_t object.
INPUTS
cg(2,mpw*nspinor*mband*nsppol)=planewave coefficients of wavefunctions at k
cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
projected with non-local projectors: cprj=<p_i|Cnk>
rf2 : the object we want to initialize (see m_rf2.F90 for more information)
dtfil <type(datafiles_type)>=variables related to files
dtset <type(dataset_type)>=all input variables for this dataset
eig0_k(mband*nsppol)=GS eigenvalues at k (hartree)
eig1_k(2*mband*mband*nsppol)=2nd-order eigenvalues at k,q (hartree)
ffnl1=nonlocal form factors
ffnl1_test=nonlocal form factors used for tests (i.e. when dtset%nonlinear_info>2)
gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
ibg=shift to be applied on the location of data in the array cprj
icg=shift to be applied on the location of data in the array cg
idir=direction of the perturbation
ikpt=number of the k-point
ipert=type of the perturbation
isppol=index of current spin component
mkmem =number of k points trated by this node (GS data).
mpi_enreg=information about MPI parallelization
mpw=maximum dimensioned size of npw or wfs at k
nband_k=number of bands at this k point for that spin polarization
ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
rf_hamk_dir2 <type(rf_hamiltonian_type)>= (used only when ipert=natom+11, so q=0)
same as rf_hamkq, but the direction of the perturbation is different
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
ddk_f<wfk_t>=struct info for DDK file.
OUTPUT
rf2%RHS_Stern
NOTES
SOURCE
102 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 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 real(dp),allocatable :: cg_jband(:,:,:),ddk_read(:,:),dudkdk(:,:),dudk_dir2(:,:) 143 real(dp),allocatable :: eig1_read(:),gvnlx1(:,:),h_cwave(:,:),s_cwave(:,:),dsusdu_loc(:,:),dsusdu_gather(:,:) 144 real(dp),allocatable,target :: dsusdu(:,:),dudk(:,:),eig1_k_stored(:) 145 real(dp), ABI_CONTIGUOUS pointer :: cwave_dudk(:,:),cwave_i(:,:),cwave_j(:,:),eig1_k_jband(:) 146 real(dp),pointer :: rhs_j(:,:) 147 type(pawcprj_type),target :: cprj_empty(0,0) 148 type(pawcprj_type),allocatable,target :: cprj_jband(:,:),dudkprj(:,:) 149 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 call timab(514,1,tsec) 158 159 !my mpi rank : 160 me=mpi_enreg%me_kpt 161 162 size_wf=gs_hamkq%npw_k*gs_hamkq%nspinor 163 size_cprj=gs_hamkq%nspinor 164 natom = gs_hamkq%natom 165 debug_mode = 0 166 if (dtset%nonlinear_info>2) debug_mode = 1 ! also active a lot of tests 167 168 !Define some attributes of the rf2 object 169 rf2%nband_k = nband_k 170 rf2%size_wf = size_wf 171 rf2%size_cprj = size_cprj 172 173 if(ipert<natom+10.or.ipert>natom+11) then 174 write(msg,'(a)') 'ipert must be equal to natom+10 or natom+11 for rf2 calculations.' 175 ABI_BUG(msg) 176 end if 177 178 !Define perturbations and idirs 179 rf2%iperts(1) = natom+1 180 rf2%iperts(2) = natom+1 181 if (ipert==natom+11) rf2%iperts(2) = natom+2 182 183 if (ipert==natom+10.and.idir<=3) then ! One perturbation, one direction 184 rf2%ndir=1 185 rf2%idirs(1)=idir ; rf2%idirs(2)=idir 186 else ! Two perturbations or/and two directions 187 rf2%ndir=2 188 call rf2_getidirs(idir,idir1,idir2) 189 rf2%idirs(1)=idir1 190 rf2%idirs(2)=idir2 191 end if 192 193 !Choose dkdk option 194 symmetric = .false. 195 antisymmetric = .false. 196 total = .false. 197 198 ! even if rf2_dkdk is 5 (nonesense), 199 ! as long as rf2_dkdk /= 0, symmetric = .true. 200 if (dtset%rf2_dkdk==2) then 201 antisymmetric = .true. 202 else if (dtset%rf2_dkdk==3) then 203 total = .true. 204 else 205 symmetric = .true. 206 endif 207 208 209 ! ************************************************************************************************** 210 ! Get info from ddk files 211 ! ************************************************************************************************** 212 213 !Allocate work spaces 214 ABI_MALLOC(eig1_read,(2*nband_k)) 215 ABI_MALLOC(ddk_read,(2,size_wf)) 216 eig1_read(:)=zero 217 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 ABI_MALLOC(eig1_k_stored,(2*rf2%ndir*nband_k**2)) 222 eig1_k_stored=zero 223 224 ! "dudk" contains du/dpert1 for every bands and ndir (=1 or 2) directions 225 ABI_MALLOC_OR_DIE(dudk,(2,rf2%ndir*nband_k*size_wf), ierr) 226 dudk=zero 227 has_dudkprj=.false. 228 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 229 ABI_MALLOC(dudkprj,(natom,rf2%ndir*nband_k*size_cprj)) 230 ncpgr_loc=1;if(ipert==natom+10.or.ipert==natom+11) ncpgr_loc=3 231 call pawcprj_alloc(dudkprj,ncpgr_loc,gs_hamkq%dimcprj) 232 choice_cprj=5 ; cpopt_cprj=0 233 has_dudkprj=.true. 234 else 235 ABI_MALLOC(dudkprj,(natom,0)) 236 end if 237 238 if (debug_mode/=0) then 239 write(msg,'(4(a,i2))') 'RF2_INIT : ipert-natom = ',ipert-natom,' , idir = ',idir,& 240 ' , ikpt = ',ikpt,' , isppol = ',isppol 241 call wrtout(std_out,msg,'COLL') 242 end if 243 244 file_index(1)=1 ! dir1 245 file_index(2)=2 ! dir2 246 if (ipert==natom+11) then ! see dfpt_looppert.F90 247 file_index(1)=3 ! dir1 248 file_index(2)=2 ! dir2 249 end if 250 251 do kdir1=1,rf2%ndir 252 idir1=rf2%idirs(kdir1) 253 ipert1=rf2%iperts(kdir1) 254 do iband=1,nband_k 255 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 shift_band1=(iband-1)*size_wf 258 shift_dir1=(kdir1-1)*nband_k*size_wf 259 dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=ddk_read(:,:) 260 ! Copy eig1_read in "eig1_k_stored" 261 shift_band1=(iband-1)*2*nband_k 262 shift_dir1_lambda=2*(kdir1-1)*nband_k**2 263 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 if (has_dudkprj.and.mpi_enreg%proc_distrb(ikpt,iband,isppol)==me) then 266 shift_cprj_band1=(iband-1)*size_cprj 267 shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj 268 cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1: & 269 & size_cprj+shift_cprj_band1+shift_cprj_dir1) 270 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,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 & gs_hamkq%ph1d,gs_hamkq%ph3d_k,gs_hamkq%ucvol,gs_hamkq%useylm) 276 end if 277 end do 278 end do 279 280 ABI_MALLOC(dudkdk,(2,0)) 281 ABI_MALLOC(dudk_dir2,(2,0)) 282 283 !Get dudkdk for ipert==natom+11 284 if(ipert==natom+11) then 285 ABI_FREE(dudkdk) 286 ABI_MALLOC(dudkdk,(2,nband_k*size_wf)) 287 if (idir>3) then 288 ABI_FREE(dudk_dir2) 289 ABI_MALLOC(dudk_dir2,(2,nband_k*size_wf)) 290 end if 291 do iband=1,nband_k 292 call ddk_f(1)%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read) 293 shift_band1=(iband-1)*size_wf 294 dudkdk(:,1+shift_band1:size_wf+shift_band1)=ddk_read(:,:) 295 ! Check that < u^(0) | u^(2) > = - Re[< u^(1) | u^(1) >] 296 if (debug_mode/=0 .and. gs_hamkq%usepaw==0) then 297 ! Compute < u^(0) | u^(2) > 298 do jband=1,nband_k 299 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 300 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 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 cwave_dudk => dudk(:,1+shift_band1:size_wf+shift_band1) 304 call sqnorm_g(dot2r,gs_hamkq%istwf_k,size_wf,cwave_dudk,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 305 dotr = dotr + dot2r 306 dotr = sqrt(dotr**2+doti**2) 307 if (dotr > tol7) then 308 write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST dudkdk iband = ',iband,& 309 ' : NOT PASSED. | < u^(0) | u^(2) > + Re[< u^(1) | u^(1) >] | = ',dotr 310 call wrtout(std_out,msg) 311 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 if (idir>3) then 321 call ddk_f(4)%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read) 322 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 ABI_FREE(ddk_read) 328 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 ABI_MALLOC(h_cwave,(2,size_wf)) 336 ABI_MALLOC(s_cwave,(2,size_wf)) 337 ABI_MALLOC(gvnlx1,(2,size_wf)) 338 h_cwave(:,:) = zero 339 s_cwave(:,:) = zero 340 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 ABI_MALLOC_OR_DIE(dsusdu,(2,rf2%ndir*nband_k*size_wf), ierr) 344 dsusdu=zero 345 346 ABI_MALLOC(rf2%amn,(2,nband_k**2)) 347 rf2%amn=zero 348 349 ABI_MALLOC(rf2%lambda_mn,(2,nband_k**2)) 350 rf2%lambda_mn(:,:)=zero 351 352 !Allocate work spaces when debug_mode is activated 353 has_cprj_jband=.false. 354 if (debug_mode/=0) then ! Only for test purposes 355 ABI_MALLOC(cg_jband,(2,size_wf*nband_k,2)) 356 cg_jband(:,:,1) = cg(:,1+icg:size_wf*nband_k+icg) 357 if (ipert==natom+11) then ! Note the multiplication by "i" 358 if (idir<=3) then 359 cg_jband(1,:,2) = -dudk(2,1:size_wf*nband_k) ! for dir1 360 cg_jband(2,:,2) = dudk(1,1:size_wf*nband_k) ! for dir1 361 else 362 cg_jband(1,:,2) = -dudk_dir2(2,1:size_wf*nband_k) ! for dir2 363 cg_jband(2,:,2) = dudk_dir2(1,1:size_wf*nband_k) ! for dir2 364 end if 365 end if 366 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 367 ABI_MALLOC(cprj_jband,(natom,size_cprj*nband_k)) 368 has_cprj_jband=.true. 369 else 370 ABI_MALLOC(cprj_jband,(natom,0)) 371 end if 372 else 373 ABI_MALLOC(cg_jband,(2,0,2)) 374 ABI_MALLOC(cprj_jband,(natom,0)) 375 end if 376 377 ! define dkdk factor 378 factor = one 379 if(ipert==natom+10 .and. idir<=3 .and. symmetric) factor=two ! in order to not compute same terms twice 380 381 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 idir1=rf2%idirs(kdir1) 389 ipert1=rf2%iperts(kdir1) 390 shift_dir1=(kdir1-1)*nband_k*size_wf 391 shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj 392 shift_dir1_lambda=(kdir1-1)*2*nband_k**2 393 if(ipert==natom+10 .and. idir<=3) then 394 shift_dir2=0 395 idir2 = idir1 396 ipert2 = ipert1 397 rf_hamk_idir => rf_hamkq 398 else 399 shift_dir2=(2-kdir1)*nband_k*size_wf 400 idir2 = rf2%idirs(3-kdir1) 401 ipert2 = rf2%iperts(3-kdir1) 402 if (kdir1==1) rf_hamk_idir => rf_hamkq 403 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 if (antisymmetric .and. rf2%ndir==1) factor = zero 408 if (total .and. rf2%ndir==2 .and. kdir1==2) factor = zero 409 410 ! define factor for antisymmetric dkdk case 411 if (rf2%ndir==2 .and. kdir1==2 .and. antisymmetric) factor = -one 412 413 ! define factor_in for accumulate_bands 414 factor_in = factor 415 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 cprj_j => cprj_empty ; cprj_dudk => cprj_empty 419 if (has_cprj_jband) then 420 call pawcprj_free(cprj_jband) 421 ncpgr_loc= 3;if(ipert1==natom+1.or.ipert1==natom+2) ncpgr_loc=1 422 icpgr_loc=-1;if(ipert1==natom+1.or.ipert1==natom+2) icpgr_loc=idir1 423 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 & 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 do jband=1,nband_k ! = band n 431 432 ! Skip bands not treated by current proc 433 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 434 435 shift_band1=(jband-1)*size_wf 436 shift_cprj_band1=(jband-1)*size_cprj 437 shift_jband_lambda=(jband-1)*2*nband_k 438 439 if (abs(occ_k(jband))>tol8) then 440 441 ! Extract first order wavefunction and eigenvalues for jband 442 eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir1_lambda:2*nband_k+shift_jband_lambda+shift_dir1_lambda) 443 cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 444 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 & mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf) 450 451 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 +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 456 457 if (debug_mode/=0) then 458 write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 1 kdir1 = ',kdir1,' jband = ',jband 459 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 do iband=1,rf2%nband_k ! = band m 466 if (abs(occ_k(iband))>tol8) then 467 shift_band2=(iband-1)*size_wf 468 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 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 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 476 if(has_cprj_jband) cprj_j => cprj_jband(:,1+shift_cprj_band1:size_cprj+shift_cprj_band1) 477 478 if (ipert1==natom+2) then 479 ! Extract ddk and multiply by i : 480 if(idir<=3) then ! in this case : idir1=idir2 481 gvnlx1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1) 482 gvnlx1(2,:) = dudk(1,1+shift_band1:size_wf+shift_band1) 483 else 484 gvnlx1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1) 485 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 & mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf) 493 494 ! Copy infos in dsusdu 495 if (gs_hamkq%usepaw==1) then 496 dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=s_cwave(:,:)& 497 +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 498 end if 499 500 if (debug_mode/=0) then 501 write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 2 kdir1 = ',kdir1,' jband = ',jband 502 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 do iband=1,rf2%nband_k ! = band m 509 if (abs(occ_k(iband))>tol8) then 510 shift_band2=(iband-1)*size_wf 511 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 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 nproc_band = xmpi_comm_size(mpi_enreg%comm_band) 523 if (nproc_band>1) then 524 525 my_nband = nband_k/nproc_band;if (mod(nband_k,nproc_band)/=0) my_nband=my_nband+1 526 ABI_MALLOC(dsusdu_loc,(2,size_wf*my_nband*rf2%ndir)) 527 ABI_MALLOC(dsusdu_gather,(2,size_wf*my_nband*rf2%ndir*nproc_band)) 528 dsusdu_loc(:,:) = zero 529 dsusdu_gather(:,:) = zero 530 531 do kdir1=1,rf2%ndir 532 indb = 1 533 shift_dir1=(kdir1-1)*size_wf*nband_k 534 shift_dir1_loc=(kdir1-1)*size_wf*my_nband 535 do jband=1,nband_k 536 ! Skip bands not treated by current proc 537 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 538 539 shift_band1=(jband-1)*size_wf 540 dsusdu_loc(:,indb+shift_dir1_loc:indb-1+size_wf+shift_dir1_loc) = & 541 dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 542 indb = indb + size_wf 543 end do 544 end do 545 546 call xmpi_allgather(dsusdu_loc,2*size_wf*my_nband*rf2%ndir,dsusdu_gather,mpi_enreg%comm_band,ierr) 547 548 do kdir1=1,rf2%ndir 549 shift_dir1=(kdir1-1)*size_wf*nband_k 550 shift_dir1_loc=(kdir1-1)*size_wf*my_nband 551 do iproc=1,nproc_band 552 shift_proc = (iproc-1)*size_wf*my_nband*rf2%ndir 553 indb = 1 554 do jband=1,my_nband 555 iband = jband+(iproc-1)*my_nband 556 if(iband<=nband_k) then 557 shift_band1=(iband-1)*size_wf 558 dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) = & 559 dsusdu_gather(:,indb+shift_dir1_loc+shift_proc:indb-1+size_wf+shift_dir1_loc+shift_proc) 560 end if 561 indb = indb + size_wf 562 end do 563 end do 564 end do 565 ABI_FREE(dsusdu_loc) 566 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 ABI_MALLOC_OR_DIE(rf2%RHS_Stern,(2,nband_k*size_wf), ierr) 575 rf2%RHS_Stern(:,:)=zero 576 577 ! Define prefactor for the H^{(2)} term 578 factor = one 579 if (total) factor = half 580 581 !Computation of terms containing H^(2) 582 if (ipert/=natom+11 .or. gs_hamkq%usepaw==1) then ! Otherwise H^(2) = 0 583 if (.not. antisymmetric) then ! check antisymmetric dkdk 584 585 !Load projected WF according to ipert and idir 586 cprj_j => cprj_empty 587 if (has_cprj_jband) then 588 call pawcprj_free(cprj_jband) 589 ncpgr_loc= 3;if(ipert==natom+1.or.ipert==natom+2) ncpgr_loc=1 590 icpgr_loc=-1;if(ipert==natom+1.or.ipert==natom+2) icpgr_loc=idir 591 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 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc) 595 end if 596 597 if (ipert==natom+10) then 598 rf_hamk_idir => rf_hamkq ! all info are in rf_hamkq 599 else if (ipert==natom+11) then 600 rf_hamk_idir => rf_hamk_dir2 ! all info are in rf_hamk_dir2 601 end if 602 603 do jband=1,nband_k 604 605 ! Skip bands not treated by current proc 606 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 607 608 if (abs(occ_k(jband))>tol8) then 609 shift_band1=(jband-1)*size_wf 610 shift_cprj_band1=(jband-1)*size_cprj 611 shift_jband_lambda=(jband-1)*2*nband_k 612 613 ! Extract GS wavefunction 614 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 615 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 eig1_k_jband => eig1_k_stored(1+shift_jband_lambda:2*nband_k+shift_jband_lambda) 619 620 if (ipert==natom+11) then 621 ! Extract ddk and multiply by i : 622 if(idir<=3) then ! in this case : idir1=idir2 623 gvnlx1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1) 624 gvnlx1(2,:) = dudk(1,1+shift_band1:size_wf+shift_band1) 625 else 626 gvnlx1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1) 627 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 & ffnl1=ffnl1,ffnl1_test=ffnl1_test) 637 638 if (debug_mode/=0) then 639 write(msg,'(a,i2)') 'RF2 TEST before accumulate_bands choice = 3 jband = ',jband 640 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 do iband=1,rf2%nband_k ! = band m 647 if (abs(occ_k(iband))>tol8) then 648 shift_band2=(iband-1)*size_wf 649 cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg) 650 if(ipert == natom+10) then 651 ipert1 = natom+1 652 ipert2 = natom+1 653 else 654 ipert1 = natom+1 655 ipert2 = natom+2 656 end if 657 call rf2_getidirs(idir,idir1,idir2) 658 call rf2_accumulate_bands(rf2,3,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,& 659 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 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 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 666 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 factor = one 675 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 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 shift_dir1=(kdir1-1)*nband_k*size_wf 686 shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj 687 shift_dir1_lambda=(kdir1-1)*2*nband_k**2 688 idir1=rf2%idirs(kdir1) 689 ipert1=rf2%iperts(kdir1) 690 if(ipert==natom+10 .and. idir<=3) then 691 idir2=idir1 692 ipert2=ipert1 693 shift_dir2=0 694 shift_dir2_lambda=0 695 rf_hamk_idir => rf_hamkq 696 else 697 idir2=rf2%idirs(2-kdir1+1) 698 ipert2=rf2%iperts(2-kdir1+1) 699 shift_dir2=(2-kdir1)*nband_k*size_wf 700 shift_dir2_lambda=(2-kdir1)*2*nband_k**2 701 if (kdir1==1) rf_hamk_idir => rf_hamk_dir2 ! dir2 702 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 if (antisymmetric .and. rf2%ndir==1) exit 707 if (total .and. rf2%ndir==2 .and. kdir1==2) cycle 708 709 ! define factor for antisymmetric dkdk case 710 if (rf2%ndir==2 .and. kdir1==2 .and. antisymmetric) factor = -one 711 712 ! define factor_in for accumulate_bands 713 factor_in = factor 714 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 cprj_j => cprj_empty ; ; cprj_dudk => cprj_empty 718 if (has_cprj_jband) then 719 call pawcprj_free(cprj_jband) 720 ncpgr_loc= 3;if(ipert2==natom+1.or.ipert2==natom+2) ncpgr_loc=1 721 icpgr_loc=-1;if(ipert2==natom+1.or.ipert2==natom+2) icpgr_loc=idir2 722 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 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc) 726 end if 727 728 do jband=1,nband_k 729 730 ! Skip bands not treated by current proc 731 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 732 733 if (abs(occ_k(jband))>tol8) then 734 shift_band1=(jband-1)*size_wf 735 shift_cprj_band1=(jband-1)*size_cprj 736 shift_jband_lambda=(jband-1)*2*nband_k 737 738 ! Extract first order wavefunction | du/dpert1 > and eigenvalues 739 eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir2_lambda:2*nband_k+shift_jband_lambda+shift_dir2_lambda) 740 cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 741 if (has_dudkprj) cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1:size_cprj+shift_cprj_band1+shift_cprj_dir1) 742 743 if (ipert2==natom+2) then 744 ! Extract dkdk and multiply by i : 745 gvnlx1(1,:) = -dudkdk(2,1+shift_band1:size_wf+shift_band1) 746 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 & mkmem,mpi_enreg,nband_k,nsppol,debug_mode,dtset%prtvol,rf_hamk_idir,size_cprj,size_wf) 753 754 if (debug_mode/=0) then 755 write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 4 kdir1 = ',kdir1,' jband = ',jband 756 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 do iband=1,rf2%nband_k ! = band m 763 if (abs(occ_k(iband))>tol8) then 764 shift_band2=(iband-1)*size_wf 765 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 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 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 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 774 call cg_zaxpy(size_wf,(/factor,zero/),h_cwave,rhs_j) 775 776 ! Compute : -factor * sum_iband ( dLambda/dpert1_{iband,jband} * dsusdu_{iband} ) 777 do iband=1,nband_k 778 if (abs(occ_k(iband))>tol8) then ! if empty band, nothing to do 779 780 ! Extract lambda_ij(iband,jband) for dir1 781 lambda_ij(1)=eig1_k_stored(2*iband-1+shift_jband_lambda+shift_dir1_lambda) 782 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 shift_band2=(iband-1)*size_wf 786 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 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 790 call cg_zaxpy(size_wf,(/-factor*lambda_ij(1), & 791 & -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 ABI_FREE(gvnlx1) 802 ABI_FREE(h_cwave) 803 ABI_FREE(s_cwave) 804 ABI_FREE(cg_jband) 805 ABI_FREE(dudk) 806 ABI_FREE(dudkdk) 807 ABI_FREE(dudk_dir2) 808 ABI_FREE(dsusdu) 809 ABI_FREE(eig1_k_stored) 810 if (has_cprj_jband) call pawcprj_free(cprj_jband) 811 ABI_FREE(cprj_jband) 812 if (has_dudkprj) call pawcprj_free(dudkprj) 813 ABI_FREE(dudkprj) 814 815 ! Compute the part of 2nd order wavefunction that belongs to the space of empty bands 816 do jband=1,nband_k 817 818 ! Skip bands not treated by current proc 819 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 820 821 shift_band1=(jband-1)*size_wf 822 if (abs(occ_k(jband))>tol8) then 823 invocc = one/occ_k(jband) 824 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 825 do iband=1,nband_k 826 if (iband /= jband) then 827 if (debug_mode/=0) then 828 if (abs(occ_k(iband) - occ_k(jband)) > tol12 .and. occ_k(iband) > tol8) then 829 write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband 830 call wrtout(std_out,msg) 831 call wrtout(ab_out,msg) 832 write(msg,'(a)') 'ERROR : occ_k(iband) /= occ_k(jband) (and both are >0)' 833 call wrtout(std_out,msg) 834 call wrtout(ab_out,msg) 835 end if 836 if (abs(eig0_k(iband) - eig0_k(jband)) < tol8 ) then 837 write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband 838 call wrtout(std_out,msg) 839 write(msg,'(a,es22.13e3)') 'WARNING : DEGENERATE BANDS Eig(jband) = Eig(jband) = ',eig0_k(jband) 840 call wrtout(std_out,msg) 841 end if 842 if ( (eig0_k(iband) - eig0_k(jband) < -tol12) .and. (jband < iband) ) then 843 write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband 844 call wrtout(std_out,msg) 845 call wrtout(ab_out,msg) 846 write(msg,'(a)') 'ERROR : Eig(jband) < Eig(iband) with jband < iband' 847 call wrtout(std_out,msg) 848 call wrtout(ab_out,msg) 849 write(msg,'(a,es22.13e3)') 'Eig(jband) = ',eig0_k(jband) 850 call wrtout(std_out,msg) 851 call wrtout(ab_out,msg) 852 write(msg,'(a,es22.13e3)') 'Eig(iband) = ',eig0_k(iband) 853 call wrtout(std_out,msg) 854 call wrtout(ab_out,msg) 855 end if 856 end if ! end tests 857 if ( abs(occ_k(iband))<tol8 ) then ! for empty bands only 858 shift_band2=(iband-1)*size_wf 859 cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg) 860 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 +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 ABI_MALLOC_OR_DIE(rf2%dcwavef,(2,nband_k*size_wf), ierr) 877 rf2%dcwavef=zero 878 879 do jband=1,nband_k 880 881 ! Skip bands not treated by current proc 882 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 883 884 shift_band1=(jband-1)*size_wf 885 if (abs(occ_k(jband))>tol8) then 886 do iband=1,nband_k 887 shift_band2=(iband-1)*size_wf 888 889 ! Extract GS wavefunction for iband 890 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 & cwave_i,rf2%dcwavef(:,1+shift_band1)) 894 895 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 -half*(eig0_k(iband)+eig0_k(jband))*rf2%amn(:,iband+(jband-1)*nband_k) 898 899 eig1_k(2*iband-1+(jband-1)*2*nband_k) = rf2%lambda_mn(1,iband+(jband-1)*nband_k) 900 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 call xmpi_barrier(mpi_enreg%comm_band) 909 910 ! ************************************************************************************************** 911 ! FINAL TEST 912 ! ************************************************************************************************** 913 914 tol_final = tol6 915 if (debug_mode/=0) then 916 do jband=1,nband_k 917 918 ! Skip bands not treated by current proc 919 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 920 921 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 shift_band1=(jband-1)*size_wf 925 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 926 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 927 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 dot2r = dotr - rf2%lambda_mn(1,jband+(jband-1)*nband_k) 929 dot2i = doti - rf2%lambda_mn(2,jband+(jband-1)*nband_k) 930 dot2r = sqrt(dot2r**2+dot2i**2) 931 if (dot2r > tol_final) then 932 write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST FINAL iband = ',jband,' : NOT PASSED dotr = ',dotr 933 call wrtout(std_out,msg) 934 call wrtout(ab_out,msg) 935 write(msg,'(2(a,es22.13E3))') ' < cwave_j | rhs_j > =',dotr,',',doti 936 call wrtout(std_out,msg) 937 call wrtout(ab_out,msg) 938 write(msg,'(2(a,es22.13E3))') ' lambda_jj =',& 939 & rf2%lambda_mn(1,jband+(jband-1)*nband_k),',',rf2%lambda_mn(2,jband+(jband-1)*nband_k) 940 call wrtout(std_out,msg) 941 call wrtout(ab_out,msg) 942 else 943 if (dot2r<tol9) dot2r = zero ! in order to hide the numerical noise 944 write(msg,'(a,i2,a,es22.13E3,a,es7.1E2)') & 945 'RF2 TEST FINAL iband = ',jband,' : OK. |test| = ',dot2r,' < ',tol_final 946 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 if (debug_mode==0) then 958 ABI_FREE(rf2%amn) 959 end if 960 961 call timab(514,2,tsec) 962 963 DBG_EXIT("COLL") 964 965 end subroutine rf2_init