TABLE OF CONTENTS
- ABINIT/m_hexc
- m_haydock/hexc_t
- m_hexc/hexc_build_hinterp
- m_hexc/hexc_compute_hinterp
- m_hexc/hexc_compute_subhinterp
- m_hexc/hexc_free
- m_hexc/hexc_init
- m_hexc/hexc_interp_free
- m_hexc/hexc_interp_init
- m_hexc/hexc_interp_matmul
- m_hexc/hexc_interp_t
- m_hexc/hexc_matmul_elphon
- m_hexc/hexc_matmul_full
- m_hexc/hexc_matmul_tda
ABINIT/m_hexc [ Modules ]
NAME
m_hexc
FUNCTION
module for excitonic hamiltonian for Haydock
COPYRIGHT
Copyright (C) 2014-2024 ABINIT group (M.Giantomassi, Y. Gillet) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_hexc 23 24 use defs_basis 25 use m_abicore 26 use m_bs_defs 27 use m_xmpi 28 use m_errors 29 use m_nctk 30 use m_haydock_io 31 use m_linalg_interfaces 32 use netcdf 33 34 use m_time, only : timab 35 use m_fstrings, only : indent, strcat, sjoin, itoa 36 use defs_datatypes, only : ebands_t, pseudopotential_type 37 use m_hide_blas, only : xdotc, xgemv 38 use m_numeric_tools, only : print_arr, symmetrize, hermitianize, wrap2_pmhalf 39 use m_crystal, only : crystal_t 40 use m_bz_mesh, only : kmesh_t, findqg0 41 use m_double_grid, only : double_grid_t, get_kpt_from_indices_coarse, compute_corresp 42 use m_wfd, only : wfdgw_t 43 use m_bse_io, only : exc_read_rcblock, exc_write_optme, exc_ham_ncwrite 44 use m_pawtab, only : pawtab_type 45 use m_vcoul, only : vcoul_t 46 use m_bseinterp, only : interpolator_t, interpolator_init, interpolator_normalize, & 47 & interpolator_free, int_alloc_work, int_free_work 48 49 implicit none 50 51 private
m_haydock/hexc_t [ Types ]
[ Top ] [ m_haydock ] [ Types ]
NAME
hexc_t
FUNCTION
Store the excitonic hamiltonian and other related information
SOURCE
63 type,public :: hexc_t 64 65 !scalars 66 integer :: comm 67 ! MPI communicator 68 69 integer :: hsize_coarse 70 ! Size of the coarse hamiltonian 71 72 integer :: hsize 73 ! Size of the hamiltonian and the kets 74 ! (= hsize_coarse without interpolation, = hsize_dense with interpolation) 75 76 integer :: nbz 77 ! Number of kpoints for the full problem 78 ! (= nbz_coarse without interpolation, = nbz_dense with interpolation) 79 80 integer :: my_t1 81 ! Lower limit of MPI paral 82 83 integer :: my_t2 84 ! Upper limit of MPI paral 85 86 integer :: my_nt 87 ! Number of transitions treat by node 88 ! = my_t2 - my_t1 + 1 89 90 integer :: nbnd_coarse 91 ! Product of number of bands conduction X valence 92 93 ! Pointers to data that are already in memory 94 type(excparam),pointer :: bsp => null() 95 ! parameters for BS 96 97 type(excfiles),pointer :: bs_files => null() 98 ! files for BSE 99 100 type(crystal_t),pointer :: crystal => null() 101 ! crystal info 102 103 type(kmesh_t),pointer :: kmesh_coarse => null() 104 ! kmesh of the coarse mesh 105 106 type(kmesh_t),pointer :: kmesh => null() 107 ! kmesh of the full problem 108 109 type(ebands_t),pointer :: ks_bst => null() 110 type(ebands_t),pointer :: qp_bst => null() 111 ! band structures of the full problem 112 113 type(wfdgw_t),pointer :: wfd_coarse => null() 114 ! Wfd of the coarse problem 115 116 type(wfdgw_t),pointer :: wfd => null() 117 ! wfd of the full problem 118 119 !arrays 120 complex(dpc),allocatable :: hreso(:,:) 121 ! Resonant part of the hamiltonian 122 123 complex(dpc),allocatable :: hcoup(:,:) 124 ! Coupling part of the hamiltonian 125 126 complex(dpc),allocatable :: diag_coarse(:) 127 ! Diagonal part of the hamiltonian with transition energies 128 129 end type hexc_t
m_hexc/hexc_build_hinterp [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_build_hinterp
FUNCTION
Pre-compute interpolated hamiltonian and store it in memory
INPUTS
OUTPUT
SIDE EFFECTS
hexc, hexc_i Pre-compute info to save CPU time when computing matmul
SOURCE
564 subroutine hexc_build_hinterp(hexc,hexc_i) 565 566 !Arguments --------------------------- 567 type(hexc_t),intent(inout) :: hexc 568 type(hexc_interp_t),intent(inout) :: hexc_i 569 570 !Local variables --------------------- 571 integer :: ierr,ncerr,ncid 572 character(len=500) :: msg 573 574 !***************************************************************************** 575 576 write(msg,"(a,f8.1,a)")"Memory needed for hinterp = ",one*(hexc_i%hsize_dense**2)*2*dpc*b2Mb," Mb" 577 call wrtout(std_out,msg,"COLL") 578 579 ABI_MALLOC_OR_DIE(hexc_i%hinterp,(hexc_i%hsize_dense,hexc_i%hsize_dense), ierr) 580 581 call hexc_compute_hinterp(hexc%BSp, hexc%hsize_coarse, hexc_i%hsize_dense, hexc_i%all_hmat, & 582 & hexc_i%interpolator%double_grid,hexc%nbnd_coarse, hexc_i%interpolator, & 583 & hexc_i%kdense2div, hexc_i%all_acoeffs,hexc_i%all_bcoeffs, hexc_i%all_ccoeffs, & 584 & hexc_i%Kmesh_dense, hexc_i%Vcp_dense, hexc%crystal%gmet, hexc_i%hinterp, hexc_i%m3_width) 585 586 ABI_SFREE(hexc_i%all_acoeffs) 587 ABI_SFREE(hexc_i%all_bcoeffs) 588 ABI_SFREE(hexc_i%all_ccoeffs) 589 590 591 if (hexc%BSp%prt_ncham) then 592 ABI_COMMENT("Printing HEXC_I.nc file") 593 ncerr = nctk_open_create(ncid, trim(hexc%BS_files%out_basename)//"_HEXC_I.nc", xmpi_comm_self) 594 NCF_CHECK_MSG(ncerr, "Creating HEXC_I file") 595 call exc_ham_ncwrite(ncid, hexc_i%Kmesh_dense, hexc%BSp, hexc_i%hsize_dense, hexc%BSp%nreh_interp, & 596 & hexc%BSp%vcks2t_interp, hexc_i%hinterp, hexc_i%diag_dense) 597 NCF_CHECK(nf90_close(ncid)) 598 end if 599 600 end subroutine hexc_build_hinterp
m_hexc/hexc_compute_hinterp [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_compute_hinterp
FUNCTION
Compute interpolated matrix elements for methods 2 and 3
INPUTS
BSp<type(excparam)=The parameter for the Bethe-Salpeter run. hsize_coarse=Size of the coarse Hamiltonian hsize_dense=Size of the dense Hamiltonian hmat(hsize_coarse,hsize_coarse,8)=Excitonic matrix grid<double_grid_t> = Correspondence between coarse and dense k-mesh. nbnd_coarse = Total number of bands interpolator<interpolator_t> = Interpolator kdense2div = Mapping kdense2div acoeffs, bcoeffs, ccoeffs = decomposition "W = a/q^2 + b/q + c" Kmesh_dense<type(kmesh_t)>=The list of k-points in the BZ, IBZ and symmetry tables. Vcp_dense<vcoul_t>=Coulomb interation in G-space on the dense Q-mesh gmet(3,3)=Metric tensor in G-space
OUTPUT
hinterp = Interpolated hamiltonian
SOURCE
749 subroutine hexc_compute_hinterp(BSp,hsize_coarse,hsize_dense,hmat,grid,nbnd_coarse,& 750 & interpolator,kdense2div,acoeffs,bcoeffs,ccoeffs,Kmesh_dense,Vcp_dense,gmet,hinterp,& 751 & m3_width) 752 753 !Arguments ------------------------------------ 754 !scalars 755 integer,intent(in) :: hsize_coarse,hsize_dense,nbnd_coarse !,ntrans 756 real(dp),intent(in) :: m3_width 757 type(excparam),intent(in) :: BSp 758 type(double_grid_t),intent(in) :: grid 759 type(vcoul_t),intent(in) :: Vcp_dense 760 type(kmesh_t),intent(in) :: Kmesh_dense 761 type(interpolator_t),target,intent(inout) :: interpolator 762 !arrays 763 integer,intent(in) :: kdense2div(grid%nbz_dense) 764 real(dp),intent(in) :: gmet(3,3) 765 complex(dpc),intent(in) :: hmat(hsize_coarse,hsize_coarse) 766 complex(dpc),intent(in) :: acoeffs(hsize_coarse,hsize_coarse) 767 complex(dpc),intent(in) :: bcoeffs(hsize_coarse,hsize_coarse) 768 complex(dpc),intent(in) :: ccoeffs(hsize_coarse,hsize_coarse) 769 complex(dpc),intent(out) :: hinterp(hsize_dense,hsize_dense) 770 771 !Local variables ------------------------------ 772 !scalars 773 integer,parameter :: spin1=1, spin2=1 774 integer :: ic,iv,iv1,ic1,ik_dense,ik_coarse,it_coarse,it_dense,idiv,ibnd_coarse,ibnd_coarse1,ineighbour 775 integer :: icp,ivp,ikp_dense,ikp_coarse,itp_coarse,itp_dense,idivp,ibndp_coarse,ibndp_coarse1,ineighbourp,itp_coarse1 776 integer :: itc,it_dense1,indwithnb, corresp_ind 777 integer :: limitnbz, inb,ierr 778 real(dp) :: factor,vc_sqrt_qbz,qnorm 779 complex(dpc) :: term 780 logical :: newway, use_herm 781 !arrays 782 real(dp) :: kmkp(3),q2(3),shift(3),qinred(3),tsec(2) 783 complex(dpc),allocatable :: Cmat(:,:,:) !Temp matrices for optimized version 784 complex(dpc),allocatable :: tmp_Cmat(:) 785 complex(dpc),allocatable :: work_coeffs(:,:) 786 integer,allocatable :: band2it(:) 787 complex(dpc),ABI_CONTIGUOUS pointer :: btemp(:),ctemp(:) 788 789 !************************************************************************ 790 791 call timab(696,1,tsec) 792 793 newway = .True. 794 use_herm = .True. 795 796 if (any(BSp%interp_mode == [2,3])) then 797 if(Vcp_dense%mode /= 'CRYSTAL' .and. Vcp_dense%mode /= 'AUXILIARY_FUNCTION') then 798 ABI_BUG('Vcp_dense%mode not implemented yet !') 799 end if 800 end if 801 802 if(BSp%nsppol > 1) then 803 ABI_BUG("nsppol > 1 not yet implemented") 804 end if 805 806 factor = one/grid%ndiv 807 808 hinterp = czero; term = czero 809 810 ABI_MALLOC_OR_DIE(Cmat,(nbnd_coarse,nbnd_coarse,interpolator%nvert), ierr) 811 Cmat = czero 812 813 ABI_MALLOC(band2it,(nbnd_coarse)) 814 call int_alloc_work(interpolator,nbnd_coarse*interpolator%nvert) 815 816 btemp => interpolator%btemp 817 ctemp => interpolator%ctemp 818 819 if(newway) then 820 ABI_MALLOC(tmp_Cmat,(nbnd_coarse)) 821 ABI_MALLOC(work_coeffs,(nbnd_coarse,interpolator%nvert)) 822 end if 823 824 do ik_dense = 1,grid%nbz_dense 825 write(std_out,*) "Kdense = ",ik_dense,"/",grid%nbz_dense 826 ik_coarse = grid%dense_to_coarse(ik_dense) 827 if(use_herm) then 828 limitnbz = ik_dense 829 else 830 limitnbz = grid%nbz_dense 831 end if 832 833 do ikp_dense = 1,limitnbz 834 ikp_coarse = grid%dense_to_coarse(ikp_dense) 835 836 do iv1 = BSp%lomo_spin(spin2),BSp%homo_spin(spin2) 837 do ic1 = BSp%lumo_spin(spin2),BSp%humo_spin(spin2) 838 itp_coarse1 = BSp%vcks2t(iv1,ic1,ikp_coarse,spin2) 839 ibndp_coarse1 = (iv1-BSp%lomo_spin(spin2))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin2)+1) 840 841 band2it(ibndp_coarse1) = itp_coarse1 842 end do 843 end do 844 845 846 if (any(BSp%interp_mode == [2,3])) then 847 ! Check if we are along the diagonal 848 kmkp = Kmesh_dense%bz(:,ik_dense) - Kmesh_dense%bz(:,ikp_dense) 849 850 call wrap2_pmhalf(kmkp(:),q2(:),shift(:)) 851 qinred = MATMUL(grid%kptrlatt_coarse,q2) 852 853 ! We are outside the diagonal 854 if (BSp%interp_mode==3 .and. ANY((ABS(qinred)-tol7) > m3_width)) cycle 855 856 qnorm = two_pi*SQRT(DOT_PRODUCT(q2,MATMUL(gmet,q2))) 857 858 if(ALL(ABS(q2(:)) < 1.e-3)) then 859 vc_sqrt_qbz = SQRT(Vcp_dense%i_sz) 860 else 861 vc_sqrt_qbz = SQRT(four_pi/qnorm**2) 862 end if 863 864 !!DEBUG CHK ! 865 !!COMPUTE Qpoint 866 !call findqg0(iq_bz,g0,kmkp,Qmesh_dense%nbz,Qmesh_dense%bz,BSp%mG0) 867 868 !! * Get iq_ibz, and symmetries from iq_bz 869 !call qmesh_dense%get_BZ_item(Qmesh_dense,iq_bz,qbz,iq_ibz,isym_q,itim_q) 870 871 !if(iq_ibz > 1 .and. ABS(vc_sqrt_qbz - Vcp_dense%vc_sqrt(1,iq_ibz)) > 1.e-3) then 872 ! write(*,*) "vc_sqrt_qbz = ",vc_sqrt_qbz 873 ! write(*,*) "Vcp_dense%vc_sqrt(1,iq_ibz) = ",Vcp_dense%vc_sqrt(1,iq_ibz) 874 ! ABI_ERROR("vcp are not the same !") 875 !else if(iq_ibz == 1 .and. ABS(vc_sqrt_qbz - SQRT(Vcp_dense%i_sz)) > 1.e-3) then 876 ! write(*,*) "vc_sqrt_qbz = ",vc_sqrt_qbz 877 ! write(*,*) "SQRT(Vcp_dense%i_sz) = ",SQRT(Vcp_dense%i_sz) 878 ! ABI_ERROR("vcp are not the same !") 879 !end if 880 !!END DEBUG CHK ! 881 end if 882 883 if(newway) then 884 885 Cmat = czero 886 887 work_coeffs = czero 888 889 do ineighbour = 1,interpolator%nvert 890 891 ! Loop over the (c, v) part of the left transition 892 do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1) 893 do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1) 894 895 it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1) 896 it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1) 897 ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1) 898 899 itc = interpolator%corresp(it_coarse,ineighbour,spin1) 900 901 if (any(BSp%interp_mode == [1,3,4])) then 902 903 !work_coeffs(:,:) = hmat(itc,band2it(:),:) 904 do inb = 1,interpolator%nvert 905 work_coeffs(:,inb) = hmat(itc,interpolator%corresp(band2it(:),inb,spin2)) 906 end do 907 908 call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,& 909 & interpolator,kdense2div,& 910 & work_coeffs,tmp_Cmat,ikp_dense,& 911 & interpolator%overlaps(:,:,:,ikp_dense,spin2)) 912 913 if(any(BSp%interp_mode == [1,4])) then 914 Cmat(ibnd_coarse,:,ineighbour) = tmp_Cmat 915 else if (BSp%interp_mode == 3) then 916 Cmat(ibnd_coarse,:,ineighbour) = -tmp_Cmat 917 end if 918 end if 919 920 921 if (any(BSp%interp_mode == [2,3])) then 922 !work_coeffs(:,:) = acoeffs(itc,band2it(:),:) 923 do inb = 1,interpolator%nvert 924 work_coeffs(:,inb) = acoeffs(itc,interpolator%corresp(band2it(:),inb,spin2)) 925 end do 926 927 call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,& 928 & interpolator,kdense2div,& 929 & work_coeffs,tmp_Cmat,ikp_dense,& 930 & interpolator%overlaps(:,:,:,ikp_dense,spin2)) 931 932 tmp_Cmat = tmp_Cmat * (vc_sqrt_qbz**2) 933 Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat 934 end if 935 936 937 if (any(BSp%interp_mode == [2,3])) then 938 !work_coeffs(:,:) = bcoeffs(itc,band2it(:),:) 939 do inb = 1,interpolator%nvert 940 work_coeffs(:,inb) = bcoeffs(itc,interpolator%corresp(band2it(:),inb,spin2)) 941 end do 942 943 call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,& 944 & interpolator,kdense2div,& 945 & work_coeffs,tmp_Cmat,ikp_dense,& 946 & interpolator%overlaps(:,:,:,ikp_dense,spin2)) 947 948 tmp_Cmat = tmp_Cmat * (vc_sqrt_qbz) 949 Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat 950 end if 951 952 953 if (any(BSp%interp_mode == [2,3])) then 954 !work_coeffs(:,:) = ccoeffs(itc,band2it(:),:) 955 do inb = 1,interpolator%nvert 956 work_coeffs(:,inb) = ccoeffs(itc,interpolator%corresp(band2it(:),inb,spin2)) 957 end do 958 959 call hexc_compute_subhinterp(BSp,grid,nbnd_coarse,& 960 & interpolator,kdense2div,& 961 & work_coeffs,tmp_Cmat,ikp_dense,& 962 & interpolator%overlaps(:,:,:,ikp_dense,spin2)) 963 964 Cmat(ibnd_coarse,:,ineighbour) = Cmat(ibnd_coarse,:,ineighbour) + tmp_Cmat 965 end if 966 end do ! ic 967 end do ! iv 968 end do ! ineighbour 969 970 else 971 ! Loop over the (c, v) part of the left transition 972 do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1) 973 do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1) 974 975 it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1) 976 it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1) 977 ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1) 978 979 ! Loop over the (c', v') part of the right transition 980 do ivp = BSp%lomo_spin(spin2),BSp%homo_spin(spin2) 981 do icp = BSp%lumo_spin(spin2),BSp%humo_spin(spin2) 982 983 itp_dense = BSp%vcks2t_interp(ivp,icp,ikp_dense,spin2) 984 itp_coarse = BSp%vcks2t(ivp,icp,ikp_coarse,spin2) 985 ibndp_coarse = (ivp-Bsp%lomo_spin(spin2))*BSp%maxnbndc+(icp-BSp%lumo_spin(spin2)+1) 986 ! Now we now it_dense, and itp_dense 987 988 idivp = kdense2div(ikp_dense) 989 990 btemp = czero; ctemp = czero 991 992 ! MG TODO: This way of looping is not optimal 993 do ineighbour = 1,interpolator%nvert 994 itc = interpolator%corresp(it_coarse,ineighbour,spin1) 995 996 do ineighbourp = 1,interpolator%nvert 997 998 do iv1 = BSp%lomo_spin(spin2),BSp%homo_spin(spin2) 999 do ic1 = BSp%lumo_spin(spin2),BSp%humo_spin(spin2) 1000 1001 ibndp_coarse1 = (iv1-BSp%lomo_spin(spin2))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin2)+1) 1002 indwithnb = (ineighbourp-1)*nbnd_coarse+ibndp_coarse1 1003 itp_coarse1 = BSp%vcks2t(iv1,ic1,ikp_coarse,spin2) 1004 corresp_ind = interpolator%corresp(itp_coarse1,ineighbourp,spin2) 1005 1006 select case (BSp%interp_mode) 1007 case (1,4) 1008 interpolator%btemp(indwithnb) = hmat(itc,corresp_ind) 1009 case (2) 1010 interpolator%btemp(indwithnb) = acoeffs(itc,corresp_ind)*(vc_sqrt_qbz**2) & 1011 & + bcoeffs(itc,corresp_ind)*(vc_sqrt_qbz) & 1012 & + ccoeffs(itc,corresp_ind) 1013 case (3) 1014 ! Diff between divergence and hmat 1015 interpolator%btemp(indwithnb) = acoeffs(itc,corresp_ind)*(vc_sqrt_qbz**2) & 1016 & + bcoeffs(itc,corresp_ind)*(vc_sqrt_qbz) & 1017 & + ccoeffs(itc,corresp_ind) & 1018 & - hmat(itc,corresp_ind) 1019 case default 1020 ABI_ERROR("Wrong Bsp%interp_mode") 1021 end select 1022 1023 ctemp(indwithnb) = & 1024 & interpolator%overlaps(iv1,ivp,ineighbourp,ikp_dense,spin2) & 1025 & * GWPC_CONJG(interpolator%overlaps(ic1,icp,ineighbourp,ikp_dense,spin2)) & 1026 & *interpolator%interp_factors(ineighbourp,idivp) 1027 end do ! ic1 1028 end do !iv1 1029 1030 end do !ineighbourp 1031 Cmat(ibnd_coarse,ibndp_coarse,ineighbour) = xdotc(interpolator%nvert*nbnd_coarse,& 1032 & ctemp,1,btemp,1) 1033 end do !ineighbour 1034 1035 end do !icp 1036 end do !ivp 1037 1038 end do !ic 1039 end do !iv 1040 1041 end if 1042 1043 do iv = BSp%lomo_spin(spin1),BSp%homo_spin(spin1) 1044 do ic = BSp%lumo_spin(spin1),BSp%humo_spin(spin1) 1045 it_dense = BSp%vcks2t_interp(iv,ic,ik_dense,spin1) 1046 it_coarse = BSp%vcks2t(iv,ic,ik_coarse,spin1) 1047 ibnd_coarse = (iv-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic-BSp%lumo_spin(spin1)+1) 1048 1049 idiv = kdense2div(ik_dense) 1050 1051 do ivp = BSp%lomo_spin(spin2),BSp%homo_spin(spin2) 1052 do icp = BSp%lumo_spin(spin2),BSp%humo_spin(spin2) 1053 itp_dense = BSp%vcks2t_interp(ivp,icp,ikp_dense,spin2) 1054 itp_coarse = BSp%vcks2t(ivp,icp,ikp_coarse,spin2) 1055 ibndp_coarse = (ivp-Bsp%lomo_spin(spin2))*BSp%maxnbndc+(icp-BSp%lumo_spin(spin2)+1) 1056 1057 ! Hermicity 1058 if(use_herm .and. it_dense < itp_dense) then 1059 continue 1060 end if 1061 1062 !btemp = czero; ctemp = czero 1063 1064 do ineighbour = 1,interpolator%nvert 1065 do iv1 = BSp%lomo_spin(spin1),BSp%homo_spin(spin1) 1066 do ic1 = BSp%lumo_spin(spin1),BSp%humo_spin(spin1) 1067 ibnd_coarse1 = (iv1-BSp%lomo_spin(spin1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(spin1)+1) 1068 it_dense1 = BSp%vcks2t_interp(iv1,ic1,ik_dense,spin1) 1069 indwithnb = (ineighbour-1)*nbnd_coarse+ibnd_coarse1 1070 1071 btemp(indwithnb) = Cmat(ibnd_coarse1,ibndp_coarse,ineighbour) 1072 1073 ctemp(indwithnb) = GWPC_CONJG(interpolator%overlaps(iv1,iv,ineighbour,ik_dense,spin1)) & 1074 & *interpolator%overlaps(ic1,ic,ineighbour,ik_dense,spin1) & 1075 & *interpolator%interp_factors(ineighbour,idiv) 1076 end do !ic1 1077 end do !iv1 1078 end do !ineighbour 1079 1080 ! Save interpolated value. 1081 hinterp(it_dense,itp_dense) = xdotc(interpolator%nvert*nbnd_coarse,ctemp,1,btemp,1) 1082 !DOT_PRODUCT(ctemp,btemp) 1083 1084 end do !icp 1085 end do !ivp 1086 1087 end do !ic 1088 end do !iv 1089 1090 end do !ikp 1091 end do !ik 1092 1093 ! Enforce hermiticity 1094 if(use_herm) then 1095 do itp_dense = 1,BSp%nreh_interp(spin2) 1096 do it_dense = itp_dense,BSp%nreh_interp(spin1) 1097 if(it_dense == itp_dense) then 1098 hinterp(itp_dense,it_dense) = DBLE(hinterp(it_dense,itp_dense)) 1099 else 1100 hinterp(itp_dense,it_dense) = CONJG(hinterp(it_dense,itp_dense)) 1101 end if 1102 end do 1103 end do 1104 end if 1105 1106 ABI_FREE(Cmat) 1107 ABI_FREE(band2it) 1108 1109 if(newway) then 1110 ABI_FREE(tmp_Cmat) 1111 ABI_FREE(work_coeffs) 1112 end if 1113 1114 nullify(ctemp) 1115 nullify(btemp) 1116 call int_free_work(interpolator) 1117 1118 hinterp = hinterp*factor 1119 1120 call timab(696,2,tsec) 1121 1122 end subroutine hexc_compute_hinterp
m_hexc/hexc_compute_subhinterp [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_compute_subhinterp
FUNCTION
Compute the interpolation for work_coeffs in dense mesh
INPUTS
BSp<excparam>=Parameters for BS run grid<double_grid_t>=Double grid info nbnd_coarse=Number of bands (= nbndv * nbndc) interpolator<interpolator_t>=Interpolation info kdense2div=Mapping between dense point and coarse point work_coeffs=Coefficients to be interpolated ikp_dense=Current kpoint overlaps=Wavefunction overlaps
OUTPUT
Cmat(nbnd_coarse) = Interpolated coefficients
SOURCE
627 subroutine hexc_compute_subhinterp(BSp,grid,nbnd_coarse,& 628 & interpolator,kdense2div,work_coeffs,Cmat,ikp_dense,overlaps) 629 630 !Arguments ------------------------------------ 631 !scalars 632 integer,intent(in) :: nbnd_coarse 633 integer,intent(in) :: ikp_dense 634 type(excparam),intent(in) :: BSp 635 type(double_grid_t),intent(in) :: grid 636 type(interpolator_t),target,intent(inout) :: interpolator 637 !arrays 638 integer,intent(in) :: kdense2div(grid%nbz_dense) 639 complex(gwpc),intent(in) :: overlaps(interpolator%mband_coarse,interpolator%mband_dense,interpolator%nvert) 640 complex(dpc),intent(in) :: work_coeffs(nbnd_coarse,interpolator%nvert) 641 complex(dpc),intent(out) :: Cmat(nbnd_coarse) 642 643 !Local variables ------------------------------ 644 !scalars 645 integer,parameter :: spin1=1,spin2=1 646 integer :: iv1,ic1 647 integer :: icp,ivp,idivp,ibndp_coarse,ibndp_coarse1,ineighbourp 648 integer :: indwithnb 649 integer :: lumo2,lomo2,humo2,homo2 650 complex(dpc) :: tmp_val, tmp2, tmp4 651 !arrays 652 complex(dpc),ABI_CONTIGUOUS pointer :: btemp(:),ctemp(:) 653 654 !********************************************************************* 655 656 btemp => interpolator%btemp 657 ctemp => interpolator%ctemp 658 659 btemp = czero 660 ctemp = czero 661 662 lumo2 = BSp%lumo_spin(spin2) 663 lomo2 = BSp%lomo_spin(spin2) 664 humo2 = BSp%humo_spin(spin2) 665 homo2 = BSp%homo_spin(spin2) 666 667 Cmat = czero 668 669 idivp = kdense2div(ikp_dense) 670 671 do ineighbourp = 1,interpolator%nvert 672 673 btemp(((ineighbourp-1)*nbnd_coarse+1):(ineighbourp*nbnd_coarse)) = & 674 & interpolator%interp_factors(ineighbourp,idivp)*work_coeffs(:,ineighbourp) 675 676 end do !ineighbourp 677 678 ! Loop over the (c', v') part of the right transition 679 do ivp = lomo2,homo2 680 do icp = lumo2,humo2 681 682 ibndp_coarse = (ivp-lomo2)*BSp%maxnbndc+(icp-lumo2+1) 683 ! Now we now it_dense, and itp_dense 684 685 do ineighbourp = 1,interpolator%nvert 686 687 do iv1 = lomo2, homo2 688 ! BSp%lumo_spin(spin2),BSp%humo_spin(spin2) 689 690 tmp4 = overlaps(iv1,ivp,ineighbourp) 691 692 do ic1 = lumo2,humo2 693 tmp2 = GWPC_CONJG(overlaps(ic1,icp,ineighbourp)) 694 ! BSp%lomo_spin(spin2),BSp%homo_spin(spin2) 695 696 ibndp_coarse1 = (iv1-lomo2)*BSp%maxnbndc+(ic1-lumo2+1) 697 indwithnb = (ineighbourp-1)*nbnd_coarse+ibndp_coarse1 698 699 ctemp(indwithnb) = & 700 & tmp4 & 701 & *tmp2 702 703 end do ! iv1 704 end do !ic1 705 706 end do !ineighbourp 707 708 tmp_val = xdotc(interpolator%nvert*nbnd_coarse,ctemp,1,btemp,1) 709 !tmp_val = DOT_PRODUCT(ctemp,btemp) 710 711 Cmat(ibndp_coarse) = tmp_val 712 end do !ivp 713 end do !icp 714 715 nullify(btemp) 716 nullify(ctemp) 717 718 end subroutine hexc_compute_subhinterp
m_hexc/hexc_free [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_free
FUNCTION
Destroy the interpolator object in memory
SOURCE
1136 subroutine hexc_free(hexc) 1137 1138 !Arguments --------------------------- 1139 type(hexc_t),intent(inout) :: hexc 1140 1141 !***************************************************************************** 1142 1143 if( associated(hexc%bsp) ) then 1144 nullify(hexc%bsp) 1145 end if 1146 1147 if( associated(hexc%crystal) ) then 1148 nullify(hexc%crystal) 1149 end if 1150 1151 if( associated(hexc%kmesh_coarse) ) then 1152 nullify(hexc%kmesh_coarse) 1153 end if 1154 1155 ABI_SFREE(hexc%hreso) 1156 ABI_SFREE(hexc%hcoup) 1157 ABI_SFREE(hexc%diag_coarse) 1158 1159 end subroutine hexc_free
m_hexc/hexc_init [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_init
FUNCTION
Construct the hexc object
INPUTS
BSp<excparam>=Parameters of BS BS_files<excparam>=Files for BS Cryst<crystal_t>=Info on the crystalline structure Kmesh_coarse<kmesh_t>=Kmesh info Wfd_coarse<wfdgw_t>=Wavefunction descriptor KS_BSt<ebands_t>=Kohn-Sham band structure QP_BSt<ebands_t>=Quasi-Particle band structure comm=communicator
OUTPUT
hexc<hexc_t>=Excitonic Hamiltonian
SOURCE
234 subroutine hexc_init(hexc, BSp, BS_files, Cryst, Kmesh_coarse, Wfd_coarse, KS_BSt, QP_BSt, comm) 235 236 !Arguments --------------------------- 237 !scalars 238 integer,intent(in) :: comm 239 type(hexc_t),intent(inout) :: hexc 240 type(excparam),intent(in),target :: BSp 241 type(excfiles),intent(in),target :: BS_files 242 type(crystal_t),intent(in),target :: Cryst 243 type(kmesh_t),intent(in),target :: Kmesh_coarse 244 type(wfdgw_t),intent(in),target :: Wfd_coarse 245 type(ebands_t),intent(in),target :: KS_BSt, QP_BSt 246 !arrays 247 248 !Local variables --------------------- 249 !scalars 250 integer :: ierr,ncid,ncerr 251 integer :: max_r, max_c 252 integer :: hsize 253 integer :: spin, spad, itt ! For diagonal ! 254 logical :: is_resonant, diago_is_real, use_mpio=.FALSE. 255 character(len=fnlen) :: hreso_fname, hcoup_fname 256 !character(len=500) :: msg 257 !arrays 258 complex(dpc),allocatable :: test(:,:) 259 260 !***************************************************************************** 261 262 hexc%bsp => BSp 263 hexc%bs_files => BS_files 264 hexc%crystal => Cryst 265 hexc%kmesh_coarse => Kmesh_coarse 266 267 268 hexc%comm = comm 269 hsize = SUM(BSp%nreh) 270 271 hexc%ks_bst => KS_BSt 272 hexc%qp_bst => QP_BSt 273 hexc%wfd => Wfd_coarse 274 hexc%wfd_coarse => Wfd_coarse 275 hexc%kmesh => Kmesh_coarse 276 hexc%hsize_coarse = hsize 277 hexc%hsize = hsize 278 hexc%nbz = Kmesh_coarse%nbz 279 280 hexc%nbnd_coarse = BSp%maxnbndv*BSp%maxnbndc 281 282 ! Divide the columns of the Hamiltonian among the nodes. 283 call xmpi_split_work(hsize,comm,hexc%my_t1,hexc%my_t2) 284 285 hexc%my_nt = hexc%my_t2 - hexc%my_t1 + 1 286 ABI_CHECK(hexc%my_nt>0,"found processor with 0 rows") 287 288 ABI_MALLOC_OR_DIE(hexc%hreso,(hsize,hexc%my_t1:hexc%my_t2), ierr) 289 290 ! Read the resonant block from file. 291 if (BS_files%in_hreso /= BSE_NOFILE) then 292 hreso_fname = BS_files%in_hreso 293 else 294 hreso_fname = BS_files%out_hreso 295 end if 296 297 is_resonant=.TRUE.; diago_is_real=(.not.BSp%have_complex_ene) 298 call exc_read_rcblock(hreso_fname,Bsp,is_resonant,diago_is_real,BSp%nsppol,BSp%nreh,hsize,& 299 & hexc%my_t1,hexc%my_t2,hexc%hreso,use_mpio,comm) 300 301 !BEGIN DEBUG 302 if (use_mpio) then 303 ABI_WARNING("Testing MPI-IO routines") 304 ABI_MALLOC_OR_DIE(test,(hsize,hexc%my_t1:hexc%my_t2), ierr) 305 diago_is_real=(.not.BSp%have_complex_ene) 306 call exc_read_rcblock(hreso_fname,Bsp,is_resonant,diago_is_real,Bsp%nsppol,Bsp%nreh,hsize,& 307 & hexc%my_t1,hexc%my_t2,test,.FALSE.,comm) 308 test = test-hexc%hreso 309 write(std_out,*)"DEBUG: Diff MPI-IO - Fortran ",MAXVAL(ABS(test)) 310 max_r=20; max_c=10 311 write(std_out,*)" **** Testing resonant block **** " 312 call print_arr(test,max_r=max_r,max_c=max_c,unit=std_out) 313 if (BSp%nsppol==2) then 314 write(std_out,*)" **** D down down ****" 315 call print_arr(test(hsize/2+1:,hsize/2+1:),max_r=max_r,max_c=max_c,unit=std_out) 316 write(std_out,*)" **** V up down ****" 317 call print_arr(test(1:hsize/2,hsize/2+1:),max_r=max_r,max_c=max_c,unit=std_out) 318 write(std_out,*)" **** V down up ****" 319 call print_arr(test(hsize/2+1:,1:hsize/2),max_r=max_r,max_c=max_c,unit=std_out) 320 end if 321 ABI_FREE(test) 322 end if 323 !END DEBUG 324 325 ! 326 ! Read coupling block. 327 if (BSp%use_coupling>0) then 328 ABI_CHECK(.not. Bsp%use_interp,"interpolation with coupling not coded!") 329 if (BS_files%in_hcoup /= BSE_NOFILE) then 330 hcoup_fname = BS_files%in_hcoup 331 else 332 hcoup_fname = BS_files%out_hcoup 333 end if 334 335 ABI_MALLOC_OR_DIE(hexc%hcoup,(hsize,hexc%my_t1:hexc%my_t2), ierr) 336 is_resonant=.FALSE.; diago_is_real=.FALSE. 337 call exc_read_rcblock(hcoup_fname,Bsp,is_resonant,diago_is_real,BSp%nsppol,BSp%nreh,hsize,& 338 & hexc%my_t1,hexc%my_t2,hexc%hcoup,use_mpio,comm) 339 !call symmetrize(hcoup,"ALL") 340 341 if (use_mpio) then 342 ABI_WARNING("Testing MPI-IO routines") 343 ABI_MALLOC_OR_DIE(test,(hsize,hexc%my_t1:hexc%my_t2), ierr) 344 diago_is_real=.FALSE. 345 call exc_read_rcblock(hcoup_fname,Bsp,is_resonant,diago_is_real,BSp%nsppol,Bsp%nreh,hsize,& 346 & hexc%my_t1,hexc%my_t2,test,.FALSE.,comm) 347 test = test-hexc%hcoup 348 write(std_out,*)"DEBUG: Diff MPI-IO - Fortran ",MAXVAL(ABS(test)) 349 max_r=20; max_c=10 350 write(std_out,*)" **** Testing coupling block **** " 351 call print_arr(test,max_r=max_r,max_c=max_c,unit=std_out) 352 if (BSp%nsppol==2) then 353 write(std_out,*)" **** D down down ****" 354 call print_arr(test(hsize/2+1:,hsize/2+1:),max_r=max_r,max_c=max_c,unit=std_out) 355 write(std_out,*)" **** V up down ****" 356 call print_arr(test(1:hsize/2,hsize/2+1:),max_r=max_r,max_c=max_c,unit=std_out) 357 write(std_out,*)" **** V down up ****" 358 call print_arr(test(hsize/2+1:,1:hsize/2),max_r=max_r,max_c=max_c,unit=std_out) 359 end if 360 ABI_FREE(test) 361 end if 362 end if 363 364 if(BSp%prt_ncham .or. BSp%use_interp) then 365 ! I want to store the diagonal part for future use (printing or interpolation) ! 366 ABI_MALLOC(hexc%diag_coarse,(hexc%hsize_coarse)) 367 spad=0 368 do spin=1,BSp%nsppol 369 if(spin==2) spad=BSp%nreh(1) 370 do itt=1,BSp%nreh(spin) ! 1 is for spin 1 371 hexc%diag_coarse(spad+itt) = Bsp%Trans(itt,spin)%en 372 end do 373 end do 374 375 if (BSp%prt_ncham) then 376 ncerr = nctk_open_create(ncid, trim(hexc%BS_files%out_basename)//"_HEXC.nc", xmpi_comm_self) 377 NCF_CHECK_MSG(ncerr, "Creating HEXC file") 378 call exc_ham_ncwrite(ncid, hexc%Kmesh_coarse, hexc%BSp, hexc%hsize_coarse, hexc%BSp%nreh, & 379 & hexc%BSp%vcks2t,hexc%hreso,hexc%diag_coarse) 380 NCF_CHECK(nf90_close(ncid)) 381 end if 382 end if 383 384 end subroutine hexc_init
m_hexc/hexc_interp_free [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_interp_free
FUNCTION
Destroy the interpolator object in memory
SOURCE
1173 subroutine hexc_interp_free(hexc_i) 1174 1175 !Arguments --------------------------- 1176 type(hexc_interp_t),intent(inout) :: hexc_i 1177 1178 !***************************************************************************** 1179 1180 ABI_SFREE(hexc_i%kdense2div) 1181 ABI_SFREE(hexc_i%div2kdense) 1182 ABI_SFREE(hexc_i%diag_dense) 1183 ABI_SFREE(hexc_i%hinterp) 1184 ABI_SFREE(hexc_i%all_hmat) 1185 ABI_SFREE(hexc_i%all_acoeffs) 1186 ABI_SFREE(hexc_i%all_bcoeffs) 1187 ABI_SFREE(hexc_i%all_ccoeffs) 1188 1189 if( associated(hexc_i%kmesh_dense) ) then 1190 nullify(hexc_i%kmesh_dense) 1191 end if 1192 1193 if( associated(hexc_i%vcp_dense) ) then 1194 nullify(hexc_i%vcp_dense) 1195 end if 1196 1197 call interpolator_free(hexc_i%interpolator) 1198 1199 end subroutine hexc_interp_free
m_hexc/hexc_interp_init [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_interp_init
FUNCTION
Construct the hexc_interp object
INPUTS
hexc<hexc_t>=Excitonic hamiltonian Kmesh_dense<kmesh_t>=Kmesh info Vcp_dense<vcoul_t>=Dense mesh info about coulomb double_grid<double_grid_t>=Link between dense and coarse mesh Wfd_dense<wfdgw_t>=Wavefunction descriptor KS_BSt_dense<ebands_t>=Kohn-Sham band structure QP_BSt_dense<ebands_t>=Quasi-Particle band structure Psps <type(pseudopotential_type)>=variables related to pseudopotentials. Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data. comm=communicator
OUTPUT
hexc_i<hexc_interp_t>=Interpolated excitonic hamiltonian
SIDE EFFECTS
hexc Will be modified so that the size of the problem is full interpolated hamiltonian Wfd, Wfd_dense The memory might be modified by computing wavefunctions
SOURCE
420 subroutine hexc_interp_init(hexc_i, hexc, m3_width, method, Kmesh_dense, Vcp_dense, & 421 & double_grid, Wfd_dense, KS_BSt_dense, QP_BSt_dense, Psps, Pawtab) 422 423 !Arguments --------------------------- 424 !scalars 425 integer,intent(in) :: method 426 real(dp),intent(in) :: m3_width 427 type(hexc_t),intent(inout) :: hexc 428 type(hexc_interp_t),intent(inout) :: hexc_i 429 type(double_grid_t),intent(in),target :: double_grid 430 type(wfdgw_t),intent(inout),target :: Wfd_dense !, Wfd 431 type(kmesh_t),intent(in),target :: Kmesh_dense 432 type(pseudopotential_type),intent(in) :: Psps 433 type(vcoul_t),intent(in),target :: Vcp_dense 434 type(ebands_t),intent(in),target :: KS_BSt_dense, QP_BSt_dense 435 !arrays 436 type(pawtab_type),intent(in) :: Pawtab(hexc%crystal%ntypat*hexc%Wfd_coarse%usepaw) 437 438 !Local variables --------------------- 439 !scalars 440 integer,parameter :: spin1 = 1 441 integer :: nsppol,ierr,ii,itt,nproc, my_rank,hsize 442 real(dp),parameter :: threshold = 0.1_dp 443 type(excparam) :: BSp 444 logical :: is_resonant, diago_is_real, use_mpio 445 !arrays 446 character(len=fnlen) :: tmpfname, hreso_fname 447 448 !***************************************************************************** 449 450 BSp = hexc%bsp 451 ABI_CHECK(BSp%nsppol == 1,"nsppol > 1 not implemented yet") 452 453 hsize = hexc%hsize_coarse 454 455 nproc = xmpi_comm_size(hexc%comm); my_rank= xmpi_comm_rank(hexc%comm) 456 nsppol = hexc%Bsp%nsppol 457 458 ABI_CHECK(nproc == 1,"Parallelization not available in interpolation") 459 460 hexc_i%m3_width = m3_width 461 462 hexc_i%kmesh_dense => Kmesh_dense 463 hexc_i%vcp_dense => Vcp_dense 464 hexc_i%hsize_dense = SUM(BSp%nreh_interp) 465 hexc%hsize = hexc_i%hsize_dense 466 hexc%nbz = Kmesh_dense%nbz 467 468 hexc%ks_bst => KS_BSt_dense 469 hexc%qp_bst => QP_BSt_dense 470 hexc%wfd => Wfd_dense 471 hexc%kmesh => Kmesh_dense 472 473 ! No parallelization ! 474 hexc%my_t1 = 1 475 hexc%my_t2 = hexc_i%hsize_dense 476 477 ! Initialize the interpolator 478 call interpolator_init(hexc_i%interpolator, double_grid, Wfd_dense, hexc%Wfd_coarse, Kmesh_dense, & 479 & hexc%Kmesh_coarse, hexc%BSp, hexc%crystal, Psps, Pawtab, method) 480 481 if(BSp%sum_overlaps) then 482 call interpolator_normalize(hexc_i%interpolator) 483 end if 484 485 ABI_MALLOC(hexc_i%kdense2div,(double_grid%nbz_dense)) 486 ABI_MALLOC(hexc_i%div2kdense,(double_grid%nbz_coarse,double_grid%ndiv)) 487 488 call compute_corresp(double_grid,hexc_i%div2kdense,hexc_i%kdense2div) 489 490 if (any(BSp%interp_mode == [2,3])) then 491 ! Read a, b, c coefficient matrices from file. 492 ! For the time being, we read the full matrix in a temporary array, and 493 ! then we store the data in a form suitable for the interpolation. 494 is_resonant=.TRUE.; diago_is_real=(.not.BSp%have_complex_ene); use_mpio=.FALSE. 495 496 if (hexc%bs_files%in_hreso /= BSE_NOFILE) then 497 hreso_fname = hexc%bs_files%in_hreso 498 else 499 hreso_fname = hexc%bs_files%out_hreso 500 end if 501 502 tmpfname = hreso_fname; ii = LEN_TRIM(hreso_fname) 503 504 ! TODO: Write new IO routines to read MPI-distributed data in a format suitable for the interpolation 505 tmpfname(ii-2:ii+1) = 'ABSR' 506 ABI_MALLOC_OR_DIE(hexc_i%all_acoeffs,(hsize,hsize), ierr) 507 call exc_read_rcblock(tmpfname,Bsp,is_resonant,diago_is_real,nsppol,BSp%nreh,hsize,1,hsize,& 508 & hexc_i%all_acoeffs,use_mpio,hexc%comm) 509 510 tmpfname(ii-2:ii+1) = 'BBSR' 511 ABI_MALLOC_OR_DIE(hexc_i%all_bcoeffs,(hsize,hsize), ierr) 512 call exc_read_rcblock(tmpfname,Bsp,is_resonant,diago_is_real,nsppol,BSp%nreh,hsize,1,hsize,& 513 & hexc_i%all_bcoeffs,use_mpio,hexc%comm) 514 515 tmpfname(ii-2:ii+1) = 'CBSR' 516 ABI_MALLOC_OR_DIE(hexc_i%all_ccoeffs,(hsize,hsize), ierr) 517 call exc_read_rcblock(tmpfname,Bsp,is_resonant,diago_is_real,nsppol,BSp%nreh,hsize,1,hsize,& 518 & hexc_i%all_ccoeffs,use_mpio,hexc%comm) 519 end if 520 521 ! Compute overlaps & compute all hmat 522 ABI_MALLOC_OR_DIE(hexc_i%all_hmat,(hsize,hsize), ierr) 523 524 hexc_i%all_hmat(:,:) = hexc%hreso(:,:) 525 526 do itt=1,hsize 527 hexc_i%all_hmat(itt,itt) = hexc_i%all_hmat(itt,itt) - hexc%diag_coarse(itt) 528 end do 529 530 ! I don't need the diag_coarse any more 531 ABI_FREE(hexc%diag_coarse) 532 533 ! Compute diagonal part of the dense Ham 534 ABI_MALLOC(hexc_i%diag_dense,(hexc_i%hsize_dense)) 535 do itt=1,BSp%nreh_interp(spin1) ! 1 is for spin 1 536 hexc_i%diag_dense(itt) = Bsp%Trans_interp(itt,spin1)%en 537 end do 538 539 end subroutine hexc_interp_init
m_hexc/hexc_interp_matmul [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_interp_matmul
FUNCTION
Compute matrix-vector product Hmat * phi by interpolating coarse Hmat
INPUTS
BSp<type(excparam)>=Parameters defining the BS calculation hsize_coarse = Size of the coarse hamiltonian hsize_dense = Size of the dense hamiltonian hmat(hsize_coarse,hsize_coarse,8) = coarse hamiltonian phi(hsize_dense) = ket on which apply the matrix grid <double_grid_t> = Correspondence between coarse and dense k-mesh. nbnd_coarse = Total number of bands interpolator<interpolator_t> = Interpolator div2kdense = Mapping from coarse and division -> dense mesh kdense2div = Mapping from dense mesh -> division
OUTPUT
hphi(hsize_dense) = Interp(hmat)*phi
SOURCE
1228 subroutine hexc_interp_matmul(BSp,hsize_coarse,hsize_dense,hmat,phi,hphi,grid,& 1229 & nbnd_coarse,interpolator,div2kdense,kdense2div) 1230 1231 !Arguments ------------------------------------ 1232 !scalars 1233 integer,intent(in) :: hsize_coarse,hsize_dense,nbnd_coarse !,ntrans 1234 type(excparam),intent(in) :: BSp 1235 type(double_grid_t),intent(in) :: grid 1236 type(interpolator_t),intent(in) :: interpolator 1237 !arrays 1238 integer,intent(in) :: div2kdense(grid%nbz_coarse,grid%ndiv), kdense2div(grid%nbz_dense) 1239 complex(dpc),intent(in) :: phi(hsize_dense) 1240 complex(dpc),intent(in) :: hmat(hsize_coarse,hsize_coarse) 1241 complex(dpc),intent(inout) :: hphi(hsize_dense) 1242 1243 !Local variables ------------------------------ 1244 !scalars 1245 integer :: itt,ik_dense,ik_coarse,it_coarse 1246 integer :: ic,iv,iv1,ic1, ibnd_coarse 1247 integer :: ibnd_coarse1 1248 integer :: ineighbour,idense,ikpt 1249 integer :: my_k1,my_k2,ind_with_nb,is, is1 1250 real(dp) :: factor 1251 complex(dpc) :: tmp 1252 logical,parameter :: use_blas=.True. 1253 !arrays 1254 integer :: allindices(nbnd_coarse) 1255 complex(dpc) :: allp(hsize_coarse,interpolator%nvert), test(hsize_coarse) 1256 complex(dpc) :: ophi(grid%nbz_dense,interpolator%nvert,nbnd_coarse) 1257 complex(dpc),allocatable :: b(:), c(:),A(:,:) 1258 complex(dpc),allocatable :: tmp_array(:), tmp_array2(:,:) 1259 1260 !************************************************************************ 1261 1262 factor = one/grid%ndiv 1263 1264 !hphi = czero 1265 1266 ! Outer index : k point in the dense zone 1267 ! Sum over vc 1268 ! Index of result : k point in the dense zone, v2,c2,neighbour 1269 1270 ! Parallelization on nbz in the coarse mesh ! 1271 my_k1 = 1 1272 my_k2 = grid%nbz_coarse 1273 1274 ABI_MALLOC(A,(interpolator%nvert*nbnd_coarse,nbnd_coarse)) 1275 ABI_MALLOC(b,(nbnd_coarse)) 1276 ABI_MALLOC(c,(interpolator%nvert*nbnd_coarse)) 1277 1278 c = czero; ophi = czero 1279 1280 do ik_dense = 1,grid%nbz_dense 1281 ! if( ik_dense is not in my set of k-points) 1282 ! ! continue 1283 ! 1284 do is1 = 1, BSp%nsppol 1285 do iv1 = BSp%lomo_spin(is1),Bsp%homo_spin(is1) 1286 do ic1 = BSp%lumo_spin(is1),Bsp%humo_spin(is1) 1287 ibnd_coarse = (iv1-BSp%lomo_spin(is1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is1)+1) 1288 itt = BSp%vcks2t_interp(iv1,ic1,ik_dense,is1) 1289 allindices(ibnd_coarse) = itt 1290 end do !ic1 1291 end do !iv1 1292 end do !is1 1293 1294 b(:) = phi(allindices(:)) 1295 1296 do is = 1, BSp%nsppol 1297 do iv = BSp%lomo_spin(is),Bsp%homo_spin(is) 1298 do ic = BSp%lumo_spin(is),Bsp%humo_spin(is) 1299 ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1) 1300 idense = Bsp%vcks2t_interp(iv,ic,ik_dense,is) 1301 1302 do ineighbour = 1,interpolator%nvert 1303 ind_with_nb = (ineighbour-1)*(nbnd_coarse)+ibnd_coarse 1304 1305 !A(ind_with_nb,:) = overlaps(allindices(:),ibnd_coarse,ineighbour) 1306 1307 ! Should be optimized !!! 1308 do iv1 = BSp%lomo_spin(is),Bsp%homo_spin(is) 1309 do ic1 = BSp%lumo_spin(is),Bsp%humo_spin(is) 1310 ibnd_coarse1 = (iv1-BSp%lomo_spin(is))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is)+1) 1311 A(ind_with_nb,ibnd_coarse1) = GWPC_CONJG(interpolator%overlaps(iv,iv1,ineighbour,ik_dense,is)) & 1312 & *interpolator%overlaps(ic,ic1,ineighbour,ik_dense,is) 1313 end do !ic1 1314 end do !iv1 1315 end do !ineighbour 1316 end do !ic 1317 end do !iv 1318 end do !is 1319 1320 if(use_blas) then 1321 call xgemv('N',interpolator%nvert*nbnd_coarse,nbnd_coarse,cone,A,interpolator%nvert*nbnd_coarse,b,1,czero,c,1) 1322 else 1323 c = MATMUL(A,b) 1324 end if 1325 1326 do is = 1, BSp%nsppol 1327 do iv = BSp%lomo_spin(is),BSp%homo_spin(is) 1328 do ic = BSp%lumo_spin(is),BSp%humo_spin(is) 1329 ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1) 1330 do ineighbour = 1,interpolator%nvert 1331 ind_with_nb = (ineighbour-1)*(nbnd_coarse)+ibnd_coarse 1332 ophi(ik_dense,ineighbour,ibnd_coarse) = c(ind_with_nb) 1333 end do !ineighbour 1334 end do !ic 1335 end do !iv 1336 end do !is 1337 1338 end do !ik_dense 1339 1340 ABI_FREE(A) 1341 ABI_FREE(b) 1342 ABI_FREE(c) 1343 1344 !call xmpi_sum_(ophi,comm,ierr) 1345 1346 ! Outer index : k,v,c in the coarse zone, ineighbour 1347 ! Sum over all k-dense relative to one coarse point 1348 ! Index of result : k,v,c in the coarse zone, ineighbour 1349 1350 ABI_MALLOC(b,(grid%ndiv)) 1351 ABI_MALLOC(c,(grid%ndiv)) 1352 1353 allp = czero 1354 1355 do is = 1, BSp%nsppol 1356 do ineighbour = 1,interpolator%nvert 1357 1358 do it_coarse = 1, BSp%nreh(is) 1359 ibnd_coarse = (Bsp%trans(it_coarse,is)%v-BSp%lomo_spin(is))*BSp%maxnbndc+& 1360 & (BSp%Trans(it_coarse,is)%c-BSp%lumo_spin(is)+1) 1361 ik_coarse = BSp%trans(it_coarse,is)%k 1362 !b(:) = interp_factors(it_coarse,ineighbour,:) 1363 b(:) = interpolator%interp_factors(ineighbour,:) 1364 !c(:) = ophi(indices(it_coarse,:),ineighbour,ibnd_coarse) 1365 c(:) = ophi(div2kdense(ik_coarse,:),ineighbour,ibnd_coarse) 1366 tmp = DOT_PRODUCT(b,c) 1367 allp(it_coarse,ineighbour) = tmp 1368 end do 1369 1370 end do 1371 end do 1372 1373 !call xmpi_sum_(allp,comm,ierr) 1374 1375 ABI_FREE(b) 1376 ABI_FREE(c) 1377 1378 ABI_MALLOC(tmp_array,(hsize_coarse)) 1379 ABI_MALLOC(tmp_array2,(hsize_coarse,hsize_coarse)) 1380 tmp_array(:) = czero 1381 tmp_array2(:,:) = czero 1382 1383 test = czero 1384 1385 ! Second step : Multiplication by hmat 1386 do ineighbour = 1,interpolator%nvert 1387 if(use_blas) then 1388 !call xgemv('N',hsize_coarse,hsize_coarse,cone,factor*(hmat(:,:,ineighbour)),hsize_coarse,allp(:,ineighbour),1,czero,tmp_array,1) 1389 !tmp_array2 = hmat(:,:,ineighbour) 1390 tmp_array2 = hmat(:,interpolator%corresp(:,ineighbour,1)) ! 1 is for spin 1391 tmp_array2 = factor*tmp_array2 1392 call xgemv('N',hsize_coarse,hsize_coarse,cone,tmp_array2,hsize_coarse,allp(:,ineighbour),1,czero,tmp_array,1) 1393 test = test + tmp_array 1394 else 1395 test = test+MATMUL(factor*(hmat(:,interpolator%corresp(:,ineighbour,1))),allp(:,ineighbour)) 1396 end if 1397 end do 1398 1399 ABI_FREE(tmp_array) 1400 ABI_FREE(tmp_array2) 1401 1402 ! Outer index : ineighbour 1403 ! Sum over all v c 1404 ! Index of result : ineighbour, k_dense, v,c 1405 ABI_MALLOC(A,(nbnd_coarse,nbnd_coarse)) 1406 ABI_MALLOC(b,(nbnd_coarse)) 1407 ABI_MALLOC(c,(nbnd_coarse)) 1408 c = czero 1409 1410 do ineighbour = 1,interpolator%nvert 1411 do ik_dense = 1,grid%nbz_dense 1412 1413 do is1 = 1, Bsp%nsppol 1414 do iv1 = Bsp%lomo_spin(is1),Bsp%homo_spin(is1) 1415 do ic1 = BSp%lumo_spin(is1), Bsp%humo_spin(is1) 1416 ibnd_coarse = (iv1-BSp%lomo_spin(is1))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is1)+1) 1417 1418 ik_coarse = grid%dense_to_coarse(ik_dense) 1419 itt = BSp%vcks2t(iv1,ic1,ik_coarse,is1) 1420 b(ibnd_coarse) = test(interpolator%corresp(itt,ineighbour,is1)) 1421 end do ! ic1 1422 end do ! iv1 1423 end do ! is1 1424 1425 do is = 1, BSp%nsppol 1426 do iv = BSp%lomo_spin(is),Bsp%homo_spin(is) 1427 do ic = BSp%lumo_spin(is),BSp%humo_spin(is) 1428 ibnd_coarse = (iv-BSp%lomo_spin(is))*Bsp%maxnbndc+(ic-BSp%lumo_spin(is)+1) 1429 idense = BSp%vcks2t_interp(iv,ic,ik_dense,is) 1430 1431 !A(ibnd_coarse,:) = CONJG(overlaps(idense,:,ineighbour)) 1432 1433 ! Should be optimized !!! 1434 do iv1 = BSp%lomo_spin(is),Bsp%homo_spin(is) 1435 do ic1 = BSp%lumo_spin(is),Bsp%humo_spin(is) 1436 ibnd_coarse1 = (iv1-BSp%lomo_spin(is))*BSp%maxnbndc+(ic1-BSp%lumo_spin(is)+1) 1437 A(ibnd_coarse,ibnd_coarse1) = (interpolator%overlaps(iv1,iv,ineighbour,ik_dense,is)) & 1438 & *GWPC_CONJG(interpolator%overlaps(ic1,ic,ineighbour,ik_dense,is)) 1439 end do !ic1 1440 end do !iv1 1441 end do ! ic 1442 end do !iv 1443 end do !is 1444 1445 if(use_blas) then 1446 call xgemv('N',nbnd_coarse,nbnd_coarse,cone,A,nbnd_coarse,b,1,czero,c,1) 1447 else 1448 c = MATMUL(A,b) 1449 end if 1450 1451 do is = 1, BSp%nsppol 1452 do iv = BSp%lomo_spin(is),Bsp%homo_spin(is) 1453 do ic = BSp%lumo_spin(is),BSp%humo_spin(is) 1454 ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1) 1455 idense = Bsp%vcks2t_interp(iv,ic,ik_dense,is) 1456 !ophi(ik_dense,ineighbour,ibnd_coarse) = c(idense) 1457 ophi(ik_dense,ineighbour,ibnd_coarse) = c(ibnd_coarse) 1458 end do 1459 end do 1460 end do 1461 1462 end do ! ik_dense 1463 end do ! ineighbour 1464 1465 !call xmpi_sum_(ophi,comm,ierr) 1466 1467 ABI_FREE(A) 1468 ABI_FREE(b) 1469 ABI_FREE(c) 1470 1471 ! Outer indices : it_dense 1472 ! Sum over neighbours 1473 ! Index of result : it_dense (ik,ic,iv dense) 1474 1475 ABI_MALLOC(b,(interpolator%nvert)) 1476 ABI_MALLOC(c,(interpolator%nvert)) 1477 1478 do is = 1, BSp%nsppol 1479 do itt = 1,BSp%nreh_interp(is) 1480 ! From itt -> ik_ibz,ic,iv 1481 ik_dense = BSp%Trans_interp(itt,is)%k 1482 ic = BSp%Trans_interp(itt,is)%c 1483 iv = BSp%Trans_interp(itt,is)%v 1484 1485 ! From ik_ibz in the dense mesh -> indices_dense 1486 ik_coarse = grid%dense_to_coarse(ik_dense) 1487 it_coarse = BSp%vcks2t(iv,ic,ik_coarse,is) 1488 1489 ibnd_coarse = (iv-BSp%lomo_spin(is))*BSp%maxnbndc+(ic-BSp%lumo_spin(is)+1) 1490 1491 ikpt = kdense2div(ik_dense) 1492 !ikpt = -1 1493 !do ix = 1,grid%ndiv 1494 ! if (indices(it_coarse,ix) == ik_dense) then 1495 ! ikpt = ix 1496 ! exit 1497 ! end if 1498 !end do 1499 !ABI_CHECK(ikpt/=-1,"Cannot find ik_dense") 1500 1501 !b = interp_factors(it_coarse,:,ikpt) 1502 b = interpolator%interp_factors(:,ikpt) 1503 c = ophi(ik_dense,:,ibnd_coarse) 1504 1505 hphi(itt) = hphi(itt) + xdotc(interpolator%nvert, b, 1, c, 1) 1506 end do 1507 end do 1508 1509 ABI_FREE(b) 1510 ABI_FREE(c) 1511 1512 end subroutine hexc_interp_matmul
m_hexc/hexc_interp_t [ Types ]
NAME
hexc_interp_t
FUNCTION
Store information about interpolation of excitonic hamiltonian
SOURCE
142 type,public :: hexc_interp_t 143 144 !scalars 145 integer :: hsize_dense 146 ! Size of the dense hamiltonian 147 148 real(dp) :: m3_width 149 ! Width of the region where M3 is applied instead of M1 150 151 type(interpolator_t) :: interpolator 152 ! Interpolator containing overlaps and interpolation info 153 154 ! Pointers to datatypes that are already in memory 155 type(kmesh_t),pointer :: kmesh_dense => null() 156 ! kmesh of the dense mesh 157 158 type(vcoul_t),pointer :: vcp_dense => null() 159 ! coulomb interaction on the dense mesh 160 161 !arrays 162 integer,allocatable :: kdense2div(:) 163 ! kdense2div(nbz_dense) 164 ! Index of kpoint -> Index of division 165 166 integer,allocatable :: div2kdense(:,:) 167 ! div2kdense(nbz_coarse,ndiv) 168 ! Index of kpoint coarse + Index of division -> Index of kdense 169 170 complex(dpc),allocatable :: diag_dense(:) 171 ! diag_dense(hsize_dense) 172 ! Diagonal part of the dense hamiltonian 173 174 complex(dpc),allocatable :: hinterp(:,:) 175 ! hinterp(hsize_dense,hsize_dense) 176 ! Interpolated hamiltonian 177 178 complex(dpc),allocatable :: all_hmat(:,:) 179 ! all_hmat,(hsize,hsize)) 180 ! Coarse excitonic matrix in a format suitable for interpolation in k-space 181 182 complex(dpc),allocatable :: all_acoeffs(:,:) 183 ! all_acoeffs(hsize,hsize)) 184 ! a coefficients in a format suitable for interpolation in k-space 185 186 complex(dpc),allocatable :: all_bcoeffs(:,:) 187 ! all_bcoeffs(hsize,hsize)) 188 ! b coefficients in a format suitable for interpolation in k-space 189 190 complex(dpc),allocatable :: all_ccoeffs(:,:) 191 ! all_ccoeffs(hsize,hsize)) 192 ! c coefficients in a format suitable for interpolation in k-space 193 194 end type hexc_interp_t
m_hexc/hexc_matmul_elphon [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_matmul_elphon
FUNCTION
Compute H | \psi > + E_{elphon} | \psi >
INPUTS
hexc<hexc_t> = Excitonic hamiltonian phi = Input ket ep_renorm = vector with electron-phonon renorms op = 'N' for H | psi >, 'C' for H^\dagger | psi >
OUTPUT
hphi = hreso * phi + ep_renorm * phi
SOURCE
1595 subroutine hexc_matmul_elphon(hexc, phi, hphi, op, ep_renorm) 1596 1597 !Arguments --------------------------- 1598 type(hexc_t),intent(in) :: hexc 1599 character,intent(in) :: op 1600 complex(dpc),intent(in) :: phi(hexc%my_nt) 1601 complex(dpc),intent(out) :: hphi(hexc%hsize) 1602 complex(dpc),intent(in) :: ep_renorm(hexc%hsize) 1603 1604 !Local variables --------------------- 1605 integer :: ierr 1606 real(dp) :: tsec(2) 1607 1608 !***************************************************************************** 1609 1610 call timab(697,1,tsec) 1611 1612 if(hexc%BSp%use_interp) then 1613 ABI_ERROR('Not yet implemented with interpolation !') 1614 else ! No interpolation 1615 ! As our matrix is hermitian (hreso), we should always use 'N' here (it is stored column-wise !) 1616 call xgemv('N',hexc%hsize,hexc%my_nt,cone,hexc%hreso,hexc%hsize,phi,1,czero,hphi,1) 1617 1618 !!! ep_renorm is stored on each cpu 1619 if (op == 'N') then 1620 hphi(hexc%my_t1:hexc%my_t2) = hphi(hexc%my_t1:hexc%my_t2) + ep_renorm(hexc%my_t1:hexc%my_t2) * phi 1621 else if(op == 'C') then 1622 hphi(hexc%my_t1:hexc%my_t2) = hphi(hexc%my_t1:hexc%my_t2) + CONJG(ep_renorm(hexc%my_t1:hexc%my_t2)) * phi 1623 end if 1624 call xmpi_sum(hphi,hexc%comm,ierr) 1625 end if 1626 1627 call timab(697,2,tsec) 1628 1629 end subroutine hexc_matmul_elphon
m_hexc/hexc_matmul_full [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_matmul_full
FUNCTION
Compute H | \psi >
INPUTS
hexc<hexc_t> = Excitonic hamiltonian hexc_i<hexc_interp_t> = Interpolated hamiltonian phi = Input ket parity = -1 or +1 parameter
OUTPUT
hphi = hreso * phi + parity * hcoup * CONJ(phi)
SOURCE
1652 subroutine hexc_matmul_full(hexc, hexc_i, phi, hphi, parity) 1653 1654 !Arguments --------------------------- 1655 integer,intent(in) :: parity 1656 type(hexc_t),intent(in) :: hexc 1657 type(hexc_interp_t),intent(in) :: hexc_i 1658 complex(dpc),intent(in) :: phi(hexc%hsize) 1659 complex(dpc),intent(out) :: hphi(hexc%hsize) 1660 1661 !Local variables --------------------- 1662 real(dp) :: tsec(2) 1663 1664 !***************************************************************************** 1665 1666 call timab(697,1,tsec) 1667 1668 ABI_UNUSED(hexc_i%hsize_dense) 1669 1670 if(hexc%BSp%use_interp) then 1671 ABI_ERROR("Coupling is not yet implemented with interpolation") 1672 else ! No interpolation 1673 hphi = MATMUL(hexc%hreso,phi) + parity * MATMUL(hexc%hcoup,CONJG(phi)) 1674 end if 1675 1676 call timab(697,2,tsec) 1677 1678 end subroutine hexc_matmul_full
m_hexc/hexc_matmul_tda [ Functions ]
[ Top ] [ m_hexc ] [ Functions ]
NAME
hexc_matmul_tda
FUNCTION
Compute H | \psi >
INPUTS
hexc<hexc_t> = Excitonic hamiltonian hexc_i<hexc_interp_t> = Interpolated excitonic hamiltonian phi = Input ket
OUTPUT
hphi = hreso * phi
SOURCE
1534 subroutine hexc_matmul_tda(hexc, hexc_i, phi, hphi) 1535 1536 !Arguments --------------------------- 1537 type(hexc_t),intent(in) :: hexc 1538 type(hexc_interp_t),intent(in) :: hexc_i 1539 complex(dpc),intent(in) :: phi(hexc%hsize) 1540 complex(dpc),intent(out) :: hphi(hexc%hsize) 1541 1542 !Local variables --------------------- 1543 integer :: ierr 1544 real(dp) :: tsec(2) 1545 1546 !***************************************************************************** 1547 1548 call timab(697,1,tsec) 1549 1550 if(hexc%BSp%use_interp) then 1551 hphi = hexc_i%diag_dense * phi 1552 1553 if (any(hexc%BSp%interp_mode == [2,3,4])) then 1554 ! hphi = hphi + MATMUL(hinterp,phi) 1555 call xgemv('N',hexc_i%hsize_dense,hexc_i%hsize_dense,cone,hexc_i%hinterp,hexc_i%hsize_dense,phi,1,cone,hphi,1) 1556 if (any(hexc%BSp%interp_mode == [2,4])) then 1557 call timab(697,2,tsec) 1558 return ! We are done 1559 end if 1560 end if 1561 1562 call hexc_interp_matmul(hexc%bsp, hexc%hsize_coarse, hexc_i%hsize_dense, hexc_i%all_hmat, phi, hphi, & 1563 & hexc_i%interpolator%double_grid,hexc%nbnd_coarse, hexc_i%interpolator, hexc_i%div2kdense, hexc_i%kdense2div) 1564 1565 else ! No interpolation 1566 call xgemv('N',hexc%hsize,hexc%my_nt,cone,hexc%hreso,hexc%hsize,phi,1,czero,hphi,1) 1567 call xmpi_sum(hphi,hexc%comm,ierr) 1568 end if 1569 1570 call timab(697,2,tsec) 1571 1572 end subroutine hexc_matmul_tda