TABLE OF CONTENTS
- ABINIT/bare_vqg
- ABINIT/m_fock
- ABINIT/strfock
- m_fock/fock_ACE_destroy
- m_fock/fock_calc_ene
- m_fock/fock_destroy
- m_fock/fock_get_getghc_call
- m_fock/fock_init
- m_fock/fock_print
- m_fock/fock_set_getghc_call
- m_fock/fock_set_ieigen
- m_fock/fock_type
- m_fock/fock_update_exc
- m_fock/fock_updatecwaveocc
- m_fock/fock_updateikpt
- m_fock/fockbz_create
ABINIT/bare_vqg [ Functions ]
NAME
bare_vqg
FUNCTION
Compute bare coulomb term in G-space on the FFT mesh i.e. 4pi/(G+q)**2
INPUTS
qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)) icutcoul=Option for the Coulomb potential cutoff technique divgq0= value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq. Used if q = Gamma gmet(3,3)=metrix tensor in G space in Bohr**-2. izero=if 1, unbalanced components of V(q,g) are set to zero # Used by the PAW library hyb_mixing=hybrid mixing coefficient for the Fock contribution hyb_mixing_sr=hybrid mixing coefficient for the short-range Fock contribution hyb_range_fock=hybrid range for separation nfft=Total number of FFT grid points. ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
OUTPUT
vqg(nfft)=4pi/(G+q)**2, G=0 component is set to divgq0/pi if q = Gamma.
NOTES
This routine operates on the full FFT mesh. DO NOT PASS MPI_TYPE One can easily implemente MPI-FFT by just calling this routine and then extracting the G-vectors treated by the node.
SOURCE
1985 subroutine bare_vqg(qphon,gsqcut,gmet,izero,hyb_mixing,hyb_mixing_sr,hyb_range_fock,nfft,nkpt_bz,ngfft,ucvol,vqg) 1986 1987 !Arguments ------------------------------------ 1988 !scalars 1989 integer,intent(in) :: izero,nfft,nkpt_bz 1990 real(dp),intent(in) :: gsqcut,hyb_mixing,hyb_mixing_sr,hyb_range_fock,ucvol 1991 !arrays 1992 integer,intent(in) :: ngfft(18) 1993 real(dp),intent(in) :: qphon(3) 1994 real(dp),intent(inout) :: gmet(3,3) 1995 real(dp),intent(out) :: vqg(nfft) 1996 1997 !Local variables------------------------------- 1998 !scalars 1999 integer,parameter :: cplex1=1 2000 integer :: i1,i2,i23,i3,id1,id2,id3 2001 integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max 2002 integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05 2003 real(dp),parameter :: tolfix=1.000000001e0_dp ! Same value as the one used in hartre 2004 real(dp) :: cutoff,gs,rcut,divgq0,gqg2p3,gqgm12,gqgm13,gqgm23,gs2,gs3 2005 !character(len=500) :: msg 2006 logical :: shortrange 2007 !arrays 2008 integer :: id(3) 2009 real(dp),allocatable :: gq(:,:) 2010 2011 ! ************************************************************************* 2012 2013 if (abs(hyb_mixing_sr)>tol8.and.abs(hyb_range_fock)<tol8) then 2014 ABI_BUG('SR mixing<>0 while range separation=0!') 2015 end if 2016 2017 !Treatment of the divergence at q+g=zero 2018 !For the time being, only Spencer-Alavi scheme... 2019 rcut= (three*nkpt_bz*ucvol/four_pi)**(one/three) 2020 divgq0= two_pi*rcut**two 2021 2022 !Initialize a few quantities 2023 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 2024 cutoff=gsqcut*tolfix 2025 vqg=zero 2026 2027 !Some peculiar values of q 2028 qeq0=0; if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1 2029 qeq05=0 2030 if (qeq0==0) then 2031 if (abs(abs(qphon(1))-half)<tol12.or.abs(abs(qphon(2))-half)<tol12.or. & 2032 & abs(abs(qphon(3))-half)<tol12) qeq05=1 2033 end if 2034 2035 !In order to speed the routine, precompute the components of g+q 2036 !Also check if the booked space was large enough... 2037 ABI_MALLOC(gq,(3,max(n1,n2,n3))) 2038 do ii=1,3 2039 id(ii)=ngfft(ii)/2+2 2040 do ing=1,ngfft(ii) 2041 ig=ing-(ing/id(ii))*ngfft(ii)-1 2042 gq(ii,ing)=ig+qphon(ii) 2043 end do 2044 end do 2045 ig1max=-1;ig2max=-1;ig3max=-1 2046 ig1min=n1;ig2min=n2;ig3min=n3 2047 2048 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 2049 2050 if (abs(hyb_mixing)>tol8) then 2051 shortrange=.false. 2052 rcut= (three*nkpt_bz*ucvol/four_pi)**(one/three) 2053 call barevcoul(rcut,qphon,gsqcut,gmet,nfft,nkpt_bz,ngfft,ucvol,vqg,shortrange) 2054 vqg=vqg*hyb_mixing 2055 if (hyb_range_fock>tol8)then 2056 vqg(1)=hyb_mixing*divgq0+hyb_mixing*(pi/hyb_range_fock**2) 2057 endif 2058 end if 2059 2060 if (abs(hyb_mixing_sr)>tol8) then 2061 shortrange=.true. 2062 rcut=hyb_range_fock 2063 call barevcoul(rcut,qphon,gsqcut,gmet,nfft,nkpt_bz,ngfft,ucvol,vqg,shortrange) 2064 vqg=vqg*hyb_mixing_sr 2065 ! if (hyb_range_fock>tol8)then 2066 ! vqg(1)=hyb_mixing*divgq0+hyb_mixing_sr*pi/hyb_range_fock**2 2067 ! endif 2068 end if 2069 2070 ! Triple loop on each dimension 2071 do i3=1,n3 2072 ig3=i3-(i3/id3)*n3-1 2073 ! Precompute some products that do not depend on i2 and i1 2074 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 2075 gqgm23=gq(3,i3)*gmet(2,3)*2 2076 gqgm13=gq(3,i3)*gmet(1,3)*2 2077 2078 do i2=1,n2 2079 ig2=i2-(i2/id2)*n2-1 2080 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 2081 gqgm12=gq(2,i2)*gmet(1,2)*2 2082 gqg2p3=gqgm13+gqgm12 2083 2084 i23=n1*(i2-1 +(n2)*(i3-1)) 2085 ! Do the test that eliminates the Gamma point outside of the inner loop 2086 ii1=1 2087 if (i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0) then 2088 ii1=2 2089 ! value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq 2090 vqg(1+i23)=hyb_mixing*divgq0 2091 2092 ! Note the combination of Spencer-Alavi and Erfc screening 2093 if (abs(hyb_range_fock)>tol8)then 2094 vqg(1+i23)=vqg(1+i23)+hyb_mixing_sr*(pi/hyb_range_fock**2) 2095 ! This would give a combination of Spencer-Alavi and Erfc screening, 2096 ! unfortunately, it modifies also the tests for pure HSE06, so was not retained. 2097 ! vqg(1+i23)=vqg(1+i23)+hyb_mixing_sr*min(divgq0,pi/(hyb_range_fock**2)) 2098 endif 2099 2100 end if 2101 2102 ! Final inner loop on the first dimension (note the lower limit) 2103 do i1=ii1,n1 2104 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 2105 2106 ii=i1+i23 2107 2108 if(gs<=cutoff)then 2109 ! Identify min/max indexes (to cancel unbalanced contributions later) 2110 ! Count (q+g)-vectors with similar norm 2111 if ((qeq05==1).and.(izero==1)) then 2112 ig1=i1-(i1/id1)*n1-1 2113 ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1) 2114 ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2) 2115 ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3) 2116 end if 2117 2118 end if ! Cut-off 2119 end do ! End loop on i1 2120 end do ! End loop on i2 2121 end do ! End loop on i3 2122 2123 2124 if (izero==1) then 2125 ! Set contribution of unbalanced components to zero 2126 if (qeq0==1) then !q=0 2127 call zerosym(vqg,cplex1,n1,n2,n3) 2128 else if (qeq05==1) then 2129 !q=1/2; this doesn't work in parallel 2130 ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2 2131 ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2 2132 ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2 2133 if (abs(abs(qphon(1))-half)<tol12) then 2134 if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max) 2135 if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min) 2136 end if 2137 if (abs(abs(qphon(2))-half)<tol12) then 2138 if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max) 2139 if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min) 2140 end if 2141 if (abs(abs(qphon(3))-half)<tol12) then 2142 if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max) 2143 if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min) 2144 end if 2145 call zerosym(vqg,cplex1,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3) 2146 end if 2147 end if 2148 2149 ABI_FREE(gq) 2150 2151 end subroutine bare_vqg
ABINIT/m_fock [ Modules ]
NAME
m_fock
FUNCTION
This module provides the definition of the fock_type used to store data for the calculation of Fock exact exchange term and the procedures to perform this calculation.
COPYRIGHT
Copyright (C) 2012-2024 ABINIT group (CMartins,FJ,FA,MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 module m_fock 25 26 use defs_basis 27 use m_abicore 28 use m_errors 29 use m_mpinfo 30 use m_xmpi 31 use libxc_functionals 32 use m_pawang 33 use m_pawtab 34 use m_pawfgr 35 use m_pawfgrtab 36 use m_pawcprj 37 use m_cgtools 38 use m_nctk 39 use m_dtset 40 41 use defs_abitypes, only : MPI_type 42 use m_time, only : timab 43 use m_fstrings, only : itoa, ftoa, sjoin 44 use m_symtk, only : mati3inv, matr3inv 45 use m_fftcore, only : sphereboundary 46 use m_fft, only : zerosym, fourwf 47 use m_kg, only : ph1d3d, getph 48 use m_kpts, only : listkk 49 use m_barevcoul, only : barevcoul 50 51 implicit none 52 53 private
ABINIT/strfock [ Functions ]
NAME
strfock
FUNCTION
Compute Fock energy contribution to stress tensor (Cartesian coordinates).
INPUTS
gsqcut=cutoff value on $G^2$ for (large) sphere inside fft box. $gsqcut=(boxcut^2)*ecut/(2._dp*(\pi^2))$ gprimd(3,3)=reciprocal space dimensional primitive translations hyb_mixing=hybrid mixing coefficient for the Fock contribution hyb_mixing_sr=hybrid mixing coefficient for the short-range Fock contribution hyb_range_fock=hybrid range for separation mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkpt_bz= number of k points in the BZ qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). rhog(2,nfft)=Fourier transform of charge density (bohr^-3) rhog2(2,nfft)= optional argument: Fourier transform of a second charge density (bohr^-3) ucvol=unit cell volume (bohr^3) vqg(nfft)=4pi/(G+q)**2
OUTPUT
fockstr(6)=components of Fock part of stress tensor (Cartesian coordinates, symmetric tensor) in hartree/bohr^3 Definition of symmetric tensor storage: store 6 unique components in the order 11, 22, 33, 32, 31, 21 (suggested by Xavier Gonze).
SOURCE
2186 subroutine strfock(gprimd,gsqcut,fockstr,hyb_mixing,hyb_mixing_sr,hyb_range_fock,mpi_enreg,nfft,ngfft,& 2187 nkpt_bz,rhog,ucvol,qphon,& 2188 rhog2) ! optional argument 2189 2190 !Arguments ------------------------------------ 2191 !scalars 2192 integer,intent(in) :: nfft,nkpt_bz 2193 real(dp),intent(in) :: gsqcut,hyb_mixing,hyb_mixing_sr,hyb_range_fock,ucvol 2194 type(MPI_type),intent(in) :: mpi_enreg 2195 !arrays 2196 integer,intent(in) :: ngfft(18) 2197 real(dp),intent(in) :: gprimd(3,3),rhog(2,nfft),qphon(3) 2198 real(dp),intent(in),optional :: rhog2(2,nfft) 2199 real(dp),intent(out) :: fockstr(6) 2200 2201 !Local variables------------------------------- 2202 !scalars 2203 integer,parameter :: im=2,re=1 2204 integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,me_fft,n1,n2,n3,nproc_fft 2205 real(dp) :: arg,cutoff,gsquar,rcut,rhogsq,tolfix=1.000000001_dp,tot,tot1,divgq0 2206 !character(len=500) :: msg 2207 !arrays 2208 real(dp) :: gcart(3),tsec(2) 2209 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 2210 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 2211 2212 ! ************************************************************************* 2213 2214 call timab(568,1,tsec) 2215 2216 if (abs(hyb_mixing_sr)>tol8.and.abs(hyb_range_fock)<tol8) then 2217 ABI_BUG('strfock: SR mixing<>0 while range separation=0!') 2218 end if 2219 2220 fockstr(:)=zero 2221 rcut= (three*nkpt_bz*ucvol/four_pi)**(one/three) 2222 irho2=0;if (present(rhog2)) irho2=1 2223 divgq0=two_pi/three*rcut**2 2224 2225 !Conduct looping over all fft grid points to find G vecs inside gsqcut 2226 !Include G**2 on surface of cutoff sphere as well as inside: 2227 cutoff=gsqcut*tolfix 2228 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 2229 me_fft=ngfft(11) 2230 nproc_fft=ngfft(10) 2231 id1=n1/2+2 2232 id2=n2/2+2 2233 id3=n3/2+2 2234 2235 2236 ii=0 2237 ! Get the distrib associated with this fft_grid 2238 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 2239 2240 do i3=1,n3 2241 ig3=i3-(i3/id3)*n3-1 2242 do i2=1,n2 2243 ig2=i2-(i2/id2)*n2-1 2244 if (fftn2_distrib(i2)==me_fft) then 2245 do i1=1,n1 2246 tot=zero; tot1=zero 2247 ig1=i1-(i1/id1)*n1-1 2248 ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1)) 2249 2250 ! Compute cartesian components of G 2251 gcart(1)=gprimd(1,1)*(dble(ig1)+qphon(1))+gprimd(1,2)*(dble(ig2)+qphon(2))+gprimd(1,3)*(dble(ig3)+qphon(3)) 2252 gcart(2)=gprimd(2,1)*(dble(ig1)+qphon(1))+gprimd(2,2)*(dble(ig2)+qphon(2))+gprimd(2,3)*(dble(ig3)+qphon(3)) 2253 gcart(3)=gprimd(3,1)*(dble(ig1)+qphon(1))+gprimd(3,2)*(dble(ig2)+qphon(2))+gprimd(3,3)*(dble(ig3)+qphon(3)) 2254 ! Compute |G+q|^2 2255 gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2 2256 ! take |rho(G)|^2 for complex rhog 2257 if (irho2==0) then 2258 rhogsq=rhog(re,ii)**2+rhog(im,ii)**2 2259 else 2260 rhogsq=rhog(re,ii)*rhog2(re,ii)+rhog(im,ii)*rhog2(im,ii) 2261 end if 2262 ! Case G=0: 2263 if(gsquar<tol10) then 2264 if (abs(hyb_mixing_sr)>tol8) cycle 2265 if (abs(hyb_mixing)>tol8) then 2266 fockstr(1)=fockstr(1)+hyb_mixing*divgq0*rhogsq 2267 fockstr(2)=fockstr(2)+hyb_mixing*divgq0*rhogsq 2268 fockstr(3)=fockstr(3)+hyb_mixing*divgq0*rhogsq 2269 cycle 2270 end if 2271 end if 2272 2273 ! Spencer-Alavi screening 2274 if (abs(hyb_mixing)>tol8) then 2275 arg=two_pi*rcut*sqrt(gsquar) 2276 tot=hyb_mixing*rhogsq*piinv/(gsquar**2)*(1-cos(arg)-arg*sin(arg)/two) 2277 tot1=hyb_mixing*rhogsq/three*rcut*sin(arg)/sqrt(gsquar) 2278 end if 2279 2280 ! Erfc screening 2281 if (abs(hyb_mixing_sr)>tol8) then 2282 arg=-gsquar*pi**2/(hyb_range_fock**2) 2283 tot=tot+hyb_mixing_sr*rhogsq*piinv/(gsquar**2)*(1.d0-exp(arg)*(1-arg)) 2284 end if 2285 fockstr(1)=fockstr(1)+tot*gcart(1)*gcart(1)+tot1 2286 fockstr(2)=fockstr(2)+tot*gcart(2)*gcart(2)+tot1 2287 fockstr(3)=fockstr(3)+tot*gcart(3)*gcart(3)+tot1 2288 fockstr(4)=fockstr(4)+tot*gcart(3)*gcart(2) 2289 fockstr(5)=fockstr(5)+tot*gcart(3)*gcart(1) 2290 fockstr(6)=fockstr(6)+tot*gcart(2)*gcart(1) 2291 end do 2292 end if 2293 end do 2294 end do 2295 2296 !Init mpi_comm 2297 if(mpi_enreg%nproc_fft>1)then 2298 call timab(48,1,tsec) 2299 call xmpi_sum(fockstr,mpi_enreg%comm_fft ,ierr) 2300 call timab(48,2,tsec) 2301 end if 2302 2303 2304 !Normalize and add term -efock/ucvol on diagonal 2305 !efock has been set to zero because it is not yet known. It will be added later. 2306 fockstr(1)=-fockstr(1) 2307 fockstr(2)=-fockstr(2) 2308 fockstr(3)=-fockstr(3) 2309 fockstr(4)=-fockstr(4) 2310 fockstr(5)=-fockstr(5) 2311 fockstr(6)=-fockstr(6) 2312 2313 call timab(568,2,tsec) 2314 2315 end subroutine strfock
m_fock/fock_ACE_destroy [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_ACE_destroy
FUNCTION
Clean and destroy fock datastructure.
INPUTS
fockACE <type(fock_ACE_type)>= all the quantities to calculate Fock exact exchange in the ACE context
SOURCE
1319 subroutine fock_ACE_destroy(fockACE) 1320 1321 !Arguments ------------------------------------ 1322 type(fock_ACE_type),pointer :: fockACE(:,:) 1323 !Local variables------------------------------- 1324 integer :: dim1,dim2,ii,jj 1325 ! ************************************************************************* 1326 DBG_ENTER("COLL") 1327 1328 dim1=size(fockACE,1) 1329 dim2=size(fockACE,2) 1330 do jj=1,dim2 1331 do ii=1,dim1 1332 if (allocated(fockACE(ii,jj)%xi)) then 1333 ABI_FREE(fockACE(ii,jj)%xi) 1334 end if 1335 end do 1336 end do 1337 DBG_EXIT("COLL") 1338 1339 end subroutine fock_ACE_destroy
m_fock/fock_calc_ene [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_calc_ene
FUNCTION
Calculate the Fock contribution to the total energy
INPUTS
fock <type(fock_type)>= all the quantities to calculate Fock exact exchange ikpt= reduced planewave coordinates.
OUTPUT
none
SIDE EFFECTS
energies <type(energies_type)>=storage for energies computed here : | e_exactX = Fock contribution to the total energy (Hartree)
NOTES
If the cgocc_bz are not updated at each iteration, be careful to calculate Fock energy at the same frequency. TO CHECK == CHANGE IN SOME DEFINTIONS
SOURCE
1367 subroutine fock_calc_ene(dtset,fock,fock_energy,ikpt,nband,occ) 1368 1369 !Arguments ------------------------------------ 1370 !scalars 1371 integer,intent(in) :: ikpt,nband 1372 real(dp),intent(inout) :: fock_energy 1373 type(dataset_type),intent(in) :: dtset 1374 type(fock_common_type),pointer :: fock 1375 !arrays 1376 real(dp),intent(in) :: occ(nband) 1377 1378 !Local variables------------------------------- 1379 integer :: iband 1380 1381 ! ************************************************************************* 1382 1383 ABI_UNUSED(fock_energy) 1384 1385 do iband=1,nband 1386 1387 ! Select only the occupied states (such that fock%occ_bz > 10^-8) 1388 if (abs(occ(iband))>tol8) then 1389 ! fock_energy=fock_energy + half*fock%eigen_ikpt(iband)*occ(iband)*dtset%wtk(ikpt) 1390 !* Sum the contribution of each occupied states at point k_i 1391 !* No need to multiply %wtk by ucvol since there is no factor 1/ucvol in the definition of %wtk 1392 1393 !* accumulate Fock contributions to the forces. 1394 ! if (fock%optfor) then 1395 fock%forces(:,:)=fock%forces(:,:)+occ(iband)*dtset%wtk(ikpt)*fock%forces_ikpt(:,:,iband) 1396 ! endif 1397 end if 1398 end do 1399 1400 end subroutine fock_calc_ene
m_fock/fock_destroy [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_destroy
FUNCTION
Clean and destroy fock datastructure.
INPUTS
fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
SOURCE
1194 subroutine fock_destroy(fock) 1195 1196 !Arguments ------------------------------------ 1197 type(fock_type),pointer :: fock 1198 1199 ! ************************************************************************* 1200 1201 if (fock%fock_common%use_ACE/=0) then 1202 ABI_FREE(fock%fockACE) 1203 end if 1204 ABI_FREE(fock%fock_common) 1205 ABI_FREE(fock%fock_BZ) 1206 ABI_FREE(fock) 1207 1208 end subroutine fock_destroy 1209 1210 subroutine fock_common_destroy(fock) 1211 1212 !Arguments ------------------------------------ 1213 type(fock_common_type),pointer :: fock 1214 1215 ! ************************************************************************* 1216 1217 DBG_ENTER("COLL") 1218 1219 ABI_SFREE(fock%atindx) 1220 ABI_SFREE(fock%typat) 1221 ! real arrays 1222 ABI_SFREE(fock%forces) 1223 ABI_SFREE(fock%nband) 1224 ABI_SFREE(fock%forces_ikpt) 1225 ABI_SFREE(fock%stress_ikpt) 1226 ABI_SFREE(fock%eigen_ikpt) 1227 ! Deallocate datatypes 1228 if (allocated(fock%pawfgrtab)) then 1229 call pawfgrtab_free(fock%pawfgrtab) 1230 ABI_FREE(fock%pawfgrtab) 1231 endif 1232 1233 ! Put the integer to 0 1234 fock%ieigen=0 1235 fock%ikpt=0 1236 fock%isppol=0 1237 1238 if (allocated(fock%symrec)) then 1239 ABI_FREE(fock%symrec) 1240 endif 1241 1242 !* [description of divergence in |q+G|=0] 1243 !* Put the real (dp) to 0 1244 fock%gsqcut=zero 1245 fock%hyb_mixing=zero 1246 fock%hyb_mixing_sr=zero 1247 fock%hyb_range_dft=zero 1248 fock%hyb_range_fock=zero 1249 1250 DBG_EXIT("COLL") 1251 end subroutine fock_common_destroy 1252 1253 1254 subroutine fock_BZ_destroy(fock) 1255 1256 !Arguments ------------------------------------ 1257 type(fock_BZ_type),pointer :: fock 1258 1259 ! ************************************************************************* 1260 1261 DBG_ENTER("COLL") 1262 1263 ABI_SFREE(fock%cwaveocc_bz) 1264 ABI_SFREE(fock%cgocc) 1265 ABI_SFREE(fock%npwarr) 1266 ABI_SFREE(fock%occ_bz) 1267 if (allocated(fock%cwaveocc_prj)) then 1268 call pawcprj_free(fock%cwaveocc_prj) 1269 ABI_FREE(fock%cwaveocc_prj) 1270 endif 1271 ! Deallocate integer arrays 1272 1273 ABI_SFREE(fock%kg_bz) 1274 ABI_SFREE(fock%nbandocc_bz) 1275 ABI_SFREE(fock%istwfk_bz) 1276 ABI_SFREE(fock%calc_phase) 1277 ABI_SFREE(fock%timerev) 1278 ABI_SFREE(fock%tab_ibg) 1279 ABI_SFREE(fock%tab_icg) 1280 ABI_SFREE(fock%tab_icp) 1281 ABI_SFREE(fock%tab_ikpt) 1282 ABI_SFREE(fock%tab_symkpt) 1283 1284 !* [description of IBZ and BZ] 1285 !* Deallocate real arrays 1286 ABI_SFREE(fock%wtk_bz) 1287 ABI_SFREE(fock%kptns_bz) 1288 ABI_SFREE(fock%phase) 1289 !* Put the integer to 0 1290 fock%nkpt_bz=0 1291 1292 !* Deallocate real arrays 1293 1294 !* Deallocate integer arrays 1295 ABI_SFREE(fock%gbound_bz) 1296 1297 !* [description of size of arrays/pointers] 1298 !* Put the integer to 0 1299 fock%mkpt=0 1300 fock%mkptband=0 1301 call destroy_mpi_enreg(fock%mpi_enreg) 1302 1303 DBG_EXIT("COLL") 1304 1305 end subroutine fock_BZ_destroy
m_fock/fock_get_getghc_call [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_get_getghc_call
FUNCTION
Returns the value of fock%getghc_call_
SOURCE
1871 pure integer function fock_get_getghc_call(fock) 1872 1873 !Arguments ------------------------------------ 1874 !scalars 1875 type(fock_common_type),intent(in) :: fock 1876 1877 ! ************************************************************************* 1878 1879 fock_get_getghc_call = fock%getghc_call_ 1880 1881 end function fock_get_getghc_call
m_fock/fock_init [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_init
FUNCTION
Init fock_t object
INPUTS
cg(2,mcg)= wavefunctions dtset <type(dataset_type)>= all input variables for this dataset gsqcut= Fourier cutoff on G^2 used to calculate charge density kg(3,mpw*mkmem)= reduced planewave coordinates. mcg= size of wave-functions array (cg) = mpw*nspinor*mband*mkmem*nsppol mpi_enreg=information about MPI parallelization npwarr_bz(nkpt)= number of planewaves in basis at this k point occ(mband*nkpt*nsppol)= occupation number for each band (often 2) at each k point
SIDE EFFECTS
fock <type(fock_type)>= all the quantities to calculate Fock exact exchange are initialized
NOTES
############################ ### Not fully tested yet ### ############################ The current version is restricted to the case nsym=1, nspinor=1 and mkmem/=0.
SOURCE
450 subroutine fock_init(atindx,cplex,dtset,fock,gsqcut,kg,mpi_enreg,nattyp,npwarr,pawang,pawfgr,pawtab,rprimd) 451 452 !Arguments ------------------------------------ 453 !scalars 454 integer, intent(in) :: cplex 455 real(dp),intent(in) :: gsqcut 456 type(dataset_type),intent(in) :: dtset 457 type(MPI_type),intent(in) :: mpi_enreg 458 type(fock_type),intent(inout),pointer :: fock 459 type(pawfgr_type),intent(in),target :: pawfgr 460 type(pawang_type),intent(in),target :: pawang 461 !arrays 462 integer, intent(in) :: atindx(dtset%natom),nattyp(dtset%ntypat), npwarr(dtset%nkpt) 463 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem) 464 real(dp), intent(in) :: rprimd(3,3) 465 type(pawtab_type), intent(in),target :: pawtab(dtset%ntypat*dtset%usepaw) 466 !Local variables------------------------------- 467 !scalars 468 integer :: iatom,ibg,icg,icp,ier,ik,ikg,ikpt,isppol,isym,itypat,jkpt,jpw,jsym,mband,mgfft,mkpt,mkptband 469 integer :: n1,n2,n3,n4,n5,n6,nband,ncpgr,nkpt_bz,nproc_hf,npwj,timrev,use_ACE,v1,v2,v3 470 integer :: my_jkpt,jkg_this_proc,my_nsppol,my_nspinor 471 real(dp) :: dksqmax,arg 472 character(len=500) :: msg 473 !arrays 474 integer :: indx(1),l_size_atm(dtset%natom),shiftg(3),symm(3,3),ident(3,3),symrec(3,3,dtset%nsym) 475 real(dp) :: gmet(3,3),gprimd(3,3),tau_nons(3),phktnons(2,1),tsec(2),Rtnons(3,dtset%nsym) 476 integer,allocatable :: dimcprj(:),indkk(:,:),kg_tmp(:),my_ikgtab(:),my_ibgtab(:,:),my_icgtab(:,:),my_icptab(:,:),invsym(:) 477 real(dp),allocatable :: kptns_hf(:,:), phase1d(:,:) 478 type(fock_common_type),pointer :: fockcommon 479 type(fock_BZ_type),pointer :: fockbz 480 481 ! ************************************************************************* 482 483 DBG_ENTER("COLL") 484 485 call timab(1501,1,tsec) 486 ABI_CHECK_IEQ(dtset%nspinor, 1, 'Hartree-Fock option can be used only with option nspinor = 1') 487 488 ! ===================================== 489 ! === Define useful local variables === 490 ! ===================================== 491 492 nkpt_bz=dtset%nkpthf 493 nproc_hf=mpi_enreg%nproc_hf 494 mband=dtset%nbandhf 495 496 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 497 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 498 499 !* Allocations 500 ABI_MALLOC(kptns_hf,(3,nkpt_bz)) 501 kptns_hf=zero 502 ABI_MALLOC(indkk,(nkpt_bz,6)) 503 indkk=0 504 ABI_MALLOC(phase1d,(2,(2*n1+1)*(2*n2+1)*(2*n3+1))) 505 phase1d=zero 506 ABI_MALLOC(kg_tmp,(3*dtset%mpw)) 507 508 !* Initialize the array my_ikgtab = shifts in arrays kg(ikg) associated to ikpt 509 ABI_MALLOC(my_ikgtab,(dtset%nkpt)) 510 ikg=0 511 do ikpt=1,dtset%nkpt 512 nband=dtset%nband(ikpt) 513 if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband,-1,mpi_enreg%me_kpt))) then 514 !* The point ikpt is treated on this processor. 515 my_ikgtab(ikpt)=ikg 516 !* The array kg is distributed, the shift ikg is incremented only on this proc. 517 ikg=ikg+npwarr(ikpt) 518 else 519 my_ikgtab(ikpt)=-1 520 !* Default value is -1. 521 end if 522 end do 523 524 !* Initialize the array my_ibgtab = shifts in arrays occ(ibg) associated to ikpt 525 !* Initialize the array my_icgtab = shifts in arrays cg(icg) associated to ikpt 526 ABI_MALLOC(my_ibgtab,(dtset%nkpt,dtset%nsppol)) 527 ABI_MALLOC(my_icgtab,(dtset%nkpt,dtset%nsppol)) 528 ABI_MALLOC(my_icptab,(dtset%nkpt,dtset%nsppol)) 529 ibg=0; icg=0 ;icp=0 530 do isppol=1,dtset%nsppol 531 do ikpt=1,dtset%nkpt 532 nband=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 533 if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband,isppol,mpi_enreg%me_kpt))) then 534 !* The states with (ikpt,isppol) are stored on this processor. 535 my_icgtab(ikpt,isppol)=icg 536 my_icptab(ikpt,isppol)=icp 537 !* The array cg is distributed, the shift icg is incremented only on this proc. 538 icg=icg+npwarr(ikpt)*nband 539 icp=icp+nband 540 else 541 my_icgtab(ikpt,isppol)=-1 542 my_icptab(ikpt,isppol)=-1 543 !* Otherwise, the states with (ikpt,isspol) are not stored on this processor and default value is -1. 544 end if 545 !* The array occ is shared among the proc, the shift ibg is always incremented. 546 my_ibgtab(ikpt,isppol)=ibg 547 ibg=ibg+nband 548 end do 549 end do 550 551 if (.not.(associated(fock))) then 552 553 ! ================================= 554 ! === Create the fock structure === 555 ! ================================= 556 ABI_MALLOC(fock,) 557 ABI_MALLOC(fock%fock_common,) 558 ABI_MALLOC(fock%fock_BZ,) 559 ! ======================================================== 560 ! === Set all the other state-dependent fields to zero === 561 ! ======================================================== 562 fockcommon=>fock%fock_common 563 fockbz=> fock%fock_BZ 564 ABI_MALLOC(fockcommon%nband,(dtset%nkpt*dtset%nsppol)) 565 do ikpt=1,dtset%nkpt*dtset%nsppol 566 fockcommon%nband(ikpt)=dtset%nband(ikpt) 567 end do 568 569 nband=dtset%mband 570 fockcommon%ikpt= 0 571 !* Will contain the k-point ikpt of the current state 572 fockcommon%isppol= 0 573 !* Will contain the spin isppol of the current state 574 fockcommon%ieigen=0 575 !* Will contain the band index of the current state 576 !* if the value is 0, the Fock contribution to the eigenvalue is not calculated. 577 ABI_MALLOC(fockcommon%eigen_ikpt,(nband)) 578 fockcommon%eigen_ikpt=0.d0 579 !* Will contain the Fock contributions to the eigenvalue of the current state 580 581 !* Compute the dimension of arrays in "spin" w.r.t parallelism 582 my_nsppol=dtset%nsppol 583 if (mpi_enreg%nproc_spkpt>1) my_nsppol=1 584 !* my_nsppol=1 when nsppol=1 or nsppol=2 and only one spin is treated by the processor. 585 !* my_nsppol=2 when nsppol=2 and no parallelization over kpt (both spins are treated by the processor). 586 587 !* Compute mkpt the size of arrays/pointers for k points w.r.t. parallelism 588 !* Compute mkptband the size of arrays/pointers for occupied states w.r.t. parallelism 589 if (nproc_hf<nkpt_bz) then 590 !* Parallelization over kpts only 591 mkpt=nkpt_bz/nproc_hf 592 if (mod(nkpt_bz,nproc_hf) /=0) mkpt=mkpt+1 593 mkptband=mkpt*mband 594 else 595 !* Parallelization over occupied states 596 if (nproc_hf<nkpt_bz*mband) then 597 mkptband=(nkpt_bz*mband)/nproc_hf 598 if (mod((nkpt_bz*mband),nproc_hf) /=0) mkptband=mkptband+1 599 mkpt=1 600 if (mod(nproc_hf,nkpt_bz) /=0) mkpt=2 601 else 602 mkptband=1 603 mkpt=1 604 end if 605 end if 606 607 ! mpi_enreg settings 608 call copy_mpi_enreg(mpi_enreg,fockbz%mpi_enreg) 609 fockbz%mpi_enreg%me_kpt=mpi_enreg%me_hf 610 fockbz%mpi_enreg%comm_kpt=mpi_enreg%comm_hf 611 fockbz%mpi_enreg%nproc_spkpt=mpi_enreg%nproc_hf 612 if (allocated(fockbz%mpi_enreg%proc_distrb)) then 613 ABI_FREE(fockbz%mpi_enreg%proc_distrb) 614 end if 615 ABI_MALLOC(fockbz%mpi_enreg%proc_distrb,(nkpt_bz,mband,1)) 616 do jkpt=1,nkpt_bz 617 fockbz%mpi_enreg%proc_distrb(jkpt,:,1)=fockbz%mpi_enreg%me_kpt 618 end do 619 620 mgfft=dtset%mgfft 621 fockcommon%usepaw=dtset%usepaw 622 if (fockcommon%usepaw==1)then 623 mgfft=dtset%mgfftdg 624 n4=dtset%ngfftdg(4) ; n5=dtset%ngfftdg(5) ; n6=dtset%ngfftdg(6) 625 end if 626 fockcommon%optfor=.false.; fockcommon%optstr=.false. 627 if(dtset%optforces==1) fockcommon%optfor=.true. 628 if (fockcommon%optfor) then 629 ABI_MALLOC(fockcommon%forces_ikpt,(3,dtset%natom,nband)) 630 ABI_MALLOC(fockcommon%forces,(3,dtset%natom)) 631 fockcommon%forces=zero 632 endif 633 use_ACE=1 ! Default. Normal users do not have access to this variable, although the next line allows experts to make tests. 634 if(dtset%userie==1729)use_ACE=0 ! Hidden possibility to disable ACE 635 636 fockcommon%use_ACE=use_ACE 637 call fockbz_create(fockbz,mgfft,dtset%mpw,mkpt,mkptband,my_nsppol,n4,n5,n6,use_ACE) 638 639 !* Initialize %mband, %mkpt, %mkptband = size of arrays 640 fockcommon%mband=mband 641 fockbz%mkpt=mkpt 642 fockbz%mkptband=mkptband 643 fockcommon%my_nsppol = my_nsppol 644 fockcommon%nsppol = dtset%nsppol 645 if (fockcommon%use_ACE/=0) then 646 ABI_MALLOC(fock%fockACE,(dtset%nkpt,dtset%nsppol)) 647 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 648 do isppol=1,dtset%nsppol 649 do ikpt=1,dtset%nkpt 650 nband=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 651 ABI_MALLOC(fock%fockACE(ikpt,isppol)%xi,(2,npwarr(ikpt)*my_nspinor,nband)) 652 end do 653 end do 654 end if 655 !========Initialze PAW data======== 656 fockcommon%ntypat=dtset%ntypat 657 fockcommon%natom=dtset%natom 658 if (fockcommon%usepaw==1) then 659 fockbz%mcprj=mkptband*my_nsppol 660 fockcommon%pawfgr => pawfgr 661 fockbz%pawang => pawang 662 fockcommon%pawtab => pawtab 663 ABI_MALLOC(fockcommon%pawfgrtab,(dtset%natom)) 664 do iatom = 1, dtset%natom 665 itypat=dtset%typat(iatom) 666 l_size_atm(iatom) = pawtab(itypat)%lcut_size 667 end do 668 call pawfgrtab_init(fockcommon%pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat) 669 call pawfgrtab_nullify(fockcommon%pawfgrtab) 670 ABI_MALLOC(fockbz%cwaveocc_prj,(dtset%natom,fockbz%mcprj)) 671 ABI_MALLOC(dimcprj,(dtset%natom)) 672 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 673 ncpgr = 0 674 if (dtset%optforces== 1) ncpgr = 3 675 ! if (dtset%optstress /= 0) ncpgr = 6 676 ! ncpgr=3*dtset%optforces+6*dtset%optstress 677 call pawcprj_alloc(fockbz%cwaveocc_prj,ncpgr,dimcprj) 678 ABI_FREE(dimcprj) 679 ABI_MALLOC(fockcommon%atindx,(dtset%natom)) 680 fockcommon%atindx=atindx 681 ABI_MALLOC(fockcommon%typat,(dtset%natom)) 682 fockcommon%typat=dtset%typat 683 end if 684 ! ========================================== 685 ! === Initialize the convergence options === 686 ! ========================================== 687 write(msg,'(2a)') ch10,'Fock_init: initialization of Fock operator parameters:' 688 call wrtout(std_out,msg) 689 690 fockcommon%fock_converged=.false. 691 fockcommon%scf_converged=.false. 692 693 !* Number of iterations with fixed occupied states when calculating the exact exchange contribution. 694 if (dtset%nnsclohf<0) then 695 ABI_ERROR('The parameter nnsclohf must be a non-negative integer.') 696 end if 697 if (dtset%nnsclohf==0) then 698 fockcommon%nnsclo_hf=1 699 msg=' - The parameter nnsclohf is set to its default value 1.' 700 call wrtout(std_out,msg) 701 !* Default value is set to 1 (updating cgocc at each step) 702 !* May be useful to put default to 3 703 else 704 fockcommon%nnsclo_hf=dtset%nnsclohf 705 write(msg,'(a,i3)') ' - The parameter nnsclohf is set to the value:', dtset%nnsclohf 706 call wrtout(std_out,msg) 707 !* value chosen by the user 708 end if 709 710 ! ========================================= 711 ! === Initialize the hybrid coefficient === 712 ! ========================================= 713 fockcommon%ixc = dtset%ixc 714 ! By convention, positive values are the default values for the ixc, 715 ! while negative values have been set by the user (and stored as negative numbers) 716 fockcommon%hyb_mixing=abs(dtset%hyb_mixing) 717 fockcommon%hyb_mixing_sr=abs(dtset%hyb_mixing_sr) 718 fockcommon%hyb_range_dft=abs(dtset%hyb_range_dft) 719 fockcommon%hyb_range_fock=abs(dtset%hyb_range_fock) 720 721 ! Set the hybrid parameters if functional from libxc for which parameters can be changed, or if the user asked to do so. 722 ! Usually, these parameters were obtained from libxc, 723 ! but the user might have possibly modified them. By the way, must define them here for the usual changeable fonctionals, 724 ! since otherwise might inherit them from the previous dataset ! 725 if(dtset%ixc<0)then 726 if (dtset%ixc==-406.or.dtset%ixc==-427.or.dtset%ixc==-428 .or. & 727 & min(dtset%hyb_mixing,dtset%hyb_mixing_sr,dtset%hyb_range_dft,dtset%hyb_range_fock)<-tol8)then 728 call libxc_functionals_set_hybridparams(hyb_mixing=fockcommon%hyb_mixing,& 729 & hyb_mixing_sr=fockcommon%hyb_mixing_sr,& 730 & hyb_range=fockcommon%hyb_range_dft) 731 end if 732 end if 733 734 735 ! ====================================================== 736 ! === Initialize the data relative to Poisson solver === 737 ! ====================================================== 738 739 !* gsqcut = cutoff value on G^2 for sphere inside the fft box (input for vhartre). 740 fockcommon%gsqcut= gsqcut 741 742 ! ======================================================= 743 ! === Initialize the properties of the k-points in BZ === 744 ! ======================================================= 745 !* Initialize %nkpt_bz = nb of k point in BZ for the calculation of exchange 746 fockbz%nkpt_bz=nkpt_bz 747 !* Initialize the array %wtk_bz = weight assigned to each k point. 748 fockbz%wtk_bz=1.0_dp/dble(nkpt_bz) 749 750 751 if (dtset%kptopt>=1 .and. dtset%kptopt<=4) then 752 ! ============================================ 753 ! === Initialize the set of k-points in BZ === 754 ! ============================================ 755 kptns_hf(:,1:nkpt_bz)=dtset%kptns_hf(:,1:nkpt_bz) 756 !* kptns_hf contains the special k points obtained by the Monkhorst & Pack method, in reduced coordinates. (output) 757 758 ! ======================================================= 759 ! === Compute the transformation to go from IBZ to BZ === 760 ! ======================================================= 761 !* Compute the reciprocal space metric. 762 call matr3inv(rprimd,gprimd) 763 gmet = MATMUL(TRANSPOSE(gprimd),gprimd) 764 765 !* Calculate the array indkk which describes how to get IBZ from BZ 766 !* dksqmax=maximal value of the norm**2 of the difference between a kpt2 vector and the closest k-point found from the kptns1 set, using symmetries. (output) 767 !* sppoldbl=1, no spin-polarisation doubling is required. 768 timrev=1 ; if (dtset%kptopt==3 .or. dtset%kptopt==4) timrev=0 769 !* timrev=1 if the use of time-reversal is allowed ; 0 otherwise 770 if (dtset%kptopt==2 .or. dtset%kptopt==3) then 771 !* No space symmetry is used, if kptopt==2 time reversal symmetry is used. 772 symm=0 ; symm(1,1)=1 ; symm(2,2)=1 ; symm(3,3)=1 773 call listkk(dksqmax,gmet,indkk(1:nkpt_bz,:),dtset%kptns,kptns_hf,dtset%nkpt, & 774 & nkpt_bz,1,1,indx,symm,timrev,xmpi_comm_self) 775 else 776 !* As in getkgrid, no use of antiferromagnetic symmetries thans to the option sppoldbl=1 777 call listkk(dksqmax,gmet,indkk(1:nkpt_bz,:),dtset%kptns,kptns_hf,dtset%nkpt, & 778 & nkpt_bz,dtset%nsym,1,dtset%symafm,dtset%symrel,timrev, xmpi_comm_self) 779 end if 780 !* indkk(nkpt_bz,6) describes the k point of IBZ that generates each k point of BZ 781 !* indkk(:,1) = k point of IBZ, kpt_ibz 782 !* indkk(:,2) = symmetry operation to apply to kpt_ibz to give the k point of BZ 783 !* (if 0, means no symmetry operation, equivalent to identity ) 784 !* indkk(:,3:5) = Umklapp vectors to apply to remain in BZ 785 !* indkk(:,6) = 1 if time-reversal was used to generate the k point of BZ, 0 otherwise 786 !* No use of symafm to generate spin down wfs from spin up wfs for the moment 787 788 else 789 if (dtset%kptopt==0) then 790 !* kptopt =0 : read directly nkpt, kpt, kptnrm and wtk in the input file 791 !* => this case is not allowed for the moment 792 ABI_ERROR('Hartree-Fock option can not be used with option kptopt=0.') 793 else 794 !* kptopt <0 : rely on kptbounds, and ndivk to set up a band structure calculation 795 !* => a band structure calculation is not yet allowed. 796 ABI_ERROR('Hartree-Fock option can not be used with option kptopt<0.') 797 end if 798 end if 799 800 !! ======================================================= 801 !! === Initialize the properties of the k-points in BZ === 802 !! ======================================================= 803 ! jkg=0 804 !!* Initialize the arrays %npwarr_bz, %kg_j, %phase_j, %gbound_j 805 ! do jkpt=1,nkpt_bz 806 ! ikpt=indkk(jkpt,1) 807 !!* ikpt = the point of IBZ that jkpt is an image of in BZ 808 ! npwj=npwarr(ikpt) 809 !!* npwj = number of planewaves in basis at point jkpt = at point ikpt 810 ! jsym=indkk(jkpt,2) 811 !!* jsym = symmetry operation to apply to get jkpt from ikpt 812 ! shiftg(:)=indkk(jkpt,3:5) 813 !!* shiftg = Bravais vector G0 to add to remain in BZ 814 ! if (jsym/=0) then 815 ! symm(:,:)=dtset%symrel(:,:,jsym) 816 ! tau_nons(:)=dtset%tnons(:,jsym) 817 !!* The symmetry operation in k-space (symm) and the non-symorphic translation (tau_nons) are now defined. 818 ! if(sum(tau_nons(:)**2)>tol8) then 819 !!* Initialize %calc_phase(jkpt) to 1 820 ! fock%calc_phase(jkpt)=1 821 !!* Compute the phase factor exp(i*2*pi*G.tau) for all G. 822 ! indx(1)=1 823 ! phase1d=zero 824 ! call getph(indx,1,n1,n2,n3,phase1d,tau_nons) 825 !!* Although the routine getph is orignally written for atomic phase factors, it does precisely what we want 826 ! arg=two_pi*(dtset%kptns(1,ikpt)*tau_nons(1) + dtset%kptns(2,ikpt)*tau_nons(2) & 827 !& + dtset%kptns(3,ikpt)*tau_nons(3)) 828 ! phktnons(1,1)=cos(arg) 829 ! phktnons(2,1)=sin(arg) 830 !! phktnons(1,1)=one 831 !! phktnons(2,1)=zero 832 !!* Convert 1D phase factors to 3D phase factors exp(i*2*pi*(k+G).tau) and store it in %phase_j 833 ! call ph1d3d(1,1,kg(:,1+tab_indikpt(1,ikpt):npwj+tab_indikpt(1,ikpt)),1,1,npwj,n1, & 834 !& n2,n3,phktnons,phase1d,fock%phase(:,1+jkg:npwj+jkg)) 835 ! end if 836 ! else 837 ! symm=0 ; symm(1,1)=1 ; symm(2,2)=1 ; symm(3,3)=1 838 ! tau_nons(:)=zero 839 ! shiftg(:)=0 840 ! end if 841 !!* Apply time-reversal symmetry if required 842 ! if(indkk(jkpt,6)/=0) then 843 !!* Initialize %timerev(jkpt) to 1 844 ! fock%timerev(jkpt)=1 845 ! symm(:,:)=-symm(:,:) 846 ! end if 847 848 !!* Initialize %istwfk_bz(jkpt) to 849 ! fock%istwfk_bz(jkpt)=dtset%istwfk(ikpt) 850 851 !!* Initialize %tab_ikpt and %tab_ibgcg 852 ! fock%tab_ikpt(jkpt)=ikpt 853 ! fock%tab_ibgcg(1:dtset%nsppol,jkpt)=tab_indikpt(2:1+dtset%nsppol,ikpt) 854 ! fock%tab_ibgcg(1+dtset%nsppol:2*dtset%nsppol,jkpt)= & 855 !& tab_indikpt(2+dtset%nsppol:2*dtset%nsppol+1,ikpt) 856 857 !!* Initialize %npwarr_bz 858 ! fock%npwarr_bz(jkpt)=npwj 859 860 !!* Initialize %kg_bz 861 ! do jpw=1,npwj 862 ! v1=kg(1,jpw+tab_indikpt(1,ikpt)) ; v2=kg(2,jpw+tab_indikpt(1,ikpt)) ; v3=kg(3,jpw+tab_indikpt(1,ikpt)) 863 ! fock%kg_bz(1,jpw+jkg)=-shiftg(1)+symm(1,1)*v1+symm(2,1)*v2+symm(3,1)*v3 864 ! fock%kg_bz(2,jpw+jkg)=-shiftg(2)+symm(1,2)*v1+symm(2,2)*v2+symm(3,2)*v3 865 ! fock%kg_bz(3,jpw+jkg)=-shiftg(3)+symm(1,3)*v1+symm(2,3)*v2+symm(3,3)*v3 866 !!* The symmetry operation symm must be transposed when used. (cf. docs about wfconv) 867 ! end do 868 869 !!* Initialize %gbound_bz 870 ! call sphereboundary(fock%gbound_bz(:,:,jkpt),fock%istwfk_bz(jkpt), & 871 !& fock%kg_bz(:,1+jkg:npwj+jkg),dtset%mgfft,npwj) 872 873 !!* Update of the shift to be applied 874 ! jkg=jkg+npwj 875 ! end do 876 877 ! ========================================================== 878 ! === Initialize the k-points in BZ and their properties === 879 ! ========================================================== 880 ! jkg=0; 881 882 do isym=1,dtset%nsym 883 call mati3inv(dtset%symrel(:,:,isym),symrec(:,:,isym)) 884 Rtnons (:,isym)= MATMUL(TRANSPOSE(symrec(:,:,isym)),dtset%tnons(:,isym)) 885 end do 886 ABI_MALLOC(fockcommon%symrec,(3,3,dtset%nsym)) 887 fockcommon%symrec=symrec 888 889 ABI_MALLOC(invsym,(dtset%nsym)) 890 invsym=0 891 ident(1,:3)=(/1,0,0/) 892 ident(2,:3)=(/0,1,0/) 893 ident(3,:3)=(/0,0,1/) 894 do isym=1,dtset%nsym 895 symm(:,:)=MATMUL(dtset%symrel(:,:,isym),dtset%symrel(:,:,isym)) 896 if (all(symm(:,:)==ident(:,:))) then 897 invsym(isym)=isym 898 else 899 do jsym=1,dtset%nsym 900 symm(:,:)=MATMUL(dtset%symrel(:,:,isym),dtset%symrel(:,:,jsym)) 901 if (all(symm(:,:)==ident(:,:))) then 902 invsym(isym)=jsym 903 cycle 904 end if 905 end do 906 end if 907 if(invsym(isym)==0) then 908 ABI_ERROR('No inverse has been found for isym') 909 end if 910 end do 911 912 jkg_this_proc=0;my_jkpt=0 913 !indkk(1:nkpt_bz,2)=(/1,1,3,1,11,7,9,1/) 914 do jkpt=1,nkpt_bz 915 916 !* If this processor does not calculate exchange with the k point jkpt, skip the rest of the k-point loop. 917 if (proc_distrb_cycle(mpi_enreg%distrb_hf,jkpt,1,mband,1,mpi_enreg%me_hf)) cycle 918 ! if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,jkpt,1,dtset%nbandhf,1,mpi_enreg%me_kpt))) then 919 !* The processor does own a copy of the array kg of ikpt ; increment the shift. 920 ! jkg=jkg+npwj 921 ! end if 922 ! Skip the rest of the k-point loop 923 ! cycle 924 ! end if 925 my_jkpt=my_jkpt+1 926 927 ikpt=indkk(jkpt,1) 928 !* ikpt = the point of IBZ that jkpt is an image of in BZ 929 npwj=npwarr(ikpt) 930 fockbz%npwarr(my_jkpt)=npwarr(ikpt) 931 !* npwj = number of planewaves in basis at point jkpt = at point ikpt 932 jsym=indkk(jkpt,2) 933 !* jsym = symmetry operation to apply to get jkpt from ikpt 934 fockbz%tab_symkpt(my_jkpt)=invsym(jsym) 935 shiftg(:)=indkk(jkpt,3:5) 936 !* shiftg = Bravais vector G0 to add to remain in BZ 937 938 !* Initialize the array %kptns_bz = the k points in full BZ 939 fockbz%kptns_bz(:,my_jkpt)=kptns_hf(:,jkpt) 940 941 !* Initialize the array %jstwfk = how is stored the wavefunction at each k point 942 if (dtset%istwfk(ikpt)/=1) then 943 fockbz%istwfk_bz(my_jkpt)=set_istwfk(kptns_hf(:,jkpt)) 944 end if 945 946 !* One can take advantage of the time-reversal symmetry in this case. 947 !* Initialize the array %wtk_bz = weight assigned to each k point. 948 ! fock%wtk_bz(my_jkpt)=dtset%wtk(jkpt)/ucvol 949 !* Caution, the definition takes into account "ucvol" ! 950 951 !* Initialize the array %npwarr_bz = number of planewaves in basis at each k point 952 ! fock%npwarr_bz(my_jkpt)=npwj 953 954 !!* Initialize the array %tab_ikpt = indices of k-point in IBZ ikpt for each k point jkpt in BZ (here,ikpt=jkpt) 955 fockbz%tab_ikpt(my_jkpt)=ikpt 956 957 958 !!* Initialize the array %tab_ibgcg = indices of cprj(ikpt)/occ(ikpt) and cg(ikpt) for each k point jkpt 959 ! if (my_nsppol==2) then 960 !!* In this case, my_nsppol=dtset%nsppol=2 961 ! fock%tab_ibgcg(1:2,my_jkpt)=tab_indikpt(2:3,ikpt) 962 ! fock%tab_ibgcg(3:4,my_jkpt)=tab_indikpt(4:5,ikpt) 963 ! else 964 ! if(mpi_enreg%my_isppoltab(1)==1) then 965 !!* In this case, my_nsppol=1 and the up spin is treated (dtset%nsppol= 1 or 2) 966 ! fock%tab_ibgcg(1,my_jkpt)=tab_indikpt(2,ikpt) 967 ! fock%tab_ibgcg(2,my_jkpt)=tab_indikpt(2+dtset%nsppol,ikpt) 968 ! else 969 !!* In this case, my_nsppol=1 and the dn spin is treated (so dtset%nsppol=2) 970 ! fock%tab_ibgcg(1,my_jkpt)=tab_indikpt(3,ikpt) 971 ! fock%tab_ibgcg(2,my_jkpt)=tab_indikpt(5,ikpt) 972 ! end if 973 ! end if 974 975 !* Initialize the array %kg_bz = reduced planewave coordinates at each k point 976 if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me_kpt))) then 977 !* We perform the test with isppol=-1 (both spins) and the occupied band (dtset%nbandhf). 978 !* We assume that paral_kgb==0 (a k-point may not be present on several proc.) 979 !* The array kg for ikpt is stored on this processor and copied in kg_tmp. 980 ikg=my_ikgtab(ikpt) 981 !* ikg = the shift in kg to get the G-vectors associated to ikpt 982 do ik=1,3 983 ! kg_tmp(1+(ik-1)*npwj:ik*npwj)=kg(ik,1+tab_indikpt(1,ikpt):npwj+tab_indikpt(1,ikpt)) 984 kg_tmp(1+(ik-1)*npwj:ik*npwj)=kg(ik,1+ikg:npwj+ikg) 985 end do 986 ! jkg=jkg+npwj 987 end if 988 !* Broadcast the array kg_tmp to all the processors of comm_kpt. 989 !* Since paral_kgb==0, all the bands of a k-point are treated on the same proc. 990 call xmpi_bcast(kg_tmp,mpi_enreg%proc_distrb(ikpt,1,1),mpi_enreg%comm_kpt,ier) 991 do ik=1,3 992 fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)=kg_tmp(1+(ik-1)*npwj:ik*npwj) 993 end do 994 995 !* Apply a symmetry operation on kg_bz if necessary 996 if (jsym/=0) then 997 symm(:,:)=dtset%symrel(:,:,jsym) 998 ! tau_nons(:)=dtset%tnons(:,jsym) 999 tau_nons(:)=-Rtnons(:,invsym(jsym)) 1000 !* The symmetry operation in k-space (symm) and the non-symorphic translation (tau_nons) are now defined. 1001 if(sum(tau_nons(:)**2)>tol8) then 1002 !* Initialize %calc_phase(jkpt) to 1 1003 fockbz%calc_phase(my_jkpt)=1 1004 !* Compute the phase factor exp(i*2*pi*G.tau) for all G. 1005 indx(1)=1 1006 phase1d=zero 1007 call getph(indx,1,n1,n2,n3,phase1d,tau_nons) 1008 !* Although the routine getph is orignally written for atomic phase factors, it does precisely what we want 1009 arg=two_pi*(dtset%kptns(1,ikpt)*tau_nons(1) + dtset%kptns(2,ikpt)*tau_nons(2) & 1010 & + dtset%kptns(3,ikpt)*tau_nons(3)) 1011 phktnons(1,1)=cos(arg) 1012 phktnons(2,1)=sin(arg) 1013 ! phktnons(1,1)=one 1014 ! phktnons(2,1)=zero 1015 !* Convert 1D phase factors to 3D phase factors exp(i*2*pi*(k+G).tau) and store it in %phase_j 1016 call ph1d3d(1,1,fockbz%kg_bz(:,1+jkg_this_proc:npwj+jkg_this_proc),1,1,npwj,n1,n2,n3, & 1017 & phktnons,phase1d,fockbz%phase(:,1+jkg_this_proc:npwj+jkg_this_proc)) 1018 end if 1019 !* Apply time-reversal symmetry if required 1020 if(indkk(jkpt,6)/=0) then 1021 !* Initialize %timerev(jkpt) to 1 1022 fockbz%timerev(my_jkpt)=1 1023 symm(:,:)=-symm(:,:) 1024 end if 1025 !* Initialize %kg_bz 1026 do jpw=1,npwj 1027 v1=fockbz%kg_bz(1,jpw+jkg_this_proc) ; v2=fockbz%kg_bz(2,jpw+jkg_this_proc) ; v3=fockbz%kg_bz(3,jpw+jkg_this_proc) 1028 fockbz%kg_bz(1,jpw+jkg_this_proc)=-shiftg(1)+symm(1,1)*v1+symm(2,1)*v2+symm(3,1)*v3 1029 fockbz%kg_bz(2,jpw+jkg_this_proc)=-shiftg(2)+symm(1,2)*v1+symm(2,2)*v2+symm(3,2)*v3 1030 fockbz%kg_bz(3,jpw+jkg_this_proc)=-shiftg(3)+symm(1,3)*v1+symm(2,3)*v2+symm(3,3)*v3 1031 !* The symmetry operation symm must be transposed when used. (cf. docs about wfconv) 1032 end do 1033 else 1034 !* Ths symmetry operation is the identity. 1035 !* Apply time-reversal symmetry if required 1036 if(indkk(jkpt,6)/=0) then 1037 !* Initialize %timerev(jkpt) to 1 1038 fockbz%timerev(my_jkpt)=1 1039 fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)=-fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc) 1040 end if 1041 end if 1042 1043 !* Initialize the array %gbound_bz = boundary of the basis sphere of G vectors at each k point 1044 call sphereboundary(fockbz%gbound_bz(:,:,my_jkpt),fockbz%istwfk_bz(my_jkpt),& 1045 & fockbz%kg_bz(:,1+jkg_this_proc:npwj+jkg_this_proc),mgfft,npwj) 1046 1047 jkg_this_proc=jkg_this_proc+npwj 1048 1049 !* Initialize the arrays %tab_ibg = shifts in arrays cprj and occ (ibg) for each k point jkpt 1050 !* Initialize the arrays %tab_icg = shifts in arrays cg(icg) for each k point jkpt 1051 if (my_nsppol==1) then 1052 fockbz%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,1+mpi_enreg%my_isppoltab(2)) 1053 fockbz%tab_icg(my_jkpt,1)=my_icgtab(ikpt,1+mpi_enreg%my_isppoltab(2)) 1054 fockbz%tab_icp(my_jkpt,1)=my_icptab(ikpt,1+mpi_enreg%my_isppoltab(2)) 1055 !* if mpy_isppoltab(2)=0, the up spin is treated (dtset%nsppol= 1 or 2) 1056 !* if mpy_isppoltab(2)=1, the dn spin is treated (so dtset%nsppol=2) 1057 1058 ! if(mpi_enreg%my_isppoltab(2)==1) then 1059 !* In this case, my_nsppol=1 and the up spin is treated (dtset%nsppol= 1 or 2) 1060 ! fock%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,1) 1061 ! fock%tab_icg(my_jkpt,1)=my_icgtab(ikpt,1) 1062 ! else 1063 !* In this case, my_nsppol=1 and the dn spin is treated (so dtset%nsppol=2) 1064 ! fock%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,2) 1065 ! fock%tab_icg(my_jkpt,1)=my_icgtab(ikpt,2) 1066 ! end if 1067 else 1068 !* In this case, my_nsppol=dtset%nsppol=2 1069 fockbz%tab_ibg(my_jkpt,:)=my_ibgtab(ikpt,:) 1070 fockbz%tab_icg(my_jkpt,:)=my_icgtab(ikpt,:) 1071 fockbz%tab_icp(my_jkpt,:)=my_icptab(ikpt,:) 1072 end if 1073 1074 enddo 1075 1076 !* Deallocation 1077 ABI_FREE(invsym) 1078 1079 end if 1080 ABI_FREE(indkk) 1081 ABI_FREE(kg_tmp) 1082 ABI_FREE(kptns_hf) 1083 ABI_FREE(my_ibgtab) 1084 ABI_FREE(my_icgtab) 1085 ABI_FREE(my_icptab) 1086 ABI_FREE(my_ikgtab) 1087 ABI_FREE(phase1d) 1088 call fock_print(fockcommon,fockbz,unit=std_out) 1089 1090 call timab(1501,2,tsec) 1091 1092 DBG_EXIT("COLL") 1093 1094 end subroutine fock_init
m_fock/fock_print [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_print
FUNCTION
Print info on the fock_type data type
INPUTS
fock<crystal_t>=The object [unit]=Unit number for output [prtvol]=Verbosity level [mode_paral]=Either "COLL" or "PERS" [header]=String to be printed as header for additional info.
OUTPUT
Only printing
SOURCE
1905 subroutine fock_print(fockcommon,fockbz,header,unit,mode_paral,prtvol) 1906 1907 !Arguments ------------------------------------ 1908 !scalars 1909 integer,optional,intent(in) :: unit,prtvol 1910 character(len=4),optional,intent(in) :: mode_paral 1911 character(len=*),optional,intent(in) :: header 1912 type(fock_common_type),intent(in) :: fockcommon 1913 type(fock_BZ_type),intent(in) :: fockbz 1914 1915 !Local variables------------------------------- 1916 integer :: my_unt,my_prtvol 1917 character(len=4) :: my_mode 1918 character(len=500) :: msg 1919 1920 ! ********************************************************************* 1921 1922 my_unt=std_out; if (PRESENT(unit)) my_unt=unit 1923 my_prtvol=0 ; if (PRESENT(prtvol)) my_prtvol=prtvol 1924 my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral 1925 1926 msg=' ==== Info on fock_type ==== ' 1927 if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== ' 1928 call wrtout(my_unt,msg,my_mode) 1929 1930 ! Important dimensions 1931 call wrtout(my_unt,sjoin(" my_nsppol ...",itoa(fockcommon%my_nsppol)),my_mode) 1932 call wrtout(my_unt,sjoin(" nkpt_bz .....",itoa(fockbz%nkpt_bz)),my_mode) 1933 1934 ! Options 1935 call wrtout(my_unt,sjoin(" nnsclo_hf .......",itoa(fockcommon%nnsclo_hf)),my_mode) 1936 call wrtout(my_unt,sjoin(" ixc .............",itoa(fockcommon%ixc)),my_mode) 1937 call wrtout(my_unt,sjoin(" hybrid mixing....",ftoa(fockcommon%hyb_mixing)),my_mode) 1938 call wrtout(my_unt,sjoin(" hybrid SR mixing ",ftoa(fockcommon%hyb_mixing_sr)),my_mode) 1939 call wrtout(my_unt,sjoin(" hybrid range DFT ",ftoa(fockcommon%hyb_range_dft)),my_mode) 1940 call wrtout(my_unt,sjoin(" hybrid range Fock",ftoa(fockcommon%hyb_range_fock)),my_mode) 1941 1942 ! write(msg,"(a,f12.1,a)")" Memory required for HF u(r) states: ",product(shape(fockbz%cwaveocc_bz)) * dp * b2Mb, " [Mb]" 1943 ! call wrtout(my_unt,msg,my_mode) 1944 1945 ! Extra info. 1946 !if (my_prtvol > 0) then 1947 ! call wrtout(my_unt,"Extra info not available",my_mode) 1948 !end if 1949 1950 end subroutine fock_print
m_fock/fock_set_getghc_call [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_set_getghc_call
FUNCTION
Set the value of fock%getghc_call, Returns the old value
SOURCE
1845 integer function fock_set_getghc_call(fock, new) result(old) 1846 1847 !Arguments ------------------------------------ 1848 !scalars 1849 integer,intent(in) :: new 1850 type(fock_common_type),intent(inout) :: fock 1851 1852 ! ************************************************************************* 1853 1854 old = fock%getghc_call_ 1855 fock%getghc_call_ = new 1856 1857 end function fock_set_getghc_call
m_fock/fock_set_ieigen [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_set_ieigen
FUNCTION
Set the value of ieigen to the value given in argument.
INPUTS
fock <type(fock_type)>= all the quantities to calculate Fock exact exchange iband= index of the band iband
OUTPUT
none
SOURCE
1160 subroutine fock_set_ieigen(fock,iband) 1161 1162 !Arguments ------------------------------------ 1163 integer, intent(in) :: iband 1164 type(fock_common_type),pointer :: fock 1165 1166 ! ************************************************************************* 1167 1168 !Nothing to do if fock pointer is not associated... 1169 1170 ! ====================================================== 1171 ! === Update the data relative to the current states === 1172 ! ====================================================== 1173 1174 !* Copy of the value iband in the field ieigen 1175 if (associated(fock)) then 1176 fock%ieigen=iband 1177 fock%iband=iband 1178 end if 1179 1180 end subroutine fock_set_ieigen
m_fock/fock_type [ Types ]
NAME
fock_type
FUNCTION
This object stores the occupied wavefunctions and other quantities needed to calculate Fock exact exchange
SOURCE
66 type, public :: fock_type 67 type(fock_common_type), pointer :: fock_common=> null() 68 type(fock_BZ_type), pointer :: fock_BZ=> null() 69 type(fock_ACE_type), pointer :: fockACE(:,:)=> null() 70 end type fock_type 71 72 73 type, public :: fock_common_type 74 75 ! Integer scalars 76 !integer :: mcgocc_bz,mkg_bz,mocc 77 !integer :: natom,ntypat 78 79 integer :: usepaw 80 ! 0 if norm-conserving psps, 1 for PAW (not implemented) 81 82 integer :: ikpt,isppol,ieigen,iband 83 ! data relative to the current states. 84 85 integer :: mband 86 ! maximum number of bands 87 88 integer :: my_nsppol 89 ! my_nsppol=1 when nsppol=1 or nsppol=2 and only one spin is treated by the processor. 90 ! my_nsppol=2 when nsppol=2 and no parallelization over kpt (both spins are treated by the processor). 91 92 integer :: natom 93 ! Number of atoms, input variable 94 95 integer :: nsppol 96 ! Number of independent spin polarizations, input variable 97 ! Note that this value does not take into account the MPI distribution of the wavefunctions. 98 99 integer :: ntypat 100 ! Number of type of atoms 101 102 integer :: nnsclo_hf 103 ! Number of iterations with fixed occupied states when calculating the exact exchange contribution. 104 105 integer :: ixc 106 ! XC option (abinit input variable) 107 108 integer :: use_ACE 109 ! option to use the ACE method of Lin Lin 110 !==0 if the normal Fock operator is to be created and/or used 111 !==1 if the ACE operator is to be created and/or used 112 113 integer ABI_PRIVATE :: getghc_call_ = 1 114 ! 1 if fock_getghc should be called in getghc, 0 otherwise 115 116 ! Logical 117 logical :: optfor 118 ! option to calculate forces 119 120 logical :: optstr 121 ! option to calculate stresses 122 123 logical :: fock_converged 124 ! .false. if the Fock cycle (with changing Fock/ACE operator) is not converged 125 ! .true. if the Fock cycle (with changing Fock/ACE operator) has converged 126 127 logical :: scf_converged 128 ! .false. if the SCF cycle (with fixed Fock/ACE operator) is not converged 129 ! .true. if the SCF cycle (with fixed Fock/ACE operator) has converged 130 131 ! Real(dp) scalars 132 133 real(dp) :: gsqcut 134 ! cutoff value on G**2 for sphere inside fft box. 135 ! (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)). Used in hartre 136 137 real(dp) :: hyb_mixing 138 ! hybrid mixing coefficient for the Fock contribution 139 140 real(dp) :: hyb_mixing_sr 141 ! hybrid mixing coefficient for the short-range Fock contribution 142 143 real(dp) :: hyb_range_dft 144 ! hybrid range for separation, used in the DFT functional 145 ! (should be equal to hyb_range_fock, but this is not true for HSE03) 146 147 real(dp) :: hyb_range_fock 148 ! hybrid range for separation, used in the fock contribution 149 150 real(dp) :: e_fock0 151 ! contribution of the Fock term to energy (computed and stored here in case of ACE) 152 153 integer, allocatable :: atindx(:) 154 ! atindx(natom)=index table for atoms (see gstate.f) 155 156 integer, allocatable :: nband(:) 157 ! nband(nkpt) 158 ! Number of bands for each k point 159 160 integer, allocatable :: symrec(:,:,:) 161 162 integer,allocatable :: typat(:) 163 ! typat(natom) 164 ! type of each atom 165 166 ! Real(dp) arrays 167 real(dp) :: stress(6) 168 ! stress(6) 169 ! contribution of the fock term to stresses 170 171 real(dp), allocatable :: stress_ikpt(:,:) 172 ! stress(6,nband) 173 ! contribution of the fock term to stresses for the current band 174 175 real(dp), allocatable :: forces_ikpt(:,:,:) 176 ! forces(3,natom,nband)) 177 ! contribution of the fock term to forces for the current band 178 179 real(dp), allocatable :: forces(:,:) 180 ! forces(3,natom)) 181 ! contribution of the fock term to forces 182 183 real(dp), allocatable :: eigen_ikpt(:) 184 ! eigen_ikpt,(nband)) 185 ! Will contain the band index of the current state 186 ! if the value is 0, the Fock contribution to the eigenvalue is not calculated. 187 188 ! Pointers to PAW-types (associated only if usepaw==1) 189 ! Note that these are references to already existing objects. 190 191 type(pawtab_type), pointer :: pawtab(:) => null() 192 type(pawfgr_type),pointer :: pawfgr => null() 193 type(pawfgrtab_type),allocatable :: pawfgrtab(:) 194 195 end type fock_common_type 196 197 type, public :: fock_BZ_type 198 199 integer :: mcprj 200 ! dimension of cwaveocc_cprj 201 202 integer :: mkpt 203 ! maximum number of k-points for Fock treated by this node 204 205 integer :: mkptband 206 ! size of occupied states stored by this node. 207 208 integer :: nkpt_bz 209 ! Number of k-points in the BZ for Fock operator 210 211 integer, allocatable :: gbound_bz(:,:,:) 212 ! gbound_bz(2*mgfft+8,2,mkpt) 213 ! Tables for zero-padded FFT of wavefunctions. 214 215 integer, allocatable :: kg_bz(:,:) 216 ! kg_bz(3,mpw*mkpt) 217 ! G-vectors for each k-point in the BZ treate by this node 218 219 integer, allocatable :: nbandocc_bz(:,:) 220 ! nbandocc_bz,(mkpt,my_nsppol)) 221 ! nb of bands at each k point 222 223 integer, allocatable :: npwarr(:) 224 ! npwarr(mkpt) 225 226 integer, allocatable :: istwfk_bz(:) 227 ! istwfk_bz,(mkpt)) 228 ! storage mode of the wavefunction at each k-point 229 230 integer, allocatable :: calc_phase(:) 231 ! calc_phase,(mkpt)) 232 ! 1 if a phase factor must be considered (0 otherwise) at each k point 233 234 integer, allocatable :: tab_symkpt(:) 235 ! tab_symkpt,(mkpt)) 236 ! indices of symmetry operation to apply to get jkpt in full BZ from ikpt in IBZ 237 238 integer, allocatable :: timerev(:) 239 ! timerev,(mkpt)) 240 ! 1 if time reversal symmetry must be used (0 otherwise) at each k point 241 242 integer, allocatable :: tab_ibg(:,:) 243 ! tab_ibg,(mkpt,my_nsppol)) 244 ! indices of cprj(ikpt)/occ(ikpt) in the arrays cprj/occ for each k-point jkpt 245 246 integer, allocatable :: tab_icg(:,:) 247 ! tab_icg,(mkpt,my_nsppol)) 248 ! indices of cg(ikpt) in the arrays cg for each k-point jkpt 249 250 integer, allocatable :: tab_icp(:,:) 251 ! tab_icg,(mkpt,my_nsppol)) 252 ! indices of cprj(ikpt) in the arrays cprj for each k-point jkpt 253 254 integer, allocatable :: tab_ikpt(:) 255 ! tab_ikpt,(mkpt)) 256 ! indices of k-point ikpt in IBZ which corresponds to each k-point jkpt in full BZ 257 258 real(dp), allocatable :: cgocc(:,:,:) 259 ! cgocc(2,npw*mkptband,my_nsppol) 260 ! wavefunction in the G-space 261 262 real(dp), allocatable :: cwaveocc_bz(:,:,:,:,:,:) 263 ! (2,n4,n5,n6,mkptband,my_nsppol)) 264 ! occupied states of each bands at each k point (used to construct Fock operator), in the real space 265 real(dp), allocatable :: occ_bz(:,:) 266 ! occ_bz(mkptband,my_nsppol)) 267 ! occupancy of each bands at each k point 268 269 real(dp), allocatable :: wtk_bz(:) 270 ! wtk_bz,(mkpt)) 271 ! weights assigned to each k point in the BZ 272 ! Caution, the definition takes into account "ucvol" ! 273 274 real(dp), allocatable :: kptns_bz(:,:) 275 ! kptns_bz(3,mkpt) 276 ! k-points in full BZ 277 278 real(dp), allocatable :: phase(:,:) 279 ! phase(2,mpw*mkpt)) 280 ! phase factor the cg array will be multiplied with at each k point 281 282 type(MPI_type) :: mpi_enreg 283 type(pawang_type),pointer :: pawang 284 type(pawcprj_type), allocatable :: cwaveocc_prj(:,:) 285 286 end type fock_BZ_type 287 !---------------------------------------------------------------------- 288 289 type,public :: fock_ACE_type 290 291 real(dp), allocatable :: xi(:,:,:) 292 293 end type fock_ACE_type 294 !---------------------------------------------------------------------- 295 296 public :: fock_init ! Initialize the object. 297 !public :: fock_from_wfk ! Initialize the object from external WFK file. 298 public :: fock_set_ieigen ! Set the value of ieigen to the value given in argument. 299 public :: fock_updateikpt ! Update the value of energies%e_xc and energies%e_xcdc with Fock contribution. 300 public :: fock_destroy ! Free memory. 301 public :: fock_ACE_destroy ! Free memory. 302 public :: fock_common_destroy ! Free memory. 303 public :: fock_bz_destroy ! Free memory. 304 public :: fock_calc_ene ! Calculate the Fock contribution to the total energy. 305 public :: fock_update_exc ! Update the value of energies%e_xc and energies%e_xcdc with Fock contribution. 306 public :: fock_updatecwaveocc ! Update in the fock datastructure the fields relative to the occupied states. 307 public :: fock_set_getghc_call ! Enable/disable the call to fock_getghc in getghc. 308 public :: fock_get_getghc_call ! Return the value of the flag used to enable/disable the call to fock_getghc in getghc. 309 public :: fock_print ! Print info on the object.
m_fock/fock_update_exc [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_update_exc
FUNCTION
Update the value of energies%e_xc and energies%e_xcdc with Fock contribution
INPUTS
OUTPUT
none energies <type(energies_type)>=storage for energies computed here : | e_fock= Fock contribution to the total energy (Hartree)
NOTES
If the cgocc_bz are not updated at each iteration, be careful to calculate Fock energy at the same frequency.
SOURCE
1423 subroutine fock_update_exc(fock_energy,xc_energy,xcdc_energy) 1424 1425 !Arguments ------------------------------------ 1426 real(dp),intent(in) :: fock_energy 1427 real(dp),intent(inout) :: xc_energy,xcdc_energy 1428 1429 ! ************************************************************************* 1430 1431 !xc_energy = fock%hyb_mixing*fock_energy 1432 !xcdc_energy = two*fock%hyb_mixing*fock_energy 1433 xc_energy = fock_energy 1434 xcdc_energy = two*fock_energy 1435 !CMartins : For an atom, ewald should be set to zero (at the beginning of the loop) and 1436 !the contribution in !|q+G|=0 should be an approximation to the missing component of Vloc in G=0 1437 !energies%e_ewald=energies%e_ewald-half*fock%divgq0*fock%wtk_bz(1)*piinv 1438 1439 end subroutine fock_update_exc
m_fock/fock_updatecwaveocc [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_updatecwaveocc
FUNCTION
Update in the fock datastructure the fields relative to the occupied states.
INPUTS
cg(2,mcg)= Input wavefunctions cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors dtset <type(dataset_type)>=all input variables for this dataset fock <type(fock_type)>= all the quantities to calculate Fock exact exchange indsym(4,nsym,natom) :: 1:3 shift, and 4 final atom, of symmetry isym operating on iatom (S^{-1}(R - t) = r0 + L, see symatm.F90 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mpi_enreg=information about MPI parallelization npwarr(nkpt)=number of planewaves in basis at this k point occ(mband*nkpt*nsppol)= occupation number for each band (often 2) at each k point ucvol= unit cell volume ($\textrm{bohr}^{3}$)
OUTPUT
SIDE EFFECTS
The field fock%cgocc_bz contains the table cg at the end. The fields kg_bz, occ_bz and fock%cwaveocc_prj are simultaneously updated.
NOTES
############################ ### Not fully tested yet ### ############################ May be improved by selecting only the occupied states with the same spin isppol.
SOURCE
1481 subroutine fock_updatecwaveocc(cg,cprj,dtset,fock,indsym,mcg,mcprj,& 1482 & mpi_enreg,nattyp,npwarr,occ,ucvol) 1483 1484 !scalars 1485 integer, intent(in) :: mcg,mcprj 1486 real(dp), intent(in) :: ucvol 1487 type(dataset_type),intent(in) :: dtset 1488 type(fock_type),intent(inout),pointer :: fock 1489 type(MPI_type),intent(in) :: mpi_enreg 1490 !arrays 1491 integer, intent(in) :: indsym(4,dtset%nsym,dtset%natom),nattyp(dtset%ntypat),npwarr(dtset%nkpt) 1492 real(dp),intent(in) :: cg(2,mcg),occ(dtset%mband*dtset%nkpt*dtset%nsppol) 1493 type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj) 1494 1495 !Local variables------------------------------- 1496 !scalars 1497 integer,parameter :: tim_fourwf0=0 1498 integer :: iatm,iatom,iband,iband0,iband_cprj,ibg,icg,icp,ier,ikpt,ilmn,isize,ispinor,isppol,itypat,jbg,jcg,jkg,jkpt,jpw,jstwfk!,ii1,ii2 1499 integer :: lmnmax,mband,mband0,mgfft,mkpt,mpw,my_jsppol,my_jband,my_jkpt 1500 integer :: nband,ncpgr,n4,n5,n6,nkpt_bz,npwj,nsppol,nspinor 1501 real(dp),parameter :: weight1=one 1502 real(dp) :: cgre,cgim,invucvol 1503 character(len=500) :: message 1504 ! arrays 1505 integer :: ngfft(18) 1506 integer, ABI_CONTIGUOUS pointer :: gbound_k(:,:),kg_k(:,:) 1507 integer,allocatable :: dimlmn(:),indlmn(:,:,:),indsym_(:,:,:),typat_srt(:) 1508 real(dp) :: tsec(2),tsec2(2),dcp(3) 1509 real(dp),allocatable :: cgocc_tmp(:),cgocc(:,:),dummytab2(:,:),dummytab3(:,:,:),phase_jkpt(:,:) 1510 type(pawcprj_type),allocatable :: cprj_tmp(:,:) 1511 type(fock_common_type),pointer :: fockcommon 1512 type(fock_BZ_type),pointer :: fockbz 1513 ! ************************************************************************* 1514 1515 call timab(1502,1,tsec) 1516 1517 ABI_CHECK(associated(fock),"fock must be associated") 1518 1519 if (associated(fock)) then 1520 1521 fockcommon=>fock%fock_common 1522 fockbz=> fock%fock_BZ 1523 1524 invucvol=1.d0/sqrt(ucvol) 1525 ! Local variables = useful dimensions 1526 mband=fockcommon%mband 1527 mkpt=fockbz%mkpt 1528 mpw=dtset%mpw 1529 mgfft=dtset%mgfft 1530 ngfft=dtset%ngfft 1531 nkpt_bz=fockbz%nkpt_bz 1532 nsppol=dtset%nsppol 1533 nspinor=1 1534 ncpgr=0 1535 1536 ! Local variables : useful arrays 1537 ABI_MALLOC(cgocc,(2,mpw)) 1538 cgocc=zero 1539 ABI_MALLOC(cgocc_tmp,(2*mpw+1)) 1540 cgocc_tmp=zero 1541 if (fockcommon%usepaw==1) then 1542 mgfft=dtset%mgfftdg 1543 ngfft=dtset%ngfftdg 1544 ABI_MALLOC(cprj_tmp,(dtset%natom,nspinor)) 1545 ABI_MALLOC(dimlmn,(dtset%natom)) 1546 call pawcprj_getdim(dimlmn,dtset%natom,nattyp,dtset%ntypat,dtset%typat,fockcommon%pawtab,"O") 1547 ncpgr = 0 1548 if (dtset%optforces== 1) ncpgr = 3 1549 ! if (dtset%optstress /= 0) ncpgr = 6 1550 call pawcprj_alloc(cprj_tmp,ncpgr,dimlmn) 1551 1552 lmnmax=maxval(fockcommon%pawtab(:)%lmn_size) 1553 ABI_MALLOC(indlmn,(6,lmnmax,dtset%ntypat)) 1554 do itypat=1,dtset%ntypat 1555 isize=size(fockcommon%pawtab(itypat)%indlmn,2) 1556 indlmn(:,1:isize,itypat)=fockcommon%pawtab(itypat)%indlmn(:,1:isize) 1557 end do 1558 ABI_MALLOC(indsym_,(4,dtset%nsym,dtset%natom)) 1559 ABI_MALLOC(typat_srt,(dtset%natom)) 1560 1561 if (dtset%nsym==1) then 1562 indsym_=0 1563 do iatom=1,dtset%natom 1564 iatm=fockcommon%atindx(iatom) 1565 typat_srt(iatm)=dtset%typat(iatom) 1566 indsym_(4,:,iatom)=iatom 1567 end do 1568 else 1569 do iatom=1,dtset%natom 1570 iatm=fockcommon%atindx(iatom) 1571 typat_srt(iatm)=dtset%typat(iatom) 1572 indsym_(1:3,:,iatm)=indsym(1:3,:,iatom) 1573 indsym_(4,:,iatm)=fockcommon%atindx(indsym(4,:,iatom)) 1574 end do 1575 end if 1576 end if 1577 1578 ! Local variables to perform FFT 1579 n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6) 1580 ABI_MALLOC(dummytab3,(n4,n5,n6)) 1581 1582 1583 if(ANY(fockbz%calc_phase(:)/=0)) then 1584 ABI_MALLOC(phase_jkpt,(2,mpw)) 1585 phase_jkpt=zero 1586 end if 1587 1588 1589 ! ======================================================= 1590 ! === Update the data relative to the occupied states === 1591 ! ======================================================= 1592 !* The arrays cgocc_bz, kg_bz, occ_bz and npwarr_bz are already allocated with the maximal size. 1593 ! if ((dtset%kptopt>=1).and.(dtset%kptopt<=4)) then 1594 ! if (dtset%kptopt/=3) then 1595 1596 do isppol=1,nsppol 1597 jbg=0 ; jcg=0 ; jkg=0 ; icp=0 1598 my_jsppol=isppol 1599 if ((isppol==2).and.(mpi_enreg%nproc_spkpt/=1)) my_jsppol=1 1600 !* Both spins are treated on the same proc., only in the case where nproc_spkpt=1; 1601 !* otherwise each proc. treats only one spin. 1602 1603 ! MG: This loop is not effient! 1604 ! Ok the number of k-points in the BZ is usually `small` when hybrids are used 1605 ! but what happens if we use a 12x12x12. 1606 ! One should loop over the IBZ, broadcast and reconstruct the star of the k-point. 1607 my_jkpt=0 1608 do jkpt=1,nkpt_bz 1609 1610 if (proc_distrb_cycle(mpi_enreg%distrb_hf,jkpt,1,mband,1,mpi_enreg%me_hf)) cycle 1611 1612 !* In this case, the processor does not calculate the exchange with any occupied state on jkpt. 1613 1614 ! if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,jkpt,1,dtset%nbandhf,jsppol,mpi_enreg%me_kpt))) then 1615 !* The state (jkpt,jband,jsppol) is stored in the array cg of this processor and copied in cgocc_tmp. 1616 ! icg=icg+dtset%nband(jkpt)*npwj 1617 ! end if 1618 ! ibg=ibg+dtset%nband(jkpt) 1619 ! cycle 1620 ! end if 1621 my_jkpt=my_jkpt+1 1622 1623 ikpt=fockbz%tab_ikpt(my_jkpt) 1624 1625 !* ikpt = the point of IBZ that jkpt is an image of in BZ 1626 npwj=npwarr(ikpt) 1627 !* npwj= number of plane wave in basis for the wavefunction 1628 jstwfk=fockbz%istwfk_bz(my_jkpt) 1629 !* jstwfk= how is stored the wavefunction 1630 ibg=fockbz%tab_ibg(my_jkpt,my_jsppol) 1631 !* ibg = shift to be applied on the location of data in the array occ 1632 icg=fockbz%tab_icg(my_jkpt,my_jsppol) 1633 !* icg = shift to be applied on the location of data in the array cg 1634 icp=fockbz%tab_icp(my_jkpt,my_jsppol) 1635 !* icp = shift to be applied on the location of data in the array cprj 1636 gbound_k => fockbz%gbound_bz(:,:,my_jkpt) 1637 !* boundary of the basis sphere of G vectors 1638 kg_k => fockbz%kg_bz(:,1+jkg:npwj+jkg) 1639 !* reduced plean wave coordinates 1640 if (fockbz%calc_phase(my_jkpt)==1) then 1641 phase_jkpt(:,1:npwj)=fockbz%phase(:,1+jkg:npwj+jkg) 1642 end if 1643 !* phase factor at k-point j 1644 1645 !* Initialize the band counter 1646 my_jband=0 1647 !FBru: Here run over all the nbandhf bands instead of just the truly occupied states 1648 ! this is a very tiny waste, but solves parallelization issues 1649 do iband=1,dtset%nbandhf 1650 1651 1652 cgocc_tmp=zero 1653 if (fockcommon%usepaw==1) then 1654 call pawcprj_set_zero(cprj_tmp) 1655 end if 1656 1657 !* To avoid segmentation fault, my_jband should not be greater than nbandhf 1658 !FBru: this error should never happen in practice since we now limit the loop to nbandhf 1659 if ((my_jband+1)>mband) then 1660 write(message,*) 'The number of occupied band',my_jband+1,' at k-point',& 1661 & ikpt,' is greater than the value of nbandhf ', mband 1662 ABI_ERROR(message) 1663 end if 1664 1665 !* If the processor does not calculate the exchange with the occupied state (jkpt,my_jband), cycle 1666 ! if (mpi_enreg%distrb_hf(jkpt,(my_jband+1),1)/=mpi_enreg%me_hf) cycle 1667 if (mpi_enreg%distrb_hf(jkpt,iband,1)/=mpi_enreg%me_hf) cycle 1668 ! if (mpi_enreg%proc_distrb(jkpt,jband,jsppol)==mpi_enreg%me_kpt) then 1669 !* The state (jkpt,jband,jsppol) is stored in the array cg of this processor ; shift are incremented. 1670 ! icg=icg+npwj 1671 ! end if 1672 ! ibg=ibg+1 1673 !* Skip the end of the loop 1674 ! cycle 1675 ! end if 1676 1677 !* increment the number of occupied bands treated on this processor 1678 my_jband = my_jband+1 1679 1680 !* In this case, the processor calculates the exchange with the occupied state (jkpt,my_jband). 1681 if (mpi_enreg%proc_distrb(ikpt,iband,isppol)==mpi_enreg%me_kpt) then 1682 !* The state (ikpt,iband,isppol) is stored in the array cg of this processor and copied in cgocc_tmp. 1683 if(icg==-1) then 1684 write(100,*) 'icg=-1',mpi_enreg%me,isppol,my_jsppol,jkpt,my_jkpt,ikpt,iband 1685 end if 1686 ! MG: Why packing re and im part? 1687 cgocc_tmp(1)=occ(iband+ibg) 1688 cgocc_tmp(2:npwj+1)=cg(1,1+(iband-1)*npwj+icg:iband*npwj+icg) 1689 cgocc_tmp(npwj+2:2*npwj+1)=cg(2,1+(iband-1)*npwj+icg:iband*npwj+icg) 1690 if (fockcommon%usepaw==1) then 1691 call pawcprj_copy(cprj(:,icp+iband:icp+iband+nspinor-1),cprj_tmp) 1692 end if 1693 end if 1694 1695 !* Broadcast the state (ikpt,iband,isppol) to all the processors of comm_kpt for cgocc 1696 call timab(1503,1,tsec2) 1697 call xmpi_bcast(cgocc_tmp,mpi_enreg%proc_distrb(ikpt,iband,isppol),mpi_enreg%comm_kpt,ier) 1698 1699 !* Broadcast the state (ikpt,iband,isppol) to all the processors of comm_kpt for cprj 1700 if (fockcommon%usepaw==1) then 1701 call pawcprj_bcast(cprj_tmp,dtset%natom,nspinor,dimlmn,ncpgr,mpi_enreg%proc_distrb(ikpt,iband,isppol),& 1702 & mpi_enreg%comm_kpt,ier) 1703 end if 1704 call timab(1503,2,tsec2) 1705 !* Keep the processors in %comm_kpt which needs the values in cgocc_tmp to build their own %cwaveocc and %occ_bz. 1706 if ((mpi_enreg%nproc_spkpt/=1).and.(nsppol==2)) then 1707 if (fockbz%timerev(my_jkpt)==mpi_enreg%my_isppoltab(isppol)) cycle 1708 !* In the case of a parallel spin-polarized calculation 1709 !* when time reversal symmetry is applied at this k-point (timrev==1), only the processors with the opposite spin (my_isppoltab==0) are kept. 1710 !* when time reversal symmetry is not applied at this k-point (timrev==0), only the processors with the same spin (my_isppoltab==1) are kept. 1711 1712 ! if (fock%timerev(my_jkpt)==1)) then 1713 ! if (mpi_enreg%my_isppoltab(isppol)==1) cycle 1714 !* In the case of a parallel spin-polarized calculation and when time reversal symmetry is applied at this k-point, 1715 !* only the processors with the opposite spin are kept. 1716 ! else 1717 ! if (mpi_enreg%my_isppoltab(isppol)==0) cycle 1718 !* only the processors with isppol are kept. 1719 ! end if 1720 end if 1721 1722 !* Copy the values of cgocc_tmp in the arrays cgocc and %occ_bz 1723 fockbz%occ_bz(my_jband+jbg,my_jsppol) = cgocc_tmp(1) 1724 cgocc(1,1:npwj) = cgocc_tmp(2:npwj+1) 1725 cgocc(2,1:npwj) = cgocc_tmp(npwj+2:2*npwj+1) 1726 1727 !* calculate cg and store it in cgocc_bz 1728 if (fockbz%calc_phase(my_jkpt)==1) then 1729 do jpw=1,npwj 1730 cgre=cgocc(1,jpw) ; cgim=cgocc(2,jpw) 1731 cgocc(1,jpw) = phase_jkpt(1,jpw)*cgre - phase_jkpt(2,jpw)*cgim 1732 cgocc(2,jpw) = phase_jkpt(1,jpw)*cgim + phase_jkpt(2,jpw)*cgre 1733 end do 1734 end if ! phase 1735 1736 !* apply time reversal symmetry if necessary 1737 if (fockbz%timerev(my_jkpt)==1) then 1738 cgocc(2,:) = - cgocc(2,:) 1739 if((mpi_enreg%nproc_spkpt==1).and.(nsppol==2)) my_jsppol=mod(my_jsppol,2)+1 1740 !* exchange spin (1 ->2 ; 2-> 1) in the sequential case. 1741 end if 1742 1743 !* apply FFT to get cwaveocc in real space 1744 1745 if (allocated(fockbz%cwaveocc_bz)) then 1746 1747 ABI_MALLOC(dummytab2,(2,npwj)) 1748 call fourwf(1,dummytab3,cgocc(:,1:npwj),dummytab2,fockbz%cwaveocc_bz(:,:,:,:,my_jband+jbg,my_jsppol), & 1749 & gbound_k,gbound_k,jstwfk,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,& 1750 & npwj,npwj,n4,n5,n6,tim_fourwf0,0,weight1,weight1,gpu_option=dtset%gpu_option) 1751 ABI_FREE(dummytab2) 1752 1753 else 1754 fockbz%cgocc(:,jcg+1+(my_jband-1)*npwj:jcg+my_jband*npwj,my_jsppol)=cgocc(:,1:npwj) 1755 end if 1756 1757 !* calculate cprj and store it in cwaveocc_prj 1758 if (fockcommon%usepaw==1) then 1759 iband_cprj=(my_jsppol-1)*fockbz%mkptband+jbg+my_jband 1760 nband=1;mband0=1;iband0=1 1761 call pawcprj_symkn(fockbz%cwaveocc_prj(:,iband_cprj:iband_cprj+nspinor-1),cprj_tmp(:,1:nspinor),& 1762 & indsym_,dimlmn,iband0,indlmn,& 1763 & fockbz%tab_symkpt(my_jkpt),fockbz%timerev(my_jkpt),dtset%kptns(:,ikpt),fockbz%pawang%l_max-1,lmnmax,& 1764 & mband0,dtset%natom,nband,nspinor,dtset%nsym,dtset%ntypat,typat_srt,fockbz%pawang%zarot) 1765 1766 if(dtset%optforces==1) then 1767 do iatom=1,dtset%natom 1768 iatm=fockcommon%atindx(iatom) 1769 do ispinor=iband_cprj,iband_cprj+nspinor-1 1770 do ilmn=1,fockcommon%pawtab(dtset%typat(iatom))%lmn_size 1771 dcp(:)= MATMUL(TRANSPOSE(fockcommon%symrec(:,:,fockbz%tab_symkpt(my_jkpt))),& 1772 & fockbz%cwaveocc_prj(iatm,ispinor)%dcp(1,:,ilmn)) 1773 fockbz%cwaveocc_prj(iatm,ispinor)%dcp(1,:,ilmn)=dcp(:) 1774 dcp(:)= MATMUL(TRANSPOSE(fockcommon%symrec(:,:,fockbz%tab_symkpt(my_jkpt))),& 1775 & fockbz%cwaveocc_prj(iatm,ispinor)%dcp(2,:,ilmn)) 1776 fockbz%cwaveocc_prj(iatm,ispinor)%dcp(2,:,ilmn)=dcp(:) 1777 end do 1778 end do 1779 end do 1780 end if 1781 1782 end if 1783 1784 !* update the shift to apply to occ in all case because this array is not distributed among the proc. 1785 ! ibg=ibg+1 1786 1787 end do ! iband 1788 1789 !* Save the true number of occupied bands in the array %nbandocc_bz 1790 fockbz%nbandocc_bz(my_jkpt,my_jsppol) = my_jband 1791 1792 !* update the shifts to apply 1793 jbg=jbg+my_jband 1794 jcg=jcg+npwj*my_jband 1795 jkg=jkg+npwj 1796 end do ! ikpt 1797 end do ! isppol 1798 if (allocated(fockbz%cwaveocc_bz)) then 1799 fockbz%cwaveocc_bz=fockbz%cwaveocc_bz*invucvol 1800 end if 1801 1802 ABI_FREE(cgocc_tmp) 1803 ABI_FREE(cgocc) 1804 if (fockcommon%usepaw==1) then 1805 ABI_FREE(indlmn) 1806 ABI_FREE(indsym_) 1807 ABI_FREE(typat_srt) 1808 ABI_FREE(dimlmn) 1809 call pawcprj_free(cprj_tmp) 1810 ABI_FREE(cprj_tmp) 1811 end if 1812 if(allocated(phase_jkpt)) then 1813 ABI_FREE(phase_jkpt) 1814 end if 1815 ABI_FREE(dummytab3) 1816 1817 1818 ! Restricted or unrestricted HF 1819 if (nsppol==1) then 1820 !* Update the array %occ_bz => May be limited to the occupied states only 1821 fockbz%occ_bz(:,:)=half*fockbz%occ_bz(:,:) 1822 1823 ! If nsppol=1, this is a restricted Hartree-Fock calculation. 1824 ! If nsppol=2, this is an unrestricted Hartree-Fock calculation. 1825 end if 1826 1827 end if 1828 1829 call timab(1502,2,tsec) 1830 1831 end subroutine fock_updatecwaveocc
m_fock/fock_updateikpt [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fock_updateikpt
FUNCTION
Update the value of ikpt,isppol for the next exact exchange calculation.
INPUTS
fock <type(fock_type)>= all the quantities to calculate Fock exact exchange ikpt= k-point index isppol= Spin index
SIDE EFFECTS
The field fock%eigen_ikpt is also set to 0.d0.
NOTES
May be improved to calculate the star of ikpt. => I think NO finally
SOURCE
1117 subroutine fock_updateikpt(fock,ikpt,isppol) 1118 1119 !Arguments ------------------------------------ 1120 integer, intent(in) :: ikpt,isppol 1121 type(fock_common_type),pointer :: fock 1122 1123 ! ************************************************************************* 1124 1125 !write (std_out,*) ' fock_updateikpt : enter' 1126 1127 ! ====================================================== 1128 ! === Update the data relative to the current states === 1129 ! ====================================================== 1130 !* Copy of the value ikpt in the field ikpt 1131 fock%ikpt=ikpt 1132 !* Copy of the value isppol in the field isppol 1133 fock%isppol=isppol 1134 !* Set all the Fock contributions to the eigenvalues to 0.d0. 1135 fock%eigen_ikpt=zero 1136 !* Set all the Fock contributions to the forces to 0.d0. 1137 if ((fock%optfor).and.(fock%use_ACE==0)) then 1138 fock%forces_ikpt=zero 1139 endif 1140 1141 end subroutine fock_updateikpt
m_fock/fockbz_create [ Functions ]
[ Top ] [ m_fock ] [ Functions ]
NAME
fockbz_create
FUNCTION
Create a fock__BZ_type structure.
INPUTS
NOTES
############################ ### Not fully tested yet ### ############################ The current version is restricted to the case nsym=1, nspinor=1 and mkmem/=0.
SOURCE
338 subroutine fockbz_create(fockbz,mgfft,mpw,mkpt,mkptband,my_nsppol,n4,n5,n6,use_ACE) 339 340 !Arguments ------------------------------------ 341 !scalars 342 integer, intent(in) :: mgfft,mpw,mkpt,mkptband,my_nsppol,n4,n5,n6,use_ACE 343 type(fock_BZ_type) , intent(inout) :: fockbz 344 345 !Local variables------------------------------- 346 347 ! ************************************************************************* 348 349 !write (std_out,*) ' fockbz_create : enter' 350 351 !* Create the array %kptns_bz = the k points in full BZ 352 ABI_MALLOC(fockbz%kptns_bz,(3,mkpt)) 353 fockbz%kptns_bz=zero 354 !* Create the array %jstwfk = how is stored the wavefunction at each k-point 355 !* By default, the table is initialized to 1 (do NOT take advantage of the time-reversal symmetry) 356 ABI_MALLOC(fockbz%istwfk_bz,(mkpt)) 357 fockbz%istwfk_bz=1 358 !* Create the array %wtk_bz = weight assigned to each k point. 359 ABI_MALLOC(fockbz%wtk_bz,(mkpt)) 360 fockbz%wtk_bz=zero 361 !* Create the array %npwarr_bz = number of planewaves in basis at each k-point 362 ! ABI_MALLOC(fockbz%npwarr_bz,(mkpt)) 363 ! fockbz%npwarr_bz=0 364 !* Create the array %kg_bz = reduced planewave coordinates at each k-point 365 ABI_MALLOC(fockbz%kg_bz,(3,mpw*mkpt)) 366 fockbz%kg_bz=0 367 !* Create the array %gbound_bz = boundary of the basis sphere of G vectors at each k-point 368 ABI_MALLOC(fockbz%gbound_bz,(2*mgfft+8,2,mkpt)) 369 fockbz%gbound_bz=0 370 371 !* Create the array %tab_ikpt = indices of k-point ikpt in IBZ which corresponds to each k-point jkpt in full BZ 372 ABI_MALLOC(fockbz%tab_ikpt,(mkpt)) 373 fockbz%tab_ikpt=0 374 !* Create the array %tab_symkpt =indices of symmetry operation to apply to get jkpt in full BZ from ikpt in IBZ 375 ABI_MALLOC(fockbz%tab_symkpt,(mkpt)) 376 fockbz%tab_symkpt=0 377 !* Create the array %tab_ibg = indices of occ(ikpt) in the arrays cprj/occ for each k-point jkpt 378 ABI_MALLOC(fockbz%tab_ibg,(mkpt,my_nsppol)) 379 fockbz%tab_ibg=0 380 !* Create the array %tab_icp = indices of cprj(ikpt) in the arrays cprj/occ for each k-point jkpt 381 ABI_MALLOC(fockbz%tab_icp,(mkpt,my_nsppol)) 382 fockbz%tab_icp=0 383 !* Create the array %tab_icg = indices of cg(ikpt) in the arrays cg for each k-point jkpt 384 ABI_MALLOC(fockbz%tab_icg,(mkpt,my_nsppol)) 385 fockbz%tab_icg=0 386 387 !* Create the array %calc_phase = 1 if a phase factor must be considered (0 otherwise) at each k point 388 ABI_MALLOC(fockbz%calc_phase,(mkpt)) 389 fockbz%calc_phase=0 390 !* Create the array %phase = phase factor the cg array will be multiplied with at each k point 391 ABI_MALLOC(fockbz%phase,(2,mpw*mkpt)) 392 fockbz%phase=zero 393 394 !* Create the array %timerev i= 1 if time reversal symmetry must be used (0 otherwise) at each k point 395 ABI_MALLOC(fockbz%timerev,(mkpt)) 396 fockbz%timerev=0 397 398 !* Create the array %cwaveocc_bz = wavefunctions of each bands at each k point 399 400 if (use_ACE==1) then 401 ABI_MALLOC(fockbz%cgocc,(2,mpw*mkptband,my_nsppol)) 402 fockbz%cgocc=zero 403 else 404 ABI_MALLOC(fockbz%cwaveocc_bz,(2,n4,n5,n6,mkptband,my_nsppol)) 405 fockbz%cwaveocc_bz=zero 406 end if 407 !* Create the array %occ_bz = occupancy of each bands at each k point => will be limited to only the occupied states 408 ABI_MALLOC(fockbz%occ_bz,(mkptband,my_nsppol)) 409 fockbz%occ_bz=zero 410 !* Create the array %nbandocc_bz = nb of bands at each k point 411 ABI_MALLOC(fockbz%nbandocc_bz,(mkpt,my_nsppol)) 412 fockbz%nbandocc_bz=0 413 414 ABI_MALLOC(fockbz%npwarr,(mkpt)) 415 fockbz%npwarr=0 416 417 end subroutine fockbz_create