TABLE OF CONTENTS


ABINIT/m_gwls_hamiltonian [ Modules ]

[ Top ] [ Modules ]

NAME

 m_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 module m_gwls_hamiltonian
29 
30 ! local modules
31 use m_gwls_utility
32 use m_gwls_wf
33 
34 ! abinit modules
35 use m_bandfft_kpt
36 use m_cgtools
37 use defs_basis
38 use defs_datatypes
39 use defs_abitypes
40 use m_abicore
41 use m_xmpi
42 use m_pawang
43 use m_errors
44 use m_ab7_mixing
45 use m_mpinfo
46 
47 use m_io_tools,         only : get_unit
48 use m_dtset,            only : dtset_copy, dtset_free
49 use m_hamiltonian,      only : gs_hamiltonian_type, copy_hamiltonian, destroy_hamiltonian, &
50 &                              load_k_hamiltonian
51 use m_paw_dmft,         only : paw_dmft_type
52 use m_pawcprj,          only : pawcprj_type
53 use m_vcoul,            only : vcoul_t, vcoul_init, vcoul_free
54 use m_crystal,          only : crystal_t, crystal_init, crystal_free, crystal_print
55 use m_io_kss,           only : make_gvec_kss
56 use m_gsphere,          only : gsphere_t, gsph_init, gsph_free, print_gsphere
57 use m_bz_mesh,          only : kmesh_t, kmesh_init, kmesh_free, kmesh_print, find_qmesh
58 use m_fft,              only : fftpac, fourwf
59 use m_getghc,           only : getghc
60 
61 implicit none
62 save
63 private

m_gwls_hamiltonian/build_H [ Functions ]

