TABLE OF CONTENTS
ABINIT/wf_mixing [ Functions ]
NAME
wf_mixing
FUNCTION
Mixing of wavefunctions in the outer loop of a double loop SCF approach. Different algorithms are implemented, depending on the value of wfmixalg.
COPYRIGHT
Copyright (C) 2017-2018 ABINIT group (XG,MT,FJ) 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
atindx1(dtset%natom)=index table for atoms, inverse of atindx dtset <type(dataset_type)>=all input variables in this dataset istep=number of call the routine (usually the outer loop in the SCF double loop) mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of cprj array mpi_enreg=information about MPI parallelization nattyp(dtset%ntypat)=number of atoms of each type in cell. npwarr(nkpt)=number of planewaves in basis at this k point pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
SIDE EFFECTS
cg(2,mcg)= plane wave wavefunction coefficient Value from previous SCF cycle is input and stored in some form Extrapolated value is output cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors Value from previous SCF cycle is input and stored in some form Extrapolated value is output scf_history_wf <type(scf_history_type)>=arrays obtained from previous SCF cycles
PARENTS
scfcv
CHILDREN
cgcprj_cholesky,dotprod_set_cgcprj,dotprodm_sumdiag_cgcprj lincom_cgcprj,pawcprj_alloc,pawcprj_axpby,pawcprj_free,pawcprj_get pawcprj_getdim,pawcprj_lincom,pawcprj_put,timab,xmpi_sum,zgesv
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine wf_mixing(atindx1,cg,cprj,dtset,istep,mcg,mcprj,mpi_enreg,& 55 & nattyp,npwarr,pawtab,scf_history_wf) 56 57 use defs_basis 58 use defs_abitypes 59 use m_scf_history 60 use m_xmpi 61 use m_profiling_abi 62 use m_errors 63 use m_cgtools 64 65 use m_pawtab, only : pawtab_type 66 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_get, pawcprj_lincom, & 67 & pawcprj_free, pawcprj_axpby, pawcprj_put, pawcprj_getdim 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'wf_mixing' 73 use interfaces_18_timing 74 use interfaces_32_util 75 use interfaces_66_wfs 76 !End of the abilint section 77 78 implicit none 79 80 !Arguments ------------------------------------ 81 !scalars 82 integer,intent(in) :: istep,mcg,mcprj 83 type(MPI_type),intent(in) :: mpi_enreg 84 type(dataset_type),intent(in) :: dtset 85 type(scf_history_type),intent(inout) :: scf_history_wf 86 !arrays 87 integer,intent(in) :: atindx1(dtset%natom),nattyp(dtset%ntypat) 88 integer,intent(in) :: npwarr(dtset%nkpt) 89 real(dp), intent(inout) :: cg(2,mcg) 90 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj) 91 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 92 93 !Local variables------------------------------- 94 !scalars 95 integer :: hermitian 96 integer :: ibdmix,ibdsp,ibg,ibg_hist,icg,icg_hist 97 integer :: ierr,ikpt,indh,ind_biorthog,ind_biorthog_eff,ind_newwf,ind_residual,inplace 98 integer :: iset2,isppol,istep_cycle,istep_new,istwf_k,kk,me_distrb,my_nspinor 99 integer :: nband_k,nbdmix,npw_k,nset1,nset2,ntypat 100 integer :: shift_set1,shift_set2,spaceComm_band,spare_mem,usepaw,wfmixalg 101 real(dp) :: alpha,beta 102 complex(dpc) :: sum_coeffs 103 !arrays 104 integer,allocatable :: ipiv(:),dimcprj(:) 105 real(dp) :: tsec(2) 106 real(dp),allocatable :: al(:,:),mmn(:,:,:) 107 real(dp),allocatable :: dotprod_res(:,:,:),dotprod_res_k(:,:,:),res_mn(:,:,:),smn(:,:,:) 108 complex(dpc),allocatable :: coeffs(:) 109 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kh(:,:) 110 111 ! ************************************************************************* 112 113 !DEBUG 114 !write(std_out,*) 115 !write(std_out,*)' wf_mixing : enter, istep= ',istep 116 !call flush(std_out) 117 !write(std_out,*)' istep,scf_history_wf%alpha=',istep,scf_history_wf%alpha 118 !write(std_out,*)' cg(1,1)=',cg(1,1) 119 !write(std_out,*)' scf_history_wf%cg(1,1,1:5)=',scf_history_wf%cg(1,1,1:5) 120 !ABI_ALLOCATE(cg_ref,(2,mcg)) 121 !cg_ref(:,:)=cg(:,:) 122 !ABI_DATATYPE_ALLOCATE(cprj_ref,(dtset%natom,mcprj)) 123 !cprj_ref(:,:)=cprj(:,:) 124 ! write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',& 125 !& scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2) 126 ! call flush(std_out) 127 !ENDDEBUG 128 129 if (istep==0) return 130 131 ntypat=dtset%ntypat 132 usepaw=dtset%usepaw 133 wfmixalg=scf_history_wf%wfmixalg 134 nbdmix=dtset%nbandhf 135 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 136 me_distrb=mpi_enreg%me_kpt 137 spaceComm_band=xmpi_comm_self 138 139 spare_mem=0 140 if(scf_history_wf%history_size==wfmixalg-1)spare_mem=1 141 142 !scf_history_wf%alpha contains dtset%wfmix 143 alpha=scf_history_wf%alpha 144 beta=one-scf_history_wf%alpha 145 146 icg=0 147 icg_hist=0 148 ibg=0 149 ibg_hist=0 150 151 !Useful array 152 ABI_ALLOCATE(dimcprj,(dtset%natom)) 153 if (usepaw==1) then 154 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,ntypat,dtset%typat,pawtab,'O') 155 end if 156 157 if(istep==1)then 158 do indh=1,scf_history_wf%history_size 159 call pawcprj_alloc(scf_history_wf%cprj(:,:,indh),0,dimcprj) 160 end do 161 end if 162 163 ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,my_nspinor*nbdmix)) 164 ABI_DATATYPE_ALLOCATE(cprj_kh,(dtset%natom,my_nspinor*nbdmix)) 165 if(usepaw==1) then 166 call pawcprj_alloc(cprj_k,0,dimcprj) 167 call pawcprj_alloc(cprj_kh,0,dimcprj) 168 end if 169 ABI_ALLOCATE(smn,(2,nbdmix,nbdmix)) 170 ABI_ALLOCATE(mmn,(2,nbdmix,nbdmix)) 171 172 if(wfmixalg>2)then 173 nset1=1 174 nset2=min(istep-1,wfmixalg-1) 175 ABI_ALLOCATE(dotprod_res_k,(2,1,nset2)) 176 ABI_ALLOCATE(dotprod_res,(2,1,nset2)) 177 ABI_ALLOCATE(res_mn,(2,wfmixalg-1,wfmixalg-1)) 178 dotprod_res=zero 179 if(istep==1)then 180 scf_history_wf%dotprod_sumdiag_cgcprj_ij=zero 181 end if 182 end if 183 184 !Explanation for the index for the wavefunction stored in scf_history_wf 185 !The reference is the cg+cprj output after the wf optimization at istep 1. 186 !It comes as input to the present routine as cgcprj input at step 2, and is usually found at indh=1. 187 188 !In the simple mixing case (wfmixalg==1), the reference is never stored, because it is used "on-the-fly" to biothogonalize the 189 !previous input (that was stored in indh=1), then generate the next input, which is stored again in indh=1 190 191 !When the storage is not spared: 192 !- the values of indh from 2 to wfmixalg store the (computed here) biorthogonalized input cgcprj, then the residual 193 !- the values of indh from wfmixalg+1 to 2*wfmixalg-1 store the biorthogonalized output cgcprj (coming as argument) 194 195 !First step 196 if (istep==1 .or. (wfmixalg==2 .and. abs(scf_history_wf%alpha-one)<tol8) ) then 197 198 indh=2 ! This input wavefunction is NOT the reference 199 if(wfmixalg==2)indh=1 ! But this does not matter in the simple mixing case that has history_size=1 200 201 ! Simply store the wavefunctions and cprj. However, nband_k might be different from nbandhf... 202 ! LOOP OVER SPINS 203 do isppol=1,dtset%nsppol 204 205 ! BIG FAT k POINT LOOP 206 do ikpt=1,dtset%nkpt 207 208 ! Select k point to be treated by this proc 209 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 210 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 211 212 npw_k=npwarr(ikpt) 213 214 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 215 if(usepaw==1) then 216 ! scf_history_wf%cprj(:,ibg_hist+1:ibg_hist+my_nspinor*nbdmix,1)=cprj(:,ibg+1:ibg+my_nspinor*nbdmix) 217 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,& 218 & dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,& 219 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 220 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,1,isppol,& 221 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 222 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 223 end if 224 225 ! Update the counters 226 ibg=ibg+my_nspinor*nband_k 227 ibg_hist=ibg_hist+my_nspinor*nbdmix 228 icg=icg+my_nspinor*nband_k*npw_k 229 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 230 231 end do 232 end do 233 234 else 235 ! From istep==2 236 237 ! First part of the computation : biorthogonalization, and computation of the residual (possibly, prediction of the next input in the case of simple mixing) 238 ! Index for the wavefunctions stored in scf_history_wf whose scalar products with the argument cgcprj will have to be computed. 239 indh=1 ! This input wavefunction is the reference 240 if(wfmixalg/=2 .and. istep==2)indh=2 ! except for istep=2 in the rmm-diis 241 242 if(wfmixalg>2)then 243 ! istep inside the cycle defined by wfmixalg, and next index. Then, indices of the wavefunction sets. 244 istep_cycle=mod((istep-2),wfmixalg-1) 245 istep_new=mod((istep-1),wfmixalg-1) 246 ind_biorthog=1+wfmixalg+istep_cycle 247 ind_residual=2+istep_cycle 248 ind_newwf=2+istep_new 249 shift_set1=ind_residual-1 250 shift_set2=1 251 end if 252 253 ! LOOP OVER SPINS 254 do isppol=1,dtset%nsppol 255 256 ! BIG FAT k POINT LOOP 257 do ikpt=1,dtset%nkpt 258 259 ! Select k point to be treated by this proc 260 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 261 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 262 263 istwf_k=dtset%istwfk(ikpt) 264 npw_k=npwarr(ikpt) 265 266 ! Biorthogonalization 267 268 if(usepaw==1) then 269 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,ibg,ikpt,1,isppol,dtset%mband,& 270 & dtset%mkmem,dtset%natom,nbdmix,nband_k,my_nspinor,dtset%nsppol,0,& 271 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 272 call pawcprj_get(atindx1,cprj_kh,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,1,isppol,& 273 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,my_nspinor,dtset%nsppol,0,& 274 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 275 end if !end usepaw=1 276 277 hermitian=0 278 if(wfmixalg==2 .or. istep==2)then 279 call dotprod_set_cgcprj(atindx1,cg,scf_history_wf%cg(:,:,indh),cprj_k,cprj_kh,dimcprj,hermitian,& 280 & 0,0,icg,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,& 281 & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw) 282 else 283 call dotprod_set_cgcprj(atindx1,scf_history_wf%cg(:,:,indh),cg,cprj_kh,cprj_k,dimcprj,hermitian,& 284 & 0,0,icg,icg_hist,ikpt,isppol,istwf_k,nbdmix,mcg,mcg,mcprj,mcprj,dtset%mkmem,& 285 & mpi_enreg,dtset%natom,nattyp,nbdmix,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,smn,usepaw) 286 end if 287 288 ! Invert S matrix, that is NOT hermitian. 289 ! Calculate M=S^-1 290 mmn=zero 291 do kk=1,nbdmix 292 mmn(1,kk,kk)=one 293 end do 294 295 ABI_ALLOCATE(ipiv,(nbdmix)) 296 ! The smn is destroyed by the following inverse call 297 call zgesv(nbdmix,nbdmix,smn,nbdmix,ipiv,mmn,nbdmix,ierr) 298 ABI_DEALLOCATE(ipiv) 299 !DEBUG 300 if(ierr/=0)then 301 MSG_ERROR(' The call to cgesv general inversion routine failed') 302 end if 303 !ENDDEBUG 304 305 ! The M matrix is used to compute the biorthogonalized set of wavefunctions, and to store it at the proper place 306 if(wfmixalg==2 .or. istep==2)then 307 inplace=1 308 call lincom_cgcprj(mmn,scf_history_wf%cg(:,:,indh),cprj_kh,dimcprj,& 309 & icg_hist,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw) 310 else 311 inplace=0 312 call lincom_cgcprj(mmn,cg,cprj_k,dimcprj,& 313 & icg,inplace,mcg,my_nspinor*nbdmix,dtset%natom,nbdmix,nbdmix,npw_k,my_nspinor,usepaw,& 314 & cgout=scf_history_wf%cg(:,:,ind_biorthog),cprjout=scf_history_wf%cprj(:,:,ind_biorthog),icgout=icg_hist) 315 end if 316 317 ! The biorthogonalised set of wavefunctions is now stored at the proper place 318 319 ! Finalize this first part of the computation, depending on the algorithm and the step. 320 321 if(wfmixalg==2)then 322 323 ! Wavefunction extrapolation, simple mixing case 324 ! alpha contains dtset%wfmix, beta contains one-alpha 325 cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)=& 326 & alpha*cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix)& 327 & +beta*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh) 328 if(usepaw==1) then 329 do ibdmix=1,nbdmix 330 call pawcprj_axpby(beta,alpha,cprj_kh(:,ibdmix:ibdmix),cprj_k(:,ibdmix:ibdmix)) 331 end do ! end loop on ibdmix 332 end if 333 334 ! Back to usual orthonormalization 335 call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,& 336 & mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw) 337 338 ! Store the newly extrapolated wavefunctions, orthonormalized, in scf_history_wf 339 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,indh)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 340 if(usepaw==1) then 341 do ibdmix=1,nbdmix 342 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,indh),dtset%natom,1,ibg_hist,ikpt,1,isppol,& 343 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 344 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 345 end do ! end loop on ibdmix 346 end if 347 348 else ! wfmixalg/=2 349 ! RMM-DIIS 350 351 if (istep==2)then 352 ! Store the argument wf as the reference for all future steps, in scf_history_wf with index 1. 353 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 354 if(usepaw==1) then 355 do ibdmix=1,nbdmix 356 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,1),dtset%natom,1,ibg_hist,ikpt,1,isppol,& 357 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 358 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 359 end do ! end loop on ibdmix 360 end if 361 end if 362 363 ind_biorthog_eff=ind_biorthog 364 if(istep==2)ind_biorthog_eff=1 ! The argument wf has not been stored in ind_biorthog 365 366 ! Compute the residual of the wavefunctions for this istep, 367 ! that replaces the previously stored set of (biorthogonalized) input wavefunctions 368 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual)=& 369 & scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)& 370 & -scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual) 371 if(usepaw==1) then 372 do ibdmix=1,nbdmix 373 call pawcprj_axpby(one,-one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff),& 374 & scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual)) 375 end do ! end loop on ibdmix 376 end if 377 378 ! Compute the new scalar products to fill the res_mn matrix 379 call dotprodm_sumdiag_cgcprj(atindx1,scf_history_wf%cg,scf_history_wf%cprj,dimcprj,& 380 & ibg,icg,ikpt,isppol,istwf_k,nbdmix,mcg,mcprj,dtset%mkmem,& 381 & mpi_enreg,scf_history_wf%history_size,dtset%natom,nattyp,nbdmix,npw_k,nset1,nset2,my_nspinor,dtset%nsppol,ntypat,& 382 & shift_set1,shift_set2,pawtab,dotprod_res_k,usepaw) 383 384 dotprod_res=dotprod_res+dotprod_res_k 385 386 ! scf_history_wf for index ind_biorthog will contain the extrapolated wavefunctions (and no more the output of the SCF loop). 387 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog)=& 388 & scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_biorthog_eff)+& 389 & (alpha-one)*scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_residual) 390 if(usepaw==1) then 391 do ibdmix=1,nbdmix 392 if(ind_biorthog/=ind_biorthog_eff)then 393 scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog)=scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog_eff) 394 end if 395 call pawcprj_axpby((alpha-one),one,scf_history_wf%cprj(:,ibdmix:ibdmix,ind_residual),& 396 & scf_history_wf%cprj(:,ibdmix:ibdmix,ind_biorthog)) 397 end do ! end loop on ibdmix 398 end if 399 400 end if 401 402 ibg=ibg+my_nspinor*nband_k 403 ibg_hist=ibg_hist+my_nspinor*nbdmix 404 icg=icg+my_nspinor*nband_k*npw_k 405 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 406 407 ! End big k point loop 408 end do 409 ! End loop over spins 410 end do 411 412 end if ! istep>=2 413 414 if(wfmixalg>2 .and. istep>1)then 415 416 !DEBUG 417 ! write(std_out,*)' ' 418 ! write(std_out,*)' Entering the residual minimisation part ' 419 ! write(std_out,*)' ' 420 ! call flush(std_out) 421 !ENDDEBUG 422 423 call timab(48,1,tsec) 424 call xmpi_sum(dotprod_res,mpi_enreg%comm_kpt,ierr) 425 call timab(48,2,tsec) 426 427 scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set1,1+shift_set2:nset2+shift_set2)=dotprod_res(:,1,1:nset2) 428 scf_history_wf%dotprod_sumdiag_cgcprj_ij(1,1+shift_set2:nset2+shift_set2,1+shift_set1)=dotprod_res(1,1,1:nset2) 429 scf_history_wf%dotprod_sumdiag_cgcprj_ij(2,1+shift_set2:nset2+shift_set2,1+shift_set1)=-dotprod_res(2,1,1:nset2) 430 431 end if ! wfmixalg>2 and istep>1 432 433 if(wfmixalg>2 .and. istep>2)then 434 435 ! Extract the relevant matrix R_mn 436 res_mn(:,1:nset2,1:nset2)=& 437 & scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,1+shift_set2:nset2+shift_set2,1+shift_set2:nset2+shift_set2) 438 439 !DEBUG 440 ! write(std_out,*)' The matrix res_mn(:,1:nset2,1:nset2) is :' 441 ! write(std_out,*)res_mn(:,1:nset2,1:nset2) 442 ! call flush(std_out) 443 !ENDDEBUG 444 445 ! Solve R_mn \alpha_n = 1_m 446 ABI_ALLOCATE(ipiv,(nset2)) 447 ABI_ALLOCATE(coeffs,(nset2)) 448 coeffs(:)=cone 449 ! The res_mn is destroyed by the following inverse call 450 call zgesv(nset2,1,res_mn,wfmixalg-1,ipiv,coeffs,nset2,ierr) 451 ABI_DEALLOCATE(ipiv) 452 ! The coefficients must sum to one 453 sum_coeffs=sum(coeffs) 454 coeffs=coeffs/sum_coeffs 455 456 !DEBUG 457 ! write(std_out,*)' The coefficients that minimize the residual have been found' 458 ! write(std_out,*)' coeffs =',coeffs 459 ! call flush(std_out) 460 !ENDDEBUG 461 end if ! wfmixalg>2 and istep>2 462 463 if(wfmixalg>2 .and. istep>1)then 464 465 ! Find the new "input" wavefunction, bi-orthogonalized, and store it replacing the adequate "old" input wavefunction. 466 467 icg=0 468 icg_hist=0 469 ibg=0 470 ibg_hist=0 471 ABI_ALLOCATE(al,(2,nset2)) 472 if(istep>2)then 473 do iset2=1,nset2 474 al(1,iset2)=real(coeffs(iset2)) ; al(2,iset2)=aimag(coeffs(iset2)) 475 end do 476 else 477 al(1,1)=one ; al(2,1)=zero 478 end if 479 480 !DEBUG 481 ! write(std_out,*)' Overload the coefficients, in order to simulate a simple mixing with wfmix ' 482 ! write(std_out,*)' Set al(1,ind_biorthog-3)=one, for ind_biorthog=',ind_biorthog 483 ! write(std_out,*)' This will feed scf_history for set ind_biorthog-3+wfmixalg=',ind_biorthog-3+wfmixalg 484 ! al(:,:)=zero 485 ! al(1,ind_biorthog-3)=one 486 ! call flush(std_out) 487 !ENDDEBUG 488 489 ! LOOP OVER SPINS 490 do isppol=1,dtset%nsppol 491 492 ! BIG FAT k POINT LOOP 493 do ikpt=1,dtset%nkpt 494 495 ! Select k point to be treated by this proc 496 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 497 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) cycle 498 499 istwf_k=dtset%istwfk(ikpt) 500 npw_k=npwarr(ikpt) 501 502 if(istep>2)then 503 ! Make the appropriate linear combination (from the extrapolated wfs) 504 cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=zero 505 do iset2=1,nset2 506 cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(1,icg+1:icg+my_nspinor*npw_k*nband_k)& 507 & +al(1,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)& 508 & -al(2,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg) 509 cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)=cg(2,icg+1:icg+my_nspinor*npw_k*nband_k)& 510 & +al(1,iset2)*scf_history_wf%cg(2,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg)& 511 & +al(2,iset2)*scf_history_wf%cg(1,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,iset2+wfmixalg) 512 end do 513 else ! One needs a simple copy from the extrapolated wavefunctions 514 cg(:,icg+1:icg+my_nspinor*npw_k*nband_k)=scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,1+wfmixalg) 515 end if 516 ! Note the storage in cprj_k. By the way, a simple copy might also be used in case istep=2. 517 if(usepaw==1) then 518 do ibdsp=1,my_nspinor*nbdmix 519 call pawcprj_lincom(al,scf_history_wf%cprj(:,ibdsp,1+wfmixalg:nset2+wfmixalg),cprj_k(:,ibdsp:ibdsp),nset2) 520 end do 521 end if 522 523 ! Store the newly extrapolated wavefunctions for this k point, still bi-orthonormalized, in scf_history_wf 524 scf_history_wf%cg(:,icg_hist+1:icg_hist+my_nspinor*npw_k*nbdmix,ind_newwf)=cg(:,icg+1:icg+my_nspinor*npw_k*nbdmix) 525 if(usepaw==1) then 526 call pawcprj_put(atindx1,cprj_k,scf_history_wf%cprj(:,:,ind_newwf),dtset%natom,1,ibg_hist,ikpt,1,isppol,& 527 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 528 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 529 end if 530 531 ! Back to usual orthonormalization for the cg and cprj_k 532 call cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf_k,mcg,my_nspinor*nband_k,dtset%mkmem,& 533 & mpi_enreg,dtset%natom,nattyp,nbdmix,npw_k,my_nspinor,dtset%nsppol,ntypat,pawtab,usepaw) 534 535 ! Need to transfer cprj_k to cprj 536 if(usepaw==1) then 537 call pawcprj_put(atindx1,cprj_k,cprj,dtset%natom,1,ibg_hist,ikpt,1,isppol,& 538 & nbdmix,dtset%mkmem,dtset%natom,nbdmix,nbdmix,dimcprj,my_nspinor,dtset%nsppol,0,& 539 & mpicomm=mpi_enreg%comm_kpt,mpi_comm_band=spaceComm_band,proc_distrb=mpi_enreg%proc_distrb) 540 end if 541 542 ibg=ibg+my_nspinor*nband_k 543 ibg_hist=ibg_hist+my_nspinor*nbdmix 544 icg=icg+my_nspinor*nband_k*npw_k 545 icg_hist=icg_hist+my_nspinor*nbdmix*npw_k 546 547 ! End big k point loop 548 end do 549 ! End loop over spins 550 end do 551 552 if(istep>2)then 553 ABI_DEALLOCATE(coeffs) 554 end if 555 ABI_DEALLOCATE(al) 556 557 end if ! wfmixalg>2 and istep>1 558 559 !DEBUG 560 ! write(std_out,*)' wf_mixing : exit ' 561 ! write(std_out,*)' scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2)=',& 562 !& scf_history_wf%dotprod_sumdiag_cgcprj_ij(:,2,2) 563 ! write(std_out,*)' cg(1:2,1:2)=',cg(1:2,1:2) 564 ! write(std_out,*)' scf_history_wf%cg(1:2,1:2,1)=',scf_history_wf%cg(1:2,1:2,1) 565 ! ABI_DEALLOCATE(cg_ref) 566 ! ABI_DATATYPE_DEALLOCATE(cprj_ref) 567 !ENDDEBUG 568 569 if(usepaw==1) then 570 call pawcprj_free(cprj_k) 571 call pawcprj_free(cprj_kh) 572 end if 573 ABI_DATATYPE_DEALLOCATE(cprj_k) 574 ABI_DATATYPE_DEALLOCATE(cprj_kh) 575 ABI_DEALLOCATE(dimcprj) 576 ABI_DEALLOCATE(mmn) 577 ABI_DEALLOCATE(smn) 578 if(wfmixalg>2)then 579 ABI_DEALLOCATE(dotprod_res_k) 580 ABI_DEALLOCATE(dotprod_res) 581 ABI_DEALLOCATE(res_mn) 582 end if 583 end subroutine wf_mixing