TABLE OF CONTENTS
- ABINIT/m_shirley
- m_shirley/ovlp_diago_and_prune
- m_shirley/ovlp_diff
- m_shirley/ovlp_free
- m_shirley/ovlp_init
- m_shirley/ovlp_t
- m_shirley/shirley_bands
- m_shirley/shirley_hks
- m_shirley/shirley_interp
- m_shirley/shirley_window
- m_shirley/shparams_t
- m_shirley/shres_free_0D
- m_shirley/shres_free_2D
- m_shirley/shres_init
- m_shirley/shres_t
- m_shirley/wfd_bloch_to_shirley
- m_shirley/wfd_shirley_to_eh
ABINIT/m_shirley [ Modules ]
NAME
m_shirley
FUNCTION
Procedures and objects for the interpolation of KS eigenvalues and eigenvectors with Shirley's method.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (MG) 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 .
PARENTS
CHILDREN
SOURCE
23 #if defined HAVE_CONFIG_H 24 #include "config.h" 25 #endif 26 27 #include "abi_common.h" 28 29 MODULE m_shirley 30 31 use defs_basis 32 use defs_datatypes 33 use defs_abitypes 34 use m_errors 35 use m_xmpi 36 use m_profiling_abi 37 use m_fft 38 use m_nctk 39 #ifdef HAVE_NETCDF 40 use netcdf 41 #endif 42 43 use m_fstrings, only : firstchar 44 use m_time, only : cwtime 45 use m_io_tools, only : flush_unit 46 use m_numeric_tools, only : print_arr, hermitianize, imax_loc, stats_t, stats_eval, wrap2_pmhalf 47 use m_blas, only : xdotc, xgemm, blas_cholesky_ortho 48 use m_abilasi, only : xheev, xhegv, xheevx, xhegvx 49 use m_fft_mesh, only : fft_check_rotrans, setmesh 50 use m_geometry, only : normv, vdotw 51 use m_fftcore, only : get_kg 52 use m_crystal, only : crystal_t 53 use m_crystal_io, only : crystal_ncwrite 54 use m_bz_mesh, only : kmesh_t, get_BZ_item, make_mesh, kmesh_free, kpath_t 55 use m_ebands, only : pack_eneocc, ebands_init, ebands_ncwrite, ebands_free 56 use m_hamiltonian, only : ddiago_ctl_type, init_ddiago_ctl, destroy_hamiltonian, init_hamiltonian, & 57 & load_spin_hamiltonian,load_k_hamiltonian, gs_hamiltonian_type 58 use m_wfd, only : wfd_get_ur, wfd_t, wfd_barrier, wfd_get_cprj, wfd_barrier, wfd_change_ngfft, wfd_print,& 59 & wfd_sym_ur, wfd_init, wfd_push_ug, wfd_reset_ur_cprj, wfd_ihave_ur, wfd_free,& 60 & wfd_test_ortho, kdata_t, kdata_init, kdata_free, wfd_set_mpicomm, wfd_itreat_spin, wfd_xdotc 61 use m_pawang, only : pawang_type 62 use m_pawrad, only : pawrad_type 63 use m_pawtab, only : pawtab_type 64 use m_paw_ij, only : paw_ij_type 65 use m_pawrhoij, only : pawrhoij_type 66 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, paw_overlap 67 use m_pawpwij, only : pawpwff_t, pawpwff_init, pawpwff_free, pawpwij_t, pawpwij_init,& 68 & pawpwij_free, paw_rho_tw_g 69 use m_pawfgr, only : pawfgr_type 70 71 implicit none 72 73 private 74 75 public :: wfd_bloch_to_shirley ! Produce new wavefunction descriptor with the optical basis set from a set of KS states, 76 public :: shirley_interp ! Interpolate the KS Hamiltonian in the Shirley basis set. 77 public :: shirley_hks ! Compute <i|H_ks|j> in the Shirley basis set for a single k-point and spin 78 public :: shirley_window ! Find an optimal basis set for the interpolation of the KS states in a energy window, 79 public :: shirley_bands ! Compute interpolated bands (no eigenvectors) 80 !public :: wfd_shirley_to_eh ! Under development
m_shirley/ovlp_diago_and_prune [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
ovlp_diago_and_prune
FUNCTION
Diagonalize the overlap matrix for a single spin, and selected the optimal basis set according the value of sh_coverage.
INPUTS
sh_coverage=Coverage factor.
OUTPUT
sh_size=Size of the optimal basis set base=Index of the first eigenvector that is excluded from the optimal basis set.
SIDE EFFECTS
O <type(ovlp_t)>= Compute the eigenvalues and store the results in O%eigene
PARENTS
m_shirley
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
1117 subroutine ovlp_diago_and_prune(O,sh_coverage,sh_size,base) 1118 1119 1120 !This section has been created automatically by the script Abilint (TD). 1121 !Do not modify the following lines by hand. 1122 #undef ABI_FUNC 1123 #define ABI_FUNC 'ovlp_diago_and_prune' 1124 use interfaces_14_hidewrite 1125 !End of the abilint section 1126 1127 implicit none 1128 1129 !Arguments ------------------------------------ 1130 !scalars 1131 integer,intent(out) :: sh_size,base 1132 real(dp),intent(in) :: sh_coverage 1133 type(ovlp_t),intent(inout) :: O 1134 1135 !Local variables ------------------------------ 1136 !scalars 1137 integer :: ii 1138 real(dp) :: sum_eige,trace,cpu,wall,gflops 1139 character(len=500) :: msg 1140 1141 !************************************************************************ 1142 1143 DBG_ENTER("COLL") 1144 1145 ! @ovlp_t 1146 call cwtime(cpu,wall,gflops,"start") 1147 ! 1148 ! Diagonalize the overlap matrix. 1149 ABI_MALLOC(O%eigene,(O%size)) 1150 1151 call xheev("Vectors","Upper",O%size,O%mat,O%eigene) 1152 O%mat_type = TYPE_EIGVEC 1153 1154 trace = SUM(O%eigene) 1155 write(msg,'(3(a,f8.2))')" Trace: ",trace," Min eig: ",MINVAL(O%eigene)," Max eig: ",MAXVAL(O%eigene) 1156 call wrtout(std_out,msg,"COLL") 1157 1158 ! Select the optimal subspace. 1159 sum_eige=zero; ii=O%size 1160 do while (sum_eige < trace * sh_coverage .and. ii/=1) 1161 sum_eige = sum_eige + O%eigene(ii) 1162 ii = ii-1 1163 end do 1164 base=ii; sh_size=O%size - base 1165 1166 ! Include degenerate states (if any) in order to have a symmetric basis set. 1167 write(std_out,'(a, i0)')"first base: ",base 1168 do ii=base,1,-1 1169 if (abs(O%eigene(ii)- O%eigene(base+1)) > tol8) exit 1170 end do 1171 base = ii; sh_size=O%size - base 1172 write(std_out,'(a, i0)')"new base: ",base 1173 1174 call cwtime(cpu,wall,gflops,"stop") 1175 write(std_out,*)"SHTIME: Ovlp_diago, cpu_time: ",cpu,", wall_time: ",wall 1176 1177 write(msg,'(2(a,i0),a)')" The Shirley optimal basis set contains ",sh_size,"/",O%size," elements." 1178 call wrtout(std_out,msg,"COLL") 1179 1180 DBG_EXIT("COLL") 1181 1182 end subroutine ovlp_diago_and_prune
m_shirley/ovlp_diff [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
ovlp_diff
FUNCTION
Debugging tool used to compare the overlap matrices calculated with and without symmetries.
INPUTS
PARENTS
m_shirley
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
1207 subroutine ovlp_diff(O1,O2,Kmesh,tol,header,unit,mode_paral,prtvol) 1208 1209 1210 !This section has been created automatically by the script Abilint (TD). 1211 !Do not modify the following lines by hand. 1212 #undef ABI_FUNC 1213 #define ABI_FUNC 'ovlp_diff' 1214 use interfaces_14_hidewrite 1215 !End of the abilint section 1216 1217 implicit none 1218 1219 !Arguments ------------------------------------ 1220 !scalars 1221 integer,optional,intent(in) :: unit,prtvol 1222 real(dp),intent(in) :: tol 1223 character(len=4),optional,intent(in) :: mode_paral 1224 character(len=*),optional,intent(in) :: header 1225 type(ovlp_t),intent(in) :: O1,O2 1226 type(kmesh_t),intent(in) :: Kmesh 1227 1228 !Local variables ------------------------------ 1229 !scalars 1230 integer :: row,col,my_unt,my_prtvol 1231 integer :: ik1_bz,ik2_bz,ik1_ibz,ik2_ibz,band1 1232 integer :: k1_sym,k2_sym,k1_tim,k2_tim 1233 complex(dpc) :: k1_eimkt,k2_eimkt 1234 character(len=4) :: my_mode 1235 character(len=500) :: msg 1236 !arrays 1237 integer :: k1_umklp(3),k2_umklp(3) 1238 real(dp) :: kk1(3),kk2(3) 1239 1240 !************************************************************************ 1241 1242 MSG_ERROR("Not tested") 1243 1244 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 1245 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 1246 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 1247 1248 msg=' ==== Diff O1-O2 ==== ' 1249 if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== ' 1250 call wrtout(my_unt,msg,my_mode) 1251 1252 ABI_CHECK(O1%size==O1%size,"Different sizes") 1253 ! 1254 ! Compare the overlap matrices 1255 do col=1,O2%size 1256 do row=1,O1%size 1257 ! 1258 if (ABS(O1%mat(row,col)-O2%mat(row,col))>tol) then 1259 ik1_bz = O1%idx2bk(1,row) 1260 band1 = O1%idx2bk(2,row) 1261 call get_BZ_item(Kmesh,ik1_bz,kk1,ik1_ibz,k1_sym,k1_tim,k1_eimkt,k1_umklp) 1262 ! 1263 ik2_bz = O1%idx2bk(1,col) 1264 band1 = O1%idx2bk(2,col) 1265 call get_BZ_item(Kmesh,ik2_bz,kk2,ik2_ibz,k2_sym,k2_tim,k2_eimkt,k2_umklp) 1266 1267 write(my_unt,'(2i3,2(2x,a,3f7.3),4f8.4)')& 1268 & row,col," bz1 ",Kmesh%bz(:,ik1_bz)," bz2 ",Kmesh%bz(:,ik2_bz),O1%mat(row,col),O2%mat(row,col) 1269 write(my_unt,'(a,i3,3i3,2f4.1,3i3)')"k2 ",ik2_bz,ik2_ibz,k2_sym,k2_tim,k2_eimkt,k2_umklp 1270 write(my_unt,'(a,i3,3i3,2f4.1,3i3)')"k1 ",ik1_bz,ik1_ibz,k1_sym,k1_tim,k1_eimkt,k1_umklp 1271 1272 end if 1273 end do 1274 end do 1275 1276 end subroutine ovlp_diff
m_shirley/ovlp_free [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
ovlp_free
FUNCTION
Free the memory.
PARENTS
m_shirley
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
415 subroutine ovlp_free(Ovlp) 416 417 418 !This section has been created automatically by the script Abilint (TD). 419 !Do not modify the following lines by hand. 420 #undef ABI_FUNC 421 #define ABI_FUNC 'ovlp_free' 422 !End of the abilint section 423 424 implicit none 425 426 !Arguments ------------------------------------ 427 !scalars 428 type(ovlp_t),intent(inout) :: Ovlp 429 430 ! ************************************************************************* 431 432 DBG_ENTER("COLL") 433 434 !@ovlp_t 435 ! integer 436 if (allocated(Ovlp%bk2idx)) then 437 ABI_FREE(Ovlp%bk2idx) 438 end if 439 if (allocated(Ovlp%idx2bk)) then 440 ABI_FREE(Ovlp%idx2bk) 441 end if 442 ! 443 ! real 444 if (allocated(Ovlp%eigene)) then 445 ABI_FREE(Ovlp%eigene) 446 end if 447 ! 448 ! complex 449 if (allocated(Ovlp%mat)) then 450 ABI_FREE(Ovlp%mat) 451 end if 452 453 DBG_EXIT("COLL") 454 455 end subroutine ovlp_free
m_shirley/ovlp_init [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
ovlp_init
FUNCTION
Calculates the upper triangle of the overlap matrix <u_i|u_j> for a given spin and for all the possible combinations (k,b|k',b') with k and k' in the full Brillouin zone. The u functions are the periodic part of the Bloch wavefunctions hence there is no selection rule in k-space for the matrix elements when k != k' The diagonal matrix elements equals one provided that the input wavefunctions are correctly normalized.
INPUTS
ewin(2,Wfd%nsppol)=Window energies for the two spins. ov_ngfft(18)=FFT mesh used for the computation of the overlap in real space. use_sym=logical flag defining whether symmetries have to be used for reconstructing the overlap matrix. Wfd<wfd_t>=Datatype gathering info on the wavefunctions. Cryst<crystal_t>= data type gathering info on symmetries and unit cell Kmesh<kmesh_t>=Datatype describing the k-point sampling used for the wavefunctions. Bands<ebands_t>=Band structure energies. Psps <type(pseudopotential_type)>=variables related to pseudopotentials Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data Pawang<pawang_type> angular mesh discretization and related data. Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data.
OUTPUT
O(Wfd%nsppol) <type(ovlp_t)>=Object with the overlap matrix elements.
PARENTS
m_shirley
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
498 subroutine ovlp_init(O,ewin,use_sym,ov_ngfft,Wfd,Cryst,Kmesh,Bands,Psps,Pawtab,Pawang,Pawrad) 499 500 501 !This section has been created automatically by the script Abilint (TD). 502 !Do not modify the following lines by hand. 503 #undef ABI_FUNC 504 #define ABI_FUNC 'ovlp_init' 505 use interfaces_14_hidewrite 506 use interfaces_41_geometry 507 use interfaces_65_paw 508 !End of the abilint section 509 510 implicit none 511 512 !Arguments ------------------------------------ 513 !scalars 514 logical,intent(in) :: use_sym 515 type(crystal_t),intent(in) :: Cryst 516 type(Pawang_type),intent(in) :: Pawang 517 type(Pseudopotential_type),intent(in) :: Psps 518 type(kmesh_t),intent(in) :: Kmesh 519 type(ebands_t),intent(in) :: Bands 520 type(wfd_t),target,intent(inout) :: Wfd 521 type(ovlp_t),intent(out) :: O(Wfd%nsppol) 522 !arrays 523 integer,intent(in) :: ov_ngfft(18) 524 real(dp),intent(in) :: ewin(2,Wfd%nsppol) 525 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 526 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*Wfd%usepaw) 527 528 !Local variables ------------------------------ 529 !scalars 530 integer,parameter :: dim_rtwg1=1,option_lob0=0 531 integer :: ik1_bz,ik2_bz,ik1_ibz,ik2_ibz,ierr,rhoxsp_method,nqpt,iqpt,rk 532 integer :: band1,band2,nband_k,nband_k2,nband_k1,npw_k1,npw_k2,spin 533 integer :: with_sym,without_sym 534 integer(i8b) :: row,col,ovlp_size,os_size 535 !MPI 536 integer(i8b) :: tot_nels,oidx,my_hsize,bsize_my_block 537 integer(i8b) :: my_cols(2),my_rows(2) 538 integer(i8b),allocatable :: t_start(:),t_stop(:),hsize_of(:) 539 !ENDMPI 540 integer :: k1_sym,k2_sym,k1_tim,k2_tim,inv_k2_sym,inv_k1_sym 541 integer :: band1_stop,nspinor,lk,rak,iq_found,comm_spin,nproc_spin,rank_spin 542 integer :: nsppol,mband,nkpt !,npw_k,istwf_k,!,ii,itypat,klmn,klmn_size 543 real(dp),parameter :: tnons_tol=tol8,tol_kdiff=tol6 544 real(dp) :: fft_fact,ksq,qsq_max,ene_bks,cpu,wall,gflops 545 complex(dpc) :: blk_ovlp,covlp,k1_eimkt,k2_eimkt,tnons_fact 546 complex(gwpc) :: paw_ovlp(1) 547 logical :: k1_isirred,k2_isirred,can_use_sym,take_cnj 548 character(len=500) :: msg 549 !arrays 550 integer,parameter :: g_gamma(3)=(/0,0,0/) 551 integer :: k1_umklp(3),k2_umklp(3) 552 integer,allocatable :: tmp_idx2bk(:,:) 553 integer,allocatable :: toinv(:,:),multable(:,:,:),nq_spl(:),k1mk2q(:,:) 554 real(dp) :: kpoint(3),fft_err(3,Cryst%nsym) !,onsite(2) 555 real(dp) :: kk1(3),kk2(3),k1mk2(3) !,ovlp_paw(2) !,r1_tau3(3) 556 real(dp),allocatable :: qpts(:,:),qmax(:) 557 complex(gwpc),allocatable :: ur1(:),ur2(:) 558 complex(gwpc),pointer :: ug2(:) !ug1(:) 559 complex(dpc),allocatable :: ovlp_ikfk(:,:,:,:) 560 logical :: k_needs_tr(2) 561 type(pawcprj_type),pointer :: Cp_k1(:,:),Cp_k2(:,:) 562 type(pawpwij_t),allocatable :: Pwij(:,:) 563 type(pawpwff_t),allocatable :: Paw_pwff(:) 564 565 !************************************************************************ 566 567 DBG_ENTER("COLL") 568 569 ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded") 570 571 ABI_UNUSED(Pawrad(1)%mesh_size) 572 573 call cwtime(cpu,wall,gflops,"start") 574 575 ! Change FFT mesh. 576 call wfd_change_ngfft(Wfd,Cryst,Psps,ov_ngfft) 577 578 ! Check if the mesh preserves the symmetries 579 if (.not.fft_check_rotrans(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,fft_err)) then 580 write(msg,'(a,3(i0,1x),a)')" Real space FFT mesh ",Wfd%ngfft(1:3)," is not symmetric. Cannot symmetrize in real space" 581 MSG_ERROR(msg) 582 end if 583 584 nsppol = Wfd%nsppol; nspinor = Wfd%nspinor 585 mband = Wfd%mband 586 587 ! Use the total number of k-points in the BZ. 588 nkpt = Kmesh%nbz 589 ! 590 ! Calculate the overlap matrix 591 ! ====================================================== 592 ! ==== <u_{ik b1}| u_{fk b2}>, ik in IBZ, fk in BZ ==== 593 ! ====================================================== 594 595 ABI_MALLOC(multable,(4,Cryst%nsym,Cryst%nsym)) 596 ABI_MALLOC(toinv,(4,Cryst%nsym)) 597 call sg_multable(Cryst%nsym,Cryst%symafm,Cryst%symrel,Cryst%tnons,tnons_tol,ierr,multable=multable,toinv=toinv) 598 ABI_CHECK(ierr==0,"Group error, cannot continue") 599 ! 600 ! For PAW we have to recalculate the projections in the IBZ setting k=0 in the exponential. 601 ! TODO: For the time being, these terms are calculated in the BZ, symmetrization will be added afterwards. 602 if (Wfd%usepaw==1) then 603 ! 604 ! <u_k2|u_k1> = <phi_k2|e^{-i(k1-k2).r}|psi_k1> 605 ! Find the list of q-points: q=k1-k2 and create table (k1, k2) --> q = k1-k2 606 ABI_MALLOC(qpts,(3,nkpt**2)) 607 ABI_MALLOC(k1mk2q,(nkpt,nkpt)) 608 609 nqpt=0; qsq_max=zero 610 do ik1_bz=1,nkpt 611 kk1 = Kmesh%bz(:,ik1_bz) 612 do ik2_bz=1,nkpt 613 kk2 = Kmesh%bz(:,ik2_bz) 614 615 k1mk2 = kk1 - kk2 616 !k1mk2 = kk2 - kk1 617 ksq = normv(k1mk2,Cryst%gmet,"G") 618 qsq_max = MAX(ksq,qsq_max) 619 620 iq_found=0 621 do iqpt=1,nqpt 622 if (ALL( ABS(k1mk2 - qpts(:,iqpt)) < tol_kdiff) ) then 623 iq_found=iqpt 624 EXIT 625 end if 626 end do 627 628 if (iq_found==0) then 629 ! Add it to the list. 630 nqpt = nqpt+1 631 !ABI_CHECK(nqpt <= size(qpts,dim=2),"too many qpts!") 632 qpts(:,nqpt) = k1mk2 633 iq_found = nqpt 634 end if 635 636 k1mk2q(ik1_bz,ik2_bz) = iq_found 637 end do 638 end do 639 640 ! TODO: This one should be correct but it's too small, why? 641 !qsq_max = qsq_max/(two_pi**2) 642 qsq_max = qsq_max/(two*pi**2) 643 ! 644 ! Set up q-grid for form factors, make qmax 20% larger than the largest expected. 645 ABI_MALLOC(nq_spl,(Cryst%ntypat)) 646 ABI_MALLOC(qmax,(Cryst%ntypat)) 647 nq_spl = Psps%mqgrid_ff 648 qmax = SQRT(qsq_max)*1.2d0 649 write(std_out,*)"Using nqpt:",nqpt," nq_spl ",nq_spl," qmax_type=",qmax 650 651 ! TODO: pass input variable to control the computation. 652 rhoxsp_method=1 ! Arnaud-Alouani 653 !rhoxsp_method=2 ! Shiskin-Kresse 654 655 ABI_DT_MALLOC(Paw_pwff,(Psps%ntypat)) 656 call pawpwff_init(Paw_pwff,rhoxsp_method,nq_spl,qmax,Cryst%gmet,Pawrad,Pawtab,Psps) 657 658 ABI_FREE(nq_spl) 659 ABI_FREE(qmax) 660 661 ABI_DT_MALLOC(Pwij,(Psps%ntypat,nqpt)) 662 do iqpt=1,nqpt 663 call pawpwij_init(Pwij(:,iqpt),1,qpts(:,iqpt),g_gamma,Cryst%rprimd,Psps,Pawtab,Paw_pwff) 664 665 ! DEBUG 666 !if (iqpt==1) then 667 ! write(std_out,*)qpts(:,iqpt) 668 ! do itypat=1,Cryst%ntypat 669 ! klmn_size = Pwij(itypat,iqpt)%lmn_size*(Pwij(itypat,iqpt)%lmn_size+1)/2 670 ! do klmn=1,klmn_size 671 ! !write(std_out,*)" mqpgij: ",Pwij(itypat,iqpt)%mqpgij(:,1,klmn),pawtab(itypat)%sij(klmn) 672 ! write(std_out,*)" diff mqpg0ij: ",Pwij(itypat,iqpt)%mqpgij(1,1,klmn)-pawtab(itypat)%sij(klmn) 673 ! Pwij(itypat,iqpt)%mqpgij(1,1,klmn) = pawtab(itypat)%sij(klmn) 674 ! Pwij(itypat,iqpt)%mqpgij(2,1,klmn) = zero 675 ! end do 676 ! end do 677 !end if 678 end do 679 680 call pawpwff_free(Paw_pwff) 681 ABI_DT_FREE(Paw_pwff) 682 ABI_FREE(qpts) 683 684 ! Allocate cprj datastructures. 685 ABI_DT_MALLOC(Cp_k1 ,(Cryst%natom,nspinor)) 686 call pawcprj_alloc(Cp_k1,0,Wfd%nlmn_atm) 687 688 ABI_DT_MALLOC(Cp_k2 ,(Cryst%natom,nspinor)) 689 call pawcprj_alloc(Cp_k2,0,Wfd%nlmn_atm) 690 end if ! PAW 691 ! 692 ! ============================================================== 693 ! ==== Reconstruct full <u_{kb}| u_{k'b'}> matrix in the BZ ==== 694 ! ============================================================== 695 ! 1) Symmetrization is done in real space. Easier, especially when k-centered G-sphere are used. 696 ! u(r,b,kbz)=e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t),b,kibz) 697 ! =e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t),b,kibz) for time-reversal 698 ! 2) Matrix is Hermitian. 699 ! 3) <u_{Sk b}| u_{Sk b'}> are obtained from the previously calculated <u_{kb}| u_{kb'}> table. 700 ! 701 ! {A,a} {B,b} = {AB, a+Ab} 702 ! 703 fft_fact = one/Wfd%nfft 704 ABI_MALLOC(ur1,(Wfd%nfft*nspinor)) 705 ABI_MALLOC(ur2,(Wfd%nfft*nspinor)) 706 707 do spin=1,nsppol 708 if (.not. wfd_itreat_spin(Wfd,spin,comm_spin,rank_spin,nproc_spin)) cycle 709 710 ! Compute table (b,k) --> indx and indx --> (b,k) where (b,k) is 711 ! the set of bands and kpoints included in the calculation of the overlap. 712 O(spin)%mband = mband 713 O(spin)%nkpt = nkpt 714 O(spin)%min_ene = ewin(1,spin) 715 O(spin)%max_ene = ewin(2,spin) 716 ! 717 ! The size of the overlap matrix and useful tables. 718 ABI_CALLOC(O(spin)%bk2idx,(mband,nkpt)) 719 ABI_CALLOC(tmp_idx2bk,(2,mband*nkpt)) 720 721 ovlp_size=0 722 do ik2_bz=1,nkpt 723 ik2_ibz = Kmesh%tab(ik2_bz) 724 nband_k2 = Wfd%nband(ik2_ibz,spin) 725 do band2=1,nband_k2 726 ene_bks = Bands%eig(band2,ik2_ibz,spin) 727 ! Is it inside the window? 728 if (ene_bks >= O(spin)%min_ene .and. ene_bks <= O(spin)%max_ene) then 729 ovlp_size = ovlp_size + 1 730 O(spin)%bk2idx(band2,ik2_bz) = ovlp_size 731 tmp_idx2bk(1,ovlp_size) = band2 732 tmp_idx2bk(2,ovlp_size) = ik2_bz 733 end if 734 end do 735 end do 736 737 O(spin)%size = ovlp_size 738 if (O(spin)%size /= int(ovlp_size,kind=i4b)) then 739 ! Stop here since we need Lapack i64! 740 MSG_ERROR("Size of overlap matrix too large for a default integer!") 741 end if 742 743 ABI_MALLOC(O(spin)%idx2bk,(2,ovlp_size)) 744 O(spin)%idx2bk = tmp_idx2bk(:,1:ovlp_size) 745 ABI_FREE(tmp_idx2bk) 746 ! 747 ! Allocate overlap matrix. Could use packed matrix to save memory, but Lapack call is slower. 748 write(msg,'(a,i0,a,f12.1,a,i0)')" Memory for the overlap matrix(spin=",spin,"): ",& 749 & two*dpc*ovlp_size**2*b2Mb," Mb; Matrix size = ",ovlp_size 750 call wrtout(std_out,msg,"COLL") 751 752 ABI_MALLOC(O(spin)%mat, (ovlp_size,ovlp_size)) 753 O(spin)%mat = czero 754 O(spin)%mat_type = TYPE_OVERLAP 755 756 ! Distribute the calculation of the matrix elements among the nodes. 757 ! * tstart and t_stop give the initial and final transition index treated by each node. 758 ! * my_hsize is the number of transitions treated by this processor 759 ! * my_cols(1:2) gives the initial and final column treated by this node. 760 ! 761 ! Distribution of the matrix elements among the nodes. 762 os_size = O(spin)%size; tot_nels = os_size * (os_size + 1_i8b) / 2 763 ABI_MALLOC(t_start,(0:nproc_spin-1)) 764 ABI_MALLOC(t_stop,(0:nproc_spin-1)) 765 !t_start = 1; t_stop = tot_nels 766 !my_rows = [1_i8b,os_size]; my_cols = [1_i8b,os_size] 767 768 call xmpi_split_work2_i8b(tot_nels,nproc_spin,t_start,t_stop,msg,ierr) 769 if (ierr==2) MSG_WARNING(msg) 770 771 ABI_MALLOC(hsize_of,(0:nproc_spin-1)) 772 hsize_of=0 773 do rak=0,nproc_spin-1 774 if (t_stop(rak)>=t_start(rak)) hsize_of(rak) = t_stop(rak)-t_start(rak)+1 775 !write(std_out,*)"tot_nels: ",tot_nels,hsize_of(rak) 776 end do 777 778 my_hsize = hsize_of(rank_spin) 779 if (my_hsize<=0) then 780 write(msg,'(a,i0)')"Wrong number of transitions: my_hsize= ",my_hsize 781 MSG_ERROR(msg) 782 end if 783 if (my_hsize /= int(my_hsize,kind=i4b)) then 784 write(msg,'(a,i0)')"Size of local block too large for a default integer, Increase the number of CPUs: my_hsize= ",my_hsize 785 MSG_ERROR(msg) 786 end if 787 ABI_FREE(hsize_of) 788 789 my_cols=0 790 do col=1,os_size 791 do row=1,col 792 oidx = row + col*(col-1_i8b)/2 793 if (oidx==t_start(rank_spin)) then 794 my_rows(1) = row 795 my_cols(1) = col 796 end if 797 if (oidx==t_stop(rank_spin)) then 798 my_rows(2) = row 799 my_cols(2) = col 800 end if 801 end do 802 end do 803 ! 804 ! * Announce the treatment of submatrix treated by each node. 805 bsize_my_block = 2*dpc*my_hsize 806 write(msg,'(4(a,i0))')& 807 & ' Treating ',my_hsize,'/',tot_nels,' matrix elements, from column ',my_cols(1),' up to column ',my_cols(2) 808 call wrtout(std_out,msg,'PERS') 809 810 ! Computation of the overlap matrix. 811 SELECT CASE (use_sym) 812 813 CASE (.FALSE.) 814 call wrtout(std_out," Version without symmetries","COLL") 815 816 do ik2_bz=1,nkpt 817 call get_BZ_item(Kmesh,ik2_bz,kk2,ik2_ibz,k2_sym,k2_tim,k2_eimkt,k2_umklp,k2_isirred) 818 nband_k2 = Wfd%nband(ik2_ibz,spin) 819 npw_k2 = Wfd%npwarr(ik2_ibz) 820 821 do band2=1,nband_k2 822 col = O(spin)%bk2idx(band2,ik2_bz) 823 824 ! Check is this col is treated by me. 825 if (my_cols(2)<col .or. my_cols(1)>col) cycle 826 827 ug2 => Wfd%Wave(band2,ik2_ibz,spin)%ug 828 ! Get the symmetrized wavefunction in real space. 829 call wfd_sym_ur(Wfd,Cryst,Kmesh,band2,ik2_bz,spin,ur2) 830 831 if (Wfd%usepaw==1) then 832 call wfd_get_cprj(Wfd,band2,ik2_ibz,spin,Cryst,Cp_k2,sorted=.FALSE.) 833 if (.not.k2_isirred) then 834 call paw_symcprj(ik2_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cp_k2) 835 endif 836 end if 837 838 do ik1_bz=1,ik2_bz 839 call get_BZ_item(Kmesh,ik1_bz,kk1,ik1_ibz,k1_sym,k1_tim,k1_eimkt,k1_umklp,k1_isirred) 840 npw_k1 = Wfd%npwarr(ik1_ibz) 841 nband_k1 = Wfd%nband(ik1_ibz,spin) 842 band1_stop = nband_k1; if (ik1_bz==ik2_bz) band1_stop = band2 843 if (Wfd%usepaw==1) iqpt = k1mk2q(ik1_bz,ik2_bz) 844 845 do band1=1,band1_stop 846 row = O(spin)%bk2idx(band1,ik1_bz) 847 848 ! Check if this element is treated by me. 849 oidx = row + col*(col-1)/2 850 if (oidx<t_start(rank_spin) .or. oidx>t_stop(rank_spin)) cycle 851 852 if (ik1_bz==ik2_bz) then 853 ! Save overlap here and skip the computation. Assume input u(g) are normalized! 854 blk_ovlp = czero; if (band1==band2) blk_ovlp = cone 855 O(spin)%mat(row,col) = blk_ovlp 856 !blk_ovlp = wfd_norm2(Wfd,Cryst,Pawtab,band,ik1_ibz,spin) result(norm2) 857 cycle 858 end if 859 860 if (ik2_bz==ik1_bz .and. k1_isirred) then 861 ! Save the overlap here and skip the symmetrization. 862 blk_ovlp = wfd_xdotc(Wfd,Cryst,Pawtab,band1,band2,ik1_ibz,spin) 863 cycle 864 end if 865 866 ! Get the symmetrized wavefunction in real space. 867 ! Here we have the SIGFAULT 868 call wfd_sym_ur(Wfd,Cryst,Kmesh,band1,ik1_bz,spin,ur1) 869 870 ! Scalar product 871 blk_ovlp = xdotc(Wfd%nfft*nspinor,ur1,1,ur2,1)/Wfd%nfft 872 873 if (Wfd%usepaw==1) then 874 ! Add onsite contribution. 875 call wfd_get_cprj(Wfd,band1,ik1_ibz,spin,Cryst,Cp_k1,sorted=.FALSE.) 876 if (.not.k1_isirred) then 877 call paw_symcprj(ik1_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cp_k1) 878 endif 879 880 !if (iqpt==1) then 881 ! onsite = paw_overlap(Cp_k1,Cp_k2,Cryst%typat,Pawtab) 882 ! blk_ovlp = blk_ovlp + DCMPLX(onsite(1),onsite(2)) 883 !else 884 paw_ovlp = czero 885 call paw_rho_tw_g(1,dim_rtwg1,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,g_gamma,& 886 & Cp_k1,Cp_k2,Pwij(:,iqpt),paw_ovlp) 887 888 ! FIXME not addinf the onsite term works just fine! 889 !paw_ovlp = czero 890 blk_ovlp = blk_ovlp + paw_ovlp(1) 891 !end if 892 end if 893 !end if 894 ! 895 ! Save overlap. 896 O(spin)%mat(row,col) = blk_ovlp 897 end do 898 end do 899 end do 900 end do 901 902 CASE (.TRUE.) 903 ! FIXME does not work yet. Problems with umklapp symmorphic operations and symmetry tables somewhere. 904 905 call wrtout(std_out," Version with symmetries","COLL") 906 MSG_WARNING(" Using version with symmetries (still under development") 907 908 ABI_CHECK(Wfd%usepaw==0,"PAW not coded yet") 909 ABI_MALLOC(ovlp_ikfk,(mband,Wfd%nkibz,mband,nkpt)) 910 ovlp_ikfk=+HUGE(one) 911 912 do ik2_bz=1,nkpt 913 call get_BZ_item(Kmesh,ik2_bz,kk2,ik2_ibz,k2_sym,k2_tim,k2_eimkt,k2_umklp) 914 nband_k2 = Wfd%nband(ik2_ibz,spin) 915 916 do band2=1,nband_k2 917 call wfd_sym_ur(Wfd,Cryst,Kmesh,band2,ik2_bz,spin,ur2) 918 919 do ik1_ibz=1,Wfd%nkibz 920 nband_k = Wfd%nband(ik1_ibz,spin) 921 do band1=1,nband_k 922 923 call wfd_get_ur(Wfd,band1,ik1_ibz,spin,ur1) 924 covlp = xdotc(Wfd%nfft*nspinor,ur1,1,ur2,1) * fft_fact 925 926 if (Wfd%usepaw==1) then ! TODO Add onsite term after IBZ-->BZ symmetrization 927 end if 928 929 ovlp_ikfk(band1,ik1_ibz,band2,ik2_bz) = covlp 930 end do 931 end do 932 ! 933 end do 934 end do 935 ! 936 ! FIXME: Equations are not completed. non-symmorphic phases are missing! 937 ! Let <r|k> indicate the periodic part of the Bloch wavefunction that transforms according to: 938 ! u_{Sk} = e^{-iSk.t} u_k(S^{-1}(r-t)) 939 ! S3 = S1^{-1} S2 940 ! 941 ! 1) <S1 k1 | S2 k2> = <k1| S1^{-1} S2 k2> e^{i (S1 k1 - S2 k2).t1} 942 ! 943 ! 2) <T S1 k1 | T S2 k2> = <S1 k1| S2 k2>* 944 945 ! <S1 k1 | T S2 k2> = <k1| S1^{-1} T S2 k2> 946 ! <T S1 k1 | S2 k2> = <k2| S2^{-1} T S1 k1>^* 947 ! <T S1 k1 | T S2 k2> = <k2| S2^{-1} S1 k1> ! Problematic if the BZ mesh is not invariant under inversion. 948 ! e.g. randomly shifted k-meshes. In this case one should use kptopt=3 949 call wrtout(std_out,"Using version with symmetries","COLL") 950 with_sym=0; without_sym=0 951 do ik2_bz=1,nkpt 952 953 call get_BZ_item(Kmesh,ik2_bz,kk2,ik2_ibz,k2_sym,k2_tim,k2_eimkt,k2_umklp) 954 955 !$ik2_ibz = Kmesh%tab(ik2_bz) 956 nband_k2 = Wfd%nband(ik2_ibz,spin) 957 958 inv_k2_sym = toinv(1,k2_sym) 959 k_needs_tr(2) = (k2_tim==2) 960 961 do ik1_bz=1,ik2_bz !nkpt 962 call get_BZ_item(Kmesh,ik1_bz,kk1,ik1_ibz,k1_sym,k1_tim,k1_eimkt,k1_umklp) 963 964 nband_k1 = Wfd%nband(ik1_ibz,spin) 965 966 inv_k1_sym = toinv(1,k1_sym) 967 968 k_needs_tr(1) = (k1_tim==2) 969 970 !sym3 = multable(1,inv_k1_sym,k2_sym) 971 972 rk= Kmesh%rottb(ik2_bz,1,inv_k1_sym) 973 974 can_use_sym = .TRUE. 975 !can_use_sym = ( ALL(ABS(Cryst%tnons(:,k2_sym))<tol6) .and. ALL(ABS(Cryst%tnons(:,k1_sym))<tol6) ) 976 !can_use_sym = can_use_sym .and. ( ALL(k2_umklp == g_gamma) .and. ALL(k1_umklp == g_gamma)) 977 978 can_use_sym = can_use_sym .and. ( ALL(k2_umklp == g_gamma) .and. ALL(k1_umklp == g_gamma) & 979 & .and. ALL(multable(2:4,inv_k1_sym,k2_sym)==g_gamma) ) 980 981 !can_use_sym = ( can_use_sym .and. & 982 !& ALL( ABS ( -Kmesh%bz(:,rk) + MATMUL(Cryst%symrec(:,:,inv_k1_sym),Kmesh%bz(:,ik2_bz)) ) < tol6 ) ) 983 984 !can_use_sym = ( can_use_sym .and. ALL(ABS(Cryst%tnons(:,k1_sym)) <tol6) .and. ALL(ABS(Cryst%tnons(:,k2_sym)) <tol6) ) 985 986 kpoint = kk1 - kk2 987 !do ii=1,3 ! Wrap in the first BZ thus enforcing traslational invariance. 988 ! call wrap2_pmhalf(kk1(ii)-kk2(ii),kpoint(ii),shift) ! TODO overloaded interface. 989 !end do 990 991 if (ANY (ABS(Cryst%tnons(:,k1_sym)) > tol6) ) then 992 !tnons_fact = EXP(j_dpc*two_pi*DOT_PRODUCT(kk1-kk2,Cryst%tnons(:,k1_sym))) 993 tnons_fact = EXP(j_dpc*two_pi*DOT_PRODUCT(kpoint,Cryst%tnons(:,k1_sym))) 994 else 995 tnons_fact = cone 996 end if 997 998 take_cnj=ALL(k_needs_tr) 999 1000 if (ALL(k_needs_tr) .or. ALL(.not.k_needs_tr) ) then 1001 lk = ik1_ibz 1002 rk= Kmesh%rottb(ik2_bz,1,inv_k1_sym) 1003 !rk= Kmesh%rottbm1(ik2_bz,1,k1_sym) 1004 else 1005 MSG_ERROR("Need TR") 1006 take_cnj=.FALSE. 1007 !if (k_needs_tr(2)) then 1008 !lk = ik1_ibz 1009 !rk=Kmesh%rottb(ik2_bz,1,inv_k1_sym) 1010 !tnons_fact = CONJG(k1_eimkt) * CONJG(k2_eimkt) * EXP(-j_dpc*two_pi*DOT_PRODUCT(kk2,r1_tau3)) 1011 !end if 1012 end if 1013 1014 if (can_use_sym) then 1015 1016 with_sym=with_sym+1 1017 do band2=1,nband_k2 1018 col=O(spin)%bk2idx(band2,ik2_bz) 1019 1020 band1_stop = nband_k1; if (ik1_bz==ik2_bz) band1_stop = band2 1021 do band1=1,band1_stop 1022 blk_ovlp = tnons_fact * ovlp_ikfk(band1,lk,band2,rk) 1023 if (take_cnj) blk_ovlp = DCONJG(blk_ovlp) 1024 row=O(spin)%bk2idx(band1,ik1_bz) 1025 if (col>=row) O(spin)%mat(row,col) = blk_ovlp 1026 end do 1027 end do 1028 1029 else 1030 without_sym=without_sym+1 1031 do band2=1,nband_k2 1032 col=O(spin)%bk2idx(band2,ik2_bz) 1033 1034 call wfd_sym_ur(Wfd,Cryst,Kmesh,band2,ik2_bz,spin,ur2) 1035 1036 band1_stop = nband_k1; if (ik1_bz==ik2_bz) band1_stop = band2 1037 1038 do band1=1,band1_stop 1039 call wfd_sym_ur(Wfd,Cryst,Kmesh,band1,ik1_bz,spin,ur1) 1040 1041 blk_ovlp = xdotc(Wfd%nfft,ur1,1,ur2,1) * fft_fact 1042 row=O(spin)%bk2idx(band1,ik1_bz) 1043 if (col>=row) O(spin)%mat(row,col) = blk_ovlp 1044 end do 1045 end do 1046 end if 1047 1048 end do 1049 end do 1050 1051 write(std_out,*)"with_sym",with_sym," without_sym",without_sym 1052 ABI_FREE(ovlp_ikfk) 1053 END SELECT 1054 1055 ABI_FREE(t_start) 1056 ABI_FREE(t_stop) 1057 ! 1058 ! Collect ovlp_mat on each node inside comm_spin. 1059 call xmpi_sum(O(spin)%mat,comm_spin,ierr) 1060 end do ! spin 1061 1062 ABI_FREE(multable) 1063 ABI_FREE(toinv) 1064 ABI_FREE(ur1) 1065 ABI_FREE(ur2) 1066 1067 ! Deallocate PAW stuff. 1068 if (Wfd%usepaw==1) then 1069 ABI_FREE(k1mk2q) 1070 call pawpwij_free(Pwij) 1071 ABI_DT_FREE(Pwij) 1072 call pawcprj_free(Cp_k1) 1073 ABI_DT_FREE(Cp_k1) 1074 call pawcprj_free(Cp_k2) 1075 ABI_DT_FREE(Cp_k2) 1076 end if 1077 1078 call cwtime(cpu,wall,gflops,"stop") 1079 write(std_out,*)"SHTIME: Ovlp_init, cpu_time: ",cpu,", wall_time: ",wall 1080 1081 DBG_EXIT("COLL") 1082 1083 end subroutine ovlp_init
m_shirley/ovlp_t [ Types ]
[ Top ] [ m_shirley ] [ Types ]
NAME
ovlp_t
FUNCTION
Store the overlap matrix <u_bk|u_{b'k'}> for a given spin.
SOURCE
195 type,public :: ovlp_t 196 197 integer :: size=0 198 ! The size of the overlap matrix. 199 200 integer :: mband=0 201 ! Maximum number of bands. 202 203 integer :: nkpt=0 204 ! Number of k-points. 205 206 integer :: mat_type = TYPE_NONE 207 ! Matrix type (overlap or eigenvectors) 208 209 real(dp) :: min_ene = -HUGE(one) 210 ! Min energy included. 211 212 real(dp) :: max_ene = +HUGE(one) 213 ! Max energy included 214 215 integer(i8b),allocatable :: bk2idx(:,:) 216 ! bk2idx(mband,nkpt) 217 ! Mapping (b,k) --> i 218 ! 0 if the state has been excluded from the energy window. 219 220 integer,allocatable :: idx2bk(:,:) 221 ! idx2bk(2,size)) 222 ! Mapping i --> (b,k) for i=1,size. k is always the in the BZ. 223 224 complex(dpc),allocatable :: mat(:,:) 225 ! mat(size,size) 226 ! Stores O_{ij} = <u_i|u_j> 227 ! The matrix is Hermitian hence one could use a packed matrix to save memory, 228 ! but the Lapack routines for the diagonalization are usually slower. 229 230 real(dp),allocatable :: eigene(:) 231 ! eigene(size) 232 ! The eigenvalues of the overlap operator. 233 234 end type ovlp_t 235 236 public :: ovlp_init ! Creation method 237 public :: ovlp_free ! Deallocate memory 238 public :: ovlp_diago_and_prune ! Diagonalize the overlap matrix and select the important subspace.
m_shirley/shirley_bands [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
shirley_bands
FUNCTION
Helper function for the interpolation of the eigenvalues with the Shirley method INPUT iWfd<wfd_t>=Input wavefunctions. Dtset=Input variables for this run (it will be removed later on) Cryst<crystal_t>= data type gathering info on symmetries and unit cell. iKmesh<kmesh_t>=K-points of the wavefunctions stored in Wfd. iBands<ebands_t>=Input energies. Psps<pseudopotential_type)>=variables related to pseudopotentials. Pawtab(Psps%ntypat)<pawtab_type>=paw tabulated starting data Pawang<pawang_type>=angular mesh discretization and related data: Pawrad<Pawrad_type>=For PAW, RADial mesh discretization and related data Pawfgr<Pawfgr_type> Pawang<pawang_type> angular mesh discretization and related data: Pawrad<Pawrad_type> Pawrhoij Paw_ij(natom)<paw_ij_type>=data structure containing PAW arrays given on (i,j) channels. ngfftc(18)=Information about the coarse 3D FFT. ngfftf(18)=Information about the dense 3D FFT used for vtrial. nfftf=Number of points in the FFT grid in vtrial. Might differ from the FFT mesh used for the wavefunctions. vtrial(nfftf,nspden)= Trial potential (Hartree) intp_mband=Number of bands to intepolate. sh_coverage=Parameter defining the coverage of the optimal basis set: Wsh contains N vectors where N is the first elements for which we have: \sum_i^N e_i >= sh_coverage * M with e_i being the eigenvalues of the overlap matrix ordered in descending order. M indicates the dimension of the overlap (related to the number of k-points in the BZ and the number of bands in the input iWfd descriptor. ewin(2,nsppol)=Energy window for the spin channels. Use ewin(1,:)=smallest_real and ewin(2,:)=greatest_real to disable the use of the energy window min_bsize=Minimum number of states in the optimal basis set. sh_fname=String with the filename of the netcdf file to be produced. No file is created if empty string. kpts(3,nk)=List of k-points for the interpolation. kweights(nk)=K-point weights (stored in oBands)
OUTPUT
oWsh<wfd_t>=New wavefunction descriptor with the optimal basis set. Note that Wfs contains a single (fake) k-point.
PARENTS
wfk_analyze
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
2500 subroutine shirley_bands(iWfd,Dtset,Cryst,iKmesh,iBands,Psps,Pawtab,Pawang,Pawrad,Pawfgr,Pawrhoij,Paw_ij,& 2501 & sh_coverage,ewin,ngfftc,ngfftf,nfftf,vtrial,intp_mband,kpts,kweights,sh_fname,oBands,oWsh) 2502 2503 2504 !This section has been created automatically by the script Abilint (TD). 2505 !Do not modify the following lines by hand. 2506 #undef ABI_FUNC 2507 #define ABI_FUNC 'shirley_bands' 2508 !End of the abilint section 2509 2510 implicit none 2511 2512 !Arguments ------------------------------------ 2513 !scalars 2514 integer,intent(in) :: nfftf,intp_mband 2515 real(dp),intent(in) :: sh_coverage 2516 character(len=fnlen),intent(in) :: sh_fname 2517 type(crystal_t),intent(in) :: Cryst 2518 type(Pawang_type),intent(in) :: Pawang 2519 type(kmesh_t),intent(in) :: iKmesh 2520 type(Pseudopotential_type),intent(in) :: Psps 2521 type(wfd_t),intent(inout) :: iWfd 2522 type(ebands_t),intent(in) :: iBands 2523 type(dataset_type),intent(in) :: Dtset 2524 type(Pawfgr_type),intent(in) :: Pawfgr 2525 type(wfd_t),intent(out) :: oWsh 2526 type(ebands_t),intent(out) :: oBands 2527 !arrays 2528 integer,intent(in) :: ngfftf(18),ngfftc(18) 2529 real(dp),intent(in) :: vtrial(nfftf,iWfd%nspden),kpts(:,:),kweights(:) 2530 real(dp),intent(in) :: ewin(2,iWfd%nsppol) 2531 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*iWfd%usepaw) 2532 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*iWfd%usepaw) 2533 type(paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*iWfd%usepaw) 2534 type(pawrhoij_type),intent(in) :: Pawrhoij(Cryst%natom*iWfd%usepaw) 2535 2536 !Local variables ------------------------------ 2537 !scalars 2538 integer :: min_bsize,npts 2539 !arrays 2540 integer,allocatable :: intp_nband(:,:) 2541 2542 !************************************************************************ 2543 2544 DBG_ENTER("COLL") 2545 2546 ! Compute shirley basis set 2547 min_bsize=intp_mband 2548 call wfd_bloch_to_shirley(iWfd,Cryst,iKmesh,iBands,Psps,Pawtab,Pawang,Pawrad,min_bsize,sh_coverage,ewin,oWsh) 2549 2550 npts = size(kpts,dim=2) 2551 ABI_MALLOC(intp_nband, (npts,iWfd%nsppol)) 2552 intp_nband=intp_mband 2553 2554 ! Energies only 2555 call shirley_interp(oWsh,"N",Dtset,Cryst,Psps,Pawtab,Pawfgr,Pawang,Pawrad,& 2556 & Pawrhoij,Paw_ij,ngfftc,ngfftf,nfftf,vtrial,& 2557 & intp_nband,intp_mband,npts,kpts,sh_fname,oBands,kweights) 2558 2559 ! Free memory 2560 ABI_FREE(intp_nband) 2561 2562 DBG_EXIT("COLL") 2563 2564 end subroutine shirley_bands
m_shirley/shirley_hks [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
shirley_hks
FUNCTION
Compute <i|H_ks|j> in the Shirley basis set for a single k-point and spin This routine is called by a single MPI node who has the entire set of Shirley wavefunctions in memory. MPI parallelization is done in the caller at the level of k-points and spins.
INPUTS
Wsh=Description with Shirley savefunction kpt(3)=K-point in reduced coordinates. Note: The k-point will be wrapped to [-1/2,1/2[ spin=Spin index. Ham_k=Ground state Hamiltonian for this k-point and spin. Cryst<crystal_t>= data type gathering info on symmetries and unit cell Psps <type(pseudopotential_type)>=variables related to pseudopotentials Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data Pawang<pawang_type> angular mesh discretization and related data: Paw_ij(natom)<type(paw_ij_type)>=data structure containing PAW arrays given on (i,j) channels. sh_size=Size of the Shirley basis set. vloc_ij(sh_size,sh_size)=Matrix element of the local part in the Shirley basis set.
OUTPUT
hk_ij(sh_size,sh_size)= <i|H_ks|j> sk_ij(sh_size,sh_size*Psps%usepaw)=PAW matrix.
PARENTS
m_shirley
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
1683 subroutine shirley_hks(Wsh,kpt,spin,Ham_k,Cryst,Psps,Pawtab,Pawang,Paw_ij,sh_size,vloc_ij,hk_ij,sk_ij) 1684 1685 1686 !This section has been created automatically by the script Abilint (TD). 1687 !Do not modify the following lines by hand. 1688 #undef ABI_FUNC 1689 #define ABI_FUNC 'shirley_hks' 1690 use interfaces_56_recipspace 1691 use interfaces_66_nonlocal 1692 use interfaces_69_wfdesc 1693 !End of the abilint section 1694 1695 implicit none 1696 1697 !Arguments ------------------------------------ 1698 !scalars 1699 integer,intent(in) :: spin,sh_size 1700 type(crystal_t),intent(in) :: Cryst 1701 type(Pawang_type),intent(in) :: Pawang 1702 type(Pseudopotential_type),intent(in) :: Psps 1703 type(wfd_t),target,intent(inout) :: Wsh 1704 type(gs_hamiltonian_type),intent(inout) :: Ham_k 1705 !arrays 1706 real(dp),intent(in) :: kpt(3) 1707 complex(dpc),intent(in) :: vloc_ij(sh_size,sh_size) 1708 complex(dpc),intent(out) :: hk_ij(sh_size,sh_size),sk_ij(sh_size,sh_size*Psps%usepaw) 1709 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw) 1710 type(paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw) 1711 1712 !Local variables ------------------------------ 1713 !scalars 1714 integer,parameter :: idir0=0,ider0=0,nnlout0=0,tim_nonlop0=0,k1=1 1715 integer :: ig,natom,npw_k,optder,matblk,mkmem_,nkpg,dimffnl,nspinortot 1716 integer :: blc1,blc2,istwf_k,useylm_,iat,nspinor 1717 complex(dpc) :: kin_ij,vnl_ij 1718 !character(len=500) :: msg 1719 type(kdata_t) :: Kdata 1720 !arrays 1721 integer :: nloalg(3) 1722 integer,pointer :: kg_k(:,:) 1723 integer,allocatable :: nlmn_sort(:) 1724 real(dp) :: k4intp(3),kptns_(3,1),ylmgr_dum(1,1,1),shifts(3) 1725 real(dp),allocatable :: kdotg(:),half_gsq(:),ylm_k(:,:),dum_ylm_gr_k(:,:,:) 1726 real(dp),pointer :: ffnl(:,:,:,:) 1727 real(dp),allocatable :: kpg_k(:,:),ph3d(:,:,:),vnl_psi(:,:),vectin(:,:),s_psi(:,:) 1728 complex(gwpc),pointer :: ug1(:) 1729 complex(gwpc),allocatable :: cvnl_psi(:),cs_psi(:),wsg(:) 1730 type(pawcprj_type),allocatable :: Cp_blc2(:,:) 1731 1732 !************************************************************************ 1733 1734 DBG_ENTER("COLL") 1735 1736 ABI_UNUSED(Paw_ij%cplex) 1737 1738 ABI_CHECK(Wsh%nspinor==1,"nspinor==2 not coded") 1739 ABI_CHECK(Wsh%nsppol==1,"Wsh%nsppol must be 1") 1740 ABI_CHECK(Wsh%paral_kgb/=1,"paral_kgb not coded") 1741 1742 natom = Cryst%natom 1743 useylm_ = Psps%useylm 1744 nloalg = Wsh%nloalg 1745 nspinor = Wsh%nspinor 1746 nspinortot = nspinor 1747 1748 ! Gamma-centered basis set. 1749 npw_k = Wsh%Kdata(k1)%npw 1750 istwf_k = Wsh%istwfk(k1) 1751 ABI_CHECK(istwf_k==1,"istwfk/=0 not coded") 1752 kg_k => Wsh%Kdata(k1)%kg_k 1753 1754 ! First wrap in the first BZ thus enforcing traslational invariance. 1755 call wrap2_pmhalf(kpt(:),k4intp(:),shifts(:)) 1756 1757 ! Continue to prepare the GS Hamiltonian. 1758 call load_spin_hamiltonian(Ham_k,spin,with_nonlocal=.true.) 1759 1760 call kdata_init(Kdata,Cryst,Psps,k4intp,istwf_k,Wsh%ngfft,Wsh%MPI_enreg,kg_k=kg_k) 1761 1762 ABI_MALLOC(nlmn_sort,(Cryst%natom)) 1763 iat=0 ! nlmn dims sorted by atom type. 1764 1765 if (Psps%usepaw==1) then 1766 nlmn_sort = Wsh%nlmn_sort 1767 else ! FIXME here lmnmax == lnmax if useylm_==0 1768 nlmn_sort = 9 1769 if (Psps%useylm==1) then 1770 nlmn_sort=Psps%lmnmax 1771 else 1772 MSG_ERROR("useylm==0 not coded") 1773 end if 1774 !do itypat=1,Cryst%ntypat 1775 ! nlmn_sort(iat+1:iat+Cryst%nattyp(itypat))=Pawtab(itypat)%lmn_size 1776 ! iat=iat+Cryst%nattyp(itypat) 1777 !end do 1778 !write(std_out,*)" hacking nlmn_sort",nlmn_sort 1779 !write(std_out,*)" Psps%lmnmax is ",Psps%lmnmax 1780 end if 1781 ! 1782 ! ============================================================================ 1783 ! ==== Evaluate <p_lmn|e^(ikr)U_i> for each k on the k-grid and each U_i ==== 1784 ! ============================================================================ 1785 ! 1786 ! Here I assume that the G-sphere is gamma-centered. 1787 ! Real Spherical Harmonics are always used to apply the non-local part even for NC pseudos. 1788 ! I did non find any easy way to extract only <p_nl|psi> from nonlop_pl. 1789 !useylm_=1 1790 1791 ABI_DT_MALLOC(Cp_blc2 ,(natom,nspinor)) 1792 call pawcprj_alloc(Cp_blc2, 0,nlmn_sort) 1793 1794 ! ======================= 1795 ! ==== Interpolation ==== 1796 ! ======================= 1797 ! 1798 ! Prepare the kinetic term. 1799 ABI_MALLOC(kdotg,(npw_k)) 1800 ABI_MALLOC(half_gsq,(npw_k)) 1801 ABI_MALLOC(wsg,(npw_k)) 1802 1803 ! TODO Add new overloaded interface. effmass_free option! 1804 do ig=1,npw_k 1805 kdotg(ig) = two_pi**2 * DOT_PRODUCT(k4intp,MATMUL(Cryst%gmet,kg_k(:,ig))) 1806 half_gsq(ig) = half * vdotw(one*kg_k(:,ig),one*kg_k(:,ig),Cryst%gmet,"G") 1807 end do 1808 1809 ABI_MALLOC(vectin,(2, npw_k*nspinor)) 1810 ABI_MALLOC(vnl_psi,(2,npw_k*nspinor)) 1811 ABI_MALLOC(cvnl_psi,(npw_k*nspinor)) 1812 1813 ABI_MALLOC(s_psi,(2,npw_k*nspinor*Psps%usepaw)) 1814 ABI_MALLOC(cs_psi,(npw_k*nspinor*Psps%usepaw)) 1815 1816 ! THIS PART IS NEEDED FOR THE CALL TO opernl although some quantities won't be used. 1817 ! Now I do things cleanly then we try to pass zero-sized arrays! 1818 ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*useylm_)) 1819 1820 if (useylm_==1) then 1821 kptns_(:,1)=k4intp; optder=0; mkmem_=1 1822 ABI_MALLOC(dum_ylm_gr_k,(npw_k,3+6*(optder/2),Psps%mpsang**2)) 1823 1824 ! Here mband is not used if paral_compil_kpt=0 1825 call initylmg(Cryst%gprimd,kg_k,kptns_,mkmem_,Wsh%MPI_enreg,Psps%mpsang,npw_k,(/1/),1,& 1826 & (/npw_k/),1,optder,Cryst%rprimd,ylm_k,dum_ylm_gr_k) 1827 1828 ABI_FREE(dum_ylm_gr_k) 1829 end if 1830 ! 1831 ! Compute (k+G) vectors (only if useylm_=1) 1832 nkpg=3*nloalg(3) 1833 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 1834 if (nkpg>0) then 1835 call mkkpg(kg_k,kpg_k,k4intp,nkpg,npw_k) 1836 endif 1837 ! 1838 ! ======================================================== 1839 ! ==== Compute nonlocal form factors ffnl at all (k+G) ==== 1840 ! ======================================================== 1841 dimffnl=1+ider0 ! Derivatives are not needed. 1842 ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat)) 1843 !ffnl => Kdata%fnl_dir0der0 1844 call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,Cryst%gmet,Cryst%gprimd,ider0,idir0,Psps%indlmn,& 1845 & kg_k,kpg_k,k4intp,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,& 1846 & Psps%ntypat,Psps%pspso,Psps%qgrid_ff,Cryst%rmet,Psps%usepaw,useylm_,ylm_k,ylmgr_dum) 1847 ABI_FREE(ylm_k) 1848 ! 1849 ! Load k-dependent part in the Hamiltonian datastructure 1850 matblk=min(NLO_MINCAT,maxval(Ham_k%nattyp)) ; if (nloalg(2)>0) matblk=natom 1851 ABI_MALLOC(ph3d,(2,npw_k,matblk)) 1852 call load_k_hamiltonian(Ham_k,kpt_k=k4intp,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,& 1853 & kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=(Wsh%paral_kgb/=1)) 1854 1855 !END BORING CODE NEEDED TO CALL opernl 1856 1857 ! 1858 ! Calculate the upper triangle of <blc2| H_k |blc1>. 1859 do blc2=1,sh_size 1860 1861 ! Calculate <G|Vnl|psi> for this k-point 1862 call wfd_vnlpsi(Wsh,blc2,k1,spin,npw_k,Cryst,Ham_k,vnl_psi,s_psi,Kext=Kdata) 1863 cvnl_psi = DCMPLX(vnl_psi(1,:),vnl_psi(2,:)) 1864 if (Psps%usepaw==1) cs_psi = DCMPLX(s_psi(1,:),s_psi(2,:)) 1865 1866 ! Upper triangle of the hk_ij matrix 1867 do blc1=1,blc2 1868 ug1 => Wsh%Wave(blc1,k1,1)%ug 1869 1870 ! Treat the matrix elements of the Vnl operator. 1871 vnl_ij = xdotc(npw_k*nspinor,ug1,1,cvnl_psi,1) 1872 ! 1873 ! =================================================== 1874 ! ==== Assemble final Hamiltonian T + Vloc + Vnl ==== 1875 ! =================================================== 1876 ! 1877 ! Kinetic energy. 1878 wsg = kdotg * Wsh%Wave(blc2,k1,1)%ug 1879 kin_ij = xdotc(npw_k,ug1,1,wsg,1) 1880 1881 wsg = half_gsq * Wsh%Wave(blc2,k1,1)%ug 1882 kin_ij = kin_ij + xdotc(npw_k,ug1,1,wsg,1) 1883 if (blc1==blc2) kin_ij = kin_ij + half * vdotw(k4intp,k4intp,Cryst%gmet,"G") 1884 ! 1885 ! Total Hamiltonian. 1886 hk_ij(blc1,blc2) = kin_ij + vloc_ij(blc1,blc2) + vnl_ij 1887 ! 1888 ! PAW Overlap operator. 1889 if (Psps%usepaw==1) sk_ij(blc1,blc2) = xdotc(npw_k*nspinor,ug1,1,cs_psi,1) 1890 end do ! blc1 1891 end do ! blc2 1892 1893 !DEALLOCATE BORING STUFF 1894 ABI_FREE(ffnl) 1895 ABI_FREE(kpg_k) 1896 ABI_FREE(ph3d) 1897 !END BORING STUFF 1898 1899 ABI_FREE(kdotg) 1900 ABI_FREE(half_gsq) 1901 ABI_FREE(wsg) 1902 ABI_FREE(vectin) 1903 ABI_FREE(vnl_psi) 1904 ABI_FREE(cvnl_psi) 1905 ABI_FREE(s_psi) 1906 ABI_FREE(cs_psi) 1907 1908 call pawcprj_free(Cp_blc2 ) 1909 ABI_DT_FREE(Cp_blc2) 1910 call kdata_free(Kdata) 1911 1912 ABI_FREE(nlmn_sort) 1913 1914 DBG_EXIT("COLL") 1915 1916 RETURN 1917 1918 ABI_UNUSED(Pawang%ngnt) 1919 ABI_UNUSED(Pawtab(1)%basis_size) 1920 1921 end subroutine shirley_hks
m_shirley/shirley_interp [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
shirley_interp
FUNCTION
Compute interpolated KS eigenvalues and optionally KS eigenvectors by direct diagonalization of the KS Hamiltonian in the Shirley basis set. Main entry point for client code.
INPUTS
jobz="V" if eigevectors are wanted. "N" if only eigenvalues are required. Dtset=Input variables for this run (it will be removed later on) Cryst<crystal_t>= data type gathering info on symmetries and unit cell Psps <pseudopotential_type>=variables related to pseudopotentials Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data Pawfgr<Pawfgr_type> Pawang<pawang_type> angular mesh discretization and related data: Pawrad<Pawrad_type> Pawrhoij Paw_ij(natom)<paw_ij_type>=data structure containing PAW arrays given on (i,j) channels. ngfftc(18)=Information about the coarse 3D FFT. ngfftf(18)=Information about the dense 3D FFT used for vtrial. nfftf=Number of points in the FFT grid in vtrial. Might differ from the FFT mesh used for the wavefunctions. vtrial(nfftf,nspden)= Trial potential (Hartree) intp_nband(intp_nk,Wsh%nsppol)=Number of interpolated bands intp_mband=Max number of interpolated bands (used for dimensioning arrays) intp_nk=Number of interpolated k-points. intp_kpt(3,intp_nk)=Reduced coordinates of the interpolated k-points. sh_fname=String with the filename of the netcdf file to be produced. No file is created if empty string. [kweights(intp_nk)]=Weights of the k-points used in oBands, Assume 0 if not present
OUTPUT
oBands=Interpolated band energies. Note: weights are initialized to zero. [oWfd]
PARENTS
m_shirley,wfk_analyze
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
1972 subroutine shirley_interp(Wsh,jobz,Dtset,Cryst,Psps,Pawtab,Pawfgr,Pawang,Pawrad,& 1973 & Pawrhoij,Paw_ij,ngfftc,ngfftf,nfftf,vtrial,& 1974 & intp_nband,intp_mband,intp_nk,intp_kpt,sh_fname,oBands,& 1975 & kweights,oWfd) ! optional 1976 1977 1978 !This section has been created automatically by the script Abilint (TD). 1979 !Do not modify the following lines by hand. 1980 #undef ABI_FUNC 1981 #define ABI_FUNC 'shirley_interp' 1982 use interfaces_14_hidewrite 1983 use interfaces_67_common 1984 !End of the abilint section 1985 1986 implicit none 1987 1988 !Arguments ------------------------------------ 1989 !scalars 1990 integer,intent(in) :: nfftf,intp_nk,intp_mband 1991 character(len=*),intent(in) :: jobz,sh_fname 1992 type(wfd_t),intent(inout) :: Wsh 1993 type(crystal_t),intent(in) :: Cryst 1994 type(Pawang_type),intent(in) :: Pawang 1995 type(Pseudopotential_type),intent(in) :: Psps 1996 type(dataset_type),intent(in) :: Dtset 1997 type(Pawfgr_type),intent(in) :: Pawfgr 1998 type(ebands_t),intent(out) :: oBands 1999 !arrays 2000 integer,intent(in) :: ngfftf(18),ngfftc(18) 2001 integer,intent(in) :: intp_nband(intp_nk,Wsh%nsppol) 2002 real(dp),intent(in) :: vtrial(nfftf,Wsh%nspden),intp_kpt(3,intp_nk) 2003 real(dp),optional,intent(in) :: kweights(intp_nk) 2004 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wsh%usepaw) 2005 type(paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Wsh%usepaw) 2006 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*Wsh%usepaw) 2007 type(pawrhoij_type),intent(in) :: Pawrhoij(Cryst%natom*Wsh%usepaw) 2008 type(wfd_t),optional,intent(out) :: oWfd 2009 2010 !Local variables ------------------------------ 2011 !scalars 2012 integer,parameter :: istwf1=1,k1=1 2013 #ifdef HAVE_NETCDF 2014 integer,parameter :: dummy_nshiftk1=1,dummy_nshiftk_orig1=1 2015 #endif 2016 integer :: ii,ib,jj,ierr,nband_k,ikpt,natom,nefound,band 2017 integer :: sh1,sh2,ldz,prtvol,pawprtvol,istwf_k 2018 integer :: mgfftc,nfftc,onband_diago,usepaw,spin 2019 integer :: comm_spin,nproc_spin,rank_spin,root 2020 integer :: npw_k,nspinor,nsppol,nspden,sh_size,intp_bantot 2021 integer :: my_shstart,my_shstop,my_kstart,my_kstop 2022 real(dp) :: fft_fact,ene_fact,cpu,wall,gflops 2023 !logical,parameter :: debug_with_diago=.TRUE. 2024 logical,parameter :: debug_with_diago=.FALSE. 2025 logical :: want_eigenvectors,do_full_diago 2026 character(len=500) :: msg,frmt1,frmt2 2027 type(ddiago_ctl_type) :: Diago_ctl 2028 type(gs_hamiltonian_type) :: Ham_k 2029 type(stats_t) :: Stats 2030 !arrays 2031 #ifdef HAVE_NETCDF 2032 integer,parameter :: dummy_ngkpt(3)=0,dummy_kptrlatt(3,3)=0 2033 real(dp),parameter :: dummy_shiftk(3,dummy_nshiftk1)=zero,dummy_shiftk_orig(3,dummy_nshiftk_orig1)=zero 2034 #endif 2035 integer :: nloalg(3),intp_npwarr(intp_nk),intp_istwfk(intp_nk) 2036 integer,allocatable :: nlmn_sort(:) 2037 real(dp) :: intp_ene(intp_mband,intp_nk,Wsh%nsppol) 2038 real(dp) :: intp_wtk(intp_nk),kpoint(3) 2039 real(dp),pointer :: diag_ene(:),diag_vec(:,:,:) 2040 real(dp),allocatable :: enek_ij(:),vnl_psi(:,:),opaw_psi(:,:) 2041 real(dp),allocatable :: copy_vtrial(:,:) 2042 real(dp),allocatable :: intp_doccde(:),intp_occ(:),ugly_ene(:) 2043 complex(gwpc),allocatable :: ur1(:),ur2(:),vloc_psi(:) 2044 complex(dpc),allocatable :: hk_ij(:,:),sk_ij(:,:),eig_vec(:,:),vloc_ij(:,:) 2045 character(len=10) :: spin_name(2) 2046 type(pawcprj_type),pointer :: diag_Cprj(:,:) 2047 !BEGIN For the output wfd. 2048 integer :: sh_npw 2049 logical,allocatable :: keep_ur(:,:,:),bks_mask(:,:,:) 2050 complex(gwpc),allocatable :: ug(:) 2051 !END For the output wfd. 2052 #ifdef HAVE_NETCDF 2053 integer :: ncid 2054 #endif 2055 2056 !************************************************************************ 2057 2058 DBG_ENTER("COLL") 2059 2060 !call wrtout(std_out,"Starting Shirley interpolation on: "asctime(),"COLL") 2061 2062 ABI_UNUSED(Pawrhoij(1)%cplex) 2063 ABI_UNUSED(Pawrad(1)%mesh_size) 2064 2065 ABI_CHECK(Wsh%nspinor==1,"nspinor==2 not coded") 2066 ABI_CHECK(Wsh%paral_kgb==0,"paral_kgb/=0 not coded") 2067 ABI_CHECK(Wsh%rfft_is_symok,"FFT not symmetric in real space") 2068 2069 if (maxval(intp_nband) > maxval(Wsh%npwarr)) then 2070 write(msg,'(a,i0,2a,i0)')& 2071 & "The number of bands to be interpolated: ",maxval(intp_nband),ch10,& 2072 & "cannot be greater than the size of the Shirley basis: ",maxval(Wsh%npwarr) 2073 MSG_ERROR(msg) 2074 end if 2075 2076 call cwtime(cpu,wall,gflops,"start") 2077 2078 ! Copy basic dimensions and parameters. 2079 nspinor = Wsh%nspinor 2080 nsppol = Wsh%nsppol 2081 nspden = Wsh%nspden 2082 nloalg = Wsh%nloalg 2083 usepaw = Wsh%usepaw 2084 natom = Wsh%natom 2085 prtvol = Wsh%prtvol 2086 pawprtvol = Wsh%pawprtvol 2087 2088 ! The coarse mesh used for the Hamiltonian. 2089 nfftc = PRODUCT(ngfftc(1:3)) 2090 mgfftc = MAXVAL(ngfftc(1:3)) 2091 2092 if (nsppol==1) spin_name=(/' ',' '/) 2093 if (nsppol==2) spin_name=(/'SPIN_UP ','SPIN_DOWN '/) 2094 2095 want_eigenvectors = firstchar(jobz,"V",csens=.FALSE.) 2096 if (want_eigenvectors) then 2097 call wrtout(std_out," Eigenvectors will be calculated." ,"COLL") 2098 else 2099 call wrtout(std_out," Interpolating energies only. Eigenvectors won't be calculated.","COLL") 2100 end if 2101 2102 !BEGIN DEBUG: Test on the orthogonalization of the input wavefunctions. 2103 ! call wfd_reset_ur_cprj(Wsh) 2104 ! call wfd_test_ortho(Wsh,Cryst,Pawtab,unit=std_out,mode_paral="COLL") 2105 !END DEBUG 2106 2107 ! vtrial might be given on a FFT mesh that is denser than the FFT used for Wsh. 2108 ! If the two FFTs differ, change the mesh for the wavefunctions. 2109 ! Another possibility would be moving vtrial from the dense to the coarse mesh. 2110 call wfd_change_ngfft(Wsh,Cryst,Psps,ngfftf) 2111 2112 ! Init a new wavefunction descriptor to store the interpolated KS states 2113 ! The optimal basis set is given on the gamma centered basis set with istwfk==1. 2114 ! Note that the basis set in oWfd is not k-centered since its u(g) are given 2115 ! in terms of a linear combination of Shirley u(g) that are Gamma-centered 2116 if (present(oWfd)) then 2117 ABI_CHECK(want_eigenvectors,"When oWfd is present, want_eigenvectors must be true") 2118 2119 ! TODO: BE careful in parallel when nsppol==2. I should recreate the communicators. 2120 ABI_MALLOC(bks_mask,(intp_mband,intp_nk,nsppol)) 2121 ABI_MALLOC(keep_ur ,(intp_mband,intp_nk,nsppol)) 2122 bks_mask=.TRUE.; keep_ur =.TRUE. 2123 2124 ! Build new wavefunction descriptor (Gamma centered) 2125 !sh_gvec => Wsh%Kdata(1)%kg_k 2126 intp_istwfk = istwf1; sh_npw = Wsh%npwarr(1) 2127 2128 call wfd_init(oWfd,Cryst,Pawtab,Psps,keep_ur,Wsh%paral_kgb,sh_npw,intp_mband,intp_nband,intp_nk,Wsh%nsppol,bks_mask,& 2129 & Wsh%nspden,Wsh%nspinor,Wsh%ecutsm,Wsh%dilatmx,intp_istwfk,intp_kpt,Wsh%ngfft,& 2130 & Wsh%Kdata(k1)%kg_k,Wsh%nloalg,Wsh%prtvol,Wsh%pawprtvol,Wsh%comm) 2131 2132 if (oWfd%prtvol > 0) then 2133 call wfd_print(oWfd,header="New wavefunction descriptor with interpolated states") 2134 end if 2135 2136 ABI_FREE(bks_mask) 2137 ABI_FREE(keep_ur) 2138 end if 2139 ! 2140 ! ======================= 2141 ! ==== Interpolation ==== 2142 ! ======================= 2143 fft_fact = one/Wsh%nfft 2144 ABI_MALLOC(ur1,(Wsh%nfft*nspinor)) 2145 ABI_MALLOC(ur2,(Wsh%nfft*nspinor)) 2146 2147 intp_ene=zero 2148 2149 ! Precompute <sh1| vloc_spin |sh2> as it is not k-dependent. 2150 do spin=1,nsppol 2151 if (.not. wfd_itreat_spin(Wsh,spin,comm_spin,rank_spin,nproc_spin)) cycle 2152 2153 ! Split the k-points inside comm_spin 2154 !my_kstart=1; my_kstop=intp_nk 2155 call xmpi_split_work(intp_nk,comm_spin,my_kstart,my_kstop,msg,ierr) 2156 if (ierr==2) MSG_WARNING(msg) 2157 write(std_out,*)"Will treat from my_kstart ",my_kstart,"to my_kstop ",my_kstop 2158 ! 2159 ! Compute the upper triangle of the <sh2|vloc|sh1> matrix. 2160 ! All MPI(spin) processors participate 2161 ! Be careful here as |\tpsi> is not normalized. 2162 ABI_MALLOC(vloc_psi,(Wsh%nfft*nspinor)) 2163 2164 sh_size = Wsh%nband(k1,spin) 2165 ABI_MALLOC(vloc_ij,(sh_size,sh_size)) 2166 vloc_ij = czero 2167 2168 !my_shstart=1; my_shstop=sh_size 2169 call xmpi_split_work(sh_size,comm_spin,my_shstart,my_shstop,msg,ierr) 2170 if (ierr==2) MSG_WARNING(msg) 2171 2172 do sh2=my_shstart,my_shstop 2173 call wfd_get_ur(Wsh,sh2,k1,spin,ur2) 2174 if (nspinor==1) then 2175 vloc_psi = ur2*vtrial(:,spin) 2176 else 2177 MSG_ERROR("vloc_psi doesn't support nspinor==2") 2178 end if 2179 ! 2180 ! Diagonal element. 2181 vloc_ij(sh2,sh2) = xdotc(Wsh%nfft,ur2,1,vloc_psi,1)*fft_fact 2182 ! 2183 ! Upper triangle. 2184 do sh1=1,sh2-1 2185 call wfd_get_ur(Wsh,sh1,k1,spin,ur1) 2186 vloc_ij(sh1,sh2) = xdotc(Wsh%nfft,ur1,1,vloc_psi,1)*fft_fact 2187 end do 2188 end do !sh2 2189 ABI_FREE(vloc_psi) 2190 2191 ! Gather results. 2192 call xmpi_sum(vloc_ij,comm_spin,ierr) 2193 ! 2194 ! ============================================= 2195 ! ==== Loop over the interpolated k-points ==== 2196 ! ============================================= 2197 2198 call init_hamiltonian(Ham_k,Psps,Pawtab,nspinor,nsppol,nspden,natom,& 2199 & Cryst%typat,Cryst%xred,Wsh%nfft,Wsh%mgfft,Wsh%ngfft,Cryst%rprimd,nloalg) 2200 2201 ABI_MALLOC(hk_ij, (sh_size,sh_size)) 2202 ABI_MALLOC(sk_ij,(sh_size,sh_size*usepaw)) 2203 2204 npw_k = Wsh%Kdata(k1)%npw 2205 2206 ABI_MALLOC(vnl_psi,(2,npw_k*nspinor)) 2207 ABI_MALLOC(opaw_psi,(2,npw_k*nspinor*usepaw)) 2208 2209 ! Loop over the K-points for the interpolation (MPI loop) 2210 ! (NB: k will be wrapped in ]-1/2,1/2]). 2211 do ikpt=my_kstart,my_kstop 2212 istwf_k = 1 ! FIXME no time-reversal tricks! 2213 kpoint = intp_kpt(:,ikpt) 2214 nband_k = intp_nband(ikpt,spin) 2215 2216 ! Compute <i|H_k|j> in the Shirley basis set. 2217 call shirley_hks(Wsh,kpoint,spin,Ham_k,Cryst,Psps,Pawtab,Pawang,Paw_ij,sh_size,vloc_ij,hk_ij,sk_ij) 2218 ! 2219 ! Diagonalize Hk_ij in the optimal Bloch supspace. 2220 ABI_MALLOC(enek_ij,(sh_size)) 2221 2222 do_full_diago = (nband_k == sh_size) 2223 if (.not.do_full_diago) then 2224 ldz=1; if (want_eigenvectors) ldz=sh_size 2225 ABI_MALLOC(eig_vec,(ldz,nband_k)) 2226 end if 2227 !do_full_diago = .True. 2228 !write(std_out,*)"full diago",do_full_diago 2229 2230 if (usepaw==0) then 2231 ! Solve H*v = e*v 2232 if (do_full_diago) then 2233 call xheev(jobz,"Upper",sh_size,hk_ij,enek_ij) 2234 else 2235 call xheevx(jobz,"Irange","Upper",sh_size,hk_ij,zero,zero,1,nband_k,-tol8,nefound,enek_ij,eig_vec,ldz) 2236 if (want_eigenvectors) hk_ij(:,1:ldz) = eig_vec 2237 end if 2238 else 2239 ! Solve H*v = e*S*v 2240 if (do_full_diago) then 2241 call xhegv(1,jobz,"Upper",sh_size,hk_ij,sk_ij,enek_ij) 2242 else 2243 call xhegvx(1,jobz,"Irange","Upper",sh_size,hk_ij,sk_ij,zero,zero,1,nband_k,-tol8,nefound,enek_ij,eig_vec,ldz) 2244 if (want_eigenvectors) hk_ij(:,1:ldz) = eig_vec 2245 end if 2246 end if 2247 ! 2248 ! Store the interpolated eigenvalues and eigenstates in the optimal basis set. 2249 if (want_eigenvectors) then 2250 ! Compute interpolated eigenvectors in G-space. 2251 ABI_CHECK(oWfd%npwarr(ikpt) == Wsh%Kdata(k1)%npw, "Mismatch oWfd, Wsh") 2252 ABI_MALLOC(ug, (nspinor*oWfd%npwarr(ikpt))) 2253 do band=1,nband_k 2254 ug = czero 2255 do ii=1,sh_size 2256 ug(:) = ug(:) + hk_ij(ii,band) * Wsh%Wave(ii,k1,spin)%ug 2257 end do 2258 call wfd_push_ug(oWfd,band,ikpt,spin,Cryst,ug,update_ur=.FALSE.,update_cprj=.FALSE.) 2259 end do 2260 ABI_FREE(ug) 2261 end if 2262 2263 if (allocated(eig_vec)) then 2264 ABI_FREE(eig_vec) 2265 end if 2266 2267 if (prtvol>0) then 2268 ! Write interpolated energies. 2269 ene_fact=Ha_eV; frmt1='(i4,4x,9(1x,f7.4))'; frmt2='(8x,9(1x,f7.4))' 2270 write(msg,'(a,3es16.8,2a)')' Eigenvalues in eV for kpt= ( ',kpoint," ), spin ",spin_name(spin) 2271 call wrtout(std_out,msg,'COLL') 2272 2273 write(msg,frmt1)ikpt,(enek_ij(ib)*ene_fact,ib=1,MIN(9,nband_k)) 2274 call wrtout(std_out,msg,'COLL') 2275 2276 if (nband_k>9) then 2277 do jj=10,nband_k,9 2278 write(msg,frmt2) (enek_ij(ib)*ene_fact,ib=jj,MIN(jj+8,nband_k)) 2279 call wrtout(std_out,msg,'COLL') 2280 end do 2281 end if 2282 call flush_unit(std_out) 2283 end if 2284 ! 2285 ! Save the interpolated energies. 2286 intp_ene(1:nband_k,ikpt,spin) = enek_ij(1:nband_k) 2287 2288 if (debug_with_diago) then ! Compare interpolated energies wrt direct diago results. 2289 2290 ! FIXME here there is a problem with Wd%ecut and Dtset%ecut. 2291 !if (ikpt==1) write(std_out,*)" CHECK: Dtset%ecut=",Dtset%ecut," Wsh%ecut= ",Wsh%ecut 2292 !call init_ddiago_ctl(Diago_ctl,"No Vectors",spin,nspinor,Wsh%ecut,kpoint,nloalg,Cryst%gmet,& 2293 !& nband_k=nband_k,effmass_free=Dtset%effmass_free,istwf_k=istwf1,prtvol=prtvol) 2294 2295 call init_ddiago_ctl(Diago_ctl,"No Vectors",spin,nspinor,Dtset%ecut,kpoint,nloalg,Cryst%gmet,& 2296 & nband_k=nband_k,effmass_free=Dtset%effmass_free,istwf_k=istwf1,prtvol=prtvol) 2297 2298 nullify(diag_ene) 2299 nullify(diag_vec) 2300 nullify(diag_Cprj) 2301 2302 ABI_MALLOC(copy_vtrial,(nfftf,nspden)) 2303 copy_vtrial = vtrial ! To avoid intent inout 2304 2305 call ks_ddiago(Diago_ctl,nband_k,nfftc,mgfftc,ngfftc,natom,& 2306 & Cryst%typat,nfftf,nspinor,nspden,nsppol,Pawtab,Pawfgr,Paw_ij,& 2307 & Psps,Cryst%rprimd,copy_vtrial,Cryst%xred,onband_diago,diag_ene,& 2308 & diag_vec,diag_Cprj,xmpi_comm_self,ierr) 2309 ABI_CHECK(ierr==0,"Fatal error. Cannot continue!") 2310 2311 ABI_FREE(copy_vtrial) 2312 2313 if (.TRUE..or.prtvol>0) then 2314 ene_fact=Ha_eV 2315 write(msg,'(a,3es16.8,2a)')' Diago Eigenvalues in eV for kpt= ( ',kpoint," ), spin ",spin_name(spin) 2316 call wrtout(std_out,msg,'COLL') 2317 2318 write(msg,'(i4,4x,9(1x,f7.4))')ikpt,(ene_fact*diag_ene(ii),ii=1,MIN(9,nband_k)) 2319 call wrtout(std_out,msg,'COLL') 2320 2321 if (nband_k>9) then 2322 do jj=10,nband_k,9 2323 write(msg,'(8x,9(1x,f7.4))')(ene_fact*diag_ene(ii),ii=jj,MIN(jj+8,nband_k)) 2324 call wrtout(std_out,msg,'COLL') 2325 end do 2326 end if 2327 end if 2328 2329 ! Compute statistical parameters 2330 Stats = stats_eval(ABS(diag_ene(1:nband_k-2) - enek_ij(1:nband_k-2))*Ha_eV*1000) 2331 2332 write(std_out,'(a,f7.2,a,i0,a,i0)')& 2333 & " MAX abs error diag_ene - intp_ene ",MAXVAL(ABS(diag_ene(1:nband_k) - enek_ij(1:nband_k)) )*Ha_eV*1000,& 2334 & " [meV] for band ",imax_loc(ABS(diag_ene(1:nband_k) - enek_ij(1:nband_k)) ),"/",nband_k 2335 write(std_out,'(4(a,es9.1),a)')& 2336 & " min: ",Stats%min,", Max: ",Stats%max,", mean: ",Stats%mean,", stdev: ",Stats%stdev," [meV]" 2337 call flush_unit(std_out) 2338 2339 if (associated(diag_ene)) then 2340 ABI_FREE(diag_ene) 2341 end if 2342 if (associated(diag_vec)) then 2343 ABI_FREE(diag_vec) 2344 end if 2345 if (associated(diag_Cprj)) then 2346 end if 2347 end if 2348 2349 ABI_FREE(enek_ij) 2350 end do ! intp_kpt 2351 ! 2352 ! Interpolation completed, deallocate memory. 2353 ABI_FREE(hk_ij) 2354 ABI_FREE(vloc_ij) 2355 ABI_FREE(sk_ij) 2356 ABI_FREE(vnl_psi) 2357 2358 ABI_FREE(opaw_psi) 2359 call destroy_hamiltonian(Ham_k) 2360 end do ! spin 2361 ! 2362 ! Collect energies on each node. 2363 call xmpi_sum(intp_ene,Wsh%comm,ierr) 2364 ! 2365 ! Free memory 2366 ABI_FREE(ur1) 2367 ABI_FREE(ur2) 2368 if (allocated(nlmn_sort)) then 2369 ABI_FREE(nlmn_sort) 2370 end if 2371 2372 ! TODO 2373 ! If vectors we have to write the unitary transformation with MPI-IO 2374 2375 ! Create ebands_t for storing the interpolated energies. 2376 intp_bantot=SUM(intp_nband) 2377 ABI_CALLOC(intp_doccde,(intp_bantot)) 2378 ABI_CALLOC(intp_occ,(intp_bantot)) 2379 intp_istwfk=1; intp_npwarr=Wsh%npwarr(k1) ! Useless but oh well 2380 2381 ! Have to reshape intp_ene --> ugly_ene due to the ugly convention used in abinit to store energies. 2382 ABI_MALLOC(ugly_ene,(intp_bantot)) 2383 ugly_ene(:)=HUGE(zero) 2384 2385 call pack_eneocc(intp_nk,nsppol,intp_mband,intp_nband,intp_bantot,intp_ene,ugly_ene) 2386 2387 ! Warning: weigths are set to zero if kweights are not specified. 2388 intp_wtk = zero; if (present(kweights)) intp_wtk = kweights 2389 2390 MSG_ERROR("Fix ebands_init") 2391 2392 !call ebands_init(intp_bantot,oBands,Dtset%nelect,intp_doccde,ugly_ene,intp_istwfk,intp_kpt,& 2393 ! intp_nband,intp_nk,intp_npwarr,nsppol,nspinor,Dtset%tphysel,Dtset%tsmear,Dtset%occopt,intp_occ,intp_wtk,& 2394 ! charge, kptopt, kptrlatt_orig, nshiftk_orig, shiftk_orig, kptrlatt, nshiftk, shiftk) 2395 2396 ABI_FREE(intp_doccde) 2397 ABI_FREE(ugly_ene) 2398 ABI_FREE(intp_occ) 2399 2400 if (present(oWfd)) then 2401 ABI_MALLOC(bks_mask, (oWfd%mband,oWfd%nkibz,oWfd%nsppol)) 2402 ! TODO 2403 ! Collect all the wavefunctions on each MPI node. 2404 do spin=1,oWfd%nsppol 2405 if (.not. wfd_itreat_spin(oWfd,spin,comm_spin,rank_spin,nproc_spin)) cycle 2406 do root=0,nproc_spin-1 2407 bks_mask=.False. !do ik_ibz=my_kstart,my_kstop 2408 bks_mask(:,my_kstart:my_kstop,spin) = .True. 2409 !call wfd_bcast_waves(Wfd,what,bks_mask,root,spin_comm,reinit_comms=.False.) 2410 end do 2411 end do 2412 ! Set the MPI communicators in the output wavefunction descriptor. 2413 ! Master nodes in the spin communicator must communicate the bks_table 2414 ! and broadcast the updated values to the nodes in the spin communicator. 2415 if (oWfd%nsppol == 2) MSG_ERROR("Not implemented error") 2416 call wfd_set_mpicomm(oWfd) 2417 ABI_FREE(bks_mask) 2418 end if 2419 2420 call cwtime(cpu,wall,gflops,"stop") 2421 write(std_out,*)"SHTIME: Interpolation, cpu_time: ",cpu,", wall_time: ",wall 2422 call flush_unit(std_out) 2423 2424 ! Write netcdf file with results. 2425 ! FIXME: k-point info are wrong. This trick is needed so that 2426 ! abipy will detect a path instead of a BZ mesh. 2427 #ifdef HAVE_NETCDF 2428 if (Wsh%my_rank == 0 .and. len_trim(sh_fname)>0) then 2429 NCF_CHECK_MSG(nctk_open_create(ncid, sh_fname, xmpi_comm_self), "Creating Shirley file") 2430 NCF_CHECK(crystal_ncwrite(Cryst, ncid)) 2431 NCF_CHECK(ebands_ncwrite(oBands, ncid)) 2432 NCF_CHECK(nf90_close(ncid)) 2433 end if 2434 #endif 2435 2436 DBG_EXIT("COLL") 2437 2438 end subroutine shirley_interp
m_shirley/shirley_window [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
shirley_window
FUNCTION
Find an optimal basis set for the interpolation of the KS states in a energy window, e.g states around Ef. The algorithm starts with a relatively coarse k-mesh and an initial energy window (sufficiently large so that the initial basis set has enough flexibility). Constructs the initial Shirley set and uses it to interpolate the KS states on a denser mesh. The energy window is progressively reduced so that the states far from Ef are progressively filtered out (this trick allows one to reduced the size of the overlap matrix and therefore the CPU time needed for its diagonalization. The iterative algorithm stops when the DOS is converged within a given tolerance. Returns a new wavefunction descriptor with the Shirley basis set. INPUT iWfd<wfd_t>=Input wavefunctions. Dtset=Input variables for this run (it will be removed later on) Cryst<crystal_t>= data type gathering info on symmetries and unit cell. iKmesh<kmesh_t>=K-points of the wavefunctions stored in Wfd. iBands<ebands_t>=Input energies. Psps<pseudopotential_type)>=variables related to pseudopotentials. Pawtab(Psps%ntypat)<pawtab_type>=paw tabulated starting data Pawang<pawang_type>=angular mesh discretization and related data: Pawrad<Pawrad_type>=For PAW, RADial mesh discretization and related data Pawfgr<Pawfgr_type> Pawang<pawang_type> angular mesh discretization and related data: Pawrad<Pawrad_type> Pawrhoij Paw_ij(natom)<paw_ij_type>=data structure containing PAW arrays given on (i,j) channels. ngfftc(18)=Information about the coarse 3D FFT. ngfftf(18)=Information about the dense 3D FFT used for vtrial. nfftf=Number of points in the FFT grid in vtrial. Might differ from the FFT mesh used for the wavefunctions. vtrial(nfftf,nspden)= Trial potential (Hartree) sh_coverage=Parameter defining the coverage of the optimal basis set: Wsh contains N vectors where N is the first elements for which we have: \sum_i^N e_i >= sh_coverage * M with e_i being the eigenvalues of the overlap matrix ordered in descending order. M indicates the dimension of the overlap (related to the number of k-points in the BZ and the number of bands in the input iWfd descriptor. ewin(2,nsppol)=Energy window for the spin channels. Use ewin(1,:)=smallest_real and ewin(2,:)=greatest_real to disable the use of the energy window min_bsize=Minimum number of states in the optimal basis set.
SIDE EFFECTS
ewin(2,nsppol)
OUTPUT
oWsh<wfd_t>=New wavefunction descriptor with the optimal basis set. Note that the descriptor contains a single (fake) k-point.
PARENTS
wfk_analyze
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
2631 subroutine shirley_window(iWfd,Dtset,Cryst,iKmesh,iBands,Psps,Pawtab,Pawang,Pawrad,Pawfgr,Pawrhoij,Paw_ij,& 2632 & sh_coverage,ewin,ngfftc,ngfftf,nfftf,vtrial,oWsh) 2633 2634 2635 !This section has been created automatically by the script Abilint (TD). 2636 !Do not modify the following lines by hand. 2637 #undef ABI_FUNC 2638 #define ABI_FUNC 'shirley_window' 2639 !End of the abilint section 2640 2641 implicit none 2642 2643 !Arguments ------------------------------------ 2644 !scalars 2645 integer,intent(in) :: nfftf 2646 real(dp),intent(in) :: sh_coverage 2647 type(crystal_t),intent(in) :: Cryst 2648 type(Pawang_type),intent(in) :: Pawang 2649 type(kmesh_t),intent(in) :: iKmesh 2650 type(Pseudopotential_type),intent(in) :: Psps 2651 type(wfd_t),intent(inout) :: iWfd 2652 type(ebands_t),intent(in) :: iBands 2653 type(dataset_type),intent(in) :: Dtset 2654 type(Pawfgr_type),intent(in) :: Pawfgr 2655 type(wfd_t),intent(out) :: oWsh 2656 !arrays 2657 integer,intent(in) :: ngfftf(18),ngfftc(18) 2658 real(dp),intent(in) :: vtrial(nfftf,iWfd%nspden) 2659 real(dp),intent(inout) :: ewin(2,iWfd%nsppol) 2660 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*iWfd%usepaw) 2661 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*iWfd%usepaw) 2662 type(paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*iWfd%usepaw) 2663 type(pawrhoij_type),intent(in) :: Pawrhoij(Cryst%natom*iWfd%usepaw) 2664 2665 !Local variables ------------------------------ 2666 !scalars 2667 integer,parameter :: max_niter=2 2668 integer :: min_bsize,intp_mband,intp_nk,iter,nsppol !,ierr 2669 logical :: converged 2670 character(len=500) :: sh_fname 2671 type(wfd_t) :: tmpWfd 2672 !type(edos_t),pointer :: newEdos,oldEdos 2673 !type(edos_t),target :: Doses(2) 2674 type(kmesh_t) :: oKmesh ! <<<< 2675 type(ebands_t) :: oBands ! <<<< 2676 !arrays 2677 integer :: kptrlatt(3,3) 2678 integer,allocatable :: intp_nband(:,:) 2679 2680 !************************************************************************ 2681 2682 DBG_ENTER("COLL") 2683 2684 ! Compute the DOS from the input bands and k-points. 2685 !call edos_init(Doses(1),iBands,iKmesh,method,step,broad) 2686 !call edos_free(Doses(1)) 2687 !call edos_free(Doses(2)) 2688 !oldEdos => Doses(1) 2689 !newEdos => Doses(2) 2690 2691 nsppol=iWfd%nsppol 2692 min_bsize=-1 2693 2694 ! Compute Shirley basis set with the initial Kmesh and the initial window energy. 2695 call wfd_bloch_to_shirley(iWfd,Cryst,iKmesh,iBands,Psps,Pawtab,Pawang,Pawrad,min_bsize,sh_coverage,ewin,oWsh) 2696 2697 do iter=1,max_niter 2698 ! Densify the k-mesh 2699 !if (iter==1) kptrlatt = reshape([5,0,0,0,5,0,0,0,5], [3,3]) 2700 !if (iter==2) kptrlatt = reshape([6,0,0,0,6,0,0,0,6], [3,3]) 2701 if (iter==1) kptrlatt = reshape([7,0,0,0,7,0,0,0,7], [3,3]) 2702 if (iter==2) kptrlatt = reshape([8,0,0,0,8,0,0,0,8], [3,3]) 2703 2704 !if (iter==1) kptrlatt = reshape([2,0,0,0,2,0,0,0,2], [3,3]) 2705 !if (iter==2) kptrlatt = reshape([3,0,0,0,3,0,0,0,3], [3,3]) 2706 !if (iter==3) kptrlatt = reshape([4,0,0,0,4,0,0,0,4], [3,3]) 2707 !if (iter==4) kptrlatt = reshape([5,0,0,0,5,0,0,0,5], [3,3]) 2708 2709 call make_mesh(oKmesh,Cryst,iKmesh%kptopt,kptrlatt,Dtset%nshiftk,Dtset%shiftk) 2710 write(std_out,*)"Moving to denser Kmesh with kptrlatt ",kptrlatt 2711 2712 ! Interpolate energies and vectors on the denser mesh. 2713 intp_nk=oKmesh%nibz 2714 2715 ABI_MALLOC(intp_nband, (intp_nk,nsppol)) 2716 ! this has a huge impact on the CPU time (full versus partial diago) 2717 ! Cannot have more eigenvectors than the number of basis elements, doh! 2718 intp_nband=10 2719 if (maxval(intp_nband) > oWsh%nband(1,1)) intp_nband = oWsh%nband(1,1) 2720 intp_mband=maxval(intp_nband) 2721 2722 ! Set the weights in oBands 2723 write(sh_fname,"(a,i0,a)")"SHWITER",iter,"_GSR.nc" 2724 call shirley_interp(oWsh,"V",Dtset,Cryst,Psps,Pawtab,Pawfgr,Pawang,Pawrad,& 2725 & Pawrhoij,Paw_ij,ngfftc,ngfftf,nfftf,vtrial,& 2726 & intp_nband,intp_mband,oKmesh%nibz,oKmesh%ibz,sh_fname,oBands,kweights=oKmesh%wt,oWfd=tmpWfd) 2727 2728 if (tmpWfd%prtvol > 0) then 2729 call wfd_print(tmpWfd,header="Shirley wavefunction descriptor") 2730 end if 2731 if (DEBUGME) then 2732 call wfd_test_ortho(tmpWfd,Cryst,Pawtab,unit=std_out,mode_paral="COLL") 2733 end if 2734 2735 ! Free memory 2736 ABI_FREE(intp_nband) 2737 2738 converged = .False. 2739 ! Compute DOS and DOS(e0) and compare wrt previous one. Use same mesh as oldEdos 2740 ! Decrease the energy window if not converged and iterate 2741 !call edos_init(newEdos,oBands,oKmesh,method,step,broad,ierr,mesh=oldEdos%mesh) 2742 !ws = bisect(newEdos%mesh, ene_min) 2743 !we = bisect(newEdos%mesh, ene_max) 2744 !if (iter>1) then 2745 !converged = (l2norm(newEdos%dos(ws:we) - oldEdos%dos(ws:we)) < dos_atol .and. & 2746 !& abs(newEdos%gef(0) - oldEdos%gef(0)) < dos_atol) 2747 !call edos_free(newEdos) 2748 if (converged) exit 2749 2750 ewin = ewin 2751 !ewin = 0.8*ewin 2752 write(std_out,*)"Finding new Shirley basis set with energy window. Iteration: ",iter 2753 write(std_out,*)"Energy window: ",ewin*Ha_eV," [eV]" 2754 2755 call wfd_free(oWsh) 2756 call wfd_bloch_to_shirley(tmpWfd,Cryst,oKmesh,oBands,Psps,Pawtab,Pawang,Pawrad,min_bsize,sh_coverage,ewin,oWsh) 2757 2758 call wfd_free(tmpWfd) 2759 2760 !if (iter==3) exit 2761 !if (.not. converged) then 2762 call kmesh_free(oKmesh) 2763 call ebands_free(oBands) 2764 !end if 2765 end do 2766 2767 DBG_EXIT("COLL") 2768 2769 end subroutine shirley_window
m_shirley/shparams_t [ Types ]
[ Top ] [ m_shirley ] [ Types ]
NAME
shparams_t
FUNCTION
Store the parameters controlling the algorith for the Shirley interpolation of the KS band energies and wavefunctions.
SOURCE
95 type,public :: shparams_t 96 97 integer :: ndivsm=20 98 ! Number of division for the smallest segment of the k-path. 99 ! Used for the interpolation of band structures 100 101 ! Options used for the generazion of the K-mesh if sampling == "mesh" 102 ! See abinit documentation 103 integer :: kptopt = 3 104 integer :: kptrlatt(3,3) = -1 105 integer,allocatable :: shiftk(:,:) 106 107 real(dp) :: ewin(2) = [smallest_real, greatest_real] 108 ! Energy window (default: unused) 109 110 real(dp) :: dos_tol = tol6 111 ! Tolerance on the converged of the DOS (used if method=="iterative") 112 113 character(len=500) :: method="one-shot" 114 ! Method used for selecting an optimal basis set inside the energy window 115 ! "one-shot" 116 ! "iterative" 117 118 character(len=500) :: sampling="path" 119 ! Type of sampling 120 ! "path" 121 ! "mesh" 122 123 character(len=500) :: jobz = "Normal" 124 ! jobz="V" if eigevectors are wanted. "N" if only eigenvalues are required. 125 126 real(dp),allocatable :: kptbounds(:,:) 127 ! List of vertex for band structure interpolation. 128 129 end type shparams_t 130 131 !public :: shparams_init
m_shirley/shres_free_0D [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
shres_free_0D
FUNCTION
Free memory.
PARENTS
m_shirley
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
316 subroutine shres_free_0D(Shres) 317 318 319 !This section has been created automatically by the script Abilint (TD). 320 !Do not modify the following lines by hand. 321 #undef ABI_FUNC 322 #define ABI_FUNC 'shres_free_0D' 323 !End of the abilint section 324 325 implicit none 326 327 !Arguments ------------------------------------ 328 !scalars 329 type(shres_t),intent(inout) :: Shres 330 331 ! ************************************************************************* 332 333 !@shres_t 334 ! real 335 if (allocated(Shres%ene)) then 336 ABI_FREE(Shres%ene) 337 end if 338 ! 339 ! complex 340 if (allocated(Shres%obloch)) then 341 ABI_FREE(Shres%obloch) 342 end if 343 344 end subroutine shres_free_0D
m_shirley/shres_free_2D [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
shres_free_2D
FUNCTION
Free memory
PARENTS
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
366 subroutine shres_free_2D(Shres) 367 368 369 !This section has been created automatically by the script Abilint (TD). 370 !Do not modify the following lines by hand. 371 #undef ABI_FUNC 372 #define ABI_FUNC 'shres_free_2D' 373 !End of the abilint section 374 375 implicit none 376 377 !Arguments ------------------------------------ 378 !scalars 379 type(shres_t),intent(inout) :: Shres(:,:) 380 381 !Local variables 382 integer :: ii,jj 383 384 ! ************************************************************************* 385 386 do jj=1,size(Shres,dim=2) 387 do ii=1,size(Shres,dim=1) 388 call shres_free_0d(Shres(ii,jj)) 389 end do 390 end do 391 392 end subroutine shres_free_2D
m_shirley/shres_init [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
shres_init
FUNCTION
PARENTS
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
262 subroutine shres_init(Shres,sh_size,nband_k,intp_ene,obloch) 263 264 265 !This section has been created automatically by the script Abilint (TD). 266 !Do not modify the following lines by hand. 267 #undef ABI_FUNC 268 #define ABI_FUNC 'shres_init' 269 !End of the abilint section 270 271 implicit none 272 273 !Arguments ------------------------------------ 274 !scalars 275 integer,intent(in) :: sh_size,nband_k 276 type(shres_t),intent(inout) :: Shres 277 !arrays 278 real(dp),intent(in) :: intp_ene(nband_k) 279 complex(dpc),intent(in) :: obloch(sh_size,nband_k) 280 281 ! ************************************************************************* 282 283 !@shres_t 284 Shres%sh_size = sh_size 285 Shres%nband_k = nband_k 286 287 ABI_MALLOC(Shres%ene,(nband_k)) 288 Shres%ene = intp_ene 289 290 ABI_MALLOC(Shres%obloch,(sh_size,nband_k)) 291 Shres%obloch = obloch 292 293 end subroutine shres_init
m_shirley/shres_t [ Types ]
[ Top ] [ m_shirley ] [ Types ]
NAME
shres_t
FUNCTION
The results of the Shirley interpolation for a given k-point and spin. TODO fix issue with the indexing as we might want to read a subset of bands for ouri nterpolation moreover one might used vl and vu in the direct diagonalization to select the energy window we are interested in.
SOURCE
157 type,public :: shres_t 158 159 integer :: sh_size 160 ! Number of Shirley optimal basis set elements. 161 162 integer :: nband_k 163 ! Number of interpolated Kohn-Sham bands. 164 165 real(dp),allocatable :: ene(:) 166 ! ene(nband_k) 167 ! Interpolated energies 168 169 complex(dpc),allocatable :: obloch(:,:) 170 ! obloch(sh_size,nband_k) 171 ! Matrix with the coeffients <O_i|Bloch_b> 172 173 end type shres_t 174 175 public :: shres_init 176 public :: shres_free 177 178 interface shres_free 179 module procedure shres_free_0D 180 module procedure shres_free_2D 181 end interface shres_free
m_shirley/wfd_bloch_to_shirley [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
wfd_bloch_to_shirley
FUNCTION
Build a new wavefunction descriptor containing the Shirley basis set. INPUT Wfd<wfd_t>=Input wavefunctions. Cryst<crystal_t>= data type gathering info on symmetries and unit cell. Kmesh<kmesh_t>=K-points of the wavefunctions stored in Wfd. Bands<ebands_t>=Input energies. Psps<pseudopotential_type)>=variables related to pseudopotentials. Pawtab(Psps%ntypat)<pawtab_type>=paw tabulated starting data Pawang<pawang_type>=angular mesh discretization and related data: Pawrad<Pawrad_type>=For PAW, RADial mesh discretization and related data min_bsize=Minimum number of states in the optimal basis set. < 0 if no constraint should be enforced. sh_coverage=Parameter defining the coverage of the optimal basis set: ewin(2,nsppol)=Energy window for the spin channels. Use ewin(1,:)=smallest_real and ewin(2,:)=greatest_real to disable the use of the energy window
OUTPUT
oWsh<wfd_t>=New wavefunction descriptor with the optimal basis set. Note that Wfs contains a single k-point, i.e. gamma oWsh contains N vectors where N is the first elements for which we have: \sum_i^N e_i >= sh_coverage * M with e_i being the eigenvalues of the overlap matrix ordered in descending order. M indicates the dimension of the overlap (related to the number of k-points in the BZ and the number of bands in the input Wfd descriptor.
TODO
Add periodic replica.
PARENTS
m_shirley
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
1327 subroutine wfd_bloch_to_shirley(Wfd,Cryst,Kmesh,Bands,Psps,Pawtab,Pawang,Pawrad,min_bsize,sh_coverage,ewin,oWsh) 1328 1329 1330 !This section has been created automatically by the script Abilint (TD). 1331 !Do not modify the following lines by hand. 1332 #undef ABI_FUNC 1333 #define ABI_FUNC 'wfd_bloch_to_shirley' 1334 use interfaces_14_hidewrite 1335 !End of the abilint section 1336 1337 implicit none 1338 1339 !Arguments ------------------------------------ 1340 !scalars 1341 integer,intent(in) :: min_bsize 1342 real(dp),intent(in) :: sh_coverage 1343 type(crystal_t),intent(in) :: Cryst 1344 type(Pawang_type),intent(in) :: Pawang 1345 type(kmesh_t),intent(in) :: Kmesh 1346 type(Pseudopotential_type),intent(in) :: Psps 1347 type(wfd_t),intent(inout) :: Wfd 1348 type(ebands_t),intent(in) :: Bands 1349 type(wfd_t),intent(out) :: oWsh 1350 !arrays 1351 real(dp) :: ewin(2,Wfd%nsppol) 1352 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 1353 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*Wfd%usepaw) 1354 1355 !Local variables ------------------------------ 1356 !scalars 1357 integer,parameter :: istwf1=1,k1=1,ndat1=1,method1=1,enforce_sym1=1 1358 integer :: ik1_bz,band1,natom,sh,midx,my_start,my_stop 1359 integer :: npw_gamma,spin !,row,col,useylm_,!fft_idx,ig 1360 integer :: npw_k,nspinor,nsppol,sh_nkibz,sh_mband,ierr 1361 integer :: ii,ik_ibz,istwf_k,tmp_nfft,comm_spin,nproc_spin,rank_spin 1362 real(dp) :: fft_fact !,norm1 !sqrt_norm1,sqrt_norm2 1363 real(dp) :: cpu,wall,gflops 1364 logical :: use_sym 1365 character(len=500) :: msg 1366 !arrays 1367 integer :: ov_ngfft(18),trial_ngfft(18) 1368 integer :: sh_istwfk(1),sh_size(Wfd%nsppol),base_spin(Wfd%nsppol) 1369 integer,allocatable :: kg_gamma(:,:) !kg_k(:,:), 1370 integer,allocatable :: sh_nband(:,:) 1371 real(dp) :: sh_kibz(3,1),gamma_point(3)=(/zero,zero,zero/) 1372 real(dp) :: pawovlp(2),fft_err(3,Cryst%nsym) 1373 complex(dpc) :: cdum 1374 !complex(dpc),allocatable :: dpc_tmp(:) 1375 complex(gwpc),allocatable :: ur1(:) !,ur2(:) 1376 complex(gwpc),allocatable :: sh_ug(:),sh_ur(:,:) 1377 logical,allocatable :: sh_keep_ur(:,:,:),sh_bks_mask(:,:,:) 1378 type(pawcprj_type),allocatable :: Cp_sh1(:,:) 1379 type(ovlp_t) :: O(Wfd%nsppol),Osym(Wfd%nsppol) 1380 1381 !************************************************************************ 1382 1383 DBG_ENTER("COLL") 1384 1385 write(msg,'(a,f9.6)')" Starting transformation Bloch ==> Shirley with coverage: ",sh_coverage 1386 call wrtout(std_out,msg,"COLL") 1387 1388 ! Consistency check. 1389 ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded.") 1390 ABI_CHECK(wfd%paral_kgb==0,"paral_kgb/=0 not coded.") 1391 ABI_CHECK(Wfd%rfft_is_symok,"Real space FFT is not symmetric.") 1392 ABI_CHECK(Wfd%nsppol==1,"nsppol==2 not coded") 1393 ABI_CHECK(all(ewin(2,:) > ewin(1,:)), "Wrong energy window") 1394 1395 if (Wfd%usepaw==1) MSG_WARNING("Shirley with PAW is still under testing") 1396 1397 ABI_UNUSED(Pawrad(1)%mesh_size) 1398 1399 nsppol = Wfd%nsppol; nspinor = Wfd%nspinor 1400 natom = Wfd%natom 1401 ! 1402 ! 1) Get the overlap matrix <u_i|u_j> for this spin. 1403 ! TODO 1404 ! *) The Overlap can be calculated using a coarse real space FFT mesh provided that it is compatible 1405 ! with the space group symmetries (enforce_sym1=1) 1406 ! We only have to make sure that ngfft encloses the k-centered G-spheres. 1407 ! *) Use method1 to get a coarse FFT mesh (no problem due the convolution here) 1408 ov_ngfft(1:3) = 0; ov_ngfft(7:) = Wfd%ngfft(7:) ! Have to preserve the previous fftalg options 1409 1410 do ik_ibz=1,Wfd%nkibz 1411 istwf_k = Wfd%istwfk(ik_ibz) 1412 ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded") 1413 npw_k = Wfd%npwarr(ik_ibz) 1414 trial_ngfft(7:) = Wfd%ngfft(7:) 1415 1416 !kg_k => Wfd%Kdata(ik_ibz)%kg_k 1417 call setmesh(Cryst%gmet,Wfd%Kdata(ik_ibz)%kg_k ,trial_ngfft,npw_k,1,npw_k,tmp_nfft,& 1418 & method1,[0,0,0],Cryst,enforce_sym1,unit=dev_null) 1419 do ii=1,3 1420 ov_ngfft(ii) = MAX(ov_ngfft(ii),trial_ngfft(ii)) 1421 end do 1422 end do 1423 !ov_ngfft(4)=2*(ov_ngfft(1)/2)+1 1424 !ov_ngfft(5)=2*(ov_ngfft(2)/2)+1 1425 !ov_ngfft(6)= ov_ngfft(3) 1426 1427 !THIS WAS THE CAUSE OF THE BUG in shirley_window!!!!!!!!!!!!!!! 1428 ov_ngfft(:) = Wfd%ngfft(:) 1429 1430 ! TODO 1431 ! 1) Be careful here in parallel, 1432 ! 2) One should always include all the degenerated states in the calculation of the 1433 ! overlap so that the optimal basis set will preserve the degeneracies of the 1434 ! interpolated eigenvalues. 1435 use_sym =.FALSE. 1436 call ovlp_init(O,ewin,use_sym,ov_ngfft,Wfd,Cryst,Kmesh,Bands,Psps,Pawtab,Pawang,Pawrad) 1437 1438 if (.FALSE.) then 1439 ! TODO Work in progress 1440 use_sym =.TRUE. 1441 call ovlp_init(Osym,ewin,use_sym,ov_ngfft,Wfd,Cryst,Kmesh,Bands,Psps,Pawtab,Pawang,Pawrad) 1442 1443 do spin=1,nsppol 1444 call ovlp_diff(O(spin),Osym(spin),Kmesh,tol6,"Error in matrix elements",unit=std_out) 1445 call ovlp_free(Osym(spin)) 1446 end do 1447 MSG_ERROR("Check done") 1448 end if 1449 ! 1450 ! 2) Diagonalize the overlap matrix selecting the optimal subspace: [base_spin(spin):ovlp_size] 1451 ! In exit O(spin)%mat stores the eigenvectors (only on those processors treating this spin). 1452 sh_size = -1; base_spin=HUGE(0) 1453 do spin=1,nsppol 1454 1455 if (.not. wfd_itreat_spin(Wfd,spin,comm_spin,rank_spin,nproc_spin)) cycle 1456 call ovlp_diago_and_prune(O(spin),sh_coverage,sh_size(spin),base_spin(spin)) 1457 ! 1458 ! Make sure we have enough states. 1459 if (sh_size(spin)<min_bsize) then 1460 if (O(spin)%size<min_bsize) then 1461 write(msg,'(2(a,i0),2a)')& 1462 & "Overlap size is ",O(spin)%size," whereas min_bsize is ",min_bsize,ch10,& 1463 & "Decrease the number of bands to be interpolated or increase the number of ab-initio input states." 1464 MSG_ERROR(msg) 1465 end if 1466 sh_size(spin) = min_bsize 1467 write(msg,'(a,2i0)')" Had to enlarge Shirley subspace since input sh_size < min_bsize: ",sh_size(spin),min_bsize 1468 MSG_COMMENT(msg) 1469 end if 1470 end do 1471 ! 1472 ! 3) Init a new wavefunction descriptor to store the Shirley basis set. 1473 ! 1474 ! *) oWsh must be allocated here since sh_size is known only after the diagonalization of the overlap. 1475 ! *) Keep the optimal wavefunctions on each node (if possible) to facilitate the interpolation over the fine k-mesh. 1476 ! *) Use Gamma-centered basis set to facilitate the operations in reciprocal space. 1477 ! *) The new basis is orthogonal, but not normalized since <U_i|U_j> = delta_ij e_i. 1478 ! 1479 ! The optimal basis set is centered on gamma and uses istwfk==1. 1480 call cwtime(cpu,wall,gflops,"start") 1481 1482 call get_kg(gamma_point,istwf1,Wfd%ecut,Cryst%gmet,npw_gamma,kg_gamma) 1483 sh_istwfk = istwf1; sh_nkibz = 1; sh_kibz(:,1) = gamma_point 1484 1485 ! Compute number of elements in the Shirley basis set and initialize the descriptor oWsh 1486 ! TODO: BE careful in parallel when nsppol==2. I should recreate the communicators. 1487 ABI_MALLOC(sh_nband,(sh_nkibz,nsppol)) 1488 do spin=1,nsppol 1489 sh_nband(:,spin)=sh_size(spin) 1490 end do 1491 sh_mband=MAXVAL(sh_nband) 1492 1493 ABI_MALLOC(sh_bks_mask,(sh_mband,sh_nkibz,nsppol)) 1494 ABI_MALLOC(sh_keep_ur ,(sh_mband,sh_nkibz,nsppol)) 1495 ! Set sh_keep_ur to False. as sh_u(r) are only needed for the computation of vloc so that we can save memory. 1496 !sh_bks_mask=.TRUE.; sh_keep_ur =.False. 1497 sh_bks_mask=.TRUE.; sh_keep_ur =.TRUE. 1498 1499 call wfd_init(oWsh,Cryst,Pawtab,Psps,sh_keep_ur,Wfd%paral_kgb,npw_gamma,sh_mband,sh_nband,sh_nkibz,nsppol,& 1500 & sh_bks_mask,Wfd%nspden,nspinor,Wfd%ecutsm,Wfd%dilatmx,sh_istwfk,sh_kibz,Wfd%ngfft,kg_gamma,Wfd%nloalg,& 1501 & Wfd%prtvol,Wfd%pawprtvol,Wfd%comm) 1502 1503 if (oWsh%prtvol > 0) then 1504 call wfd_print(oWsh,header="Shirley wavefunction descriptor") 1505 end if 1506 1507 ABI_FREE(sh_keep_ur) 1508 ABI_FREE(sh_bks_mask) 1509 ABI_FREE(sh_nband) 1510 1511 !DEBUG 1512 ! call wfd_change_ngfft(Wfd,Cryst,Psps,ov_ngfft) 1513 ! call wfd_change_ngfft(oWsh,Cryst,Psps,ov_ngfft) 1514 !DEBUG 1515 ! 1516 ! ===================================================================== 1517 ! ==== Rotate the input wavefunctions to get the optimal basis set ==== 1518 ! ===================================================================== 1519 ABI_CHECK(Wfd%nfft == oWsh%nfft, "Different nfft for Wfd and Wsh!") 1520 fft_fact = one/Wfd%nfft 1521 ABI_MALLOC(ur1,(Wfd%nfft*nspinor)) 1522 1523 if (.not.fft_check_rotrans(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,fft_err)) then 1524 write(msg,'(a,3(i0,1x),a)')" Real space FFT mesh ",Wfd%ngfft(1:3)," is not symmetric. Cannot symmetrize in real space" 1525 MSG_ERROR(msg) 1526 end if 1527 1528 do spin=1,nsppol 1529 if (.not. wfd_itreat_spin(Wfd,spin,comm_spin,rank_spin,nproc_spin)) cycle 1530 write(msg,'(a,f12.1,a)')' Memory needed for Shirley u(r) = ',two*gwpc*Wfd%nfft*nspinor*sh_size(spin)*b2Mb,' [Mb]' 1531 call wrtout(std_out,msg,'PERS') 1532 1533 ! Allocate space for Shirley wavefunctions in real space. 1534 ABI_MALLOC(sh_ur, (Wfd%nfft*nspinor,sh_size(spin))) 1535 sh_ur = czero 1536 1537 ! MPI loop over the index of the overlap operator. 1538 call xmpi_split_work(O(spin)%size,comm_spin,my_start,my_stop,msg,ierr) 1539 1540 do midx=my_start,my_stop 1541 ! Retrieve the k-point index in the BZ and the band 1542 band1 = O(spin)%idx2bk(1,midx) 1543 ik1_bz = O(spin)%idx2bk(2,midx) 1544 1545 ! Symmetrize the wavefunction in real space. 1546 call wfd_sym_ur(Wfd,Cryst,Kmesh,band1,ik1_bz,spin,ur1) 1547 1548 ! Construct the new optimal basis set. 1549 do sh=1,sh_size(spin) 1550 sh_ur(:,sh) = sh_ur(:,sh) + O(spin)%mat(midx,base_spin(spin)+sh) * ur1 1551 end do 1552 end do 1553 1554 ! Gather results 1555 call xmpi_sum(sh_ur,comm_spin,ierr) 1556 1557 ! NC pseudos: Normalize the basis set using the eigenvalues of the overlap matrix. 1558 do sh=1,sh_size(spin) 1559 sh_ur(:,sh) = sh_ur(:,sh)/SQRT(O(spin)%eigene(base_spin(spin)+sh)) 1560 !norm1 = xdotc(Wfd%nfft*nspinor,sh_ur(:,sh),1,sh_ur(:,sh),1) * fft_fact 1561 !sh_ur(:,sh) = sh_ur(:,sh)/SQRT(norm1) 1562 !write(std_out,*)" sh_ur integrates to: ", xdotc(Wfd%nfft*nspinor,sh_ur(:,sh),1,sh_ur(:,sh),1) * fft_fact 1563 end do 1564 ! 1565 ! Construct new optimal basis set in G-space and save results in oWsh. 1566 ABI_MALLOC(sh_ug,(npw_gamma*nspinor)) 1567 1568 do sh=1,sh_size(spin) 1569 ! Transform sh_u(r) from the FFT mesh to the G-sphere. 1570 ! Don't try to update u(r) and PAW cpcrj since they are not used here. 1571 1572 ! gbound_k => oWsh%Kdata(k1)%gbound 1573 call fft_ur(npw_gamma,oWsh%nfft,nspinor,ndat1,oWsh%mgfft,oWsh%ngfft,istwf1,& 1574 & kg_gamma,oWsh%Kdata(k1)%gbound,sh_ur(:,sh),sh_ug) 1575 ! NC: Normalize the basis set using the eigenvalues of the overlap matrix. 1576 !if (Wfd%usepaw==0) sh_ug = sh_ug/SQRT(O(spin)%eigene(base_spin(spin)+sh)) 1577 1578 ! Store Shirley u(G). 1579 call wfd_push_ug(oWsh,sh,k1,spin,Cryst,sh_ug,update_ur=.FALSE.,update_cprj=.FALSE.) 1580 end do 1581 ! 1582 ! PAW: Normalize the basis set using the eigenvalues of the overlap matrix. 1583 !if (.FALSE. .and. Wfd%usepaw==1) then 1584 if (Wfd%usepaw==1) then 1585 call wfd_reset_ur_cprj(oWsh) 1586 ABI_DT_MALLOC(Cp_sh1,(natom,nspinor)) 1587 call pawcprj_alloc(Cp_sh1,0,Wfd%nlmn_atm) 1588 1589 do sh=1,sh_size(spin) 1590 ! TODO: Write function to compute the scalar product at the same k, taking into account istwfk. 1591 !ug1 => oWsh%Wave(sh,k1,spin)%ug 1592 !cdum = xdotc(npw_gamma*nspinor,ug1,1,ug1,1) * Cryst%ucvol 1593 cdum = xdotc(npw_gamma*nspinor,oWsh%Wave(sh,k1,spin)%ug,1,oWsh%Wave(sh,k1,spin)%ug,1) 1594 !if (istwf_k>1) then 1595 ! cdum=two*DBLE(cdum) 1596 ! if (istwf_k==2) cdum=cdum-CONJG(ug1(1))*ug1(1) 1597 !end if 1598 1599 call wfd_get_cprj(oWsh,sh,k1,spin,Cryst,Cp_sh1,sorted=.FALSE.) 1600 1601 pawovlp = paw_overlap(Cp_sh1,Cp_sh1,Cryst%typat,Pawtab) 1602 !pawovlp = pawovlp * Cryst%ucvol 1603 cdum = cdum + CMPLX(pawovlp(1),pawovlp(2)) 1604 1605 !norm1 = xdotc(Wfd%nfft*nspinor,sh_ur(:,sh),1,sh_ur(:,sh),1) * fft_fact * Cryst%ucvol 1606 !sh_ur(:,sh) = sh_ur(:,sh)/SQRT(O(spin)%eigene(base_spin(spin)+sh)) 1607 !write(std_out,*)" sh_ur integrates to: ", xdotc(Wfd%nfft*nspinor,sh_ur(:,sh),1,sh_ur(:,sh),1) * fft_fact 1608 !write(std_out,*)" PAW test: ",DBLE(cdum),O(spin)%eigene(base_spin(spin)+sh) 1609 1610 sh_ug = oWsh%Wave(sh,k1,spin)%ug/SQRT(DBLE(cdum)) 1611 call wfd_push_ug(oWsh,sh,k1,spin,Cryst,sh_ug,update_ur=.FALSE.,update_cprj=.FALSE.) 1612 end do 1613 1614 call pawcprj_free(Cp_sh1) 1615 ABI_DT_FREE(Cp_sh1) 1616 end if 1617 1618 ABI_FREE(sh_ur) 1619 ABI_FREE(sh_ug) 1620 ! 1621 ! Free the overlap eigenvectors. 1622 call ovlp_free(O(spin)) 1623 end do ! spin 1624 1625 do spin=1,nsppol ! Just to be on the safe side. 1626 call ovlp_free(O(spin)) 1627 end do 1628 1629 ABI_FREE(kg_gamma) 1630 ABI_FREE(ur1) 1631 1632 ! Set the MPI communicators. 1633 call wfd_set_mpicomm(oWsh) 1634 1635 call cwtime(cpu,wall,gflops,"stop") 1636 write(std_out,*)"SHTIME: Rotation, cpu_time: ",cpu,", wall_time: ",wall 1637 1638 DBG_EXIT("COLL") 1639 1640 end subroutine wfd_bloch_to_shirley
m_shirley/wfd_shirley_to_eh [ Functions ]
[ Top ] [ m_shirley ] [ Functions ]
NAME
wfd_shirley_to_eh
FUNCTION
Return a new wavefunction descriptor containing the basis set for the e-h manifold. INPUT Wfd<wfd_t> Cryst<crystal_t>= data type gathering info on symmetries and unit cell Psps<pseudopotential_type)>=variables related to pseudopotentials Pawtab(Psps%ntypat)<pawtab_type>=paw tabulated starting data Pawang <pawang_type>=angular mesh discretization and related data: eh_coverage
OUTPUT
Weh<wfd_t>
PARENTS
CHILDREN
blas_cholesky_ortho,cwtime,fft_ur,fftbox_execute,fftbox_plan3_many flush_unit,get_kg,kgindex,ovlp_diago_and_prune,ovlp_free wfd_change_ngfft,wfd_get_ur,wfd_init,wfd_print,wfd_push_ug wfd_test_ortho,wrtout,xgemm
SOURCE
2802 subroutine wfd_shirley_to_eh(Wsh,Cryst,Psps,Pawtab,Pawang,Pawrad,min_bsize,eh_coverage,Weh,sh2eh) 2803 2804 2805 !This section has been created automatically by the script Abilint (TD). 2806 !Do not modify the following lines by hand. 2807 #undef ABI_FUNC 2808 #define ABI_FUNC 'wfd_shirley_to_eh' 2809 use interfaces_14_hidewrite 2810 use interfaces_53_ffts 2811 !End of the abilint section 2812 2813 implicit none 2814 2815 !Arguments ------------------------------------ 2816 !scalars 2817 integer,intent(in) :: min_bsize 2818 real(dp),intent(in) :: eh_coverage 2819 type(crystal_t),intent(in) :: Cryst 2820 type(Pawang_type),intent(in) :: Pawang 2821 type(Pseudopotential_type),intent(in) :: Psps 2822 type(wfd_t),target,intent(inout) :: Wsh 2823 type(wfd_t),target,intent(out) :: Weh 2824 !arrays 2825 !integer,intent(in) :: ov_ngfft(18) 2826 complex(gwpc),pointer :: sh2eh(:,:) 2827 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wsh%usepaw) 2828 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*Wsh%usepaw) 2829 2830 !Local variables ------------------------------ 2831 !scalars 2832 integer,parameter :: istwf1=1,k1=1,s1=1,nkpt1=1,ndat1=1 2833 integer :: ij,band1,band2,natom,eh,midx,ig,mband,nstates,ovlp_size 2834 integer :: npw_gamma,usepaw !,spin !,row,col,useylm_,!,ierr 2835 integer :: fft_idx,nspinor,eh_nkibz,eh_mband,npw,col,band3,band4,nsppol 2836 real(dp) :: fft_fact,norm1 !sqrt_norm1,sqrt_norm2 2837 real(dp) :: cpu,wall,gflops 2838 character(len=500) :: msg 2839 type(fftbox_plan3_t) :: plan 2840 type(ovlp_t) :: Oeh 2841 !arrays 2842 integer :: ov_ngfft(18) 2843 integer :: eh_istwfk(1),eh_size,base 2844 integer,allocatable :: kg_gamma(:,:) !,gbound_k(:,:),kg_k(:,:), 2845 integer,allocatable :: eh_nband(:,:) 2846 integer,pointer :: igfft0(:) 2847 real(dp) :: eh_kibz(3,1),gamma_point(3)=(/zero,zero,zero/) 2848 !real(dp) :: pawovlp(2) 2849 complex(dpc) :: cfft_fact 2850 complex(dpc),allocatable :: dpc_tmp(:) 2851 complex(gwpc),allocatable :: ur1(:),ur2(:),ur12(:),cf_ovlp(:,:) 2852 complex(gwpc),target,allocatable :: ur3(:),ur4(:) !ur1(:),ur2(:), 2853 !complex(gwpc),pointer :: ug1(:),pt_ur1(:),pt_ur2(:) 2854 complex(gwpc),pointer :: pt_ur3(:),pt_ur4(:) 2855 complex(dpc),allocatable :: ur34_big(:,:) 2856 complex(gwpc),allocatable :: ehg(:),eh_ur(:,:),eh_ug(:,:),gwpc_tmp(:) 2857 logical,allocatable :: eh_keep_ur(:,:,:),eh_bks_mask(:,:,:),kg_mask(:) 2858 2859 !************************************************************************ 2860 2861 DBG_ENTER("COLL") 2862 2863 ABI_CHECK(Wsh%nspinor==1,"nspinor==2 not coded.") 2864 ABI_CHECK(Wsh%paral_kgb==0,"paral_kgb/=0 not coded.") 2865 ABI_CHECK(Wsh%rfft_is_symok,"Real space FFT is not symmetric.") 2866 ABI_CHECK(Wsh%nsppol==1,"nsppol must be 1!") 2867 2868 ABI_UNUSED(Pawang%l_max) 2869 ABI_UNUSED(Pawrad(1)%mesh_size) 2870 2871 nspinor = Wsh%nspinor 2872 usepaw = Wsh%usepaw 2873 natom = Wsh%natom 2874 nsppol = Wsh%nsppol 2875 mband = Wsh%mband 2876 nstates = Wsh%nband(k1,s1) 2877 ABI_CHECK(xmpi_comm_size(Wsh%comm)==1,"ovlp2_init is not parallelized") 2878 2879 write(msg,'(a,f9.6)')" Transformation Shirley --> e-h basis set with eh_coverage : ",eh_coverage 2880 call wrtout(std_out,msg,"COLL") 2881 2882 ! 1) Get the overlap matrix <S_i^* S_j|S_k^* S_l> for this spin. 2883 ! The Overlap is calculated using the coarse real space FFT ov_ngfft 2884 ov_ngfft = Wsh%ngfft 2885 call wfd_change_ngfft(Wsh,Cryst,Psps,ov_ngfft) 2886 2887 ovlp_size = nstates**2 2888 Oeh%mband = mband 2889 Oeh%nkpt = nkpt1 2890 Oeh%min_ene = smallest_real 2891 Oeh%max_ene = greatest_real 2892 ! 2893 ! The size of the overlap matrix and useful tables. 2894 ABI_MALLOC(Oeh%bk2idx,(mband,nkpt1)) 2895 Oeh%bk2idx=0 2896 2897 Oeh%size = ovlp_size 2898 ABI_MALLOC(Oeh%idx2bk,(2,ovlp_size)) 2899 Oeh%idx2bk = 0 2900 ! 2901 ! Allocate overlap matrix. Could use packed matrix to save memory, but Lapack call is slower. 2902 write(msg,'(a,f8.2,a)')" out of memory in Oeh%mat, requiring :",two*dpc*ovlp_size**2*b2Gb," Gb" 2903 ABI_MALLOC(Oeh%mat, (ovlp_size,ovlp_size)) 2904 2905 !Oeh%mat = -HUGE(one) 2906 write(msg,'(a,f12.1,a,i0)')" Memory required for the overlap matrix: ",two*dpc*ovlp_size**2*b2Mb," Mb; Matrix size= ",ovlp_size 2907 call wrtout(std_out,msg,"COLL") 2908 call flush_unit(std_out) 2909 ! 2910 ! Calculate the overlap matrix -------------------------------------------------------------------------- 2911 ! 1) Symmetrization in k-space is not needed here 2912 ! 2) Matrix is Hermitian. 2913 ! 2914 fft_fact = one/Wsh%nfft 2915 ABI_MALLOC(ur3,(Wsh%nfft*nspinor)) 2916 ABI_MALLOC(ur4,(Wsh%nfft*nspinor)) 2917 2918 ! TODO Temporary implementation used to speed up this part. 2919 ABI_MALLOC(ur34_big,(Wsh%nfft,nstates**2)) 2920 2921 Oeh%mat_type = TYPE_OVERLAP 2922 npw = Wsh%npwarr(k1) 2923 2924 call cwtime(cpu,wall,gflops,"start") 2925 col = 0 2926 do band4=1,nstates 2927 ! 2928 if (wfd_ihave_ur(Wsh,band4,k1,s1,how="Stored")) then 2929 pt_ur4 => Wsh%Wave(band4,k1,s1)%ur 2930 else 2931 call wfd_get_ur(Wsh,band4,k1,s1,ur4) 2932 pt_ur4 => ur4 2933 end if 2934 ! 2935 do band3=1,nstates 2936 col = col+1 2937 if (wfd_ihave_ur(Wsh,band3,k1,s1,how="Stored")) then 2938 pt_ur3 => Wsh%Wave(band3,k1,s1)%ur 2939 else 2940 call wfd_get_ur(Wsh,band3,k1,s1,ur3) 2941 pt_ur3 => ur3 2942 end if 2943 !if (Wsh%usepaw==1) then 2944 !call wfd_get_cprj(Wsh,band3,k1,s1,Cryst,Cp_k3,sorted=.FALSE.) 2945 !end if 2946 ur34_big(:,col) = CONJG(pt_ur3) * pt_ur4 2947 Oeh%idx2bk(1,col) = band3 2948 Oeh%idx2bk(2,col) = band4 2949 end do 2950 end do 2951 2952 cfft_fact = cone/Wsh%nfft 2953 !Oeh%mat = cfft_fact * MATMUL( CONJG(TRANSPOSE(ur34_big)), ur34_big) 2954 call xgemm("C","N",nstates**2,nstates**2,Wsh%nfft,cfft_fact,ur34_big,Wsh%nfft,ur34_big,Wsh%nfft,czero,Oeh%mat,nstates**2) 2955 2956 !call print_arr(Oeh%mat,max_r=10,max_c=12,unit=std_out,mode_paral="COLL") 2957 2958 ABI_FREE(ur3) 2959 ABI_FREE(ur4) 2960 2961 call cwtime(cpu,wall,gflops,"stop") 2962 write(std_out,*)"SHTIME: Ovlp2_build, cpu_time: ",cpu,", wall_time: ",wall 2963 call flush_unit(std_out) 2964 2965 ! 2) Diagonalize the overlap matrix selecting the optimal subspace: [ base(spin):ovlp_size ] 2966 call ovlp_diago_and_prune(Oeh,eh_coverage,eh_size,base) ! In exit Oeh%mat stores the eigenvectors. 2967 ! 2968 ! Make sure we have enough states. 2969 if (eh_size < min_bsize) then 2970 if (Oeh%size<min_bsize) then 2971 write(msg,'(2(a,i0),2a)')& 2972 & " Overlap size is ",Oeh%size," whereas min_bsize is ",min_bsize,ch10,& 2973 & " Decrease the number of bands to be interpolated or increase the number of ab-initio input states." 2974 MSG_ERROR(msg) 2975 end if 2976 eh_size = min_bsize 2977 write(msg,'(a,2i0)')" Had to enlarge Shirley subspace since input eh_size < min_bsize: ",eh_size,min_bsize 2978 MSG_COMMENT(msg) 2979 end if 2980 2981 call cwtime(cpu,wall,gflops,"start") 2982 ! 2983 ! 3) Init a new wavefunction descriptor to store the optimal basis set. 2984 ! *) Weh must be allocated here since eh_size is know only after the diagonalization of the overlap. 2985 ! *) Keep the optimal wavefunctions on each node (if possible) to facilitate the interpolation over the fine k-mesh. 2986 ! *) Use Gamma-centered basis set to facilitate the operations in reciprocal space. 2987 ! *) The new basis is orthogonal, but not normalized since <U_i|U_j> = delta_ij e_i. 2988 ! 2989 ! The optimal basis set is given on the gamma centered basis set with istwfk==1. 2990 ! FIXME temporary hacking. There is a bug somewhere in kpgsph 2991 call get_kg(gamma_point,istwf1,Wsh%ecut,Cryst%gmet,npw_gamma,kg_gamma) 2992 !call get_kg(gamma_point,istwf1,14.0_dp,Cryst%gmet,npw_gamma,kg_gamma) 2993 ! 2994 ! * Index of the G-sphere in the FFT box. 2995 ABI_MALLOC(igfft0,(npw_gamma)) 2996 ABI_MALLOC(kg_mask,(npw_gamma)) 2997 call kgindex(igfft0,kg_gamma,kg_mask,Wsh%MPI_enreg,Wsh%ngfft,npw_gamma) 2998 2999 ABI_CHECK(ALL(kg_mask),"FFT para not yet implemented") 3000 ABI_FREE(kg_mask) 3001 3002 eh_istwfk=istwf1; eh_nkibz= 1 3003 eh_kibz(:,1) = gamma_point 3004 3005 ! TODO: BE careful in parallel when nsppol==2. I should recreate the communicators. 3006 ABI_MALLOC(eh_nband,(eh_nkibz,nsppol)) 3007 eh_nband = eh_size 3008 !call wfd_change_ngfft(Wsh,Cryst,Psps,ov_ngfft) 3009 3010 eh_mband=MAXVAL(eh_nband) 3011 ABI_MALLOC(eh_bks_mask,(eh_mband,eh_nkibz,nsppol)) 3012 ABI_MALLOC(eh_keep_ur ,(eh_mband,eh_nkibz,nsppol)) 3013 eh_bks_mask=.TRUE.; eh_keep_ur =.TRUE. 3014 3015 call wfd_init(Weh,Cryst,Pawtab,Psps,eh_keep_ur,Wsh%paral_kgb,npw_gamma,eh_mband,eh_nband,eh_nkibz,nsppol,& 3016 & eh_bks_mask,Wsh%nspden,nspinor,Wsh%ecutsm,Wsh%dilatmx,eh_istwfk,eh_kibz,Wsh%ngfft,kg_gamma,Wsh%nloalg,& 3017 & Wsh%prtvol,Wsh%pawprtvol,Wsh%comm) 3018 3019 if (Weh%prtvol > 0) then 3020 call wfd_print(Weh,header="Shirley wavefunction descriptor") 3021 end if 3022 3023 ABI_FREE(eh_keep_ur) 3024 ABI_FREE(eh_bks_mask) 3025 ABI_FREE(eh_nband) 3026 ! 3027 ! ===================================================================== 3028 ! ==== Rotate the input wavefunctions to get the optimal basis set ==== 3029 ! ===================================================================== 3030 3031 fft_fact = one/Wsh%nfft 3032 ABI_MALLOC(ur1,(Wsh%nfft*nspinor)) 3033 ABI_MALLOC(ur2,(Wsh%nfft*nspinor)) 3034 ABI_MALLOC(ur12,(Wsh%nfft*nspinor)) 3035 ! 3036 write(msg,'(a,f12.1,a)')' Memory needed for storing eh_ur= ',two*gwpc*Wsh%nfft*nspinor*eh_size*b2Mb,' [Mb]' 3037 call wrtout(std_out,msg,'PERS') 3038 3039 ABI_MALLOC(eh_ur,(Wsh%nfft*nspinor,eh_size)) 3040 eh_ur = czero 3041 3042 ABI_MALLOC(eh_ug,(npw_gamma*nspinor,eh_size)) 3043 eh_ug = czero 3044 3045 do midx=1,ovlp_size ! Loop over the single particle orbitals. 3046 band1 = Oeh%idx2bk(1,midx) ! TODO to be removed. 3047 band2 = Oeh%idx2bk(2,midx) 3048 3049 call wfd_get_ur(Wsh,band1,k1,s1,ur1) 3050 call wfd_get_ur(Wsh,band2,k1,s1,ur2) 3051 ur12 = CONJG(ur1) * ur2 ! Use same convention as the one used in ovlp2_init. 3052 ! 3053 ! Construct the new optimal basis set. 3054 do eh=1,eh_size 3055 eh_ur(:,eh) = eh_ur(:,eh) + Oeh%mat(midx,base+eh) * ur12 3056 end do 3057 ! 3058 end do 3059 ! 3060 ! NC: Normalize the basis set. 3061 do eh=1,eh_size 3062 norm1 = xdotc(Weh%nfft*nspinor,eh_ur(:,eh),1,eh_ur(:,eh),1) * fft_fact 3063 eh_ur(:,eh) = eh_ur(:,eh)/SQRT(norm1) 3064 !write(std_out,*)" eh_ur integrates to: ", xdotc(Weh%nfft*nspinor,eh_ur(:,eh),1,eh_ur(:,eh),1) * fft_fact 3065 end do 3066 ! 3067 ! From the FFT mesh to the G-sphere. 3068 ABI_MALLOC(ehg,(npw_gamma*nspinor)) 3069 ABI_MALLOC(gwpc_tmp,(Wsh%nfft*nspinor)) 3070 ! 3071 ! ============================================================================ 3072 ! ==== Construct new optimal basis set in G-space and save results in Weh ==== 3073 ! ============================================================================ 3074 do eh=1,eh_size 3075 gwpc_tmp = eh_ur(:,eh) 3076 3077 #if 1 3078 ABI_MALLOC(dpc_tmp,(Wsh%nfft*nspinor)) 3079 dpc_tmp = eh_ur(:,eh) 3080 3081 call fftbox_plan3_many(plan,Wsh%nspinor*ndat1,Wsh%ngfft(1:3),Wsh%ngfft(1:3),Wsh%ngfft(7),-1) 3082 call fftbox_execute(plan,dpc_tmp) 3083 ! 3084 do ig=1,npw_gamma 3085 fft_idx = igfft0(ig) 3086 if (fft_idx/=0) then ! G-G0 belong to the FFT mesh. 3087 if (fft_idx>Wsh%nfft .or.fft_idx<0) then 3088 MSG_ERROR("fft_idx bug") 3089 end if 3090 ehg(ig)=dpc_tmp(fft_idx) 3091 else ! Set this component to zero. 3092 MSG_ERROR("fft_idx bug") 3093 ehg(ig)=czero 3094 end if 3095 end do 3096 ABI_FREE(dpc_tmp) 3097 #else 3098 ! FIXME does not work anymore. 3099 ! TODO Check with the new version 3100 gbound_k => Weh%Kdata(k1)%gbound 3101 call fft_ur(npw_gamma,Wsh%nfft,nspinor,ndat1,Wsh%mgfft,Wsh%ngfft,istwf1,kg_gamma,gbound_k,gwpc_tmp,ehg) 3102 #endif 3103 ! NC: Normalize the basis set using the eigenvalues of the overlap matrix. 3104 !if (usepaw==0) then 3105 ! ehg = ehg/SQRT(Oeh%eigene(base+eh)) 3106 !end if 3107 !%call wfd_push_ug(Weh,eh,k1,s1,Cryst,ehg,update_ur=.FALSE.,update_cprj=.FALSE.) 3108 eh_ug(:,eh) = ehg 3109 end do 3110 ! 3111 ! ====================================== 3112 ! ==== Orthonormalize the basis set ==== 3113 ! ====================================== 3114 ABI_MALLOC(cf_ovlp,(eh_size,eh_size)) 3115 3116 call blas_cholesky_ortho(npw_gamma,eh_size,eh_ug,cf_ovlp) 3117 3118 ABI_FREE(cf_ovlp) 3119 ! 3120 ! Push data. 3121 do eh=1,eh_size 3122 call wfd_push_ug(Weh,eh,k1,s1,Cryst,eh_ug(:,eh),update_ur=.FALSE.,update_cprj=.FALSE.) 3123 end do 3124 3125 ABI_FREE(gwpc_tmp) 3126 ABI_FREE(ehg) 3127 ABI_FREE(eh_ur) 3128 ABI_FREE(eh_ug) 3129 ABI_FREE(kg_gamma) 3130 ABI_FREE(igfft0) 3131 ABI_FREE(ur1) 3132 ABI_FREE(ur2) 3133 ABI_FREE(ur12) 3134 3135 call cwtime(cpu,wall,gflops,"stop") 3136 write(std_out,*)"SHTIME: EH_Rotation, cpu_time: ",cpu,", wall_time:",wall 3137 ! 3138 ! ================================== 3139 ! ==== sh2eh = <S_i* S_j| EH_a> ==== 3140 ! ================================== 3141 call cwtime(cpu,wall,gflops,"start") 3142 ABI_MALLOC(sh2eh,(ovlp_size,eh_size)) 3143 sh2eh=czero 3144 3145 call wfd_change_ngfft(Weh,Cryst,Psps,Wsh%ngfft) ! Make sure the two set of wave are on the same mesh. 3146 ABI_MALLOC(ur1,(Weh%nfft*nspinor)) 3147 3148 do eh=1,eh_size 3149 call wfd_get_ur(Weh,eh,k1,s1,ur1) 3150 do ij=1,ovlp_size 3151 sh2eh(ij,eh) = DOT_PRODUCT(ur34_big(:,ij),ur1) 3152 end do 3153 end do 3154 ABI_FREE(ur1) 3155 3156 ABI_FREE(ur34_big) 3157 3158 call cwtime(cpu,wall,gflops,"stop") 3159 write(std_out,*)"SHTIME: EH_Projection, cpu_time: ",cpu,",wall_time:",wall 3160 ! 3161 if (DEBUGME) then 3162 call wfd_test_ortho(Weh,Cryst,Pawtab,unit=std_out,mode_paral="COLL") 3163 end if 3164 ! 3165 ! Deallocate memory. 3166 call ovlp_free(Oeh) 3167 3168 DBG_EXIT("COLL") 3169 3170 end subroutine wfd_shirley_to_eh