[ Top ] [ m_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

1875 subroutine build_H(dtset2,mpi_enreg2,cpopt2,cg2,gs_hamk2,kg_k2,kinpw2)
1876 
1877 !use m_bandfft_kpt
1878 use m_cgtools
1879 use m_wfutils
1880 
1881 !This section has been created automatically by the script Abilint (TD).
1882 !Do not modify the following lines by hand.
1883 #undef ABI_FUNC
1884 #define ABI_FUNC 'build_H'
1885 !End of the abilint section
1886 
1887 implicit none
1888 
1889 !Arguments of gw_sternheimer, reveived as argument by build_H-------------------------
1890 type(dataset_type),  intent(in) :: dtset2
1891 type(MPI_type),   intent(inout) :: mpi_enreg2
1892 type(gs_hamiltonian_type), intent(inout) :: gs_hamk2
1893 
1894 integer, intent(in) :: cpopt2
1895 integer, intent(in) :: kg_k2(3,gs_hamk2%npw_k)
1896 !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.
1897 real(dp), intent(in) :: cg2(2,dtset2%mpw*dtset2%nspinor*dtset2%mband*dtset2%mkmem*dtset2%nsppol)
1898 real(dp), intent(in) :: kinpw2(gs_hamk2%npw_k)
1899 !Local variables : most of them are in the module for now.
1900 real(dp) :: commutation_error
1901 integer  :: io_unit
1902 integer  :: io_unit_debug
1903 integer  :: ierr
1904 integer  :: cplx
1905 integer  :: j, k
1906 integer  :: dimph3d
1907 integer  :: mb
1908 
1909 integer  :: mpi_communicator
1910 
1911 character(128):: filename_debug
1912 
1913 
1914 real(dp), allocatable :: wfk_tmp1(:,:) ,wfk_tmp2(:,:)
1915 
1916 ! *************************************************************************
1917 
1918 ! Hartree-Fock cannot be used with GWLS.
1919 if(dtset2%usefock==1 .and. associated(gs_hamk2%fockcommon)) then
1920   MSG_ERROR(' build_H : Hartree-Fock option can not be used with optdriver==66 (GWLS calculations).')
1921 end if
1922 
1923 !First we copy the data structure types
1924 call dtset_copy(dtset,dtset2)  !dtset = dtset2
1925 
1926 call copy_mpi_enreg(mpi_enreg2,mpi_enreg)
1927 
1928 call copy_hamiltonian(gs_hamk,gs_hamk2)
1929 
1930 !Then we copy the standard types
1931 cpopt   = cpopt2
1932 npw_k = gs_hamk2%npw_k
1933 
1934 ABI_ALLOCATE(cg,(2,dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol))
1935 cg = cg2
1936 
1937 ABI_ALLOCATE(vlocal,(gs_hamk2%n4,gs_hamk2%n5,gs_hamk2%n6,gs_hamk2%nvloc))
1938 vlocal = gs_hamk2%vlocal
1939 gs_hamk%vlocal => vlocal
1940 
1941 ABI_ALLOCATE(kg_k,(3,npw_k))
1942 kg_k = kg_k2
1943 ABI_ALLOCATE(kinpw,(npw_k))
1944 kinpw = kinpw2
1945 
1946 dimffnl=0; if (blocksize<=1) dimffnl=size(gs_hamk2%ffnl_k,2)
1947 ABI_ALLOCATE(ffnl,(npw_k,dimffnl,gs_hamk%lmnmax,gs_hamk%ntypat))
1948 dimph3d=0; if (blocksize<=1) dimph3d=gs_hamk2%matblk
1949 ABI_ALLOCATE(ph3d,(2,npw_k,dimph3d))
1950 
1951 !Initializing variables from dataset
1952 nfft  =  dtset%nfft
1953 nline = dtset%nline
1954 tolwfr = dtset%tolwfr
1955 ngfft  = dtset%ngfft
1956 mcg = dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
1957 nspinor=dtset%nspinor
1958 n1=ngfft(1)
1959 n2=ngfft(2)
1960 n3=ngfft(3)
1961 n4=ngfft(4)
1962 n5=ngfft(5)
1963 n6=ngfft(6)
1964 mgfft=maxval(ngfft(1:3))
1965 nbandv = int(dtset%nelect)/2
1966 ABI_ALLOCATE(istwfk,(dtset%nkpt))
1967 istwfk(:)=dtset%istwfk
1968 
1969 !Initializing variables from gs_hamk
1970 ucvol = gs_hamk%ucvol
1971 ABI_ALLOCATE(gbound,(2*mgfft+8,2))
1972 !gbound = gs_hamk%gbound_k !Must be done later for bandft paralelism
1973 
1974 !Parameters which need to be set by hand for now...
1975 weight           = 1           ! The weight of the k-pts, which sum to 1.
1976 tim_fourwf       = 0
1977 ckpt             = 1           ! The kpt of the wf to FFT. We only consider gamma point for now.
1978 
1979 ! ndat = 1 used to be hard coded, which works fine for FFT only parallelism.
1980 ndat              = 1                                                ! # of FFT to do in parallel in fourwf
1981 
1982 eshift           = 0.0         ! For PAW, opt_sij=-1 gives <G|H-lambda.S|C> in ghc instead of <G|H|C>
1983 sij_opt          = 0           ! Option in PAW to tell getghc what to compute...
1984 tim_getghc       = 0           ! timing code of the calling subroutine(can be set to 0 if not attributed)
1985 type_calc        = 0           ! option governing which part of Hamitonian is to be applied: 0: whole Hamiltonian
1986 nband            = dtset%mband
1987 ispden           = 1           !When required, the spin index to be used. We don't support spin polarized calculations yet...
1988 
1989 
1990 
1991 !========================================================================================================================
1992 ! Trying to implement band parallelism (Bruno, 14/10/2013)
1993 !------------------------------------------------------------------------------------------------------------------------
1994 !
1995 ! Band parallelism is described in the article on the LOBPCG method by F. Bottin et al,
1996 ! "Large scale ab initio calculations based on three levels of parallelisation", Computational Materials Science 42 (2008) 329-336
1997 !
1998 ! The article describes the ABINIT implementation specifically, suggesting that the information there is directly related
1999 ! to the code.
2000 !
2001 ! The main points in order to understand band + FFT parallelism are:
2002 !
2003 !     1) The processors are arranged in a 2D Cartesian MPI topology, or dimensions M x P, with
2004 !                              M  = nproc_bands
2005 !                              P  = nproc_fft
2006 !        There are of course M x P processors
2007 !
2008 !     2) the wavefunction coefficients (namely the array cg) are distributed on these processors. Unfortunately,
2009 !        two different distribution schemes are necessary: one for linear algebra (matrix products, etc), and one
2010 !        consistent with the Goedecker FFTs. The code must thus change the data distribution back and forth.
2011 !
2012 !     3) the "linear algebra" distribution describes the array cg directly. The G vectors are distributed among  all
2013 !        M x P processors, such that npw = npw_tot / (nproc_bands x nproc_fft). Every processor has information about
2014 !        all bands. Schematically:
2015 !
2016 !
2017 !                                        <- P = nproc_fft ->
2018 !                             ----------------------------------
2019 !                             |  n    |  n     |  n    |  n    |       ( for all band indices n)
2020 !                     ^       |  G_11 | ...    | ...   |  G_P1 |
2021 !                     |       ----------------------------------
2022 !             M = nproc_bands |  n    |  n     |  n    |  n    |
2023 !                     |       |  G_1. | ...    | ...   |  G_P. |
2024 !                     v       ----------------------------------
2025 !                             |  n    |  n     |  n    |  n    |
2026 !                             |  G_1M | ...    | ...   |  G_PM |
2027 !                             ----------------------------------
2028 !
2029 !
2030 !        the LOBPCG algorithm acts on blocks of bands. Thus, although each processor has information about all bands,
2031 !        we must conceptually view the bands grouped in blocks of size M, as {c_{1}(G),...,c_{M}(G)},
2032 !        {c_{M+1}(G)...,c_{2M}(G)}, ...
2033 !
2034 !        With this conceptual grouping, for a given block each processor has M x NG/(M x P) coefficients, where NG is the
2035 !        total number of G-vectors.
2036 !
2037 !     4) The distribution of information above is not appropriate for parallel FFTs. The data for one block is thus
2038 !        transposed and takes the form
2039 !
2040 !                                        <- P = nproc_fft ->
2041 !                             ----------------------------------
2042 !                             |  1    |  1     | ...   |  1    |
2043 !                     ^       |  G_1  |  G_2   | ...   |  G_P  |
2044 !                     |       ----------------------------------
2045 !             M = nproc_bands |  ..   | ...    | ...   |  ...  |
2046 !                     |       |  G_1  |  G_2   | ...   |  G_P  |
2047 !                     v       ----------------------------------
2048 !                             |  M    |  M     | ...   |  M    |
2049 !                             |  G_1  |  G_2   | ...   |  G_P  |
2050 !                             ----------------------------------
2051 !
2052 !     where the set of G vectors G_i = (G_{i1}, G_{i2}, ..., G_{iM}). Thus, a given processor has information about a single
2053 !     band, but more G vectors. Each processor has
2054 !                             NG/(M x P) x M = NG/P coefficients, just like before.
2055 !     Also, it is clear that the information is only communicated in the columns of the diagram above, avoiding
2056 !     "all processors to all processors" communication.
2057 !
2058 !     The data is now distributed properly to do parallel FFTs! Each row in the diagram above corresponds to FFTs done
2059 !     in parallel over nproc_fft processors, on a given band. There are M rows running thus in parallel!
2060 !
2061 !     The underlying ABINIT routines are equiped to handle FFT parallelism, not band distributed parallelism.
2062 !     prep_getghc.F90 and lobpgcwf.F90 show how to rearange information to be able to apply basic ABINIT routines.
2063 !
2064 !     The code below is inspired / guessed from lobpcgwf and prep_getghc. WE ASSUME THERE IS NO SPINORS!  CODE SHOULD
2065 !     BE HEAVILY REVIEWED for k-points and spinors, if and when they come to be useful.
2066 !
2067 !========================================================================================================================
2068 
2069 ! variables for band parallelism, following the language in lobpcgwf
2070 nspinor          = 1                                      ! how many spinors are present. Hardcoded to 1 for now.
2071 blocksize        = mpi_enreg%nproc_band*mpi_enreg%bandpp  ! how many bands treated in a block; # of FFT done in parallel in getghc
2072 nbdblock         = nband/blocksize                        ! how many blocks of bands are there
2073 ikpt_this_proc   = 1                                      ! Assuming only 1 kpoint, for molecules.
2074 ! This will have to be reviewed for crystals!
2075 
2076 npw_kb           = npw_k*blocksize
2077 npw_g            = gs_hamk2%npw_fft_k
2078 
2079 if (blocksize>1) then
2080   kg_k_gather  => bandfft_kpt(ikpt_this_proc)%kg_k_gather
2081   ffnl_gather  => bandfft_kpt(ikpt_this_proc)%ffnl_gather
2082   ph3d_gather  => bandfft_kpt(ikpt_this_proc)%ph3d_gather
2083   kinpw_gather => bandfft_kpt(ikpt_this_proc)%kinpw_gather
2084 else
2085   ffnl  = gs_hamk2%ffnl_k
2086   ph3d  = gs_hamk2%ph3d_k
2087   kg_k_gather  => kg_k
2088   ffnl_gather  => ffnl
2089   ph3d_gather  => ph3d
2090   kinpw_gather => kinpw
2091 endif
2092 call load_k_hamiltonian(gs_hamk,kinpw_k=kinpw_gather,kg_k=kg_k_gather,ffnl_k=ffnl_gather,ph3d_k=ph3d_gather)
2093 
2094 gbound = gs_hamk%gbound_k
2095 
2096 !Set up wfk routines
2097 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)
2098 
2099 ! prepare the valence wavefunctions and projection operator to work in band+FFT parallel
2100 call DistributeValenceWavefunctions()
2101 call DistributeValenceKernel()
2102 
2103 cplx             = 2                                      ! wavefunctions have complex coefficients
2104 
2105 ! Initialize the dummy variables
2106 ABI_DATATYPE_ALLOCATE(conjgrprj,(0,0))
2107 ABI_ALLOCATE(dummy2,(0,0))
2108 ABI_ALLOCATE(dummy3,(0,0,0))
2109 
2110 !Initialisation of the total counter for iterations of SQMR :
2111 ktot = 0
2112 
2113 !Allocations of working arrays for the module
2114 !- for pc_k function
2115 ABI_ALLOCATE(scprod2,(2,nband))
2116 
2117 !- for precondition function (and set_precondition subroutine)
2118 ABI_ALLOCATE(pcon,(npw_g))
2119 pcon = one
2120 
2121 !- for (private) working wf
2122 ABI_ALLOCATE(psik1,(2,npw_k))
2123 ABI_ALLOCATE(psik2,(2,npw_k))
2124 ABI_ALLOCATE(psik3,(2,npw_k))
2125 ABI_ALLOCATE(psik4,(2,npw_k))
2126 ABI_ALLOCATE(psikb1,(2,npw_kb))
2127 ABI_ALLOCATE(psikb2,(2,npw_kb))
2128 ABI_ALLOCATE(psikb3,(2,npw_kb))
2129 ABI_ALLOCATE(psikb4,(2,npw_kb))
2130 ABI_ALLOCATE(psig1,(2,npw_g))
2131 ABI_ALLOCATE(psig2,(2,npw_g))
2132 ABI_ALLOCATE(psig3,(2,npw_g))
2133 ABI_ALLOCATE(psig4,(2,npw_g))
2134 ABI_ALLOCATE(psir1,(2,n4,n5,n6))
2135 ABI_ALLOCATE(psir2,(2,n4,n5,n6))
2136 ABI_ALLOCATE(psir3,(2,n4,n5,n6))
2137 ABI_ALLOCATE(psidg,(2,nfft))
2138 ABI_ALLOCATE(denpot,(2*n4,n5,n6))
2139 
2140 psir1 = zero
2141 psir2 = zero
2142 psir3 = zero
2143 denpot = zero
2144 
2145 !Construct the vector of eigenvalues and write then to std output
2146 ABI_ALLOCATE(eig,(nband))
2147 eig=zero
2148 
2149 write(std_out,*) ch10,"Eigenvalues computation check, routine build_H:",ch10
2150 io_unit_debug = get_unit()
2151 write(filename_debug,'(A,I0.4,A)') "DEBUG_PROC=",mpi_enreg%me,".dat"
2152 
2153 
2154 
2155 mpi_communicator =  mpi_enreg%comm_bandfft
2156 
2157 open(file=filename_debug,status=files_status_old,unit=io_unit_debug)
2158 
2159 write(io_unit_debug,'(A)') " Parameters:"
2160 write(io_unit_debug,'(A,I5)') "                 nband = ", nband
2161 write(io_unit_debug,'(A,I5)') "             blocksize = ", blocksize
2162 write(io_unit_debug,'(A,I5)') "                 npw_k = ", npw_k
2163 write(io_unit_debug,'(A,I5)') "              nbdblock = ", nbdblock
2164 
2165 ! temporary wavefunction array, for data in the "linear algebra" distribution
2166 ABI_ALLOCATE( wfk_tmp1, (2,npw_k))
2167 ABI_ALLOCATE( wfk_tmp2, (2,npw_k))
2168 
2169 do n=1, nband
2170 ! Extract i^t/h wavefunction
2171 wfk_tmp1(:,1:npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k)
2172 
2173 ! Are the wavefunctions normalized?
2174 !tmpc = cg_zdotc(npw_k,wfk_tmp1,wfk_tmp1)
2175 !call xmpi_sum(tmpc,mpi_enreg%comm_bandfft,ierr) ! sum on all processors
2176 
2177 ! DEBUGGING CODE
2178 write(std_out,'(A,I5,A,2F24.16)') "band ", n, ", <psi_n | psi_n > =", norm_k(wfk_tmp1)
2179 write(io_unit_debug,'(A,I5,A,2F24.16)') "band ", n, ", <psi_n | psi_n > =", norm_k(wfk_tmp1)
2180 flush(io_unit_debug)
2181 end do
2182 
2183 
2184 ! loop on blocks of bands. All bands in a given block will be done in parallel!
2185 
2186 ! loop on block of bands
2187 do iblock = 1, nbdblock
2188 
2189 ! loop on states for this block
2190 do iband = 1, blocksize
2191 
2192 n = (iblock-1)*blocksize+iband
2193 
2194 psikb1(:,(iband-1)*npw_k+1:iband*npw_k)   = cg(:,(n-1)*npw_k+1:n*npw_k)
2195 
2196 end do
2197 
2198 ! change configuration of the data
2199 call wf_block_distribute(psikb1,  psig1,1) ! LA -> FFT
2200 
2201 ! Apply hamiltonian on wavefunction. In bandFFT parallel, Hpsik calls prep_getghc, which knows how to transform the
2202 ! distribution of the data from "linear algebra"-like to "FFT"-like. The output is *probably* in the "linear algebra"-like
2203 ! distribution.
2204 call Hpsik(psig2,psig1)
2205 
2206 ! change configuration of the data
2207 call wf_block_distribute(psikb2,  psig2,2) ! LA -> FFT
2208 
2209 ! extract the wavefunctions band by band
2210 do iband=1, blocksize
2211 
2212 n = blocksize*(iblock-1) + iband ! band index
2213 
2214 wfk_tmp1(:,1:npw_k) = psikb1(:,(iband-1)*npw_k+1:iband*npw_k)
2215 wfk_tmp2(:,1:npw_k) = psikb2(:,(iband-1)*npw_k+1:iband*npw_k)
2216 
2217 tmpc   = cg_zdotc(npw_k,wfk_tmp1,wfk_tmp2)
2218 
2219 call xmpi_sum(tmpc,mpi_enreg%comm_bandfft,ierr) ! sum on all processors
2220 eig(n) = tmpc(1)
2221 
2222 write(std_out,'(A,I5,A,F24.16,A)') "(build_H) band ", n, ", eig =", eig(n), " Ha."
2223 
2224 ! DEBUGGING CODE
2225 write(io_unit_debug,'(A,I5,A,F24.16,A)') "band ", n, ", eig =", eig(n), " Ha."
2226 flush(io_unit_debug)
2227 flush(io_unit_debug)
2228 end do
2229 
2230 end do
2231 
2232 ABI_DEALLOCATE( wfk_tmp1 )
2233 ABI_DEALLOCATE( wfk_tmp2 )
2234 close(io_unit_debug)
2235 
2236 
2237 
2238 ! Finishing the construction of vxc (transcribing it from the double real grid (for the density)
2239 ! to the single real grid (for the wfs).
2240 ! Assummes we only need one spin component; one transcription per spin being needed.
2241 if(allocated(vxc_dg)) then
2242   ABI_ALLOCATE(vxc,(n4,n5,n6,dtset%nspden))
2243   vxc = zero
2244   call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vxc_dg(:,ispden),vxc(:,:,:,ispden),2)
2245 end if
2246 
2247 timrev = 1 !Assumes istwfk *1. 1:time reversal symmetry NOT present | 2:" " " present
2248 ABI_ALLOCATE(title,(dtset%ntypat))
2249 do i=1,dtset%ntypat
2250 title(i) = "Bloup" ! The clean way would be to get the psps structure in this module
2251 ! (build_vxc is called from a place in GS calculations where it is available;
2252 ! should be the easiest way). For now, this allows the code to run.
2253 end do
2254 call crystal_init(dtset%amu_orig(:,1),Cryst,dtset%spgroup,dtset%natom,dtset%npsp,&
2255 &                 dtset%ntypat,dtset%nsym,dtset%rprimd_orig(:,:,1),dtset%typat,&
2256 &                 dtset%xred_orig(:,:,1),dtset%ziontypat,dtset%znucl,timrev,.false.,.false.,title,&
2257 &                 dtset%symrel,dtset%tnons,dtset%symafm)
2258 ABI_DEALLOCATE(title)
2259 call crystal_print(Cryst)
2260 
2261 !TODO : Should be put in a separate build_vc constructor, and should be called right after build_H in the context of optdriver 66.
2262 if(dtset%optdriver==66) then
2263 
2264   !Set up of the k-points and tables in the whole BZ
2265   call kmesh_init(Kmesh,Cryst,dtset%nkpt,dtset%kptns,dtset%kptopt,wrap_1zone=.false.)
2266   call kmesh_print(Kmesh,"K-mesh for the wavefunctions",std_out)
2267   call find_qmesh(Qmesh,Cryst,Kmesh)
2268   call kmesh_print(Qmesh,"Q-mesh for the screening function",std_out)
2269 
2270 
2271   !------------------------------
2272   !Building the vc_sqrt structure
2273   !------------------------------
2274   ! The sqrt_vc code relies on a different parallelism scheme than Vanilla ABINIT.
2275   ! The Gsphere object must be built on a single processor having all the G-vectors
2276   ! available.
2277 
2278   !gsph_init need gvec built according to the KSS convention, so that kg_k is not suitable (not properly sorted).
2279   npw_serial=npw_k
2280   call xmpi_sum(npw_serial,mpi_enreg%comm_bandfft,ierr)
2281 
2282   ecut_eff = dtset%ecut*(dtset%dilatmx)**2
2283   call make_gvec_kss(dtset%nkpt,dtset%kptns,ecut_eff,dtset%symmorphi,dtset%nsym,dtset%symrel,dtset%tnons,Cryst%gprimd,&
2284   &                      dtset%prtvol,npw_serial,gvec,ierr)
2285 
2286   call gsph_init(Gsphere,Cryst,npw_serial,gvec=gvec)
2287 
2288   call print_gsphere(Gsphere)
2289 
2290   call vcoul_init(Vcp,Gsphere,Cryst,Qmesh,Kmesh,dtset%rcut,dtset%icutcoul,dtset%vcutgeo,dtset%ecutsigx,npw_serial,&
2291   &               dtset%nkpt,dtset%kptns,dtset%ngfft,mpi_enreg%comm_world)
2292 
2293   ! Since Vcp%vc_sqrt is sorted according to the KSS convention for G vectors
2294   ! BUT will be applied to GS wavefunctions (where G vectors are sorted otherwise)
2295   ! It is necessary to construct a vc_sqrt vector with the GS sorting.
2296   ! Moreover, if we are in parallel (over band / FFTs), Vcp%vc_sqrt is the serial version
2297   ! and need to be distributed according to the GS scheme (Linear Algebra configuration).
2298   ABI_ALLOCATE(vc_sqrt,(npw_k))
2299   vc_sqrt=zero
2300   k=0
2301   do i=1,npw_k
2302   do j=1,npw_serial
2303   if(all(kg_k(:,i)==gvec(:,j))) k=j
2304   end do
2305   vc_sqrt(i)=Vcp%vc_sqrt(k,1)
2306   end do
2307 
2308 end if
2309 
2310 !--------------------------------------------------------------------------------
2311 ! Security check : The eigenstates of the projector and the hamiltonian need to
2312 ! agree down to the precision requested (tolwfr). Otherwise, SQMR is doomed.
2313 ! Now that the density is read and the eigenstates are calculated in the GW run,
2314 ! it is less useful. A test on tolwfr could be sufficient.
2315 !
2316 ! This remains a good tool to debug band+fft parallelism
2317 !--------------------------------------------------------------------------------
2318 
2319 ! only write on the head node!
2320 if (mpi_enreg%me == 0) then
2321   io_unit = get_unit()
2322   open(file='build_H.log',status=files_status_old,unit=io_unit)
2323   write(io_unit,10) '#----------------------------------------------------------------------------'
2324   write(io_unit,10) '#'
2325   write(io_unit,10) '#               This file presents the results of a small test to check      '
2326   write(io_unit,10) '#               how well the Hamiltonian commutes with the projection        '
2327   write(io_unit,10) '#               operator.                                                    '
2328   write(io_unit,10) '#'
2329   write(io_unit,10) '# Definitions:'
2330   write(io_unit,10) '#'
2331   write(io_unit,10) '#                    P     : projections on conduction states'
2332   write(io_unit,10) '#                    H     : Hamiltonian operator'
2333   write(io_unit,10) '#                  | psi > : eigenstate'
2334   write(io_unit,10) '#'
2335   flush(io_unit)
2336 end if
2337 
2338 psikb1 = zero
2339 
2340 ! sum all valence states, on copy in every band block.
2341 do i=1,nbandv
2342 
2343 do mb = 1, blocksize
2344 
2345 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)
2346 
2347 end do
2348 
2349 end do
2350 
2351 ! normalize to one! Only sum the fist band block
2352 tmpc = cg_zdotc(npw_k,psikb1,psikb1)
2353 call xmpi_sum(tmpc,mpi_communicator ,ierr) ! sum on all processors
2354 psikb1 = psikb1/sqrt(tmpc(1))
2355 
2356 
2357 ! change data distribution
2358 call wf_block_distribute(psikb1,  psig1,1) ! LA -> FFT
2359 
2360 ! Apply P.H operator
2361 call Hpsik(psig2 ,psig1)
2362 call pc_k_valence_kernel(psig2)
2363 
2364 ! Apply H.P operator
2365 call pc_k_valence_kernel(psig1)
2366 call Hpsik(psig1)
2367 
2368 ! compute error
2369 psig3 = psig1 - psig2 ! explicitly write the difference in an array
2370 tmpc   = cg_zdotc(npw_g,psig3,psig3)
2371 commutation_error = tmpc(1)
2372 
2373 call xmpi_sum(commutation_error , mpi_enreg%comm_fft, ierr) ! sum on all processors working on FFT!
2374 
2375 ! only write on the head node!
2376 if (mpi_enreg%me == 0) then
2377   write(io_unit,20) '   || (PH -HP) |b> ||^2 =  ',commutation_error
2378   write(io_unit,20) '           tolwfr       =  ',tolwfr
2379 end if
2380 
2381 if(commutation_error > tolwfr) then
2382   !write(io_unit,10) '# || (PH -HP) |b> ||^2 > tolwfr ==> Decision taken exit!'
2383 
2384   if (mpi_enreg%me == 0) write(io_unit,10) '# || (PH -HP) |b> ||^2 > tolwfr ==> This must be fixed!'
2385 
2386   write(std_out,20) "WARNING-build_H: The tolerance tolwfr=",tolwfr
2387   write(std_out,20) "                 is smaller than the commutation error ||(PH-HP)|b>||^2=",commutation_error
2388   write(std_out,10) "                 Either set tolwfr to a less stringent value in this calculation"
2389   write(std_out,10) "                 or to a more stringent value in the wavefunction calculation."
2390 end if
2391 
2392 if (mpi_enreg%me == 0) close(io_unit)
2393 
2394 10 format(A)
2395 20 format(A,ES12.3)
2396 end subroutine build_H

