TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 2015-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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)
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<wfk_t>=struct info for DDK file.
OUTPUT
rf2%RHS_Stern
NOTES
PARENTS
dfpt_vtowfk
CHILDREN
cg_zaxpy,dotprod_g,getcprj,pawcprj_alloc,pawcprj_free,pawcprj_get rf2_accumulate_bands,rf2_apply_hamiltonian,rf2_getidirs,sqnorm_g wfk_read_bks,wrtout,xmpi_allgather,xmpi_barrier
SOURCE
63 #if defined HAVE_CONFIG_H 64 #include "config.h" 65 #endif 66 67 #include "abi_common.h" 68 69 subroutine rf2_init(cg,cprj,rf2,dtset,dtfil,eig0_k,eig1_k,gs_hamkq,ibg,icg,idir,ikpt,ipert,isppol,mkmem,& 70 mpi_enreg,mpw,nband_k,nsppol,rf_hamkq,rf_hamk_dir2,occ_k,rocceig,ddk_f) 71 72 use defs_basis 73 use defs_datatypes 74 use defs_abitypes 75 use m_xmpi 76 use m_errors 77 use m_wfk 78 use m_hamiltonian 79 use m_cgtools 80 use m_rf2 81 82 use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_copy,pawcprj_get,pawcprj_free 83 84 !This section has been created automatically by the script Abilint (TD). 85 !Do not modify the following lines by hand. 86 #undef ABI_FUNC 87 #define ABI_FUNC 'rf2_init' 88 use interfaces_14_hidewrite 89 use interfaces_66_nonlocal 90 !End of the abilint section 91 92 implicit none 93 94 ! ************************************************************************* 95 !Arguments ------------------------------- 96 !scalars 97 98 integer,intent(in) :: ibg,icg,idir,ipert,isppol,ikpt 99 integer,intent(in) :: mkmem,mpw,nband_k,nsppol 100 type(datafiles_type),intent(in) :: dtfil 101 type(dataset_type),intent(in) :: dtset 102 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 103 type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq,rf_hamk_dir2 104 type(MPI_type),intent(in) :: mpi_enreg 105 106 !arrays 107 real(dp),intent(in),target :: cg(2,mpw*gs_hamkq%nspinor*dtset%mband*mkmem*nsppol) 108 real(dp),intent(in) :: eig0_k(dtset%mband) 109 real(dp),intent(inout) :: eig1_k(2*dtset%mband**2) ! Here eig1_k contains 2nd order eigenvalues... 110 real(dp),intent(in) :: occ_k(nband_k),rocceig(nband_k,nband_k) 111 type(pawcprj_type),intent(in) :: cprj(gs_hamkq%natom,gs_hamkq%nspinor*dtset%mband*mkmem*nsppol*gs_hamkq%usecprj) 112 type(rf2_t),intent(inout) :: rf2 113 type(wfk_t),intent(inout) :: ddk_f(4) 114 ! 115 !Local variables------------------------------- 116 !scalars 117 integer,parameter :: berryopt=0,iorder_cprj=0,level=19,tim_getghc=1,tim_getgh1c=1,tim_getgh2c=1 118 integer :: choice_cprj,cpopt_cprj,iband,icpgr_loc,idir1,idir2,idir_cprj,ierr 119 integer :: indb,ipert1,ipert2,iproc,jband,kdir1 120 integer :: me,my_nband,natom,ncpgr_loc,nproc_band,print_info 121 integer :: size_cprj,size_wf,shift_band1,shift_band2,shift_cprj_band1,shift_cprj_dir1,shift_proc 122 integer :: shift_dir1_lambda,shift_dir2_lambda,shift_dir1,shift_dir1_loc,shift_dir2,shift_jband_lambda 123 logical :: has_cprj_jband,has_dudkprj 124 real(dp) :: doti,dotr,dot2r,invocc,tol_final,factor 125 character(len=500) :: msg 126 !arrays 127 integer :: file_index(2) 128 real(dp) :: lambda_ij(2) 129 real(dp),allocatable :: cg_jband(:,:,:),ddk_read(:,:),dudkdk(:,:),dudk_dir2(:,:) 130 real(dp),allocatable :: eig1_read(:),gvnl1(:,:),h_cwave(:,:),s_cwave(:,:),dsusdu_loc(:,:),dsusdu_gather(:,:) 131 real(dp),allocatable,target :: dsusdu(:,:),dudk(:,:),eig1_k_stored(:) 132 real(dp), ABI_CONTIGUOUS pointer :: cwave_dudk(:,:),cwave_i(:,:),cwave_j(:,:),eig1_k_jband(:) 133 real(dp),pointer :: rhs_j(:,:) 134 type(pawcprj_type),target :: cprj_empty(0,0) 135 type(pawcprj_type),allocatable,target :: cprj_jband(:,:),dudkprj(:,:) 136 type(pawcprj_type),pointer :: cprj_dudk(:,:),cprj_j(:,:) 137 type(rf_hamiltonian_type),pointer :: rf_hamk_idir 138 139 ! ********************************************************************* 140 141 DBG_ENTER("COLL") 142 143 !my mpi rank : 144 me=mpi_enreg%me_kpt 145 146 size_wf=gs_hamkq%npw_k*gs_hamkq%nspinor 147 size_cprj=gs_hamkq%nspinor 148 natom = gs_hamkq%natom 149 print_info = 0 150 if (dtset%prtvol == -level.or.dtset%prtvol == -20.or.dtset%prtvol == -21) print_info = 1 ! also active a lot of tests 151 152 !Define some attributes of the rf2 object 153 rf2%nband_k = nband_k 154 rf2%size_wf = size_wf 155 rf2%size_cprj = size_cprj 156 157 if(ipert<natom+10.or.ipert>natom+11) then 158 write(msg,'(a)') 'ipert must be equal to natom+10 or natom+11 for rf2 calculations.' 159 MSG_BUG(msg) 160 end if 161 162 !Define perturbations and idirs 163 rf2%iperts(1) = natom+1 164 rf2%iperts(2) = natom+1 165 if (ipert==natom+11) rf2%iperts(2) = natom+2 166 167 if (ipert==natom+10.and.idir<=3) then ! One perturbation, one direction 168 rf2%ndir=1 169 rf2%idirs(1)=idir ; rf2%idirs(2)=idir 170 else ! Two perturbations or/and two directions 171 rf2%ndir=2 172 call rf2_getidirs(idir,idir1,idir2) 173 rf2%idirs(1)=idir1 174 rf2%idirs(2)=idir2 175 end if 176 177 ! ************************************************************************************************** 178 ! Get info from ddk files 179 ! ************************************************************************************************** 180 181 !Allocate work spaces 182 ABI_ALLOCATE(eig1_read,(2*nband_k)) 183 ABI_ALLOCATE(ddk_read,(2,size_wf)) 184 eig1_read(:)=zero 185 ddk_read(:,:)=zero 186 187 ! "eig1_k_stored" contains dLambda_{nm}/dpert every bands n and m and ndir (=1 or 2) directions 188 ! pert = k_dir (wavevector) or E_dir (electric field) 189 ABI_ALLOCATE(eig1_k_stored,(2*rf2%ndir*nband_k**2)) 190 eig1_k_stored=zero 191 192 ! "dudk" contains du/dpert1 for every bands and ndir (=1 or 2) directions 193 ABI_STAT_ALLOCATE(dudk,(2,rf2%ndir*nband_k*size_wf), ierr) 194 ABI_CHECK(ierr==0, "out of memory in m_rf2 : dudk") 195 dudk=zero 196 has_dudkprj=.false. 197 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 198 ABI_DATATYPE_ALLOCATE(dudkprj,(natom,rf2%ndir*nband_k*size_cprj)) 199 ncpgr_loc=1;if(ipert==natom+10.or.ipert==natom+11) ncpgr_loc=3 200 call pawcprj_alloc(dudkprj,ncpgr_loc,gs_hamkq%dimcprj) 201 choice_cprj=5 ; cpopt_cprj=0 202 has_dudkprj=.true. 203 else 204 ABI_DATATYPE_ALLOCATE(dudkprj,(natom,0)) 205 end if 206 207 if (print_info/=0) then 208 write(msg,'(4(a,i2))') 'RF2_INIT : ipert-natom = ',ipert-natom,' , idir = ',idir,& 209 ' , ikpt = ',ikpt,' , isppol = ',isppol 210 call wrtout(std_out,msg,'COLL') 211 end if 212 213 file_index(1)=1 ! dir1 214 file_index(2)=2 ! dir2 215 if (ipert==natom+11) then ! see dfpt_looppert.F90 216 file_index(1)=3 ! dir1 217 file_index(2)=2 ! dir2 218 end if 219 220 do kdir1=1,rf2%ndir 221 idir1=rf2%idirs(kdir1) 222 ipert1=rf2%iperts(kdir1) 223 do iband=1,nband_k 224 call wfk_read_bks(ddk_f(file_index(kdir1)),iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read) 225 ! Copy ddk_read in "dudk" 226 shift_band1=(iband-1)*size_wf 227 shift_dir1=(kdir1-1)*nband_k*size_wf 228 dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=ddk_read(:,:) 229 ! Copy eig1_read in "eig1_k_stored" 230 shift_band1=(iband-1)*2*nband_k 231 shift_dir1_lambda=2*(kdir1-1)*nband_k**2 232 eig1_k_stored(1+shift_band1+shift_dir1_lambda:2*nband_k+shift_band1+shift_dir1_lambda)=eig1_read(:) 233 ! Get this dudk projected on NL projectors 234 if (has_dudkprj.and.mpi_enreg%proc_distrb(ikpt,iband,isppol)==me) then 235 shift_cprj_band1=(iband-1)*size_cprj 236 shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj 237 cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1: & 238 & size_cprj+shift_cprj_band1+shift_cprj_dir1) 239 idir_cprj=0;if (dudkprj(1,1)%ncpgr/=3) idir_cprj=idir1 240 call getcprj(choice_cprj,cpopt_cprj,ddk_read,cprj_dudk,gs_hamkq%ffnl_k,idir_cprj,& 241 & gs_hamkq%indlmn,gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kpg_k,gs_hamkq%kpt_k,& 242 & gs_hamkq%lmnmax,gs_hamkq%mgfft,mpi_enreg,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,& 243 & gs_hamkq%nloalg,gs_hamkq%npw_k,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkxred,& 244 & gs_hamkq%ph1d,gs_hamkq%ph3d_k,gs_hamkq%ucvol,gs_hamkq%useylm) 245 end if 246 end do 247 end do 248 249 ABI_ALLOCATE(dudkdk,(2,0)) 250 ABI_ALLOCATE(dudk_dir2,(2,0)) 251 252 !Get dudkdk for ipert==natom+11 253 if(ipert==natom+11) then 254 ABI_DEALLOCATE(dudkdk) 255 ABI_ALLOCATE(dudkdk,(2,nband_k*size_wf)) 256 if (idir>3) then 257 ABI_DEALLOCATE(dudk_dir2) 258 ABI_ALLOCATE(dudk_dir2,(2,nband_k*size_wf)) 259 end if 260 do iband=1,nband_k 261 call wfk_read_bks(ddk_f(1),iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read) 262 shift_band1=(iband-1)*size_wf 263 dudkdk(:,1+shift_band1:size_wf+shift_band1)=ddk_read(:,:) 264 ! Check that < u^(0) | u^(2) > = - Re[< u^(1) | u^(1) >] 265 if (print_info/=0 .and. gs_hamkq%usepaw==0) then 266 ! Compute < u^(0) | u^(2) > 267 do jband=1,nband_k 268 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 269 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,ddk_read,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 270 if (idir<=3 .and. iband==jband .and. abs(occ_k(iband))>tol8) then 271 ! Compute < u^(1) | u^(1) > = Re[< u^(1) | u^(1) >] 272 cwave_dudk => dudk(:,1+shift_band1:size_wf+shift_band1) 273 call sqnorm_g(dot2r,gs_hamkq%istwf_k,size_wf,cwave_dudk,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 274 dotr = dotr + dot2r 275 dotr = sqrt(dotr**2+doti**2) 276 if (dotr > tol7) then 277 write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST dudkdk iband = ',iband,& 278 ' : NOT PASSED. | < u^(0) | u^(2) > + Re[< u^(1) | u^(1) >] | = ',dotr 279 call wrtout(std_out,msg) 280 ! else 281 ! write(msg,'(a,i2,a)') 'RF2 TEST dudkdk iband = ',iband,' : OK.' 282 ! call wrtout(std_out,msg) 283 end if 284 end if ! idir<=3 285 end do ! jband 286 end if ! print_info 287 !Read ddk for idir2 288 if (idir>3) then 289 call wfk_read_bks(ddk_f(4),iband,ikpt,isppol,xmpio_single,cg_bks=ddk_read,eig1_bks=eig1_read) 290 dudk_dir2(:,1+shift_band1:size_wf+shift_band1)=ddk_read(:,:) 291 end if 292 end do !iband 293 end if ! ipert=natom+11 294 295 ABI_DEALLOCATE(ddk_read) 296 ABI_DEALLOCATE(eig1_read) 297 298 ! ************************************************************************************************** 299 ! COMPUTATION OF "dsusdu", A PART OF "A_mn" AND A PART OF "Lambda_mn" (see defs in m_rf2) 300 ! ************************************************************************************************** 301 302 !Allocate work spaces for one band 303 ABI_ALLOCATE(h_cwave,(2,size_wf)) 304 ABI_ALLOCATE(s_cwave,(2,size_wf)) 305 ABI_ALLOCATE(gvnl1,(2,size_wf)) 306 h_cwave(:,:) = zero 307 s_cwave(:,:) = zero 308 gvnl1(:,:) = zero 309 310 ! "dsusdu" contains dS/dpert_dir |u_band> + S|du_band/dpert1> for every bands and ndir (=1 or 2) directions 311 ABI_STAT_ALLOCATE(dsusdu,(2,rf2%ndir*nband_k*size_wf), ierr) 312 ABI_CHECK(ierr==0, "out of memory in rf2_init : dsusdu") 313 dsusdu=zero 314 315 ABI_ALLOCATE(rf2%amn,(2,nband_k**2)) 316 rf2%amn=zero 317 318 ABI_ALLOCATE(rf2%lambda_mn,(2,nband_k**2)) 319 rf2%lambda_mn(:,:)=zero 320 321 !Allocate work spaces when print_info is activated 322 has_cprj_jband=.false. 323 if (print_info/=0) then ! Only for test purposes 324 ABI_ALLOCATE(cg_jband,(2,size_wf*nband_k,2)) 325 cg_jband(:,:,1) = cg(:,1+icg:size_wf*nband_k+icg) 326 if (ipert==natom+11) then ! Note the multiplication by "i" 327 if (idir<=3) then 328 cg_jband(1,:,2) = -dudk(2,1:size_wf*nband_k) ! for dir1 329 cg_jband(2,:,2) = dudk(1,1:size_wf*nband_k) ! for dir1 330 else 331 cg_jband(1,:,2) = -dudk_dir2(2,1:size_wf*nband_k) ! for dir2 332 cg_jband(2,:,2) = dudk_dir2(1,1:size_wf*nband_k) ! for dir2 333 end if 334 end if 335 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 336 ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,size_cprj*nband_k)) 337 has_cprj_jband=.true. 338 else 339 ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,0)) 340 end if 341 else 342 ABI_ALLOCATE(cg_jband,(2,0,2)) 343 ABI_DATATYPE_ALLOCATE(cprj_jband,(natom,0)) 344 end if 345 346 factor=one 347 if(ipert==natom+10 .and. idir<=3) factor=two ! in order to not compute same terms twice 348 349 do kdir1=1,rf2%ndir 350 ! First iteration (kdir1=1) : 351 ! pert1 = rf2%iperts(1) along rf2%idirs(1) 352 ! pert2 = rf2%iperts(2) along rf2%idirs(2) 353 ! Second iteration (kdir1=2) : 354 ! pert1 = rf2%iperts(2) along rf2%idirs(2) 355 ! pert2 = rf2%iperts(1) along rf2%idirs(1) 356 idir1=rf2%idirs(kdir1) 357 ipert1=rf2%iperts(kdir1) 358 shift_dir1=(kdir1-1)*nband_k*size_wf 359 shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj 360 shift_dir1_lambda=(kdir1-1)*2*nband_k**2 361 if(ipert==natom+10 .and. idir<=3) then 362 shift_dir2=0 363 idir2 = idir1 364 ipert2 = ipert1 365 rf_hamk_idir => rf_hamkq 366 else 367 shift_dir2=(2-kdir1)*nband_k*size_wf 368 idir2 = rf2%idirs(3-kdir1) 369 ipert2 = rf2%iperts(3-kdir1) 370 if (kdir1==1) rf_hamk_idir => rf_hamkq 371 if (kdir1==2) rf_hamk_idir => rf_hamk_dir2 372 end if 373 374 ! Load projected WF according to ipert1 and idir1 375 cprj_j => cprj_empty ; cprj_dudk => cprj_empty 376 if (has_cprj_jband) then 377 call pawcprj_free(cprj_jband) 378 ncpgr_loc= 3;if(ipert1==natom+1.or.ipert1==natom+2) ncpgr_loc=1 379 icpgr_loc=-1;if(ipert1==natom+1.or.ipert1==natom+2) icpgr_loc=idir1 380 call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj) 381 call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,& 382 & isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,& 383 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc) 384 end if 385 386 ! LOOP OVER BANDS 387 do jband=1,nband_k ! = band n 388 389 ! Skip bands not treated by current proc 390 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 391 392 shift_band1=(jband-1)*size_wf 393 shift_cprj_band1=(jband-1)*size_cprj 394 shift_jband_lambda=(jband-1)*2*nband_k 395 396 if (abs(occ_k(jband))>tol8) then 397 398 ! Extract first order wavefunction and eigenvalues for jband 399 eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir1_lambda:2*nband_k+shift_jband_lambda+shift_dir1_lambda) 400 cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 401 if (has_dudkprj) cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1:size_cprj+shift_cprj_band1+shift_cprj_dir1) 402 403 ! Compute H^(0) | du/dpert1 > (in h_cwave) and S^(0) | du/dpert1 > (in s_cwave) 404 call rf2_apply_hamiltonian(rf2,cg_jband,cprj_jband,cwave_dudk,cprj_dudk,h_cwave,s_cwave,dtfil,& 405 & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,0,0,ikpt,isppol,dtset%mband,mkmem,& 406 & mpi_enreg,nsppol,print_info,dtset%prtvol,rf_hamk_idir) 407 408 if (gs_hamkq%usepaw==0) s_cwave(:,:)=cwave_dudk(:,:) ! Store | du/dpert1 > in s_cwave 409 410 ! Copy infos in dsusdu 411 dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=s_cwave(:,:)& 412 +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 413 414 if (print_info/=0) then 415 write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 1 kdir1 = ',kdir1,' jband = ',jband 416 call wrtout(std_out,msg) 417 end if 418 419 ! For every occupied iband, we compute : 420 ! < du/dpert2(iband) | H^(0) | du/dpert1(jband) > and add it to lambda_mn 421 ! < du/dpert2(iband) | S^(0) | du/dpert1(jband) > and add it to amn 422 do iband=1,rf2%nband_k ! = band m 423 if (abs(occ_k(iband))>tol8) then 424 shift_band2=(iband-1)*size_wf 425 cwave_dudk => dudk(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2) 426 call rf2_accumulate_bands(rf2,1,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,& 427 jband,print_info,cwave_dudk,h_cwave,s_cwave) 428 end if 429 end do 430 431 ! Extract GS wavefunction for jband 432 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 433 if(has_cprj_jband) cprj_j => cprj_jband(:,1+shift_cprj_band1:size_cprj+shift_cprj_band1) 434 435 if (ipert1==natom+2) then 436 ! Extract ddk and multiply by i : 437 if(idir<=3) then ! in this case : idir1=idir2 438 gvnl1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1) 439 gvnl1(2,:) = dudk(1,1+shift_band1:size_wf+shift_band1) 440 else 441 gvnl1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1) 442 gvnl1(2,:) = dudk_dir2(1,1+shift_band1:size_wf+shift_band1) 443 end if 444 end if 445 446 ! Compute dH/dpert1 | u^(0) > (in h_cwave) and dS/dpert1 | u^(0) > (in s_cwave) 447 call rf2_apply_hamiltonian(rf2,cg_jband,cprj_jband,cwave_j,cprj_j,h_cwave,s_cwave,dtfil,& 448 & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,idir1,ipert1,ikpt,isppol,& 449 & dtset%mband,mkmem,mpi_enreg,nsppol,print_info,dtset%prtvol,rf_hamk_idir) 450 451 ! Copy infos in dsusdu 452 if (gs_hamkq%usepaw==1) then 453 dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1)=s_cwave(:,:)& 454 +dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 455 end if 456 457 if (print_info/=0) then 458 write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 2 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) | dH/dpert1 | u^(0)(jband) > and add it to lambda_mn 464 ! < du/dpert2(iband) | dS/dpert1 | u^(0)(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,2,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,& 470 jband,print_info,cwave_dudk,h_cwave,s_cwave) 471 end if 472 end do 473 474 end if ! empty band test 475 end do ! jband 476 end do ! idir1 477 478 ! Allgather dsusdu 479 nproc_band = xmpi_comm_size(mpi_enreg%comm_band) 480 if (nproc_band>1) then 481 482 my_nband = nband_k/nproc_band;if (mod(nband_k,nproc_band)/=0) my_nband=my_nband+1 483 ABI_ALLOCATE(dsusdu_loc,(2,size_wf*my_nband*rf2%ndir)) 484 ABI_ALLOCATE(dsusdu_gather,(2,size_wf*my_nband*rf2%ndir*nproc_band)) 485 dsusdu_loc(:,:) = zero 486 dsusdu_gather(:,:) = zero 487 488 do kdir1=1,rf2%ndir 489 indb = 1 490 shift_dir1=(kdir1-1)*size_wf*nband_k 491 shift_dir1_loc=(kdir1-1)*size_wf*my_nband 492 do jband=1,nband_k 493 ! Skip bands not treated by current proc 494 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 495 496 shift_band1=(jband-1)*size_wf 497 dsusdu_loc(:,indb+shift_dir1_loc:indb-1+size_wf+shift_dir1_loc) = & 498 dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 499 indb = indb + size_wf 500 end do 501 end do 502 503 call xmpi_allgather(dsusdu_loc,2*size_wf*my_nband*rf2%ndir,dsusdu_gather,mpi_enreg%comm_band,ierr) 504 505 do kdir1=1,rf2%ndir 506 shift_dir1=(kdir1-1)*size_wf*nband_k 507 shift_dir1_loc=(kdir1-1)*size_wf*my_nband 508 do iproc=1,nproc_band 509 shift_proc = (iproc-1)*size_wf*my_nband*rf2%ndir 510 indb = 1 511 do jband=1,my_nband 512 iband = jband+(iproc-1)*my_nband 513 if(iband<=nband_k) then 514 shift_band1=(iband-1)*size_wf 515 dsusdu(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) = & 516 dsusdu_gather(:,indb+shift_dir1_loc+shift_proc:indb-1+size_wf+shift_dir1_loc+shift_proc) 517 end if 518 indb = indb + size_wf 519 end do 520 end do 521 end do 522 ABI_DEALLOCATE(dsusdu_loc) 523 ABI_DEALLOCATE(dsusdu_gather) 524 525 end if 526 527 ! ************************************************************************************************** 528 ! COMPUTATION OF "RHS_Stern", THE LAST PART OF "A_mn" AND A PART OF "Lambda_mn" 529 ! ************************************************************************************************** 530 531 ABI_STAT_ALLOCATE(rf2%RHS_Stern,(2,nband_k*size_wf), ierr) 532 ABI_CHECK(ierr==0, "out of memory in m_rf2 : RHS_Stern") 533 rf2%RHS_Stern(:,:)=zero 534 535 !Computation of terms containing H^(2) 536 if (ipert/=natom+11 .or. gs_hamkq%usepaw==1) then ! Otherwise H^(2) = 0 537 538 !Load projected WF according to ipert and idir 539 cprj_j => cprj_empty 540 if (has_cprj_jband) then 541 call pawcprj_free(cprj_jband) 542 ncpgr_loc= 3;if(ipert==natom+1.or.ipert==natom+2) ncpgr_loc=1 543 icpgr_loc=-1;if(ipert==natom+1.or.ipert==natom+2) icpgr_loc=idir 544 call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj) 545 call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,& 546 & isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,& 547 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc) 548 end if 549 550 if (ipert==natom+10) then 551 rf_hamk_idir => rf_hamkq ! all info are in rf_hamkq 552 else if (ipert==natom+11) then 553 rf_hamk_idir => rf_hamk_dir2 ! all info are in rf_hamk_dir2 554 end if 555 556 do jband=1,nband_k 557 558 ! Skip bands not treated by current proc 559 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 560 561 if (abs(occ_k(jband))>tol8) then 562 shift_band1=(jband-1)*size_wf 563 shift_cprj_band1=(jband-1)*size_cprj 564 shift_jband_lambda=(jband-1)*2*nband_k 565 566 ! Extract GS wavefunction 567 eig1_k_jband => null() 568 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 569 if(has_cprj_jband) cprj_j => cprj_jband(:,1+shift_cprj_band1:size_cprj+shift_cprj_band1) 570 571 if (ipert==natom+11) then 572 ! Extract ddk and multiply by i : 573 if(idir<=3) then ! in this case : idir1=idir2 574 gvnl1(1,:) = -dudk(2,1+shift_band1:size_wf+shift_band1) 575 gvnl1(2,:) = dudk(1,1+shift_band1:size_wf+shift_band1) 576 else 577 gvnl1(1,:) = -dudk_dir2(2,1+shift_band1:size_wf+shift_band1) 578 gvnl1(2,:) = dudk_dir2(1,1+shift_band1:size_wf+shift_band1) 579 end if 580 end if 581 582 ! Compute : d^2H/(dpert1 dpert2)|u^(0)> (in h_cwave) 583 ! and : d^2S/(dpert1 dpert2)|u^(0)> (in s_cwave) 584 call rf2_apply_hamiltonian(rf2,cg_jband,cprj_jband,cwave_j,cprj_j,h_cwave,s_cwave,dtfil,& 585 & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,idir,ipert,ikpt,isppol,dtset%mband,& 586 & mkmem,mpi_enreg,nsppol,print_info,dtset%prtvol,rf_hamk_idir) 587 588 if (print_info/=0) then 589 write(msg,'(a,i2)') 'RF2 TEST before accumulate_bands choice = 3 jband = ',jband 590 call wrtout(std_out,msg) 591 end if 592 593 ! For every occupied iband, we compute : 594 ! < u^(0)(iband) | d^2H/(dpert1 dpert2) | u^(0)(jband) > and add it to lambda_mn 595 ! < u^(0)(iband) | d^2S/(dpert1 dpert2) | u^(0)(jband) > and add it to amn 596 do iband=1,rf2%nband_k ! = band m 597 if (abs(occ_k(iband))>tol8) then 598 shift_band2=(iband-1)*size_wf 599 cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg) 600 if(ipert == natom+10) then 601 ipert1 = natom+1 602 ipert2 = natom+1 603 else 604 ipert1 = natom+1 605 ipert2 = natom+2 606 end if 607 call rf2_getidirs(idir,idir1,idir2) 608 call rf2_accumulate_bands(rf2,3,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,& 609 jband,print_info,cwave_i,h_cwave,s_cwave) 610 end if 611 end do 612 613 ! Add d^2H/(dk_dir1 dk_dir2)|u^(0)> to RHS_Stern : 614 if (gs_hamkq%usepaw==1) h_cwave(:,:)=h_cwave(:,:)-eig0_k(jband)*s_cwave(:,:) ! if PAW : we add H^(2)-eps^(0) S^(2) 615 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 616 call cg_zaxpy(size_wf,(/one,zero/),h_cwave,rhs_j) 617 618 end if ! empty band test 619 end do ! jband 620 end if ! H^(2) exists 621 622 !Computation of terms containing H^(1) 623 do kdir1=1,rf2%ndir 624 ! First iteration (kdir1=1) : 625 ! pert1 = rf2%iperts(1) along rf2%idirs(1) 626 ! pert2 = rf2%iperts(2) along rf2%idirs(2) 627 ! Second iteration (kdir1=2) : 628 ! pert1 = rf2%iperts(2) along rf2%idirs(2) 629 ! pert2 = rf2%iperts(1) along rf2%idirs(1) 630 shift_dir1=(kdir1-1)*nband_k*size_wf 631 shift_cprj_dir1=(kdir1-1)*nband_k*size_cprj 632 shift_dir1_lambda=(kdir1-1)*2*nband_k**2 633 idir1=rf2%idirs(kdir1) 634 ipert1=rf2%iperts(kdir1) 635 if(ipert==natom+10 .and. idir<=3) then 636 idir2=idir1 637 ipert2=ipert1 638 shift_dir2=0 639 shift_dir2_lambda=0 640 rf_hamk_idir => rf_hamkq 641 else 642 idir2=rf2%idirs(2-kdir1+1) 643 ipert2=rf2%iperts(2-kdir1+1) 644 shift_dir2=(2-kdir1)*nband_k*size_wf 645 shift_dir2_lambda=(2-kdir1)*2*nband_k**2 646 if (kdir1==1) rf_hamk_idir => rf_hamk_dir2 ! dir2 647 if (kdir1==2) rf_hamk_idir => rf_hamkq ! dir1 648 end if 649 650 ! Load projected WF according to ipert2 and idir2 651 cprj_j => cprj_empty ; ; cprj_dudk => cprj_empty 652 if (has_cprj_jband) then 653 call pawcprj_free(cprj_jband) 654 ncpgr_loc= 3;if(ipert2==natom+1.or.ipert2==natom+2) ncpgr_loc=1 655 icpgr_loc=-1;if(ipert2==natom+1.or.ipert2==natom+2) icpgr_loc=idir2 656 call pawcprj_alloc(cprj_jband,ncpgr_loc,gs_hamkq%dimcprj) 657 call pawcprj_get(gs_hamkq%atindx1,cprj_jband,cprj,natom,1,ibg,ikpt,iorder_cprj,& 658 & isppol,dtset%mband,mkmem,natom,nband_k,nband_k,gs_hamkq%nspinor,nsppol,dtfil%unpaw,& 659 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,ncpgr=3,icpgr=icpgr_loc) 660 end if 661 662 do jband=1,nband_k 663 664 ! Skip bands not treated by current proc 665 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 666 667 if (abs(occ_k(jband))>tol8) then 668 shift_band1=(jband-1)*size_wf 669 shift_cprj_band1=(jband-1)*size_cprj 670 shift_jband_lambda=(jband-1)*2*nband_k 671 672 ! Extract first order wavefunction | du/dpert1 > and eigenvalues 673 eig1_k_jband => eig1_k_stored(1+shift_jband_lambda+shift_dir2_lambda:2*nband_k+shift_jband_lambda+shift_dir2_lambda) 674 cwave_dudk => dudk(:,1+shift_band1+shift_dir1:size_wf+shift_band1+shift_dir1) 675 if (has_dudkprj) cprj_dudk => dudkprj(:,1+shift_cprj_band1+shift_cprj_dir1:size_cprj+shift_cprj_band1+shift_cprj_dir1) 676 677 if (ipert2==natom+2) then 678 ! Extract dkdk and multiply by i : 679 gvnl1(1,:) = -dudkdk(2,1+shift_band1:size_wf+shift_band1) 680 gvnl1(2,:) = dudkdk(1,1+shift_band1:size_wf+shift_band1) 681 end if 682 683 ! Compute dH/dpert2 | du/dpert1 > (in h_cwave) and dS/dpert2 | du/dpert1 > (in s_cwave) 684 call rf2_apply_hamiltonian(rf2,cg_jband,cprj_jband,cwave_dudk,cprj_dudk,h_cwave,s_cwave,dtfil,& 685 & eig0_k,eig1_k_jband,jband,gs_hamkq,gvnl1,idir2,ipert2,ikpt,isppol,& 686 & dtset%mband,mkmem,mpi_enreg,nsppol,print_info,dtset%prtvol,rf_hamk_idir) 687 688 if (print_info/=0) then 689 write(msg,'(2(a,i2))') 'RF2 TEST before accumulate_bands choice = 4 kdir1 = ',kdir1,' jband = ',jband 690 call wrtout(std_out,msg) 691 end if 692 693 ! For every occupied iband, we compute : 694 ! < u^(0)(iband) | dH/dpert2 | du/dpert1(jband) > and add it to lambda_mn 695 ! < u^(0)(iband) | dS/dpert2 | du/dpert1(jband) > and add it to amn 696 do iband=1,rf2%nband_k ! = band m 697 if (abs(occ_k(iband))>tol8) then 698 shift_band2=(iband-1)*size_wf 699 cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg) 700 call rf2_accumulate_bands(rf2,4,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,& 701 jband,print_info,cwave_i,h_cwave,s_cwave) 702 end if 703 end do 704 705 ! Add dH/dpert2 | du/dpert1 > to RHS_Stern : 706 if (gs_hamkq%usepaw==1) h_cwave(:,:)=h_cwave(:,:)-eig0_k(jband)*s_cwave(:,:) ! if PAW : we add H^(1)-eps^(0) S^(1) 707 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 708 call cg_zaxpy(size_wf,(/factor*one,zero/),h_cwave,rhs_j) 709 710 ! Compute : -factor * sum_iband ( dLambda/dpert1_{iband,jband} * dsusdu_{iband} ) 711 do iband=1,nband_k 712 if (abs(occ_k(iband))>tol8) then ! if empty band, nothing to do 713 714 ! Extract lambda_ij(iband,jband) for dir1 715 lambda_ij(1)=eig1_k_stored(2*iband-1+shift_jband_lambda+shift_dir1_lambda) 716 lambda_ij(2)=eig1_k_stored(2*iband +shift_jband_lambda+shift_dir1_lambda) 717 718 ! Extract dsusdu for iband and pert2 (in cwave_i) 719 shift_band2=(iband-1)*size_wf 720 cwave_i => dsusdu(:,1+shift_band2+shift_dir2:size_wf+shift_band2+shift_dir2) 721 722 ! Compute Lambda_{iband,jband} * dsusdu_{iband} and add it to RHS_Stern 723 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 724 call cg_zaxpy(size_wf,(/-factor*lambda_ij(1),-factor*lambda_ij(2)/),cwave_i,rhs_j) !do not forget the minus sign! 725 726 end if ! empty iband test 727 end do ! iband 728 729 end if ! empty jband test 730 end do ! jband 731 732 end do ! kdir1 733 734 ABI_DEALLOCATE(gvnl1) 735 ABI_DEALLOCATE(h_cwave) 736 ABI_DEALLOCATE(s_cwave) 737 ABI_DEALLOCATE(cg_jband) 738 ABI_DEALLOCATE(dudk) 739 ABI_DEALLOCATE(dudkdk) 740 ABI_DEALLOCATE(dudk_dir2) 741 ABI_DEALLOCATE(dsusdu) 742 ABI_DEALLOCATE(eig1_k_stored) 743 if (has_cprj_jband) call pawcprj_free(cprj_jband) 744 ABI_DATATYPE_DEALLOCATE(cprj_jband) 745 if (has_dudkprj) call pawcprj_free(dudkprj) 746 ABI_DATATYPE_DEALLOCATE(dudkprj) 747 748 ! Compute the part of 2nd order wavefunction that belongs to the space of empty bands 749 do jband=1,nband_k 750 751 ! Skip bands not treated by current proc 752 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 753 754 shift_band1=(jband-1)*size_wf 755 if (abs(occ_k(jband))>tol8) then 756 invocc = one/occ_k(jband) 757 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 758 do iband=1,nband_k 759 if (iband /= jband) then 760 if (print_info/=0) then 761 if (abs(occ_k(iband) - occ_k(jband)) > tol12 .and. occ_k(iband) > tol8) then 762 write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband 763 call wrtout(std_out,msg) 764 write(msg,'(a)') 'ERROR : occ_k(iband) /= occ_k(jband) (and both are >0)' 765 call wrtout(std_out,msg) 766 end if 767 if (abs(eig0_k(iband) - eig0_k(jband)) < tol8 ) then 768 write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband 769 call wrtout(std_out,msg) 770 write(msg,'(a,es22.13e3)') 'WARNING : DEGENERATE BANDS Eig(jband) = Eig(jband) = ',eig0_k(jband) 771 call wrtout(std_out,msg) 772 end if 773 if ( (eig0_k(iband) - eig0_k(jband) < -tol12) .and. (jband < iband) ) then 774 write(msg,'(a,i2,a,i2)') 'RF2 TEST ACTIVE SPACE : jband = ',jband,' iband = ',iband 775 call wrtout(std_out,msg) 776 write(msg,'(a)') 'ERROR : Eig(jband) < Eig(iband) with jband < iband' 777 call wrtout(std_out,msg) 778 write(msg,'(a,es22.13e3)') 'Eig(jband) = ',eig0_k(jband) 779 call wrtout(std_out,msg) 780 write(msg,'(a,es22.13e3)') 'Eig(iband) = ',eig0_k(iband) 781 call wrtout(std_out,msg) 782 end if 783 end if ! end tests 784 if ( abs(occ_k(iband))<tol8 ) then ! for empty bands only 785 shift_band2=(iband-1)*size_wf 786 cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg) 787 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_i,rhs_j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 788 ! Store it in a_mn 789 ! /!\ There is a factor "-2" to simplify the use of amn in the following. 790 ! /!\ Occupied and empty bands will be treated in a same way. 791 rf2%amn(:,iband+(jband-1)*nband_k)=-two*rocceig(iband,jband)*invocc*(/dotr,doti/)& 792 +rf2%amn(:,iband+(jband-1)*nband_k) 793 end if ! empty band test 794 end if ! iband \= jband 795 end do ! iband 796 end if ! empty band test 797 end do ! jband 798 799 ! ************************************************************************************************** 800 ! COMPUTATION OF "dcwavef" AND "Lambda_mn" FROM "A_mn" 801 ! ************************************************************************************************** 802 803 ABI_STAT_ALLOCATE(rf2%dcwavef,(2,nband_k*size_wf), ierr) 804 ABI_CHECK(ierr==0, "out of memory in m_rf2 : dcwavef") 805 rf2%dcwavef=zero 806 807 do jband=1,nband_k 808 809 ! Skip bands not treated by current proc 810 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 811 812 shift_band1=(jband-1)*size_wf 813 if (abs(occ_k(jband))>tol8) then 814 do iband=1,nband_k 815 shift_band2=(iband-1)*size_wf 816 817 ! Extract GS wavefunction for iband 818 cwave_i => cg(:,1+shift_band2+icg:size_wf+shift_band2+icg) 819 820 call cg_zaxpy(size_wf,-half*rf2%amn(:,iband+(jband-1)*nband_k), & 821 & cwave_i,rf2%dcwavef(:,1+shift_band1)) 822 823 if (abs(occ_k(iband))>tol8 .and. abs(occ_k(jband))>tol8) then 824 rf2%lambda_mn(:,iband+(jband-1)*nband_k) = rf2%lambda_mn(:,iband+(jband-1)*nband_k) & 825 -half*(eig0_k(iband)+eig0_k(jband))*rf2%amn(:,iband+(jband-1)*nband_k) 826 827 eig1_k(2*iband-1+(jband-1)*2*nband_k) = rf2%lambda_mn(1,iband+(jband-1)*nband_k) 828 eig1_k(2*iband +(jband-1)*2*nband_k) = rf2%lambda_mn(2,iband+(jband-1)*nband_k) 829 830 end if ! empty band test 831 end do ! iband 832 end if ! empty band test 833 end do ! jband 834 835 ! For the following, "rf2%lambda_mn" and "rf2%RHS_Stern" must be computed for every bands 836 call xmpi_barrier(mpi_enreg%comm_band) 837 838 ! ************************************************************************************************** 839 ! FINAL TEST 840 ! ************************************************************************************************** 841 842 tol_final = tol6 843 if (print_info/=0) then 844 do jband=1,nband_k 845 846 ! Skip bands not treated by current proc 847 if(mpi_enreg%proc_distrb(ikpt,jband,isppol)/=me) cycle 848 849 if (abs(occ_k(jband))>tol8) then 850 ! write(msg,'(3(a,i2))') 'RF2 TEST FINAL : ipert=',ipert-natom,' idir=',idir,' jband=',jband 851 ! call wrtout(std_out,msg) 852 shift_band1=(jband-1)*size_wf 853 rhs_j => rf2%RHS_Stern(:,1+shift_band1:size_wf+shift_band1) 854 cwave_j => cg(:,1+shift_band1+icg:size_wf+shift_band1+icg) 855 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,rhs_j,cwave_j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 856 ! write(msg,'(2(a,es22.13E3))') 'RF2 TEST FINAL : dot =',dotr,',',doti 857 ! call wrtout(std_out,msg) 858 ! write(msg,'(2(a,es22.13E3))') 'RF2 TEST FINAL : lambda_jj =',rf2%lambda_mn(1,jband+(jband-1)*nband_k),& 859 ! ',',rf2%lambda_mn(2,jband+(jband-1)*nband_k) 860 ! call wrtout(std_out,msg) 861 dotr = dotr - rf2%lambda_mn(1,jband+(jband-1)*nband_k) 862 doti = doti - (-rf2%lambda_mn(2,jband+(jband-1)*nband_k)) ! be careful : complex conjugate of lambda_mn 863 ! NOTE : 864 ! If lambda_nn^(2) can be comlex (possible for ipert==natom+11, as H^(1) and H^(2) are not hermitian), 865 ! the test works if we take the conjugate here. 866 ! For real systems, lambda_nn^(2) is always real (empirical assumption...). 867 dotr = sqrt(dotr**2+doti**2) 868 if (dotr > tol_final) then 869 write(msg,'(a,i2,a,es22.13E3)') 'RF2 TEST FINAL iband = ',jband,' : NOT PASSED dotr = ',dotr 870 call wrtout(std_out,msg) 871 else 872 write(msg,'(a,i2,a,es22.13E3,a,es7.1E2)') & 873 'RF2 TEST FINAL iband = ',jband,' : OK. |test| = ',dotr,' < ',tol_final 874 call wrtout(std_out,msg) 875 end if 876 end if 877 end do 878 end if 879 880 ! ************************************************************************************************** 881 ! JOB FINISHED 882 ! ************************************************************************************************** 883 884 ! Deallocations of arrays 885 if (print_info==0) ABI_DEALLOCATE(rf2%amn) 886 887 ! call timab(566,2,tsec) 888 889 DBG_EXIT("COLL") 890 891 end subroutine rf2_init