TABLE OF CONTENTS
- ABINIT/gwls_hamiltonian
- gwls_hamiltonian/build_H
- gwls_hamiltonian/build_vxc
- gwls_hamiltonian/destroy_H
- gwls_hamiltonian/dft_xc_energy
- gwls_hamiltonian/DistributeValenceKernel
- gwls_hamiltonian/DistributeValenceWavefunctions
- gwls_hamiltonian/exchange
- gwls_hamiltonian/g_to_r
- gwls_hamiltonian/gr_to_g
- gwls_hamiltonian/Hpsik
- gwls_hamiltonian/Hpsikc
- gwls_hamiltonian/kbkb_to_kb
- gwls_hamiltonian/pc_k_valence_kernel
- gwls_hamiltonian/precondition
- gwls_hamiltonian/precondition_cplx
- gwls_hamiltonian/set_precondition
- gwls_hamiltonian/sqrt_vc_k
- gwls_hamiltonian/unset_precondition
- gwls_hamiltonian/wf_block_distribute
ABINIT/gwls_hamiltonian [ Modules ]
NAME
gwls_hamiltonian
FUNCTION
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (JLJ, BR, MC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 28 29 module gwls_hamiltonian 30 31 ! local modules 32 use m_gwls_utility 33 use gwls_wf 34 35 use m_bandfft_kpt 36 use m_cgtools 37 38 39 40 ! abinit modules 41 use defs_basis 42 use defs_datatypes 43 use m_pawcprj, only : pawcprj_type 44 use defs_abitypes 45 use m_profiling_abi 46 use m_xmpi 47 use m_pawang 48 use m_errors 49 use m_ab7_mixing 50 use m_mpinfo 51 52 use m_io_tools, only : get_unit 53 use m_dtset, only : dtset_copy, dtset_free 54 use m_hamiltonian, only : gs_hamiltonian_type, copy_hamiltonian, destroy_hamiltonian, & 55 & load_k_hamiltonian 56 use m_paw_dmft, only : paw_dmft_type 57 58 use m_vcoul, only : vcoul_t, vcoul_init, vcoul_free 59 use m_crystal, only : crystal_t, crystal_init, crystal_free, crystal_print 60 use m_io_kss, only : make_gvec_kss 61 use m_gsphere, only : gsphere_t, gsph_init, gsph_free, print_gsphere 62 use m_bz_mesh, only : kmesh_t, kmesh_init, kmesh_free, kmesh_print, find_qmesh 63 64 implicit none 65 save 66 private
gwls_hamiltonian/build_H [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
build_H
FUNCTION
.
INPUTS
OUTPUT
PARENTS
vtowfk
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1886 subroutine build_H(dtset2,mpi_enreg2,cpopt2,cg2,gs_hamk2,kg_k2,kinpw2) 1887 1888 !use m_bandfft_kpt 1889 use m_cgtools 1890 use m_wfutils 1891 1892 !This section has been created automatically by the script Abilint (TD). 1893 !Do not modify the following lines by hand. 1894 #undef ABI_FUNC 1895 #define ABI_FUNC 'build_H' 1896 use interfaces_53_ffts 1897 !End of the abilint section 1898 1899 implicit none 1900 1901 !Arguments of gw_sternheimer, reveived as argument by build_H------------------------- 1902 type(dataset_type), intent(in) :: dtset2 1903 type(MPI_type), intent(inout) :: mpi_enreg2 1904 type(gs_hamiltonian_type), intent(inout) :: gs_hamk2 1905 1906 integer, intent(in) :: cpopt2 1907 integer, intent(in) :: kg_k2(3,gs_hamk2%npw_k) 1908 !Since mcg = npw_k*nband when there is only gamma, the next line works. If there is not only Gamma, mcg is the proper size of cg. 1909 real(dp), intent(in) :: cg2(2,dtset2%mpw*dtset2%nspinor*dtset2%mband*dtset2%mkmem*dtset2%nsppol) 1910 real(dp), intent(in) :: kinpw2(gs_hamk2%npw_k) 1911 !Local variables : most of them are in the module for now. 1912 real(dp) :: commutation_error 1913 integer :: io_unit 1914 integer :: io_unit_debug 1915 integer :: ierr 1916 integer :: cplx 1917 integer :: j, k 1918 integer :: dimph3d 1919 integer :: mb 1920 1921 integer :: mpi_communicator 1922 1923 character(128):: filename_debug 1924 1925 1926 real(dp), allocatable :: wfk_tmp1(:,:) ,wfk_tmp2(:,:) 1927 1928 ! ************************************************************************* 1929 1930 ! Hartree-Fock cannot be used with GWLS. 1931 if(dtset2%usefock==1 .and. associated(gs_hamk2%fockcommon)) then 1932 MSG_ERROR(' build_H : Hartree-Fock option can not be used with optdriver==66 (GWLS calculations).') 1933 end if 1934 1935 !First we copy the data structure types 1936 call dtset_copy(dtset,dtset2) !dtset = dtset2 1937 1938 call copy_mpi_enreg(mpi_enreg2,mpi_enreg) 1939 1940 call copy_hamiltonian(gs_hamk,gs_hamk2) 1941 1942 !Then we copy the standard types 1943 cpopt = cpopt2 1944 npw_k = gs_hamk2%npw_k 1945 1946 ABI_ALLOCATE(cg,(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol)) 1947 cg = cg2 1948 1949 ABI_ALLOCATE(vlocal,(gs_hamk2%n4,gs_hamk2%n5,gs_hamk2%n6,gs_hamk2%nvloc)) 1950 vlocal = gs_hamk2%vlocal 1951 gs_hamk%vlocal => vlocal 1952 1953 ABI_ALLOCATE(kg_k,(3,npw_k)) 1954 kg_k = kg_k2 1955 ABI_ALLOCATE(kinpw,(npw_k)) 1956 kinpw = kinpw2 1957 1958 dimffnl=0; if (blocksize<=1) dimffnl=size(gs_hamk2%ffnl_k,2) 1959 ABI_ALLOCATE(ffnl,(npw_k,dimffnl,gs_hamk%lmnmax,gs_hamk%ntypat)) 1960 dimph3d=0; if (blocksize<=1) dimph3d=gs_hamk2%matblk 1961 ABI_ALLOCATE(ph3d,(2,npw_k,dimph3d)) 1962 1963 !Initializing variables from dataset 1964 nfft = dtset%nfft 1965 nline = dtset%nline 1966 tolwfr = dtset%tolwfr 1967 ngfft = dtset%ngfft 1968 mcg = dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 1969 nspinor=dtset%nspinor 1970 n1=ngfft(1) 1971 n2=ngfft(2) 1972 n3=ngfft(3) 1973 n4=ngfft(4) 1974 n5=ngfft(5) 1975 n6=ngfft(6) 1976 mgfft=maxval(ngfft(1:3)) 1977 nbandv = int(dtset%nelect)/2 1978 ABI_ALLOCATE(istwfk,(dtset%nkpt)) 1979 istwfk(:)=dtset%istwfk 1980 1981 !Initializing variables from gs_hamk 1982 ucvol = gs_hamk%ucvol 1983 ABI_ALLOCATE(gbound,(2*mgfft+8,2)) 1984 !gbound = gs_hamk%gbound_k !Must be done later for bandft paralelism 1985 1986 !Parameters which need to be set by hand for now... 1987 weight = 1 ! The weight of the k-pts, which sum to 1. 1988 tim_fourwf = 0 1989 ckpt = 1 ! The kpt of the wf to FFT. We only consider gamma point for now. 1990 1991 ! ndat = 1 used to be hard coded, which works fine for FFT only parallelism. 1992 ndat = 1 ! # of FFT to do in parallel in fourwf 1993 1994 eshift = 0.0 ! For PAW, opt_sij=-1 gives <G|H-lambda.S|C> in ghc instead of <G|H|C> 1995 sij_opt = 0 ! Option in PAW to tell getghc what to compute... 1996 tim_getghc = 0 ! timing code of the calling subroutine(can be set to 0 if not attributed) 1997 type_calc = 0 ! option governing which part of Hamitonian is to be applied: 0: whole Hamiltonian 1998 nband = dtset%mband 1999 ispden = 1 !When required, the spin index to be used. We don't support spin polarized calculations yet... 2000 2001 2002 2003 !======================================================================================================================== 2004 ! Trying to implement band parallelism (Bruno, 14/10/2013) 2005 !------------------------------------------------------------------------------------------------------------------------ 2006 ! 2007 ! Band parallelism is described in the article on the LOBPCG method by F. Bottin et al, 2008 ! "Large scale ab initio calculations based on three levels of parallelisation", Computational Materials Science 42 (2008) 329-336 2009 ! 2010 ! The article describes the ABINIT implementation specifically, suggesting that the information there is directly related 2011 ! to the code. 2012 ! 2013 ! The main points in order to understand band + FFT parallelism are: 2014 ! 2015 ! 1) The processors are arranged in a 2D Cartesian MPI topology, or dimensions M x P, with 2016 ! M = nproc_bands 2017 ! P = nproc_fft 2018 ! There are of course M x P processors 2019 ! 2020 ! 2) the wavefunction coefficients (namely the array cg) are distributed on these processors. Unfortunately, 2021 ! two different distribution schemes are necessary: one for linear algebra (matrix products, etc), and one 2022 ! consistent with the Goedecker FFTs. The code must thus change the data distribution back and forth. 2023 ! 2024 ! 3) the "linear algebra" distribution describes the array cg directly. The G vectors are distributed among all 2025 ! M x P processors, such that npw = npw_tot / (nproc_bands x nproc_fft). Every processor has information about 2026 ! all bands. Schematically: 2027 ! 2028 ! 2029 ! <- P = nproc_fft -> 2030 ! ---------------------------------- 2031 ! | n | n | n | n | ( for all band indices n) 2032 ! ^ | G_11 | ... | ... | G_P1 | 2033 ! | ---------------------------------- 2034 ! M = nproc_bands | n | n | n | n | 2035 ! | | G_1. | ... | ... | G_P. | 2036 ! v ---------------------------------- 2037 ! | n | n | n | n | 2038 ! | G_1M | ... | ... | G_PM | 2039 ! ---------------------------------- 2040 ! 2041 ! 2042 ! the LOBPCG algorithm acts on blocks of bands. Thus, although each processor has information about all bands, 2043 ! we must conceptually view the bands grouped in blocks of size M, as {c_{1}(G),...,c_{M}(G)}, 2044 ! {c_{M+1}(G)...,c_{2M}(G)}, ... 2045 ! 2046 ! With this conceptual grouping, for a given block each processor has M x NG/(M x P) coefficients, where NG is the 2047 ! total number of G-vectors. 2048 ! 2049 ! 4) The distribution of information above is not appropriate for parallel FFTs. The data for one block is thus 2050 ! transposed and takes the form 2051 ! 2052 ! <- P = nproc_fft -> 2053 ! ---------------------------------- 2054 ! | 1 | 1 | ... | 1 | 2055 ! ^ | G_1 | G_2 | ... | G_P | 2056 ! | ---------------------------------- 2057 ! M = nproc_bands | .. | ... | ... | ... | 2058 ! | | G_1 | G_2 | ... | G_P | 2059 ! v ---------------------------------- 2060 ! | M | M | ... | M | 2061 ! | G_1 | G_2 | ... | G_P | 2062 ! ---------------------------------- 2063 ! 2064 ! where the set of G vectors G_i = (G_{i1}, G_{i2}, ..., G_{iM}). Thus, a given processor has information about a single 2065 ! band, but more G vectors. Each processor has 2066 ! NG/(M x P) x M = NG/P coefficients, just like before. 2067 ! Also, it is clear that the information is only communicated in the columns of the diagram above, avoiding 2068 ! "all processors to all processors" communication. 2069 ! 2070 ! The data is now distributed properly to do parallel FFTs! Each row in the diagram above corresponds to FFTs done 2071 ! in parallel over nproc_fft processors, on a given band. There are M rows running thus in parallel! 2072 ! 2073 ! The underlying ABINIT routines are equiped to handle FFT parallelism, not band distributed parallelism. 2074 ! prep_getghc.F90 and lobpgcwf.F90 show how to rearange information to be able to apply basic ABINIT routines. 2075 ! 2076 ! The code below is inspired / guessed from lobpcgwf and prep_getghc. WE ASSUME THERE IS NO SPINORS! CODE SHOULD 2077 ! BE HEAVILY REVIEWED for k-points and spinors, if and when they come to be useful. 2078 ! 2079 !======================================================================================================================== 2080 2081 ! variables for band parallelism, following the language in lobpcgwf 2082 nspinor = 1 ! how many spinors are present. Hardcoded to 1 for now. 2083 blocksize = mpi_enreg%nproc_band*mpi_enreg%bandpp ! how many bands treated in a block; # of FFT done in parallel in getghc 2084 nbdblock = nband/blocksize ! how many blocks of bands are there 2085 ikpt_this_proc = 1 ! Assuming only 1 kpoint, for molecules. 2086 ! This will have to be reviewed for crystals! 2087 2088 npw_kb = npw_k*blocksize 2089 npw_g = gs_hamk2%npw_fft_k 2090 2091 if (blocksize>1) then 2092 kg_k_gather => bandfft_kpt(ikpt_this_proc)%kg_k_gather 2093 ffnl_gather => bandfft_kpt(ikpt_this_proc)%ffnl_gather 2094 ph3d_gather => bandfft_kpt(ikpt_this_proc)%ph3d_gather 2095 kinpw_gather => bandfft_kpt(ikpt_this_proc)%kinpw_gather 2096 else 2097 ffnl = gs_hamk2%ffnl_k 2098 ph3d = gs_hamk2%ph3d_k 2099 kg_k_gather => kg_k 2100 ffnl_gather => ffnl 2101 ph3d_gather => ph3d 2102 kinpw_gather => kinpw 2103 endif 2104 call load_k_hamiltonian(gs_hamk,kinpw_k=kinpw_gather,kg_k=kg_k_gather,ffnl_k=ffnl_gather,ph3d_k=ph3d_gather) 2105 2106 gbound = gs_hamk%gbound_k 2107 2108 !Set up wfk routines 2109 call set_wf(ucvol/nfft,n4,n5,n6,npw_k,blocksize,npw_g,mpi_enreg%comm_bandfft,mpi_enreg%comm_fft,mpi_enreg%comm_band) 2110 2111 ! prepare the valence wavefunctions and projection operator to work in band+FFT parallel 2112 call DistributeValenceWavefunctions() 2113 call DistributeValenceKernel() 2114 2115 cplx = 2 ! wavefunctions have complex coefficients 2116 2117 ! Initialize the dummy variables 2118 ABI_DATATYPE_ALLOCATE(conjgrprj,(0,0)) 2119 ABI_ALLOCATE(dummy2,(0,0)) 2120 ABI_ALLOCATE(dummy3,(0,0,0)) 2121 2122 !Initialisation of the total counter for iterations of SQMR : 2123 ktot = 0 2124 2125 !Allocations of working arrays for the module 2126 !- for pc_k function 2127 ABI_ALLOCATE(scprod2,(2,nband)) 2128 2129 !- for precondition function (and set_precondition subroutine) 2130 ABI_ALLOCATE(pcon,(npw_g)) 2131 pcon = one 2132 2133 !- for (private) working wf 2134 ABI_ALLOCATE(psik1,(2,npw_k)) 2135 ABI_ALLOCATE(psik2,(2,npw_k)) 2136 ABI_ALLOCATE(psik3,(2,npw_k)) 2137 ABI_ALLOCATE(psik4,(2,npw_k)) 2138 ABI_ALLOCATE(psikb1,(2,npw_kb)) 2139 ABI_ALLOCATE(psikb2,(2,npw_kb)) 2140 ABI_ALLOCATE(psikb3,(2,npw_kb)) 2141 ABI_ALLOCATE(psikb4,(2,npw_kb)) 2142 ABI_ALLOCATE(psig1,(2,npw_g)) 2143 ABI_ALLOCATE(psig2,(2,npw_g)) 2144 ABI_ALLOCATE(psig3,(2,npw_g)) 2145 ABI_ALLOCATE(psig4,(2,npw_g)) 2146 ABI_ALLOCATE(psir1,(2,n4,n5,n6)) 2147 ABI_ALLOCATE(psir2,(2,n4,n5,n6)) 2148 ABI_ALLOCATE(psir3,(2,n4,n5,n6)) 2149 ABI_ALLOCATE(psidg,(2,nfft)) 2150 ABI_ALLOCATE(denpot,(2*n4,n5,n6)) 2151 2152 psir1 = zero 2153 psir2 = zero 2154 psir3 = zero 2155 denpot = zero 2156 2157 !Construct the vector of eigenvalues and write then to std output 2158 ABI_ALLOCATE(eig,(nband)) 2159 eig=zero 2160 2161 write(std_out,*) ch10,"Eigenvalues computation check, routine build_H:",ch10 2162 io_unit_debug = get_unit() 2163 write(filename_debug,'(A,I0.4,A)') "DEBUG_PROC=",mpi_enreg%me,".dat" 2164 2165 2166 2167 mpi_communicator = mpi_enreg%comm_bandfft 2168 2169 open(file=filename_debug,status=files_status_old,unit=io_unit_debug) 2170 2171 write(io_unit_debug,'(A)') " Parameters:" 2172 write(io_unit_debug,'(A,I5)') " nband = ", nband 2173 write(io_unit_debug,'(A,I5)') " blocksize = ", blocksize 2174 write(io_unit_debug,'(A,I5)') " npw_k = ", npw_k 2175 write(io_unit_debug,'(A,I5)') " nbdblock = ", nbdblock 2176 2177 ! temporary wavefunction array, for data in the "linear algebra" distribution 2178 ABI_ALLOCATE( wfk_tmp1, (2,npw_k)) 2179 ABI_ALLOCATE( wfk_tmp2, (2,npw_k)) 2180 2181 do n=1, nband 2182 ! Extract i^t/h wavefunction 2183 wfk_tmp1(:,1:npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k) 2184 2185 ! Are the wavefunctions normalized? 2186 !tmpc = cg_zdotc(npw_k,wfk_tmp1,wfk_tmp1) 2187 !call xmpi_sum(tmpc,mpi_enreg%comm_bandfft,ierr) ! sum on all processors 2188 2189 ! DEBUGGING CODE 2190 write(std_out,'(A,I5,A,2F24.16)') "band ", n, ", <psi_n | psi_n > =", norm_k(wfk_tmp1) 2191 write(io_unit_debug,'(A,I5,A,2F24.16)') "band ", n, ", <psi_n | psi_n > =", norm_k(wfk_tmp1) 2192 flush(io_unit_debug) 2193 end do 2194 2195 2196 ! loop on blocks of bands. All bands in a given block will be done in parallel! 2197 2198 ! loop on block of bands 2199 do iblock = 1, nbdblock 2200 2201 ! loop on states for this block 2202 do iband = 1, blocksize 2203 2204 n = (iblock-1)*blocksize+iband 2205 2206 psikb1(:,(iband-1)*npw_k+1:iband*npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k) 2207 2208 end do 2209 2210 ! change configuration of the data 2211 call wf_block_distribute(psikb1, psig1,1) ! LA -> FFT 2212 2213 ! Apply hamiltonian on wavefunction. In bandFFT parallel, Hpsik calls prep_getghc, which knows how to transform the 2214 ! distribution of the data from "linear algebra"-like to "FFT"-like. The output is *probably* in the "linear algebra"-like 2215 ! distribution. 2216 call Hpsik(psig2,psig1) 2217 2218 ! change configuration of the data 2219 call wf_block_distribute(psikb2, psig2,2) ! LA -> FFT 2220 2221 ! extract the wavefunctions band by band 2222 do iband=1, blocksize 2223 2224 n = blocksize*(iblock-1) + iband ! band index 2225 2226 wfk_tmp1(:,1:npw_k) = psikb1(:,(iband-1)*npw_k+1:iband*npw_k) 2227 wfk_tmp2(:,1:npw_k) = psikb2(:,(iband-1)*npw_k+1:iband*npw_k) 2228 2229 tmpc = cg_zdotc(npw_k,wfk_tmp1,wfk_tmp2) 2230 2231 call xmpi_sum(tmpc,mpi_enreg%comm_bandfft,ierr) ! sum on all processors 2232 eig(n) = tmpc(1) 2233 2234 write(std_out,'(A,I5,A,F24.16,A)') "(build_H) band ", n, ", eig =", eig(n), " Ha." 2235 2236 ! DEBUGGING CODE 2237 write(io_unit_debug,'(A,I5,A,F24.16,A)') "band ", n, ", eig =", eig(n), " Ha." 2238 flush(io_unit_debug) 2239 flush(io_unit_debug) 2240 end do 2241 2242 end do 2243 2244 ABI_DEALLOCATE( wfk_tmp1 ) 2245 ABI_DEALLOCATE( wfk_tmp2 ) 2246 close(io_unit_debug) 2247 2248 2249 2250 ! Finishing the construction of vxc (transcribing it from the double real grid (for the density) 2251 ! to the single real grid (for the wfs). 2252 ! Assummes we only need one spin component; one transcription per spin being needed. 2253 if(allocated(vxc_dg)) then 2254 ABI_ALLOCATE(vxc,(n4,n5,n6,dtset%nspden)) 2255 vxc = zero 2256 call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vxc_dg(:,ispden),vxc(:,:,:,ispden),2) 2257 end if 2258 2259 timrev = 1 !Assumes istwfk *1. 1:time reversal symmetry NOT present | 2:" " " present 2260 ABI_ALLOCATE(title,(dtset%ntypat)) 2261 do i=1,dtset%ntypat 2262 title(i) = "Bloup" ! The clean way would be to get the psps structure in this module 2263 ! (build_vxc is called from a place in GS calculations where it is available; 2264 ! should be the easiest way). For now, this allows the code to run. 2265 end do 2266 call crystal_init(dtset%amu_orig(:,1),Cryst,dtset%spgroup,dtset%natom,dtset%npsp,& 2267 & dtset%ntypat,dtset%nsym,dtset%rprimd_orig(:,:,1),dtset%typat,& 2268 & dtset%xred_orig(:,:,1),dtset%ziontypat,dtset%znucl,timrev,.false.,.false.,title,& 2269 & dtset%symrel,dtset%tnons,dtset%symafm) 2270 ABI_DEALLOCATE(title) 2271 call crystal_print(Cryst) 2272 2273 !TODO : Should be put in a separate build_vc constructor, and should be called right after build_H in the context of optdriver 66. 2274 if(dtset%optdriver==66) then 2275 2276 !Set up of the k-points and tables in the whole BZ 2277 call kmesh_init(Kmesh,Cryst,dtset%nkpt,dtset%kptns,dtset%kptopt,wrap_1zone=.false.) 2278 call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out) 2279 call find_qmesh(Qmesh,Cryst,Kmesh) 2280 call kmesh_print(Qmesh,"Q-mesh for the screening function",std_out) 2281 2282 2283 !------------------------------ 2284 !Building the vc_sqrt structure 2285 !------------------------------ 2286 ! The sqrt_vc code relies on a different parallelism scheme than Vanilla ABINIT. 2287 ! The Gsphere object must be built on a single processor having all the G-vectors 2288 ! available. 2289 2290 !gsph_init need gvec built according to the KSS convention, so that kg_k is not suitable (not properly sorted). 2291 npw_serial=npw_k 2292 call xmpi_sum(npw_serial,mpi_enreg%comm_bandfft,ierr) 2293 2294 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 2295 call make_gvec_kss(dtset%nkpt,dtset%kptns,ecut_eff,dtset%symmorphi,dtset%nsym,dtset%symrel,dtset%tnons,Cryst%gprimd,& 2296 & dtset%prtvol,npw_serial,gvec,ierr) 2297 2298 call gsph_init(Gsphere,Cryst,npw_serial,gvec=gvec) 2299 2300 call print_gsphere(Gsphere) 2301 2302 call vcoul_init(Vcp,Gsphere,Cryst,Qmesh,Kmesh,dtset%rcut,dtset%icutcoul,dtset%vcutgeo,dtset%ecutsigx,npw_serial,& 2303 & dtset%nkpt,dtset%kptns,dtset%ngfft,mpi_enreg%comm_world) 2304 2305 ! Since Vcp%vc_sqrt is sorted according to the KSS convention for G vectors 2306 ! BUT will be applied to GS wavefunctions (where G vectors are sorted otherwise) 2307 ! It is necessary to construct a vc_sqrt vector with the GS sorting. 2308 ! Moreover, if we are in parallel (over band / FFTs), Vcp%vc_sqrt is the serial version 2309 ! and need to be distributed according to the GS scheme (Linear Algebra configuration). 2310 ABI_ALLOCATE(vc_sqrt,(npw_k)) 2311 vc_sqrt=zero 2312 k=0 2313 do i=1,npw_k 2314 do j=1,npw_serial 2315 if(all(kg_k(:,i)==gvec(:,j))) k=j 2316 end do 2317 vc_sqrt(i)=Vcp%vc_sqrt(k,1) 2318 end do 2319 2320 end if 2321 2322 !-------------------------------------------------------------------------------- 2323 ! Security check : The eigenstates of the projector and the hamiltonian need to 2324 ! agree down to the precision requested (tolwfr). Otherwise, SQMR is doomed. 2325 ! Now that the density is read and the eigenstates are calculated in the GW run, 2326 ! it is less useful. A test on tolwfr could be sufficient. 2327 ! 2328 ! This remains a good tool to debug band+fft parallelism 2329 !-------------------------------------------------------------------------------- 2330 2331 ! only write on the head node! 2332 if (mpi_enreg%me == 0) then 2333 io_unit = get_unit() 2334 open(file='build_H.log',status=files_status_old,unit=io_unit) 2335 write(io_unit,10) '#----------------------------------------------------------------------------' 2336 write(io_unit,10) '#' 2337 write(io_unit,10) '# This file presents the results of a small test to check ' 2338 write(io_unit,10) '# how well the Hamiltonian commutes with the projection ' 2339 write(io_unit,10) '# operator. ' 2340 write(io_unit,10) '#' 2341 write(io_unit,10) '# Definitions:' 2342 write(io_unit,10) '#' 2343 write(io_unit,10) '# P : projections on conduction states' 2344 write(io_unit,10) '# H : Hamiltonian operator' 2345 write(io_unit,10) '# | psi > : eigenstate' 2346 write(io_unit,10) '#' 2347 flush(io_unit) 2348 end if 2349 2350 psikb1 = zero 2351 2352 ! sum all valence states, on copy in every band block. 2353 do i=1,nbandv 2354 2355 do mb = 1, blocksize 2356 2357 psikb1(:,(mb-1)*npw_k+1:mb*npw_k) = psikb1(:,(mb-1)*npw_k+1:mb*npw_k) + cg(:,(i-1)*npw_k+1:i*npw_k) 2358 2359 end do 2360 2361 end do 2362 2363 ! normalize to one! Only sum the fist band block 2364 tmpc = cg_zdotc(npw_k,psikb1,psikb1) 2365 call xmpi_sum(tmpc,mpi_communicator ,ierr) ! sum on all processors 2366 psikb1 = psikb1/sqrt(tmpc(1)) 2367 2368 2369 ! change data distribution 2370 call wf_block_distribute(psikb1, psig1,1) ! LA -> FFT 2371 2372 ! Apply P.H operator 2373 call Hpsik(psig2 ,psig1) 2374 call pc_k_valence_kernel(psig2) 2375 2376 ! Apply H.P operator 2377 call pc_k_valence_kernel(psig1) 2378 call Hpsik(psig1) 2379 2380 ! compute error 2381 psig3 = psig1 - psig2 ! explicitly write the difference in an array 2382 tmpc = cg_zdotc(npw_g,psig3,psig3) 2383 commutation_error = tmpc(1) 2384 2385 call xmpi_sum(commutation_error , mpi_enreg%comm_fft, ierr) ! sum on all processors working on FFT! 2386 2387 ! only write on the head node! 2388 if (mpi_enreg%me == 0) then 2389 write(io_unit,20) ' || (PH -HP) |b> ||^2 = ',commutation_error 2390 write(io_unit,20) ' tolwfr = ',tolwfr 2391 end if 2392 2393 if(commutation_error > tolwfr) then 2394 !write(io_unit,10) '# || (PH -HP) |b> ||^2 > tolwfr ==> Decision taken exit!' 2395 2396 if (mpi_enreg%me == 0) write(io_unit,10) '# || (PH -HP) |b> ||^2 > tolwfr ==> This must be fixed!' 2397 2398 write(std_out,20) "WARNING-build_H: The tolerance tolwfr=",tolwfr 2399 write(std_out,20) " is smaller than the commutation error ||(PH-HP)|b>||^2=",commutation_error 2400 write(std_out,10) " Either set tolwfr to a less stringent value in this calculation" 2401 write(std_out,10) " or to a more stringent value in the wavefunction calculation." 2402 end if 2403 2404 if (mpi_enreg%me == 0) close(io_unit) 2405 2406 10 format(A) 2407 20 format(A,ES12.3) 2408 end subroutine build_H
gwls_hamiltonian/build_vxc [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
build_vxc
FUNCTION
.
INPUTS
OUTPUT
PARENTS
scfcv
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1676 subroutine build_vxc(vxc2,nfft2,nspden2) 1677 !Only transcribe the argument vxc2 in the module; the change from dg to sg is done in build_H (and stored in vxc), since the 1678 !arguments of fftpac are built in build_H. 1679 1680 !This section has been created automatically by the script Abilint (TD). 1681 !Do not modify the following lines by hand. 1682 #undef ABI_FUNC 1683 #define ABI_FUNC 'build_vxc' 1684 !End of the abilint section 1685 1686 implicit none 1687 1688 !We need the dimensions of vxc since they don't exist yet in the module; build_vxc being called before build_H. 1689 integer, intent(in) :: nfft2, nspden2 1690 real(dp), intent(in) :: vxc2(nfft2,nspden2) 1691 1692 ! ************************************************************************* 1693 1694 ABI_ALLOCATE(vxc_dg,(nfft2,nspden2)) 1695 vxc_dg = vxc2 1696 1697 end subroutine build_vxc
gwls_hamiltonian/destroy_H [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
destroy_H
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_sternheimer
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1723 subroutine destroy_H 1724 1725 1726 !This section has been created automatically by the script Abilint (TD). 1727 !Do not modify the following lines by hand. 1728 #undef ABI_FUNC 1729 #define ABI_FUNC 'destroy_H' 1730 !End of the abilint section 1731 1732 implicit none 1733 1734 call dtset_free(dtset) 1735 call destroy_hamiltonian(gs_hamk) 1736 1737 call crystal_free(Cryst) 1738 call kmesh_free(Kmesh) 1739 call kmesh_free(Qmesh) 1740 call gsph_free(Gsphere) 1741 call vcoul_free(Vcp) 1742 1743 call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg) 1744 call destroy_mpi_enreg(mpi_enreg) 1745 1746 !NOTE : the syntax if(allocated(a)) ABI_DEALLOCATE(a) result in an error if "a" is not allocated; since the macro replace 1747 !ABI_ALLOCATE by more than one line of text, the second lines and up get outside the if... if() then syntax is equired. 1748 if(allocated(cg)) then 1749 ABI_DEALLOCATE(cg) 1750 end if 1751 if(allocated(gbound)) then 1752 ABI_DEALLOCATE(gbound) 1753 end if 1754 if(allocated(kg_k)) then 1755 ABI_DEALLOCATE(kg_k) 1756 end if 1757 if(allocated(ffnl)) then 1758 ABI_DEALLOCATE(ffnl) 1759 end if 1760 if(allocated(ph3d)) then 1761 ABI_DEALLOCATE(ph3d) 1762 end if 1763 if(allocated(kinpw)) then 1764 ABI_DEALLOCATE(kinpw) 1765 end if 1766 if(allocated(vxc)) then 1767 ABI_DEALLOCATE(vxc) 1768 end if 1769 if(allocated(vlocal)) then 1770 ABI_DEALLOCATE(vlocal) 1771 end if 1772 if(allocated(conjgrprj)) then 1773 ABI_DATATYPE_DEALLOCATE(conjgrprj) 1774 end if 1775 if(allocated(istwfk)) then 1776 ABI_DEALLOCATE(istwfk) 1777 end if 1778 if(allocated(dummy2)) then 1779 ABI_DEALLOCATE(dummy2) 1780 end if 1781 if(allocated(dummy3)) then 1782 ABI_DEALLOCATE(dummy3) 1783 end if 1784 if(allocated(eig)) then 1785 ABI_DEALLOCATE(eig) 1786 end if 1787 if(allocated(scprod2)) then 1788 ABI_DEALLOCATE(scprod2) 1789 end if 1790 if(allocated(pcon)) then 1791 ABI_DEALLOCATE(pcon) 1792 end if 1793 if(allocated(psik1)) then 1794 ABI_DEALLOCATE(psik1) 1795 end if 1796 if(allocated(psik2)) then 1797 ABI_DEALLOCATE(psik2) 1798 end if 1799 if(allocated(psik3)) then 1800 ABI_DEALLOCATE(psik3) 1801 end if 1802 if(allocated(psik4)) then 1803 ABI_DEALLOCATE(psik4) 1804 end if 1805 if(allocated(psikb1)) then 1806 ABI_DEALLOCATE(psikb1) 1807 end if 1808 if(allocated(psikb2)) then 1809 ABI_DEALLOCATE(psikb2) 1810 end if 1811 if(allocated(psikb3)) then 1812 ABI_DEALLOCATE(psikb3) 1813 end if 1814 if(allocated(psikb4)) then 1815 ABI_DEALLOCATE(psikb4) 1816 end if 1817 if(allocated(psig1)) then 1818 ABI_DEALLOCATE(psig1) 1819 end if 1820 if(allocated(psig2)) then 1821 ABI_DEALLOCATE(psig2) 1822 end if 1823 if(allocated(psig3)) then 1824 ABI_DEALLOCATE(psig3) 1825 end if 1826 if(allocated(psig4)) then 1827 ABI_DEALLOCATE(psig4) 1828 end if 1829 if(allocated(psir1)) then 1830 ABI_DEALLOCATE(psir1) 1831 end if 1832 if(allocated(psir2)) then 1833 ABI_DEALLOCATE(psir2) 1834 end if 1835 if(allocated(psir3)) then 1836 ABI_DEALLOCATE(psir3) 1837 end if 1838 if(allocated(psidg)) then 1839 ABI_DEALLOCATE(psidg) 1840 end if 1841 if(allocated(vxc_dg)) then 1842 ABI_DEALLOCATE(vxc_dg) 1843 end if 1844 if(allocated(denpot)) then 1845 ABI_DEALLOCATE(denpot) 1846 end if 1847 if(allocated(kernel_wavefunctions_FFT)) then 1848 ABI_DEALLOCATE(kernel_wavefunctions_FFT) 1849 end if 1850 if(allocated(valence_wavefunctions_FFT)) then 1851 ABI_DEALLOCATE(valence_wavefunctions_FFT) 1852 end if 1853 if(associated(gvec)) then 1854 ABI_DEALLOCATE(gvec) 1855 end if 1856 if(allocated(vc_sqrt)) then 1857 ABI_DEALLOCATE(vc_sqrt) 1858 end if 1859 1860 end subroutine destroy_H
gwls_hamiltonian/dft_xc_energy [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
dft_xc_energy
FUNCTION
.
INPUTS
OUTPUT
PARENTS
CHILDREN
SOURCE
805 function dft_xc_energy(e) 806 807 808 !This section has been created automatically by the script Abilint (TD). 809 !Do not modify the following lines by hand. 810 #undef ABI_FUNC 811 #define ABI_FUNC 'dft_xc_energy' 812 use interfaces_53_ffts 813 !End of the abilint section 814 815 implicit none 816 real(dp) :: dft_xc_energy 817 integer, intent(in) :: e 818 819 820 integer :: cplex, option, ierr 821 real(dp), allocatable :: psik_e(:,:) !Working array to store the wavefunction 822 real(dp), allocatable :: psik_e_alltoall(:,:) !Working array to store the wavefunction 823 824 real(dp), allocatable :: psik_out(:,:) !Working array to store the wavefunction 825 826 real(dp) :: tmpc(2) 827 integer :: mpi_communicator 828 real(dp) :: dft_xc_energy_tmp 829 830 ! ************************************************************************* 831 832 !-------------------------------------------------------------------------------- 833 ! The goal of this routine is to compute the xc energy using 834 ! the DFT Vxc array. 835 !-------------------------------------------------------------------------------- 836 837 ! Parallel-MPI code. This is inspired from code within the getghc routine, which 838 ! applies a complex potential to a wavefunction. 839 840 841 cplex = 1 ! real potential 842 option = 2 ! multiply wavefunction by potential 843 nspinor= 1 844 ! Allocate wavefunction arrays which contains the coefficients of the e-state 845 ABI_ALLOCATE(psik_e, (2,npw_kb)) 846 ABI_ALLOCATE(psik_e_alltoall,(2,npw_g)) 847 ABI_ALLOCATE(psik_out, (3,npw_g)) 848 849 850 ! Only fetch the e-state, setting other states in block to zero! 851 psik_e(:,:) = zero 852 psik_e(:,1:npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k) 853 psik_out(:,:) = zero 854 855 856 ! change the configuration of the data 857 call wf_block_distribute(psik_e, psik_e_alltoall,1) ! LA -> FFT 858 859 860 ! Call fourwf to generate the product, in k-space 861 ! Computation: 862 ! psik_e_alltoall(k) -> psi(r) 863 ! res(r) = (vxc(r) x psi(r)) 864 ! psik_out(k) <- res(r) 865 ! psir3 is a dummy, not used here. 866 867 call fourwf(cplex,vxc(:,:,:,ispden),psik_e_alltoall,psik_out,psir3,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,& 868 & mpi_enreg,1,ngfft,npw_g,npw_g,n4,n5,n6,option,dtset%paral_kgb,tim_fourwf,weight,weight) 869 870 871 tmpc = cg_zdotc(npw_g, psik_e_alltoall,psik_out) 872 873 mpi_communicator = mpi_enreg%comm_fft 874 call xmpi_sum(tmpc,mpi_communicator, ierr) ! sum on all processors working on FFT! 875 876 dft_xc_energy_tmp = tmpc(1) 877 878 mpi_communicator = mpi_enreg%comm_band 879 call xmpi_sum(dft_xc_energy_tmp,mpi_communicator, ierr) ! sum on all processors working on FFT! 880 881 dft_xc_energy = dft_xc_energy_tmp 882 883 ABI_DEALLOCATE(psik_e) 884 ABI_DEALLOCATE(psik_e_alltoall) 885 ABI_DEALLOCATE(psik_out) 886 887 end function dft_xc_energy
gwls_hamiltonian/DistributeValenceKernel [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
DistributeValenceKernel
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_hamiltonian
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
321 subroutine DistributeValenceKernel() 322 !-------------------------------------------------------------------------------- 323 ! 324 ! This subroutine distributes, once and for all, the kernel of the static 325 ! SQMR operator. 326 ! 327 ! In this first, quick and dirty implementation, we simply distribute 328 ! ALL the valence bands on ALL the FFT groups. This is not efficient, but 329 ! kernel projections are not a bottleneck of the computation. 330 ! 331 ! A better (forthcoming) algorithm would only distribute the actual kernel, 332 ! not all valence bands. 333 !-------------------------------------------------------------------------------- 334 335 !This section has been created automatically by the script Abilint (TD). 336 !Do not modify the following lines by hand. 337 #undef ABI_FUNC 338 #define ABI_FUNC 'DistributeValenceKernel' 339 !End of the abilint section 340 341 implicit none 342 343 integer :: mb, n 344 345 346 real(dp), allocatable :: psik_n(:,:) !wavefunctions in LA format 347 real(dp), allocatable :: psik_n_alltoall(:,:) !wavefunctions in FFT format 348 349 ! ************************************************************************* 350 351 352 !=================================================================== 353 ! Allocate the global array which will contain the valence states 354 !=================================================================== 355 ABI_ALLOCATE( kernel_wavefunctions_FFT, (2,npw_g, nband)) 356 357 kernel_wavefunctions_FFT = zero 358 359 !==================================================== 360 ! Allocate working arrays 361 !==================================================== 362 363 ABI_ALLOCATE(psik_n, (2,npw_kb)) 364 ABI_ALLOCATE(psik_n_alltoall,(2,npw_g)) 365 366 367 ! loop on all valence states 368 do n = 1, nband 369 370 ! Copy multiple instances of this valence state in the array; 371 ! this way, each FFT group will have ALL the valence states! 372 do mb = 1, blocksize 373 psik_n(:,(mb-1)*npw_k+1:mb*npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k) 374 end do 375 376 ! change configuration of the data 377 call wf_block_distribute(psik_n, psik_n_alltoall,1) ! LA -> FFT 378 379 ! copy data in global array 380 kernel_wavefunctions_FFT(:,:,n) = psik_n_alltoall(:,:) 381 382 end do 383 384 !==================================================== 385 ! cleanup 386 !==================================================== 387 388 ABI_DEALLOCATE(psik_n) 389 ABI_DEALLOCATE(psik_n_alltoall) 390 391 end subroutine DistributeValenceKernel
gwls_hamiltonian/DistributeValenceWavefunctions [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
DistributeValenceWavefunctions
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_hamiltonian
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
218 subroutine DistributeValenceWavefunctions() 219 !-------------------------------------------------------------------------------- 220 ! 221 ! This subroutine distributes, once and for all, the valence wavefunctions to 222 ! the FFT configuration. They will thus be ready to be used within the 223 ! susceptibility operator. 224 ! 225 !-------------------------------------------------------------------------------- 226 227 !This section has been created automatically by the script Abilint (TD). 228 !Do not modify the following lines by hand. 229 #undef ABI_FUNC 230 #define ABI_FUNC 'DistributeValenceWavefunctions' 231 !End of the abilint section 232 233 implicit none 234 235 integer :: iblk, mb, v 236 237 238 real(dp), allocatable :: psik_v(:,:) !wavefunctions in LA format 239 real(dp), allocatable :: psik_v_alltoall(:,:) !wavefunctions in FFT format 240 241 ! ************************************************************************* 242 243 244 245 !=================================================================== 246 ! Allocate the global array which will contain the valence states 247 !=================================================================== 248 ABI_ALLOCATE( valence_wavefunctions_FFT, (2,npw_g,nbdblock)) 249 250 valence_wavefunctions_FFT = zero 251 252 253 !==================================================== 254 ! Allocate working arrays 255 !==================================================== 256 257 ABI_ALLOCATE(psik_v, (2,npw_kb)) 258 ABI_ALLOCATE(psik_v_alltoall,(2,npw_g)) 259 260 261 262 ! loop on all blocks of states, 263 do iblk = 1, nbdblock 264 265 ! loop on valence states for this block; if the state is conduction, fill with zeros 266 do mb = 1, blocksize 267 268 v = (iblk-1)*blocksize+mb 269 270 if (v <= nbandv) then 271 psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = cg(:,(v-1)*npw_k+1:v*npw_k) 272 else 273 psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = zero 274 end if 275 276 end do 277 278 ! change configuration of the data 279 call wf_block_distribute(psik_v, psik_v_alltoall,1) ! LA -> FFT 280 281 ! copy data in global array 282 283 valence_wavefunctions_FFT(:,:,iblk) = psik_v_alltoall(:,:) 284 285 end do 286 287 !==================================================== 288 ! cleanup 289 !==================================================== 290 291 ABI_DEALLOCATE(psik_v) 292 ABI_DEALLOCATE(psik_v_alltoall) 293 294 295 end subroutine DistributeValenceWavefunctions
gwls_hamiltonian/exchange [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
exchange
FUNCTION
.
INPUTS
OUTPUT
PARENTS
CHILDREN
SOURCE
656 function exchange(e, Lbasis_lanczos) 657 658 !use m_bandfft_kpt 659 use m_cgtools 660 !================================================================================ 661 ! This subroutine computes the exchange energy in band+FFT parallel 662 ! 663 !================================================================================ 664 665 !This section has been created automatically by the script Abilint (TD). 666 !Do not modify the following lines by hand. 667 #undef ABI_FUNC 668 #define ABI_FUNC 'exchange' 669 !End of the abilint section 670 671 implicit none 672 real(dp) :: exchange 673 674 integer, intent(in) :: e 675 676 ! If these arguments are provided, the exchange energy is to be projected on this subspace 677 complex(dpc), optional, intent(in) :: Lbasis_lanczos(:,:) ! complex array which contains the Lanczos basis 678 679 real(dp), allocatable :: psik_e(:,:) !Working array to store the wavefunction 680 681 real(dp), allocatable :: psik_v(:,:) !Working array to store the wavefunction 682 683 real(dp), allocatable :: psik_out(:,:) !Working array to store the wavefunction 684 685 686 integer :: iblk, mb 687 688 integer :: l, lmax 689 690 real(dp) :: tmpc(2) 691 692 logical :: project 693 694 ! ************************************************************************* 695 696 !-------------------------------------------------------------------------------- 697 ! Determine if the exhcange energy must be projected on the Lanczos basis 698 ! a truncated Coulomb potential 699 !-------------------------------------------------------------------------------- 700 701 project = .false. 702 if (present(Lbasis_lanczos)) then 703 project = .true. 704 lmax = size(Lbasis_lanczos, 2) 705 end if 706 707 !-------------------------------------------------------------------------------- 708 ! The goal of this routine is to compute the exchange energy using 709 ! a truncated Coulomb potential 710 !-------------------------------------------------------------------------------- 711 712 exchange = 0.0_dp 713 714 !==================================================== 715 ! build the block of wavefunctions which will 716 ! contain copies the e-state 717 !==================================================== 718 719 ABI_ALLOCATE(psik_e,(2,npw_kb)) 720 721 ! fill psik_e with as many copies of the e-state as there are band processors; 722 ! that way, upon LA -> FFT, each row of fft processors will have the e-state 723 do v = 1, blocksize 724 psik_e(:,(v-1)*npw_k+1:v*npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k) 725 end do 726 727 !==================================================== 728 ! Allocate the block of wavefunctions which will 729 ! contain the valence states 730 !==================================================== 731 732 ABI_ALLOCATE(psik_v, (2,npw_kb)) 733 ABI_ALLOCATE(psik_out, (2,npw_kb)) 734 735 ! loop on all blocks of states, 736 do iblk = 1, nbdblock 737 738 ! loop on valence states for this block; if the state is conduction, fill with zeros 739 do mb = 1, blocksize 740 741 v = (iblk-1)*blocksize+mb 742 743 if (v <= nbandv) then 744 psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = cg(:,(v-1)*npw_k+1:v*npw_k) 745 else 746 psik_v(:,(mb-1)*npw_k+1:mb*npw_k) = zero 747 end if 748 749 end do 750 751 call kbkb_to_kb(psik_out,psik_v,psik_e) 752 753 ! apply Coulomb potential, and take norm: cumulate the exchange energy 754 do mb = 1, blocksize 755 756 psik1 = psik_out(:,(mb-1)*npw_k+1:mb*npw_k) 757 call sqrt_vc_k(psik1) 758 759 if (project) then 760 ! project on the Lanczos basis 761 do l = 1, lmax 762 psik2(1,:) = dble(Lbasis_lanczos(:,l)) 763 psik2(2,:) = dimag(Lbasis_lanczos(:,l)) 764 tmpc = scprod_k(psik2,psik1) 765 exchange = exchange - (tmpc(1)**2+tmpc(2)**2) 766 end do 767 else 768 ! compute the "exact" exchange energy 769 exchange = exchange - norm_k(psik1)**2 770 end if 771 772 end do 773 774 end do 775 776 ABI_DEALLOCATE(psik_e) 777 778 ABI_DEALLOCATE(psik_v) 779 780 ABI_DEALLOCATE(psik_out) 781 782 end function exchange
gwls_hamiltonian/g_to_r [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
g_to_r
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_DielectricArray,gwls_GWlanczos,gwls_GenerateEpsilon gwls_LanczosBasis,gwls_hamiltonian,gwls_model_polarisability gwls_polarisability
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1459 subroutine g_to_r(psi_out,psi_in) 1460 1461 1462 !This section has been created automatically by the script Abilint (TD). 1463 !Do not modify the following lines by hand. 1464 #undef ABI_FUNC 1465 #define ABI_FUNC 'g_to_r' 1466 use interfaces_53_ffts 1467 !End of the abilint section 1468 1469 implicit none 1470 real(dp), intent(out) :: psi_out(2,n4,n5,n6) 1471 real(dp), intent(in) :: psi_in(2,npw_g) 1472 integer :: option, cplex 1473 integer :: nproc_fft, me_fft, nd3 !, ierr 1474 1475 ! ************************************************************************* 1476 1477 option = 0 ! fft wavefunction to real space 1478 cplex = 2 ! complex potential 1479 1480 psig4 = psi_in 1481 psi_out = zero 1482 call fourwf(cplex,dummy3,psig4,dummy2,psi_out,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,mpi_enreg, & 1483 1,ngfft,npw_g,npw_g,n4,n5,n6,option,dtset%paral_kgb,tim_fourwf,weight,weight) 1484 psi_out = psi_out/sqrt(ucvol) 1485 1486 !! This comes from prep_fourwf 1487 !nproc_fft=mpi_enreg%nproc_fft 1488 !if (nproc_fft>1) then 1489 ! me_fft=mpi_enreg%me_fft 1490 ! if (me_fft>0) then 1491 ! nd3=(ngfft(3)-1)/nproc_fft+1 1492 ! psi_out(:,:,:,me_fft*nd3+1:me_fft*nd3+nd3)=psi_out(:,:,:,1:nd3) 1493 ! psi_out(:,:,:,1:nd3)=zero 1494 ! end if 1495 ! call xmpi_sum(psi_out,mpi_enreg%comm_fft,ierr) 1496 !end if 1497 1498 !Instead of communications the whole real space vector on all comm_fft CPUs, 1499 !it is possible to let each CPU keep it's part only and communicate only the sums. 1500 !Then, however, we need to clean the trash in the part of the real space 1501 !vector that's not attributed to the given comm_fft CPU. 1502 nproc_fft=mpi_enreg%nproc_fft 1503 if (nproc_fft>1) then 1504 me_fft=mpi_enreg%me_fft 1505 if (me_fft>0) then 1506 nd3=(ngfft(3)-1)/nproc_fft+1 1507 psi_out(:,:,:,me_fft*nd3+1:me_fft*nd3+nd3)=psi_out(:,:,:,1:nd3) 1508 psi_out(:,:,:,1:nd3)=zero 1509 end if 1510 !!! call xmpi_sum(psi_out,mpi_enreg%comm_fft,ierr) 1511 end if 1512 1513 end subroutine g_to_r
gwls_hamiltonian/gr_to_g [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
gr_to_g
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_DielectricArray,gwls_GWlanczos,gwls_LanczosBasis,gwls_hamiltonian gwls_model_polarisability,gwls_polarisability
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1540 subroutine gr_to_g(psig_out,psir_in,psig_in) 1541 1542 1543 !This section has been created automatically by the script Abilint (TD). 1544 !Do not modify the following lines by hand. 1545 #undef ABI_FUNC 1546 #define ABI_FUNC 'gr_to_g' 1547 use interfaces_53_ffts 1548 !End of the abilint section 1549 1550 implicit none 1551 real(dp), intent(in) :: psir_in(2,n4,n5,n6) 1552 real(dp), intent(in), optional :: psig_in(2,npw_g) 1553 real(dp), intent(out) :: psig_out(2,npw_g) 1554 1555 integer :: i1, i2, i3 1556 integer :: cplex, option 1557 1558 ! ************************************************************************* 1559 1560 cplex = 2 ! complex potential 1561 option= 2 ! multiply wavefunction by potential 1562 1563 if(.not. present(psig_in)) then 1564 psig4(:,:) = psig_out(:,:) 1565 else 1566 psig4(:,:) = psig_in(:,:) 1567 end if 1568 1569 do i3=1,n6 1570 do i2=1,n5 1571 do i1=1,n4 1572 denpot(2*i1-1,i2,i3)= psir_in(1,i1,i2,i3) 1573 denpot(2*i1 ,i2,i3)= psir_in(2,i1,i2,i3) 1574 end do 1575 end do 1576 end do 1577 1578 call fourwf( cplex, & ! complex potential 1579 denpot, & ! real space wavefunction, in denpot format 1580 psig4, & ! fourier space wavefunction 1581 psig_out, & ! result, in FFT configuration 1582 psir3,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,mpi_enreg,1, & ! Various other arguments 1583 ngfft,npw_g,npw_g,n4,n5,n6,option,dtset%paral_kgb,tim_fourwf,weight,weight) 1584 1585 end subroutine gr_to_g
gwls_hamiltonian/Hpsik [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
Hpsik
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_DielectricArray,gwls_LanczosResolvents,gwls_hamiltonian gwls_lineqsolver,gwls_polarisability
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1224 subroutine Hpsik(psi_out,psi_in,cte) 1225 1226 1227 !This section has been created automatically by the script Abilint (TD). 1228 !Do not modify the following lines by hand. 1229 #undef ABI_FUNC 1230 #define ABI_FUNC 'Hpsik' 1231 use interfaces_66_wfs 1232 !End of the abilint section 1233 1234 implicit none 1235 1236 !External variables 1237 real(dp), intent(inout) :: psi_out(2,npw_g) 1238 real(dp), intent(inout), optional :: psi_in(2,npw_g) 1239 real(dp), intent(in), optional :: cte 1240 1241 ! ************************************************************************* 1242 1243 !psikb1 is allocated in init_hamiltonian and destroyed in destroy_hamiltonian, to avoid doing it at each application of 1244 !the hamiltonian... 1245 1246 !We are keeping track of how the module treat each argument of the hamiltonian (subroutine getghc) here. 1247 !Legend : T = transmitted from 79_seqpar_mpi/vtowfk.F90 through the call gwls_sternheimer 1248 ! I = input argument (allocated and filled) 1249 ! O = output argument (allocated) 1250 ! D = dummy argument (declared but not allocated. Read/Write attempts should trigger a segfault.) 1251 !call getghc(cpopt, T 1252 !psi, I 1253 !conjgrprj, D 1254 !Hpsik, O 1255 !gsc, D (output for PAW : <G|S|C>) 1256 !gs_hamk, T 1257 !gvnlc, D (<G|Vnonlocal|C>) 1258 !eshift, I (<G|H-eshift.S|C>) 1259 !mpi_enreg, T 1260 !ndat, Fixed to 1 (# of FFTs to do in //) 1261 !dtset%prtvol, T 1262 !sij_opt, Fixed to 0 (PAW dependant : 0-><G|H|C> ; 1-><G|H|C> & <G|S|C> ; -1-><G|H-cte.S|C>) 1263 !tim_getghc, Fixed to 0 (identity of the timer of the calling subroutine. 1:cgwf 5:lobpcg) 1264 !type_calc) Fixed to 0 (0:whole H N:some part only) 1265 1266 ! It is necessary to pass this optional argument to getghc if paral_kgb=1, or else the code exits cleanly with a "BUG" message... 1267 ! It will be important to properly understand what this means when we start using kpoints... 1268 1269 if(present(psi_in)) then 1270 1271 call getghc(cpopt,psi_in,conjgrprj,psi_out,dummy2,gs_hamk,psig4,eshift,& 1272 & mpi_enreg,ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc) 1273 1274 else 1275 psig3 = psi_out 1276 1277 call getghc(cpopt,psig3,conjgrprj,psi_out,dummy2,gs_hamk,psig4,eshift,& 1278 & mpi_enreg,ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc) 1279 1280 end if 1281 if(present(cte)) then 1282 if(present(psi_in)) then 1283 psi_out = psi_out - cte*psi_in 1284 else 1285 psi_out = psi_out - cte*psig3 1286 end if 1287 end if 1288 1289 end subroutine Hpsik
gwls_hamiltonian/Hpsikc [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
Hpsikc
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_lineqsolver
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1315 subroutine Hpsikc(psi_out,psi_in,cte) 1316 1317 1318 !This section has been created automatically by the script Abilint (TD). 1319 !Do not modify the following lines by hand. 1320 #undef ABI_FUNC 1321 #define ABI_FUNC 'Hpsikc' 1322 !End of the abilint section 1323 1324 implicit none 1325 1326 !External variables 1327 complex(dpc), intent(out) :: psi_out(npw_g) 1328 complex(dpc), intent(in) :: psi_in(npw_g) 1329 complex(dpc), intent(in), optional :: cte 1330 1331 ! ************************************************************************* 1332 1333 1334 psig1(1,:) = dble(psi_in) 1335 psig1(2,:) = dimag(psi_in) 1336 1337 call Hpsik(psig1) 1338 1339 psi_out = dcmplx(psig1(1,:),psig1(2,:)) 1340 1341 if(present(cte)) psi_out = psi_out - cte*psi_in 1342 end subroutine Hpsikc
gwls_hamiltonian/kbkb_to_kb [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
kbkb_to_kb
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_hamiltonian
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1611 subroutine kbkb_to_kb(psik_out,psik_in_1,psik_in_2) 1612 !---------------------------------------------------------------------------------------------------- 1613 ! This function computes the direct product of two wavefunctions in real space, 1614 ! psi_out(r) = psi_in_1^*(r)*psi_in_2(r) (without complex conjugating any of the 2 wavefunctions) 1615 ! and returns the result in fourier space, psi_out(k). 1616 ! 1617 ! The two input wavefunctions are in k space. 1618 ! 1619 ! 1620 !---------------------------------------------------------------------------------------------------- 1621 1622 !This section has been created automatically by the script Abilint (TD). 1623 !Do not modify the following lines by hand. 1624 #undef ABI_FUNC 1625 #define ABI_FUNC 'kbkb_to_kb' 1626 !End of the abilint section 1627 1628 implicit none 1629 real(dp), intent(out) :: psik_out(2,npw_kb) 1630 real(dp), intent(inout) :: psik_in_1(2,npw_kb), psik_in_2(2,npw_kb) 1631 1632 ! ************************************************************************* 1633 1634 ! change configuration of the data : LA -> FFT 1635 call wf_block_distribute(psik_in_1, psig1,1) 1636 call wf_block_distribute(psik_in_2, psig2,1) 1637 1638 ! Fourier transform the first input 1639 call g_to_r(psir1, psig1) 1640 1641 ! Complex conjugate in real space the first input 1642 psir1(2,:,:,:) = -psir1(2,:,:,:) 1643 1644 ! Fourrier transform the second input, multiply component-wise 1645 call gr_to_g(psig2,psir1) 1646 1647 ! change configuration of the output : FFT -> LA 1648 call wf_block_distribute(psik_out, psig2,2) 1649 1650 end subroutine kbkb_to_kb
gwls_hamiltonian/pc_k_valence_kernel [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
pc_k_valence_kernel
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_DielectricArray,gwls_Projected_AT,gwls_hamiltonian gwls_lineqsolver,gwls_model_polarisability,gwls_polarisability
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
418 subroutine pc_k_valence_kernel(psi_inout,n) 419 420 !================================================================================ 421 ! This routine projects out of the valence kernel. It assumes the input/output 422 ! array is distributed in the FFT configuration, and that the global 423 ! array containing the kernel (defined in this module) is already prepared 424 ! and ready to be used. 425 !================================================================================ 426 427 !This section has been created automatically by the script Abilint (TD). 428 !Do not modify the following lines by hand. 429 #undef ABI_FUNC 430 #define ABI_FUNC 'pc_k_valence_kernel' 431 !End of the abilint section 432 433 implicit none 434 435 real(dp), intent(inout) :: psi_inout(2,npw_g) 436 integer , intent(in), optional :: n 437 438 439 real(dp), allocatable :: psi_projected(:,:) 440 441 real(dp) :: tmpc(2) 442 443 integer :: v, n_max 444 445 integer :: mpi_communicator, ierr 446 447 ! ************************************************************************* 448 449 450 ABI_ALLOCATE( psi_projected, (2,npw_g)) 451 452 mpi_communicator = mpi_enreg%comm_fft 453 454 if (present(n)) then 455 n_max = n 456 else 457 n_max = nbandv 458 end if 459 460 !==================================================== 461 ! Compute projection 462 !==================================================== 463 464 psi_projected(:,:) = zero 465 466 do v = 1, n_max 467 468 ! compute overlap of kernel member with function 469 tmpc = cg_zdotc(npw_g ,kernel_wavefunctions_FFT(:,:,v),psi_inout(:,:)) 470 471 ! Communicate results 472 call xmpi_sum(tmpc, mpi_communicator, ierr) ! sum on all processors working on FFT! 473 474 ! add overlap with array to project 475 psi_projected(1,:) = psi_projected(1,:) + (tmpc(1)*kernel_wavefunctions_FFT(1,:,v)-tmpc(2)*kernel_wavefunctions_FFT(2,:,v)) 476 psi_projected(2,:) = psi_projected(2,:) + (tmpc(1)*kernel_wavefunctions_FFT(2,:,v)+tmpc(2)*kernel_wavefunctions_FFT(1,:,v)) 477 478 !psi_inout(1,:) = psi_inout(1,:) -( tmpc(1)*kernel_wavefunctions_FFT(1,:,v)-tmpc(2)*kernel_wavefunctions_FFT(2,:,v) ) 479 !psi_inout(2,:) = psi_inout(2,:) -( tmpc(1)*kernel_wavefunctions_FFT(2,:,v)+tmpc(2)*kernel_wavefunctions_FFT(1,:,v) ) 480 481 482 end do 483 484 if (present(n)) then 485 psi_inout(:,:) = psi_projected(:,:) 486 else 487 psi_inout(:,:) = psi_inout(:,:) - psi_projected(:,:) 488 end if 489 490 491 492 ABI_DEALLOCATE( psi_projected) 493 494 end subroutine pc_k_valence_kernel
gwls_hamiltonian/precondition [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
precondition
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_lineqsolver
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1081 subroutine precondition(psi_out,psi_in) 1082 1083 1084 !This section has been created automatically by the script Abilint (TD). 1085 !Do not modify the following lines by hand. 1086 #undef ABI_FUNC 1087 #define ABI_FUNC 'precondition' 1088 !End of the abilint section 1089 1090 implicit none 1091 1092 real(dp), intent(out) :: psi_out(2,npw_g) 1093 real(dp), intent(in) :: psi_in(2,npw_g) 1094 1095 ! ************************************************************************* 1096 1097 do i=1,npw_g 1098 psi_out(:,i) = psi_in(:,i)*pcon(i) 1099 end do 1100 1101 end subroutine precondition
gwls_hamiltonian/precondition_cplx [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
precondition_cplx
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_lineqsolver
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1127 subroutine precondition_cplx(psi_out,psi_in) 1128 1129 1130 !This section has been created automatically by the script Abilint (TD). 1131 !Do not modify the following lines by hand. 1132 #undef ABI_FUNC 1133 #define ABI_FUNC 'precondition_cplx' 1134 !End of the abilint section 1135 1136 implicit none 1137 1138 complex(dpc), intent(out) :: psi_out(npw_g) 1139 complex(dpc), intent(in) :: psi_in(npw_g) 1140 1141 ! ************************************************************************* 1142 1143 do i=1,npw_g 1144 psi_out(i) = psi_in(i)*pcon(i) 1145 end do 1146 end subroutine precondition_cplx
gwls_hamiltonian/set_precondition [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
set_precondition
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_LanczosResolvents,gwls_lineqsolver
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
913 subroutine set_precondition(lambda,omega) 914 !-------------------------------------------------------------------------------- 915 ! This subroutine preconditions the problem 916 ! 917 ! A x = b, 918 ! 919 ! with: 920 ! 921 ! omega lambda Operator 922 ! ------------------------------------------------------ 923 ! absent absent A = (H - lambda_0) (value of lambda_0 not important) 924 ! present present A = (H - lambda)^2 + omega^2 925 ! other cases not implemented 926 ! 927 ! 928 ! In the above, b = psik. 929 ! 930 !-------------------------------------------------------------------------------- 931 932 ! TODO : 933 ! - eliminate the 2 "if(kinpw(i) < huge(0.0_dp)*1.0d-11)" 934 ! since ecutsm = 0.0 always (check if that's true in this gw_sternheimer subroutine). 935 936 !This section has been created automatically by the script Abilint (TD). 937 !Do not modify the following lines by hand. 938 #undef ABI_FUNC 939 #define ABI_FUNC 'set_precondition' 940 !End of the abilint section 941 942 implicit none 943 944 real(dp), intent(in), optional :: lambda, omega 945 946 real(dp) :: poly, x 947 logical :: omega_imaginary 948 949 950 !integer,save :: counter = 0 951 !integer :: io_unit 952 !character(18):: filename 953 !logical :: file_exists 954 955 ! ************************************************************************* 956 957 958 if (present(lambda) .and. present(omega)) then 959 omega_imaginary =.true. 960 else 961 omega_imaginary =.false. 962 end if 963 964 !io_unit = get_unit() 965 !filename = "PRECONDITIONER.log" 966 !inquire(file=filename,exist=file_exists) 967 968 !if (file_exists) then 969 ! open(io_unit,file=filename,position='append',status=files_status_old) 970 !else 971 ! open(io_unit,file=filename,status=files_status_new) 972 ! write(io_unit,10) "#=======================================================================================" 973 ! write(io_unit,10) "# " 974 ! write(io_unit,10) "# This file contains information regarding the preconditioning scheme for SQMR. " 975 ! write(io_unit,10) "# " 976 ! write(io_unit,10) "#=======================================================================================" 977 !end if 978 979 !counter = counter + 1 980 !write(io_unit,10) "# " 981 !write(io_unit,15) "# Call #:", counter 982 !if (present(lambda) .and. present(omega)) then 983 ! write(io_unit,10) "# lambda and omega are present: imaginary frequency case" 984 !else 985 ! write(io_unit,10) "# lambda and omega are absent: omega = 0 case" 986 !end if 987 988 !do i=1,npw_k 989 990 do i = 1, npw_g 991 992 if(omega_imaginary) then 993 x = (kinpw_gather(i)-lambda)**2 + omega**2 994 else 995 x = kinpw_gather(i) 996 end if 997 998 if(x < huge(0.0_dp)*1.0d-11) then 999 poly = 27.0 + x*(18.0 + x*(12.0 + 8.0*x)) 1000 pcon(i) = poly/(poly + 16.0*(x**4)) 1001 !pcon(i) = 1.0/(1.0+x) !I don't know why, it gives better results for Silane than the above polynomial. 1002 else 1003 pcon(i) = zero 1004 end if 1005 end do 1006 1007 !write(io_unit,30) " prec( 1: 10) = ",pcon(1:10) 1008 !write(io_unit,30) " prec(npw_k-10:npw_k) = ",pcon(npw_k-10:npw_k) 1009 1010 !close(io_unit) 1011 1012 !10 format(A) 1013 !15 format(A,I8) 1014 !30 format(A,1000(ES24.12)) 1015 1016 end subroutine set_precondition
gwls_hamiltonian/sqrt_vc_k [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
sqrt_vc_k
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_DielectricArray,gwls_GWlanczos,gwls_LanczosBasis,gwls_hamiltonian gwls_model_polarisability,gwls_polarisability
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1173 subroutine sqrt_vc_k(psi_inout) 1174 1175 1176 !This section has been created automatically by the script Abilint (TD). 1177 !Do not modify the following lines by hand. 1178 #undef ABI_FUNC 1179 #define ABI_FUNC 'sqrt_vc_k' 1180 !End of the abilint section 1181 1182 implicit none 1183 1184 !External variables 1185 real(dp), intent(inout) :: psi_inout(2,npw_k) 1186 1187 !Internal variable 1188 complex(dpc) :: c 1189 1190 ! ************************************************************************* 1191 1192 do i=1,npw_k 1193 c = vc_sqrt(i) * cmplx(psi_inout(1,i),psi_inout(2,i),dpc) 1194 psi_inout(1,i) = dble (c) 1195 psi_inout(2,i) = dimag(c) 1196 end do 1197 end subroutine sqrt_vc_k
gwls_hamiltonian/unset_precondition [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
unset_precondition
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_lineqsolver
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
1042 subroutine unset_precondition() 1043 1044 1045 !This section has been created automatically by the script Abilint (TD). 1046 !Do not modify the following lines by hand. 1047 #undef ABI_FUNC 1048 #define ABI_FUNC 'unset_precondition' 1049 !End of the abilint section 1050 1051 implicit none 1052 ! ************************************************************************* 1053 1054 pcon = one 1055 end subroutine unset_precondition
gwls_hamiltonian/wf_block_distribute [ Functions ]
[ Top ] [ gwls_hamiltonian ] [ Functions ]
NAME
wf_block_distribute
FUNCTION
.
INPUTS
OUTPUT
PARENTS
gwls_DielectricArray,gwls_GWlanczos,gwls_LanczosBasis,gwls_Projected_AT gwls_Projected_BT,gwls_hamiltonian,gwls_model_polarisability gwls_polarisability
CHILDREN
copy_hamiltonian,copy_mpi_enreg,crystal_init,crystal_print distributevalencekernel,distributevalencewavefunctions,dtset_copy fftpac,find_qmesh,gsph_init,hpsik,kmesh_init,kmesh_print load_k_hamiltonian,make_gvec_kss,pc_k_valence_kernel,print_gsphere set_wf,vcoul_init,wf_block_distribute,xmpi_sum
SOURCE
522 subroutine wf_block_distribute(psik, psik_alltoall, direction) 523 524 525 !This section has been created automatically by the script Abilint (TD). 526 !Do not modify the following lines by hand. 527 #undef ABI_FUNC 528 #define ABI_FUNC 'wf_block_distribute' 529 !End of the abilint section 530 531 implicit none 532 533 !================================================================================ 534 ! 535 ! This subroutine distributes a block of wavefunctions in the "linear algebra" 536 ! configuration, to a configuration appropriate to perform FFT (or apply Hamiltonian). 537 ! 538 ! The code below is inspired by the subroutine prep_getghc, which rearranges 539 ! data over MPI processors prior to applying the Hamiltonian. 540 ! 541 ! input: 542 ! npw_kb : dimension of wavefunctions in the "linear algebra" configuration, 543 ! which means the G-vectors are distributed over all 544 ! processors, and every processor has information about all the bands. 545 ! 546 ! 547 ! npw_g : dimension of wavefunctions in the "FFT" configuration, 548 ! which means the G-vectors are distributed along rows of the 549 ! MPI topology, and a given row only has one band. 550 ! 551 ! direction : Are we going from LA to FFT, or from FFT to LA? 552 ! 1 : LA -> FFT 553 ! 2 : LA <- FFT 554 ! 555 ! input/output: 556 ! psik : block of wavefunctions, in LA configuration 557 ! psik_alltoall : wavefunction, in FFT configuration 558 !================================================================================ 559 560 integer , intent(in) :: direction ! flag which determines the direction of the transfer 561 real(dp), intent(inout) :: psik(2,npw_kb) ! block of wavefunctions, in "linear algebra" configuration 562 real(dp), intent(inout) :: psik_alltoall(2,npw_g) ! wavefunction in "FFT" configuration; a single band, but more G-vectors 563 564 565 integer :: ier, spaceComm 566 567 integer,allocatable :: rdisplsloc(:) 568 integer,allocatable :: recvcountsloc(:) 569 integer,allocatable :: sdisplsloc(:) 570 integer,allocatable :: sendcountsloc(:) 571 572 integer :: nproc_band, bandpp 573 574 integer,pointer :: rdispls(:) 575 integer,pointer :: recvcounts(:) 576 integer,pointer :: sdispls(:) 577 integer,pointer :: sendcounts(:) 578 579 ! ************************************************************************* 580 581 582 ! extract information in order to perform MPI communication. 583 ! This code comes from prep_getghc 584 nproc_band = mpi_enreg%nproc_band 585 bandpp = mpi_enreg%bandpp 586 587 if(mpi_enreg%nproc_band*mpi_enreg%bandpp > 1) then 588 589 spaceComm=mpi_enreg%comm_fft 590 if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_band 591 592 ABI_ALLOCATE(sendcountsloc,(nproc_band)) 593 ABI_ALLOCATE(sdisplsloc ,(nproc_band)) 594 ABI_ALLOCATE(recvcountsloc,(nproc_band)) 595 ABI_ALLOCATE(rdisplsloc ,(nproc_band)) 596 597 recvcounts =>bandfft_kpt(ikpt_this_proc)%recvcounts(:) 598 sendcounts =>bandfft_kpt(ikpt_this_proc)%sendcounts(:) 599 rdispls =>bandfft_kpt(ikpt_this_proc)%rdispls (:) 600 sdispls =>bandfft_kpt(ikpt_this_proc)%sdispls (:) 601 602 recvcountsloc(:)= recvcounts(:)*2*nspinor*bandpp 603 rdisplsloc(:) = rdispls(:)*2*nspinor*bandpp 604 sendcountsloc(:)= sendcounts(:)*2*nspinor 605 sdisplsloc(:) = sdispls(:)*2*nspinor 606 607 ! use MPI to communicate information! 608 if (direction == 1) then 609 ! LA -> FFT 610 call xmpi_alltoallv(psik, sendcountsloc, sdisplsloc, psik_alltoall, recvcountsloc,rdisplsloc, spaceComm, ier) 611 612 else if (direction == 2) then 613 ! FFT -> LA 614 615 call xmpi_alltoallv(psik_alltoall,recvcountsloc,rdisplsloc, psik, sendcountsloc,sdisplsloc,spaceComm,ier) 616 617 end if 618 619 ABI_DEALLOCATE(sendcountsloc) 620 ABI_DEALLOCATE(sdisplsloc ) 621 ABI_DEALLOCATE(recvcountsloc) 622 ABI_DEALLOCATE(rdisplsloc ) 623 624 else 625 626 if(direction == 1) psik_alltoall = psik 627 628 if(direction == 2) psik = psik_alltoall 629 630 end if 631 632 end subroutine wf_block_distribute