m_gwls_hamiltonian/build_vxc [ Functions ]

[ Top ] [ m_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

1665 subroutine build_vxc(vxc2,nfft2,nspden2)
1666 !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
1667 !arguments of fftpac are built in build_H.
1668 
1669 !This section has been created automatically by the script Abilint (TD).
1670 !Do not modify the following lines by hand.
1671 #undef ABI_FUNC
1672 #define ABI_FUNC 'build_vxc'
1673 !End of the abilint section
1674 
1675 implicit none
1676 
1677 !We need the dimensions of vxc since they don't exist yet in the module; build_vxc being called before build_H.
1678 integer, intent(in) :: nfft2, nspden2
1679 real(dp), intent(in) :: vxc2(nfft2,nspden2)
1680 
1681 ! *************************************************************************
1682 
1683 ABI_ALLOCATE(vxc_dg,(nfft2,nspden2))
1684 vxc_dg = vxc2
1685 
1686 end subroutine build_vxc

m_gwls_hamiltonian/destroy_H [ Functions ]

[ Top ] [ m_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

1712 subroutine destroy_H
1713 
1714 
1715 !This section has been created automatically by the script Abilint (TD).
1716 !Do not modify the following lines by hand.
1717 #undef ABI_FUNC
1718 #define ABI_FUNC 'destroy_H'
1719 !End of the abilint section
1720 
1721 implicit none
1722 
1723 call dtset_free(dtset)
1724 call destroy_hamiltonian(gs_hamk)
1725 
1726 call crystal_free(Cryst)
1727 call kmesh_free(Kmesh)
1728 call kmesh_free(Qmesh)
1729 call gsph_free(Gsphere)
1730 call vcoul_free(Vcp)
1731 
1732 call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg)
1733 call destroy_mpi_enreg(mpi_enreg)
1734 
1735 !NOTE : the syntax if(allocated(a)) ABI_DEALLOCATE(a) result in an error if "a" is not allocated; since the macro replace
1736 !ABI_ALLOCATE by more than one line of text, the second lines and up get outside the if... if() then syntax is equired.
1737 if(allocated(cg)) then
1738   ABI_DEALLOCATE(cg)
1739 end if
1740 if(allocated(gbound)) then
1741   ABI_DEALLOCATE(gbound)
1742 end if
1743 if(allocated(kg_k)) then
1744   ABI_DEALLOCATE(kg_k)
1745 end if
1746 if(allocated(ffnl)) then
1747   ABI_DEALLOCATE(ffnl)
1748 end if
1749 if(allocated(ph3d)) then
1750   ABI_DEALLOCATE(ph3d)
1751 end if
1752 if(allocated(kinpw)) then
1753   ABI_DEALLOCATE(kinpw)
1754 end if
1755 if(allocated(vxc)) then
1756   ABI_DEALLOCATE(vxc)
1757 end if
1758 if(allocated(vlocal)) then
1759   ABI_DEALLOCATE(vlocal)
1760 end if
1761 if(allocated(conjgrprj)) then
1762   ABI_DATATYPE_DEALLOCATE(conjgrprj)
1763 end if
1764 if(allocated(istwfk)) then
1765   ABI_DEALLOCATE(istwfk)
1766 end if
1767 if(allocated(dummy2)) then
1768   ABI_DEALLOCATE(dummy2)
1769 end if
1770 if(allocated(dummy3)) then
1771   ABI_DEALLOCATE(dummy3)
1772 end if
1773 if(allocated(eig)) then
1774   ABI_DEALLOCATE(eig)
1775 end if
1776 if(allocated(scprod2)) then
1777   ABI_DEALLOCATE(scprod2)
1778 end if
1779 if(allocated(pcon)) then
1780   ABI_DEALLOCATE(pcon)
1781 end if
1782 if(allocated(psik1)) then
1783   ABI_DEALLOCATE(psik1)
1784 end if
1785 if(allocated(psik2)) then
1786   ABI_DEALLOCATE(psik2)
1787 end if
1788 if(allocated(psik3)) then
1789   ABI_DEALLOCATE(psik3)
1790 end if
1791 if(allocated(psik4)) then
1792   ABI_DEALLOCATE(psik4)
1793 end if
1794 if(allocated(psikb1)) then
1795   ABI_DEALLOCATE(psikb1)
1796 end if
1797 if(allocated(psikb2)) then
1798   ABI_DEALLOCATE(psikb2)
1799 end if
1800 if(allocated(psikb3)) then
1801   ABI_DEALLOCATE(psikb3)
1802 end if
1803 if(allocated(psikb4)) then
1804   ABI_DEALLOCATE(psikb4)
1805 end if
1806 if(allocated(psig1)) then
1807   ABI_DEALLOCATE(psig1)
1808 end if
1809 if(allocated(psig2)) then
1810   ABI_DEALLOCATE(psig2)
1811 end if
1812 if(allocated(psig3)) then
1813   ABI_DEALLOCATE(psig3)
1814 end if
1815 if(allocated(psig4)) then
1816   ABI_DEALLOCATE(psig4)
1817 end if
1818 if(allocated(psir1)) then
1819   ABI_DEALLOCATE(psir1)
1820 end if
1821 if(allocated(psir2)) then
1822   ABI_DEALLOCATE(psir2)
1823 end if
1824 if(allocated(psir3)) then
1825   ABI_DEALLOCATE(psir3)
1826 end if
1827 if(allocated(psidg)) then
1828   ABI_DEALLOCATE(psidg)
1829 end if
1830 if(allocated(vxc_dg)) then
1831   ABI_DEALLOCATE(vxc_dg)
1832 end if
1833 if(allocated(denpot)) then
1834   ABI_DEALLOCATE(denpot)
1835 end if
1836 if(allocated(kernel_wavefunctions_FFT)) then
1837   ABI_DEALLOCATE(kernel_wavefunctions_FFT)
1838 end if
1839 if(allocated(valence_wavefunctions_FFT)) then
1840   ABI_DEALLOCATE(valence_wavefunctions_FFT)
1841 end if
1842 if(associated(gvec)) then
1843   ABI_DEALLOCATE(gvec)
1844 end if
1845 if(allocated(vc_sqrt)) then
1846   ABI_DEALLOCATE(vc_sqrt)
1847 end if
1848 
1849 end subroutine destroy_H

m_gwls_hamiltonian/dft_xc_energy [ Functions ]

[ Top ] [ m_gwls_hamiltonian ] [ Functions ]

NAME

  dft_xc_energy

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

798 function dft_xc_energy(e)
799 
800 
801 !This section has been created automatically by the script Abilint (TD).
802 !Do not modify the following lines by hand.
803 #undef ABI_FUNC
804 #define ABI_FUNC 'dft_xc_energy'
805 !End of the abilint section
806 
807 implicit none
808 real(dp) :: dft_xc_energy
809 integer, intent(in) :: e
810 
811 
812 integer :: cplex, option, ierr
813 real(dp), allocatable :: psik_e(:,:)             !Working array to store the wavefunction
814 real(dp), allocatable :: psik_e_alltoall(:,:)    !Working array to store the wavefunction
815 
816 real(dp), allocatable :: psik_out(:,:)  !Working array to store the wavefunction
817 
818 real(dp) :: tmpc(2)
819 integer  :: mpi_communicator
820 real(dp) :: dft_xc_energy_tmp
821 
822 ! *************************************************************************
823 
824 !--------------------------------------------------------------------------------
825 ! The goal of this routine is to compute the xc energy using
826 ! the DFT Vxc array.
827 !--------------------------------------------------------------------------------
828 
829 ! Parallel-MPI code. This is inspired from code within the getghc routine, which
830 ! applies a complex potential to a wavefunction.
831 
832 
833 cplex  = 1 ! real potential
834 option = 2 ! multiply wavefunction by potential
835 nspinor= 1
836 ! Allocate wavefunction arrays which contains the coefficients of the e-state
837 ABI_ALLOCATE(psik_e,         (2,npw_kb))
838 ABI_ALLOCATE(psik_e_alltoall,(2,npw_g))
839 ABI_ALLOCATE(psik_out,       (3,npw_g))
840 
841 
842 ! Only fetch the e-state, setting other states in block to zero!
843 psik_e(:,:)       = zero
844 psik_e(:,1:npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k)
845 psik_out(:,:)     = zero
846 
847 
848 ! change the configuration of the data
849 call wf_block_distribute(psik_e,  psik_e_alltoall,1) ! LA -> FFT
850 
851 
852 ! Call fourwf to generate the product, in k-space
853 ! Computation:
854 !                       psik_e_alltoall(k) -> psi(r)
855 !                       res(r)   =  (vxc(r) x psi(r))
856 !                       psik_out(k) <- res(r)
857 !  psir3 is a dummy, not used here.
858 
859 call fourwf(cplex,vxc(:,:,:,ispden),psik_e_alltoall,psik_out,psir3,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,&
860 &             mpi_enreg,1,ngfft,npw_g,npw_g,n4,n5,n6,option,dtset%paral_kgb,tim_fourwf,weight,weight)
861 
862 
863 tmpc = cg_zdotc(npw_g, psik_e_alltoall,psik_out)
864 
865 mpi_communicator =  mpi_enreg%comm_fft
866 call xmpi_sum(tmpc,mpi_communicator, ierr) ! sum on all processors working on FFT!
867 
868 dft_xc_energy_tmp = tmpc(1)
869 
870 mpi_communicator =  mpi_enreg%comm_band
871 call xmpi_sum(dft_xc_energy_tmp,mpi_communicator, ierr) ! sum on all processors working on FFT!
872 
873 dft_xc_energy = dft_xc_energy_tmp
874 
875 ABI_DEALLOCATE(psik_e)
876 ABI_DEALLOCATE(psik_e_alltoall)
877 ABI_DEALLOCATE(psik_out)
878 
879 end function dft_xc_energy

m_gwls_hamiltonian/DistributeValenceKernel [ Functions ]

[ Top ] [ m_gwls_hamiltonian ] [ Functions ]

NAME

  DistributeValenceKernel

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_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

314 subroutine DistributeValenceKernel()
315 !--------------------------------------------------------------------------------
316 !
317 ! This subroutine distributes, once and for all, the kernel of the static
318 ! SQMR operator.
319 !
320 ! In this first, quick and dirty implementation, we simply distribute
321 ! ALL the valence bands on ALL the FFT groups. This is not efficient, but
322 ! kernel projections are not a bottleneck of the computation.
323 !
324 ! A better (forthcoming) algorithm would only distribute the actual kernel,
325 ! not all valence bands.
326 !--------------------------------------------------------------------------------
327 
328 !This section has been created automatically by the script Abilint (TD).
329 !Do not modify the following lines by hand.
330 #undef ABI_FUNC
331 #define ABI_FUNC 'DistributeValenceKernel'
332 !End of the abilint section
333 
334 implicit none
335 
336 integer  :: mb, n
337 
338 
339 real(dp), allocatable :: psik_n(:,:)             !wavefunctions in LA format
340 real(dp), allocatable :: psik_n_alltoall(:,:)    !wavefunctions in FFT format
341 
342 ! *************************************************************************
343 
344 
345 !===================================================================
346 ! Allocate the global array which will contain the valence states
347 !===================================================================
348 ABI_ALLOCATE( kernel_wavefunctions_FFT, (2,npw_g, nband))
349 
350 kernel_wavefunctions_FFT = zero
351 
352 !====================================================
353 ! Allocate working arrays
354 !====================================================
355 
356 ABI_ALLOCATE(psik_n,         (2,npw_kb))
357 ABI_ALLOCATE(psik_n_alltoall,(2,npw_g))
358 
359 
360 ! loop on all valence states
361 do n = 1, nband
362 
363 ! Copy multiple instances of this valence state in the array;
364 ! this way, each FFT group will have ALL the valence states!
365 do mb = 1, blocksize
366 psik_n(:,(mb-1)*npw_k+1:mb*npw_k)   = cg(:,(n-1)*npw_k+1:n*npw_k)
367 end do
368 
369 ! change configuration of the data
370 call wf_block_distribute(psik_n,  psik_n_alltoall,1) ! LA -> FFT
371 
372 ! copy data in global array
373 kernel_wavefunctions_FFT(:,:,n) = psik_n_alltoall(:,:)
374 
375 end do
376 
377 !====================================================
378 ! cleanup
379 !====================================================
380 
381 ABI_DEALLOCATE(psik_n)
382 ABI_DEALLOCATE(psik_n_alltoall)
383 
384 end subroutine DistributeValenceKernel

m_gwls_hamiltonian/DistributeValenceWavefunctions [ Functions ]

[ Top ] [ m_gwls_hamiltonian ] [ Functions ]

NAME

  DistributeValenceWavefunctions

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_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

211 subroutine DistributeValenceWavefunctions()
212 !--------------------------------------------------------------------------------
213 !
214 ! This subroutine distributes, once and for all, the valence wavefunctions to
215 ! the FFT configuration. They will thus be ready to be used within the
216 ! susceptibility operator.
217 !
218 !--------------------------------------------------------------------------------
219 
220 !This section has been created automatically by the script Abilint (TD).
221 !Do not modify the following lines by hand.
222 #undef ABI_FUNC
223 #define ABI_FUNC 'DistributeValenceWavefunctions'
224 !End of the abilint section
225 
226 implicit none
227 
228 integer  :: iblk, mb, v
229 
230 
231 real(dp), allocatable :: psik_v(:,:)             !wavefunctions in LA format
232 real(dp), allocatable :: psik_v_alltoall(:,:)    !wavefunctions in FFT format
233 
234 ! *************************************************************************
235 
236 
237 
238 !===================================================================
239 ! Allocate the global array which will contain the valence states
240 !===================================================================
241 ABI_ALLOCATE( valence_wavefunctions_FFT, (2,npw_g,nbdblock))
242 
243 valence_wavefunctions_FFT = zero
244 
245 
246 !====================================================
247 ! Allocate working arrays
248 !====================================================
249 
250 ABI_ALLOCATE(psik_v,         (2,npw_kb))
251 ABI_ALLOCATE(psik_v_alltoall,(2,npw_g))
252 
253 
254 
255 ! loop on all blocks of states,
256 do iblk = 1, nbdblock
257 
258 ! loop on valence states for this block; if the state is conduction, fill with zeros
259 do mb = 1, blocksize
260 
261 v = (iblk-1)*blocksize+mb
262 
263 if (v <= nbandv) then
264   psik_v(:,(mb-1)*npw_k+1:mb*npw_k)   = cg(:,(v-1)*npw_k+1:v*npw_k)
265 else
266   psik_v(:,(mb-1)*npw_k+1:mb*npw_k)   = zero
267 end if
268 
269 end do
270 
271 ! change configuration of the data
272 call wf_block_distribute(psik_v,  psik_v_alltoall,1) ! LA -> FFT
273 
274 ! copy data in global array
275 
276 valence_wavefunctions_FFT(:,:,iblk) = psik_v_alltoall(:,:)
277 
278 end do
279 
280 !====================================================
281 ! cleanup
282 !====================================================
283 
284 ABI_DEALLOCATE(psik_v)
285 ABI_DEALLOCATE(psik_v_alltoall)
286 
287 
288 end subroutine DistributeValenceWavefunctions

m_gwls_hamiltonian/exchange [ Functions ]

[ Top ] [ m_gwls_hamiltonian ] [ Functions ]

NAME

  exchange

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

649 function exchange(e, Lbasis_lanczos)
650 
651 !use m_bandfft_kpt
652 use m_cgtools
653 !================================================================================
654 ! This subroutine computes the exchange energy in band+FFT parallel
655 !
656 !================================================================================
657 
658 !This section has been created automatically by the script Abilint (TD).
659 !Do not modify the following lines by hand.
660 #undef ABI_FUNC
661 #define ABI_FUNC 'exchange'
662 !End of the abilint section
663 
664 implicit none
665 real(dp) :: exchange
666 
667 integer, intent(in) :: e
668 
669 ! If these arguments are provided, the exchange energy is to be projected on this subspace
670 complex(dpc), optional, intent(in) :: Lbasis_lanczos(:,:)  ! complex array which contains the Lanczos basis
671 
672 real(dp), allocatable :: psik_e(:,:)             !Working array to store the wavefunction
673 
674 real(dp), allocatable :: psik_v(:,:)             !Working array to store the wavefunction
675 
676 real(dp), allocatable :: psik_out(:,:)           !Working array to store the wavefunction
677 
678 
679 integer :: iblk, mb
680 
681 integer :: l, lmax
682 
683 real(dp) :: tmpc(2)
684 
685 logical  :: project
686 
687 ! *************************************************************************
688 
689 !--------------------------------------------------------------------------------
690 ! Determine if the exhcange energy must be projected on the Lanczos basis
691 ! a truncated Coulomb potential
692 !--------------------------------------------------------------------------------
693 
694 project = .false.
695 if (present(Lbasis_lanczos)) then
696   project = .true.
697   lmax = size(Lbasis_lanczos, 2)
698 end if
699 
700 !--------------------------------------------------------------------------------
701 ! The goal of this routine is to compute the exchange energy using
702 ! a truncated Coulomb potential
703 !--------------------------------------------------------------------------------
704 
705 exchange = 0.0_dp
706 
707 !====================================================
708 ! build the block of wavefunctions which will
709 ! contain copies the e-state
710 !====================================================
711 
712 ABI_ALLOCATE(psik_e,(2,npw_kb))
713 
714 ! fill psik_e with as many copies of the e-state as there are band processors;
715 ! that way, upon LA -> FFT, each row of fft processors will have the e-state
716 do v = 1, blocksize
717 psik_e(:,(v-1)*npw_k+1:v*npw_k)   = cg(:,(e-1)*npw_k+1:e*npw_k)
718 end do
719 
720 !====================================================
721 ! Allocate the block of wavefunctions which will
722 ! contain the valence states
723 !====================================================
724 
725 ABI_ALLOCATE(psik_v,         (2,npw_kb))
726 ABI_ALLOCATE(psik_out,         (2,npw_kb))
727 
728 ! loop on all blocks of states,
729 do iblk = 1, nbdblock
730 
731 ! loop on valence states for this block; if the state is conduction, fill with zeros
732 do mb = 1, blocksize
733 
734 v = (iblk-1)*blocksize+mb
735 
736 if (v <= nbandv) then
737   psik_v(:,(mb-1)*npw_k+1:mb*npw_k)   = cg(:,(v-1)*npw_k+1:v*npw_k)
738 else
739   psik_v(:,(mb-1)*npw_k+1:mb*npw_k)   = zero
740 end if
741 
742 end do
743 
744 call kbkb_to_kb(psik_out,psik_v,psik_e)
745 
746 ! apply Coulomb potential, and take norm: cumulate the exchange energy
747 do mb = 1, blocksize
748 
749 psik1 = psik_out(:,(mb-1)*npw_k+1:mb*npw_k)
750 call sqrt_vc_k(psik1)
751 
752 if (project) then
753   ! project on the Lanczos basis
754   do l = 1, lmax
755   psik2(1,:) = dble(Lbasis_lanczos(:,l))
756   psik2(2,:) = dimag(Lbasis_lanczos(:,l))
757   tmpc = scprod_k(psik2,psik1)
758   exchange = exchange - (tmpc(1)**2+tmpc(2)**2)
759   end do
760 else
761   ! compute the "exact" exchange energy
762   exchange = exchange - norm_k(psik1)**2
763 end if
764 
765 end do
766 
767 end do
768 
769 ABI_DEALLOCATE(psik_e)
770 
771 ABI_DEALLOCATE(psik_v)
772 
773 ABI_DEALLOCATE(psik_out)
774 
775 end function exchange

m_gwls_hamiltonian/g_to_r [ Functions ]

[ Top ] [ m_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

1450 subroutine g_to_r(psi_out,psi_in)
1451 
1452 
1453 !This section has been created automatically by the script Abilint (TD).
1454 !Do not modify the following lines by hand.
1455 #undef ABI_FUNC
1456 #define ABI_FUNC 'g_to_r'
1457 !End of the abilint section
1458 
1459 implicit none
1460 real(dp), intent(out) :: psi_out(2,n4,n5,n6)
1461 real(dp), intent(in)  :: psi_in(2,npw_g)
1462 integer :: option, cplex
1463 integer :: nproc_fft, me_fft, nd3 !, ierr
1464 
1465 ! *************************************************************************
1466 
1467 option = 0 ! fft wavefunction to real space
1468 cplex  = 2 ! complex potential
1469 
1470 psig4 = psi_in
1471 psi_out = zero
1472 call fourwf(cplex,dummy3,psig4,dummy2,psi_out,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,mpi_enreg, &
1473 1,ngfft,npw_g,npw_g,n4,n5,n6,option,dtset%paral_kgb,tim_fourwf,weight,weight)
1474 psi_out = psi_out/sqrt(ucvol)
1475 
1476 !! This comes from prep_fourwf
1477 !nproc_fft=mpi_enreg%nproc_fft
1478 !if (nproc_fft>1) then
1479 !  me_fft=mpi_enreg%me_fft
1480 !  if (me_fft>0) then
1481 !    nd3=(ngfft(3)-1)/nproc_fft+1
1482 !    psi_out(:,:,:,me_fft*nd3+1:me_fft*nd3+nd3)=psi_out(:,:,:,1:nd3)
1483 !    psi_out(:,:,:,1:nd3)=zero
1484 !  end if
1485 !  call xmpi_sum(psi_out,mpi_enreg%comm_fft,ierr)
1486 !end if
1487 
1488 !Instead of communications the whole real space vector on all comm_fft CPUs,
1489 !it is possible to let each CPU keep it's part only and communicate only the sums.
1490 !Then, however, we need to clean the trash in the part of the real space
1491 !vector that's not attributed to the given comm_fft CPU.
1492 nproc_fft=mpi_enreg%nproc_fft
1493 if (nproc_fft>1) then
1494   me_fft=mpi_enreg%me_fft
1495   if (me_fft>0) then
1496     nd3=(ngfft(3)-1)/nproc_fft+1
1497     psi_out(:,:,:,me_fft*nd3+1:me_fft*nd3+nd3)=psi_out(:,:,:,1:nd3)
1498     psi_out(:,:,:,1:nd3)=zero
1499   end if
1500   !!!  call xmpi_sum(psi_out,mpi_enreg%comm_fft,ierr)
1501 end if
1502 
1503 end subroutine g_to_r

m_gwls_hamiltonian/gr_to_g [ Functions ]

[ Top ] [ m_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

1530 subroutine gr_to_g(psig_out,psir_in,psig_in)
1531 
1532 
1533 !This section has been created automatically by the script Abilint (TD).
1534 !Do not modify the following lines by hand.
1535 #undef ABI_FUNC
1536 #define ABI_FUNC 'gr_to_g'
1537 !End of the abilint section
1538 
1539 implicit none
1540 real(dp), intent(in)  :: psir_in(2,n4,n5,n6)
1541 real(dp), intent(in), optional :: psig_in(2,npw_g)
1542 real(dp), intent(out) :: psig_out(2,npw_g)
1543 
1544 integer :: i1, i2, i3
1545 integer :: cplex, option
1546 
1547 ! *************************************************************************
1548 
1549 cplex = 2 ! complex potential
1550 option= 2 ! multiply wavefunction by potential
1551 
1552 if(.not. present(psig_in)) then
1553   psig4(:,:) = psig_out(:,:)
1554 else
1555   psig4(:,:) = psig_in(:,:)
1556 end if
1557 
1558 do i3=1,n6
1559 do i2=1,n5
1560 do i1=1,n4
1561 denpot(2*i1-1,i2,i3)= psir_in(1,i1,i2,i3)
1562 denpot(2*i1  ,i2,i3)= psir_in(2,i1,i2,i3)
1563 end do
1564 end do
1565 end do
1566 
1567 call fourwf(          cplex, & ! complex potential
1568 denpot, & ! real space wavefunction, in denpot format
1569 psig4, & ! fourier space wavefunction
1570 psig_out, & ! result, in FFT configuration
1571 psir3,gbound,gbound,istwfk(ckpt),kg_k_gather,kg_k_gather,mgfft,mpi_enreg,1, & ! Various other arguments
1572 ngfft,npw_g,npw_g,n4,n5,n6,option,dtset%paral_kgb,tim_fourwf,weight,weight)
1573 
1574 end subroutine gr_to_g

m_gwls_hamiltonian/Hpsik [ Functions ]

[ Top ] [ m_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

1216 subroutine Hpsik(psi_out,psi_in,cte)
1217 
1218 
1219 !This section has been created automatically by the script Abilint (TD).
1220 !Do not modify the following lines by hand.
1221 #undef ABI_FUNC
1222 #define ABI_FUNC 'Hpsik'
1223 !End of the abilint section
1224 
1225 implicit none
1226 
1227 !External variables
1228 real(dp), intent(inout) :: psi_out(2,npw_g)
1229 real(dp), intent(inout), optional :: psi_in(2,npw_g)
1230 real(dp), intent(in), optional :: cte
1231 
1232 ! *************************************************************************
1233 
1234 !psikb1 is allocated in init_hamiltonian and destroyed in destroy_hamiltonian, to avoid doing it at each application of
1235 !the hamiltonian...
1236 
1237 !We are keeping track of how the module treat each argument of the hamiltonian (subroutine getghc) here.
1238 !Legend : T = transmitted from 79_seqpar_mpi/vtowfk.F90 through the call gwls_sternheimer
1239 !         I = input argument (allocated and filled)
1240 !         O = output argument (allocated)
1241 !         D = dummy argument (declared but not allocated. Read/Write attempts should trigger a segfault.)
1242 !call getghc(cpopt,             T
1243 !psi,              I
1244 !conjgrprj,        D
1245 !Hpsik,            O
1246 !gsc,              D              (output for PAW : <G|S|C>)
1247 !gs_hamk,          T
1248 !gvnlc,            D              (<G|Vnonlocal|C>)
1249 !eshift,           I              (<G|H-eshift.S|C>)
1250 !mpi_enreg,        T
1251 !ndat,             Fixed to 1     (# of FFTs to do in //)
1252 !dtset%prtvol,     T
1253 !sij_opt,          Fixed to 0     (PAW dependant : 0-><G|H|C> ; 1-><G|H|C> & <G|S|C> ; -1-><G|H-cte.S|C>)
1254 !tim_getghc,       Fixed to 0     (identity of the timer of the calling subroutine. 1:cgwf 5:lobpcg)
1255 !type_calc)        Fixed to 0     (0:whole H N:some part only)
1256 
1257 ! It is necessary to pass this optional argument to getghc if paral_kgb=1, or else the code exits cleanly with a "BUG" message...
1258 ! It will be important to properly understand what this means when we start using kpoints...
1259 
1260 if(present(psi_in)) then
1261 
1262   call getghc(cpopt,psi_in,conjgrprj,psi_out,dummy2,gs_hamk,psig4,eshift,&
1263 &             mpi_enreg,ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc)
1264 
1265 else
1266   psig3 = psi_out
1267 
1268   call getghc(cpopt,psig3,conjgrprj,psi_out,dummy2,gs_hamk,psig4,eshift,&
1269 &             mpi_enreg,ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc)
1270 
1271 end if
1272 if(present(cte)) then
1273   if(present(psi_in)) then
1274     psi_out = psi_out - cte*psi_in
1275   else
1276     psi_out = psi_out - cte*psig3
1277   end if
1278 end if
1279 
1280 end subroutine Hpsik

m_gwls_hamiltonian/Hpsikc [ Functions ]

[ Top ] [ m_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

1306 subroutine Hpsikc(psi_out,psi_in,cte)
1307 
1308 
1309 !This section has been created automatically by the script Abilint (TD).
1310 !Do not modify the following lines by hand.
1311 #undef ABI_FUNC
1312 #define ABI_FUNC 'Hpsikc'
1313 !End of the abilint section
1314 
1315 implicit none
1316 
1317 !External variables
1318 complex(dpc), intent(out) :: psi_out(npw_g)
1319 complex(dpc), intent(in)  :: psi_in(npw_g)
1320 complex(dpc), intent(in), optional :: cte
1321 
1322 ! *************************************************************************
1323 
1324 
1325 psig1(1,:) = dble(psi_in)
1326 psig1(2,:) = dimag(psi_in)
1327 
1328 call Hpsik(psig1)
1329 
1330 psi_out = dcmplx(psig1(1,:),psig1(2,:))
1331 
1332 if(present(cte)) psi_out = psi_out - cte*psi_in
1333 end subroutine Hpsikc

m_gwls_hamiltonian/kbkb_to_kb [ Functions ]

[ Top ] [ m_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

1600 subroutine kbkb_to_kb(psik_out,psik_in_1,psik_in_2)
1601 !----------------------------------------------------------------------------------------------------
1602 ! This function computes the direct product of two wavefunctions in real space,
1603 !         psi_out(r) = psi_in_1^*(r)*psi_in_2(r)  (without complex conjugating any of the 2 wavefunctions)
1604 ! and returns the result in fourier space, psi_out(k).
1605 !
1606 ! The two input wavefunctions are in k space.
1607 !
1608 !
1609 !----------------------------------------------------------------------------------------------------
1610 
1611 !This section has been created automatically by the script Abilint (TD).
1612 !Do not modify the following lines by hand.
1613 #undef ABI_FUNC
1614 #define ABI_FUNC 'kbkb_to_kb'
1615 !End of the abilint section
1616 
1617 implicit none
1618 real(dp), intent(out) :: psik_out(2,npw_kb)
1619 real(dp), intent(inout)  :: psik_in_1(2,npw_kb), psik_in_2(2,npw_kb)
1620 
1621 ! *************************************************************************
1622 
1623 ! change configuration of the data : LA -> FFT
1624 call wf_block_distribute(psik_in_1, psig1,1)
1625 call wf_block_distribute(psik_in_2, psig2,1)
1626 
1627 ! Fourier transform the first input
1628 call g_to_r(psir1, psig1)
1629 
1630 ! Complex conjugate in real space the first input
1631 psir1(2,:,:,:) = -psir1(2,:,:,:)
1632 
1633 ! Fourrier transform the second input, multiply component-wise
1634 call gr_to_g(psig2,psir1)
1635 
1636 ! change configuration of the output : FFT -> LA
1637 call wf_block_distribute(psik_out, psig2,2)
1638 
1639 end subroutine kbkb_to_kb

m_gwls_hamiltonian/pc_k_valence_kernel [ Functions ]

[ Top ] [ m_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

411 subroutine pc_k_valence_kernel(psi_inout,n)
412 
413 !================================================================================
414 ! This routine projects out of the valence kernel. It assumes the input/output
415 ! array is distributed in the FFT configuration, and that the global
416 ! array containing the kernel (defined in this module) is already prepared
417 ! and ready to be used.
418 !================================================================================
419 
420 !This section has been created automatically by the script Abilint (TD).
421 !Do not modify the following lines by hand.
422 #undef ABI_FUNC
423 #define ABI_FUNC 'pc_k_valence_kernel'
424 !End of the abilint section
425 
426 implicit none
427 
428 real(dp), intent(inout) :: psi_inout(2,npw_g)
429 integer , intent(in), optional :: n
430 
431 
432 real(dp), allocatable :: psi_projected(:,:)
433 
434 real(dp) :: tmpc(2)
435 
436 integer :: v, n_max
437 
438 integer :: mpi_communicator, ierr
439 
440 ! *************************************************************************
441 
442 
443 ABI_ALLOCATE( psi_projected, (2,npw_g))
444 
445 mpi_communicator =  mpi_enreg%comm_fft
446 
447 if (present(n)) then
448   n_max = n
449 else
450   n_max = nbandv
451 end if
452 
453 !====================================================
454 ! Compute projection
455 !====================================================
456 
457 psi_projected(:,:) = zero
458 
459 do v = 1, n_max
460 
461 ! compute overlap of kernel member with function
462 tmpc  = cg_zdotc(npw_g ,kernel_wavefunctions_FFT(:,:,v),psi_inout(:,:))
463 
464 ! Communicate results
465 call xmpi_sum(tmpc, mpi_communicator, ierr) ! sum on all processors working on FFT!
466 
467 ! add overlap with array to project
468 psi_projected(1,:) = psi_projected(1,:) + (tmpc(1)*kernel_wavefunctions_FFT(1,:,v)-tmpc(2)*kernel_wavefunctions_FFT(2,:,v))
469 psi_projected(2,:) = psi_projected(2,:) + (tmpc(1)*kernel_wavefunctions_FFT(2,:,v)+tmpc(2)*kernel_wavefunctions_FFT(1,:,v))
470 
471 !psi_inout(1,:) = psi_inout(1,:) -( tmpc(1)*kernel_wavefunctions_FFT(1,:,v)-tmpc(2)*kernel_wavefunctions_FFT(2,:,v) )
472 !psi_inout(2,:) = psi_inout(2,:) -( tmpc(1)*kernel_wavefunctions_FFT(2,:,v)+tmpc(2)*kernel_wavefunctions_FFT(1,:,v) )
473 
474 
475 end do
476 
477 if (present(n)) then
478   psi_inout(:,:) = psi_projected(:,:)
479 else
480   psi_inout(:,:) = psi_inout(:,:) - psi_projected(:,:)
481 end if
482 
483 
484 
485 ABI_DEALLOCATE( psi_projected)
486 
487 end subroutine pc_k_valence_kernel

m_gwls_hamiltonian/precondition [ Functions ]

[ Top ] [ m_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

1073 subroutine precondition(psi_out,psi_in)
1074 
1075 
1076 !This section has been created automatically by the script Abilint (TD).
1077 !Do not modify the following lines by hand.
1078 #undef ABI_FUNC
1079 #define ABI_FUNC 'precondition'
1080 !End of the abilint section
1081 
1082 implicit none
1083 
1084 real(dp), intent(out) :: psi_out(2,npw_g)
1085 real(dp), intent(in)  :: psi_in(2,npw_g)
1086 
1087 ! *************************************************************************
1088 
1089 do i=1,npw_g
1090 psi_out(:,i) = psi_in(:,i)*pcon(i)
1091 end do
1092 
1093 end subroutine precondition

m_gwls_hamiltonian/precondition_cplx [ Functions ]

[ Top ] [ m_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

1119 subroutine precondition_cplx(psi_out,psi_in)
1120 
1121 
1122 !This section has been created automatically by the script Abilint (TD).
1123 !Do not modify the following lines by hand.
1124 #undef ABI_FUNC
1125 #define ABI_FUNC 'precondition_cplx'
1126 !End of the abilint section
1127 
1128 implicit none
1129 
1130 complex(dpc), intent(out) :: psi_out(npw_g)
1131 complex(dpc), intent(in)  :: psi_in(npw_g)
1132 
1133 ! *************************************************************************
1134 
1135 do i=1,npw_g
1136 psi_out(i) = psi_in(i)*pcon(i)
1137 end do
1138 end subroutine precondition_cplx

m_gwls_hamiltonian/set_precondition [ Functions ]

[ Top ] [ m_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

 905 subroutine set_precondition(lambda,omega)
 906 !--------------------------------------------------------------------------------
 907 ! This subroutine preconditions the problem
 908 !
 909 !                                 A x = b,
 910 !
 911 ! with:
 912 !
 913 !        omega         lambda              Operator
 914 !        ------------------------------------------------------
 915 !     absent        absent         A =   (H - lambda_0)  (value of lambda_0 not important)
 916 !     present       present        A =   (H - lambda)^2 + omega^2
 917 !                         other cases not implemented
 918 !
 919 !
 920 ! In the above, b = psik.
 921 !
 922 !--------------------------------------------------------------------------------
 923 
 924 ! TODO :
 925 ! - eliminate the 2 "if(kinpw(i) < huge(0.0_dp)*1.0d-11)"
 926 !   since ecutsm = 0.0 always (check if that's true in this gw_sternheimer subroutine).
 927 
 928 !This section has been created automatically by the script Abilint (TD).
 929 !Do not modify the following lines by hand.
 930 #undef ABI_FUNC
 931 #define ABI_FUNC 'set_precondition'
 932 !End of the abilint section
 933 
 934 implicit none
 935 
 936 real(dp), intent(in), optional :: lambda, omega
 937 
 938 real(dp) :: poly, x
 939 logical :: omega_imaginary
 940 
 941 
 942 !integer,save :: counter = 0
 943 !integer      :: io_unit
 944 !character(18):: filename
 945 !logical      :: file_exists
 946 
 947 ! *************************************************************************
 948 
 949 
 950 if (present(lambda) .and. present(omega))  then
 951   omega_imaginary        =.true.
 952 else
 953   omega_imaginary        =.false.
 954 end if
 955 
 956 !io_unit  = get_unit()
 957 !filename = "PRECONDITIONER.log"
 958 !inquire(file=filename,exist=file_exists)
 959 
 960 !if (file_exists) then
 961 !      open(io_unit,file=filename,position='append',status=files_status_old)
 962 !else
 963 !      open(io_unit,file=filename,status=files_status_new)
 964 !        write(io_unit,10) "#======================================================================================="
 965 !        write(io_unit,10) "#                                                                                       "
 966 !        write(io_unit,10) "#   This file contains information regarding the preconditioning scheme for SQMR.       "
 967 !        write(io_unit,10) "#                                                                                       "
 968 !        write(io_unit,10) "#======================================================================================="
 969 !end if
 970 
 971 !counter = counter + 1
 972 !write(io_unit,10) "#                                                                                       "
 973 !write(io_unit,15) "#   Call #:", counter
 974 !if (present(lambda) .and. present(omega))  then
 975 !        write(io_unit,10) "#    lambda and omega are present: imaginary frequency case"
 976 !else
 977 !        write(io_unit,10) "#    lambda and omega are absent: omega = 0 case"
 978 !end if
 979 
 980 !do i=1,npw_k
 981 
 982 do i = 1, npw_g
 983 
 984 if(omega_imaginary) then
 985   x = (kinpw_gather(i)-lambda)**2 + omega**2
 986 else
 987   x = kinpw_gather(i)
 988 end if
 989 
 990 if(x < huge(0.0_dp)*1.0d-11) then
 991   poly    = 27.0 + x*(18.0 + x*(12.0 + 8.0*x))
 992   pcon(i) = poly/(poly + 16.0*(x**4))
 993   !pcon(i) = 1.0/(1.0+x) !I don't know why, it gives better results for Silane than the above polynomial.
 994 else
 995   pcon(i) = zero
 996 end if
 997 end do
 998 
 999 !write(io_unit,30) "                prec(       1:   10) =  ",pcon(1:10)
1000 !write(io_unit,30) "                prec(npw_k-10:npw_k) =  ",pcon(npw_k-10:npw_k)
1001 
1002 !close(io_unit)
1003 
1004 !10 format(A)
1005 !15 format(A,I8)
1006 !30 format(A,1000(ES24.12))
1007 
1008 end subroutine set_precondition

m_gwls_hamiltonian/sqrt_vc_k [ Functions ]

[ Top ] [ m_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

1165 subroutine sqrt_vc_k(psi_inout)
1166 
1167 
1168 !This section has been created automatically by the script Abilint (TD).
1169 !Do not modify the following lines by hand.
1170 #undef ABI_FUNC
1171 #define ABI_FUNC 'sqrt_vc_k'
1172 !End of the abilint section
1173 
1174 implicit none
1175 
1176 !External variables
1177 real(dp), intent(inout) :: psi_inout(2,npw_k)
1178 
1179 !Internal variable
1180 complex(dpc) :: c
1181 
1182 ! *************************************************************************
1183 
1184 do i=1,npw_k
1185 c = vc_sqrt(i) * cmplx(psi_inout(1,i),psi_inout(2,i),dpc)
1186 psi_inout(1,i) = dble (c)
1187 psi_inout(2,i) = dimag(c)
1188 end do
1189 end subroutine sqrt_vc_k

m_gwls_hamiltonian/unset_precondition [ Functions ]

[ Top ] [ m_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

1034 subroutine unset_precondition()
1035 
1036 
1037 !This section has been created automatically by the script Abilint (TD).
1038 !Do not modify the following lines by hand.
1039 #undef ABI_FUNC
1040 #define ABI_FUNC 'unset_precondition'
1041 !End of the abilint section
1042 
1043 implicit none
1044 ! *************************************************************************
1045 
1046 pcon = one
1047 end subroutine unset_precondition

m_gwls_hamiltonian/wf_block_distribute [ Functions ]

[ Top ] [ m_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

515 subroutine wf_block_distribute(psik, psik_alltoall, direction)
516 
517 
518 !This section has been created automatically by the script Abilint (TD).
519 !Do not modify the following lines by hand.
520 #undef ABI_FUNC
521 #define ABI_FUNC 'wf_block_distribute'
522 !End of the abilint section
523 
524 implicit none
525 
526 !================================================================================
527 !
528 ! This subroutine distributes a block of wavefunctions in the "linear algebra"
529 ! configuration, to a configuration appropriate to perform FFT (or apply Hamiltonian).
530 !
531 ! The code below is inspired by the subroutine prep_getghc, which rearranges
532 ! data over MPI processors prior to applying the Hamiltonian.
533 !
534 ! input:
535 !             npw_kb :   dimension of wavefunctions in the "linear algebra" configuration,
536 !                               which means the G-vectors are distributed over all
537 !                               processors, and every processor has information about all the bands.
538 !
539 !
540 !             npw_g       :   dimension of wavefunctions in the "FFT" configuration,
541 !                               which means the G-vectors are distributed along rows of the
542 !                               MPI topology, and a given row only has one band.
543 !
544 !              direction    :   Are we going from LA to FFT, or from FFT to LA?
545 !                              1 : LA  -> FFT
546 !                              2 : LA <-  FFT
547 !
548 ! input/output:
549 !              psik          : block of wavefunctions, in LA configuration
550 !              psik_alltoall : wavefunction, in FFT configuration
551 !================================================================================
552 
553 integer , intent(in)    :: direction                   ! flag which determines the direction of the transfer
554 real(dp), intent(inout) :: psik(2,npw_kb)       ! block of wavefunctions, in "linear algebra" configuration
555 real(dp), intent(inout) :: psik_alltoall(2,npw_g)    ! wavefunction in "FFT" configuration; a single band, but more G-vectors
556 
557 
558 integer  :: ier, spaceComm
559 
560 integer,allocatable :: rdisplsloc(:)
561 integer,allocatable :: recvcountsloc(:)
562 integer,allocatable :: sdisplsloc(:)
563 integer,allocatable :: sendcountsloc(:)
564 
565 integer :: nproc_band,  bandpp
566 
567 integer,pointer :: rdispls(:)
568 integer,pointer :: recvcounts(:)
569 integer,pointer :: sdispls(:)
570 integer,pointer :: sendcounts(:)
571 
572 ! *************************************************************************
573 
574 
575 ! extract information in order to perform MPI communication.
576 ! This code comes from prep_getghc
577 nproc_band = mpi_enreg%nproc_band
578 bandpp     = mpi_enreg%bandpp
579 
580 if(mpi_enreg%nproc_band*mpi_enreg%bandpp > 1) then
581 
582   spaceComm=mpi_enreg%comm_fft
583   if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_band
584 
585   ABI_ALLOCATE(sendcountsloc,(nproc_band))
586   ABI_ALLOCATE(sdisplsloc   ,(nproc_band))
587   ABI_ALLOCATE(recvcountsloc,(nproc_band))
588   ABI_ALLOCATE(rdisplsloc   ,(nproc_band))
589 
590   recvcounts   =>bandfft_kpt(ikpt_this_proc)%recvcounts(:)
591   sendcounts   =>bandfft_kpt(ikpt_this_proc)%sendcounts(:)
592   rdispls      =>bandfft_kpt(ikpt_this_proc)%rdispls   (:)
593   sdispls      =>bandfft_kpt(ikpt_this_proc)%sdispls   (:)
594 
595   recvcountsloc(:)= recvcounts(:)*2*nspinor*bandpp
596   rdisplsloc(:)   = rdispls(:)*2*nspinor*bandpp
597   sendcountsloc(:)= sendcounts(:)*2*nspinor
598   sdisplsloc(:)   = sdispls(:)*2*nspinor
599 
600   ! use MPI to communicate information!
601   if (direction == 1) then
602     ! LA -> FFT
603     call xmpi_alltoallv(psik, sendcountsloc, sdisplsloc, psik_alltoall, recvcountsloc,rdisplsloc, spaceComm, ier)
604 
605   else if (direction == 2) then
606     ! FFT -> LA
607 
608     call xmpi_alltoallv(psik_alltoall,recvcountsloc,rdisplsloc, psik, sendcountsloc,sdisplsloc,spaceComm,ier)
609 
610   end if
611 
612   ABI_DEALLOCATE(sendcountsloc)
613   ABI_DEALLOCATE(sdisplsloc   )
614   ABI_DEALLOCATE(recvcountsloc)
615   ABI_DEALLOCATE(rdisplsloc   )
616 
617 else
618 
619   if(direction == 1) psik_alltoall = psik
620 
621   if(direction == 2) psik = psik_alltoall
622 
623 end if
624 
625 end subroutine wf_block_distribute