From c477b848acaf20fcfbbb44fb179e57ec644f248d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 9 Jul 2019 16:17:10 +0200 Subject: [PATCH] basis correction --- GoXC | 2 +- examples/basis.Be.STO-3G | 13 + examples/molecule.H2 | 4 +- input/basis | 34 +- input/dft | 6 +- input/methods | 2 +- input/molecule | 4 +- input/options | 2 +- input/weight | 34 +- src/QuAcK/AO_values_grid.f90 | 101 +++ src/QuAcK/MO_values_grid.f90 | 44 + src/QuAcK/QuAcK.f90 | 44 + src/QuAcK/density.f90 | 53 +- src/{eDFT/dft_grid.f90 => QuAcK/dft_grid.f} | 0 src/QuAcK/ec_srlda.f90 | 46 + src/QuAcK/excitation_density_SOSEX_RI.f90 | 65 ++ src/QuAcK/f_grid.f90 | 44 + src/QuAcK/mu_grid.f90 | 46 + src/QuAcK/quadrature_grid.f90 | 77 ++ src/QuAcK/read_grid.f90 | 47 + src/QuAcK/srlda.f90 | 125 +++ src/eDFT/DIIS_extrapolation.f90 | 54 -- src/eDFT/Makefile | 15 +- src/eDFT/eDFT.f90 | 49 +- src/eDFT/orthogonalization_matrix.f90 | 91 +- src/utils/Makefile | 4 +- src/utils/map_module.f90 | 902 -------------------- src/utils/map_module.mod | Bin 4548 -> 0 bytes 28 files changed, 849 insertions(+), 1059 deletions(-) create mode 100644 examples/basis.Be.STO-3G create mode 100644 src/QuAcK/AO_values_grid.f90 create mode 100644 src/QuAcK/MO_values_grid.f90 rename src/{eDFT/dft_grid.f90 => QuAcK/dft_grid.f} (100%) create mode 100644 src/QuAcK/ec_srlda.f90 create mode 100644 src/QuAcK/excitation_density_SOSEX_RI.f90 create mode 100644 src/QuAcK/f_grid.f90 create mode 100644 src/QuAcK/mu_grid.f90 create mode 100644 src/QuAcK/quadrature_grid.f90 create mode 100644 src/QuAcK/read_grid.f90 create mode 100644 src/QuAcK/srlda.f90 delete mode 100644 src/eDFT/DIIS_extrapolation.f90 delete mode 100644 src/utils/map_module.f90 delete mode 100644 src/utils/map_module.mod diff --git a/GoXC b/GoXC index e91cb5a..52f4e57 100755 --- a/GoXC +++ b/GoXC @@ -13,6 +13,6 @@ then cp examples/basis."$1"."$2" input/basis cp examples/basis."$1"."$2" input/weight ./bin/IntPak - ./bin/xcDFT + ./bin/eDFT fi diff --git a/examples/basis.Be.STO-3G b/examples/basis.Be.STO-3G new file mode 100644 index 0000000..de39473 --- /dev/null +++ b/examples/basis.Be.STO-3G @@ -0,0 +1,13 @@ +1 3 +S 3 1.00 + 0.3016787069D+02 0.1543289673D+00 + 0.5495115306D+01 0.5353281423D+00 + 0.1487192653D+01 0.4446345422D+00 +S 3 1.00 + 0.1314833110D+01 -0.9996722919D-01 + 0.3055389383D+00 0.3995128261D+00 + 0.9937074560D-01 0.7001154689D+00 +P 3 1.00 + 0.1314833110D+01 0.1559162750D+00 + 0.3055389383D+00 0.6076837186D+00 + 0.9937074560D-01 0.3919573931D+00 diff --git a/examples/molecule.H2 b/examples/molecule.H2 index a38d8a2..39795fb 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -1,5 +1,5 @@ -# nAt nEl nEla nElb nCore nRyd - 2 2 1 1 0 0 +# nAt nEla nElb nCore nRyd + 2 1 1 0 0 # Znuc x y z H 0. 0. -0.7 H 0. 0. +0.7 diff --git a/input/basis b/input/basis index b9ca7b5..eb58543 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,29 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 8 1.00 + 17880.0000000 0.0007380 + 2683.0000000 0.0056770 + 611.5000000 0.0288830 + 173.5000000 0.1085400 + 56.6400000 0.2909070 + 20.4200000 0.4483240 + 7.8100000 0.2580260 + 1.6530000 0.0150630 +S 8 1.00 + 17880.0000000 -0.0001720 + 2683.0000000 -0.0013570 + 611.5000000 -0.0067370 + 173.5000000 -0.0276630 + 56.6400000 -0.0762080 + 20.4200000 -0.1752270 + 7.8100000 -0.1070380 + 1.6530000 0.5670500 S 1 1.00 - 0.2976000 1.0000000 + 0.4869000 1.0000000 +P 3 1.00 + 28.3900000 0.0460870 + 6.2700000 0.2401810 + 1.6950000 0.5087440 P 1 1.00 - 1.2750000 1.0000000 + 0.4317000 1.0000000 +D 1 1.00 + 2.2020000 1.0000000 diff --git a/input/dft b/input/dft index ca9be70..90787c0 100644 --- a/input/dft +++ b/input/dft @@ -1,9 +1,9 @@ # exchange rung: Hartree = 0, LDA = 1 (S51), GGA = 2(G96,B88), hybrid = 4, Hartree-Fock = 666 - 1 S51 + 666 HF # correlation rung: Hartree = 0, LDA = 1 (W38,VWN5,C16,LF19), GGA = 2(LYP), hybrid = 4(B3LYP), Hartree-Fock = 666 - 1 VWN5 + 0 HF # quadrature grid SG-n - 3 + 1 # Number of states in ensemble (nEns) 1 # Ensemble weights: wEns(1),...,wEns(nEns-1) diff --git a/input/methods b/input/methods index 8204ed8..341e057 100644 --- a/input/methods +++ b/input/methods @@ -5,7 +5,7 @@ # CCD CCSD CCSD(T) F F F # CIS TDHF ADC - F T F + F F F # GF2 GF3 F F # G0W0 evGW qsGW diff --git a/input/molecule b/input/molecule index c78e87e..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 5 5 0 0 # Znuc x y z - He 0.0 0.0 0.0 + Ne 0.0 0.0 0.0 diff --git a/input/options b/input/options index 86d7547..d3d1e53 100644 --- a/input/options +++ b/input/options @@ -9,6 +9,6 @@ # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 10 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize - 64 0.00001 T 20 F F T F F F F + 64 0.00001 T 20 F F F F F F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/input/weight b/input/weight index b9ca7b5..eb58543 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,29 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 8 1.00 + 17880.0000000 0.0007380 + 2683.0000000 0.0056770 + 611.5000000 0.0288830 + 173.5000000 0.1085400 + 56.6400000 0.2909070 + 20.4200000 0.4483240 + 7.8100000 0.2580260 + 1.6530000 0.0150630 +S 8 1.00 + 17880.0000000 -0.0001720 + 2683.0000000 -0.0013570 + 611.5000000 -0.0067370 + 173.5000000 -0.0276630 + 56.6400000 -0.0762080 + 20.4200000 -0.1752270 + 7.8100000 -0.1070380 + 1.6530000 0.5670500 S 1 1.00 - 0.2976000 1.0000000 + 0.4869000 1.0000000 +P 3 1.00 + 28.3900000 0.0460870 + 6.2700000 0.2401810 + 1.6950000 0.5087440 P 1 1.00 - 1.2750000 1.0000000 + 0.4317000 1.0000000 +D 1 1.00 + 2.2020000 1.0000000 diff --git a/src/QuAcK/AO_values_grid.f90 b/src/QuAcK/AO_values_grid.f90 new file mode 100644 index 0000000..a5a72a6 --- /dev/null +++ b/src/QuAcK/AO_values_grid.f90 @@ -0,0 +1,101 @@ +subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & + nGrid,root,AO,dAO) + +! Compute values of the AOs and their derivatives with respect to the cartesian coordinates + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nShell + double precision,intent(in) :: CenterShell(maxShell,3) + integer,intent(in) :: TotAngMomShell(maxShell) + integer,intent(in) :: KShell(maxShell) + double precision,intent(in) :: DShell(maxShell,maxK) + double precision,intent(in) :: ExpShell(maxShell,maxK) + double precision,intent(in) :: root(3,nGrid) + integer,intent(in) :: nGrid + +! Local variables + + integer :: atot,nShellFunction,a(3) + integer,allocatable :: ShellFunction(:,:) + double precision :: rASq,xA,yA,zA,norm_coeff,prim + + integer :: iSh,iShF,iK,iG,iBas + +! Output variables + + double precision,intent(out) :: AO(nBas,nGrid) + double precision,intent(out) :: dAO(3,nBas,nGrid) + +! Initialization + + iBas = 0 + AO(:,:) = 0d0 + dAO(:,:,:) = 0d0 + +!------------------------------------------------------------------------ +! Loops over shells +!------------------------------------------------------------------------ + do iSh=1,nShell + + atot = TotAngMomShell(iSh) + nShellFunction = (atot*atot + 3*atot + 2)/2 + allocate(ShellFunction(1:nShellFunction,1:3)) + call generate_shell(atot,nShellFunction,ShellFunction) + + do iShF=1,nShellFunction + + iBas = iBas + 1 + a(:) = ShellFunction(iShF,:) + + do iG=1,nGrid + + xA = root(1,iG) - CenterShell(iSh,1) + yA = root(2,iG) - CenterShell(iSh,2) + zA = root(3,iG) - CenterShell(iSh,3) + +! Calculate distance for exponential + + rASq = xA**2 + yA**2 + zA**2 + +!------------------------------------------------------------------------ +! Loops over contraction degrees +!------------------------------------------------------------------------- + do iK=1,KShell(iSh) + +! Calculate the exponential part + + prim = DShell(iSh,iK)*norm_coeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq) + AO(iBas,iG) = AO(iBas,iG) + prim + + prim = -2d0*ExpShell(iSh,iK)*prim + dAO(:,iBas,iG) = dAO(:,iBas,iG) + prim + + enddo + + dAO(1,iBas,iG) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(1,iBas,iG) + if(a(1) > 0) dAO(1,iBas,iG) = dAO(1,iBas,iG) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iBas,iG) + + dAO(2,iBas,iG) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(2,iBas,iG) + if(a(2) > 0) dAO(2,iBas,iG) = dAO(2,iBas,iG) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iBas,iG) + + dAO(3,iBas,iG) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(3,iBas,iG) + if(a(3) > 0) dAO(3,iBas,iG) = dAO(3,iBas,iG) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iBas,iG) + +! Calculate polynmial part + + AO(iBas,iG) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iBas,iG) + + enddo + + enddo + deallocate(ShellFunction) + enddo +!------------------------------------------------------------------------ +! End loops over shells +!------------------------------------------------------------------------ + +end subroutine AO_values_grid diff --git a/src/QuAcK/MO_values_grid.f90 b/src/QuAcK/MO_values_grid.f90 new file mode 100644 index 0000000..a0eeb35 --- /dev/null +++ b/src/QuAcK/MO_values_grid.f90 @@ -0,0 +1,44 @@ +subroutine MO_values_grid(nBas,nGrid,c,AO,dAO,MO,dMO) + +! Compute values of the MOs and their derivatives with respect to the cartesian coordinates + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nGrid + double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: dAO(3,nBas,nGrid) + +! Local variables + + integer :: p,mu,iG + +! Output variables + + double precision,intent(out) :: MO(nBas,nGrid) + double precision,intent(out) :: dMO(3,nBas,nGrid) + +! Initialization + + MO(:,:) = 0d0 + dMO(:,:,:) = 0d0 + + do p=1,nBas + do mu=1,nBas + do iG=1,ngrid + + MO(p,iG) = MO(p,iG) + c(mu,p)*AO(mu,iG) + + dMO(1,p,iG) = dMO(1,p,iG) + c(mu,p)*dAO(1,mu,iG) + dMO(2,p,iG) = dMO(2,p,iG) + c(mu,p)*dAO(2,mu,iG) + dMO(3,p,iG) = dMO(3,p,iG) + c(mu,p)*dAO(3,mu,iG) + + end do + end do + end do + +end subroutine MO_values_grid diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 01e0bae..2db13b5 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -76,6 +76,20 @@ program QuAcK double precision :: dt logical :: doDrift + integer :: SGn + integer :: nRad + integer :: nAng + integer :: nGrid + double precision,allocatable :: root(:,:) + double precision,allocatable :: weight(:) + double precision,allocatable :: AO(:,:) + double precision,allocatable :: dAO(:,:,:) + double precision,allocatable :: MO(:,:) + double precision,allocatable :: dMO(:,:,:) + double precision,allocatable :: rho(:) + double precision,allocatable :: f(:) + double precision,allocatable :: mu(:) + ! Hello World write(*,*) @@ -583,6 +597,36 @@ program QuAcK endif +!------------------------------------------------------------------------ +! Basis set correction +!------------------------------------------------------------------------ + +!------------------------------------------------------------------------ +! Construct quadrature grid +!------------------------------------------------------------------------ + + SGn = 1 + + call read_grid(SGn,nRad,nAng,nGrid) + + allocate(root(ncart,nGrid),weight(nGrid)) + + call quadrature_grid(nRad,nAng,nGrid,root,weight) + +!------------------------------------------------------------------------ +! Calculate AO values at grid points +!------------------------------------------------------------------------ + + allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid),MO(nBas,nGrid),dMO(ncart,nBas,nGrid)) + allocate(rho(nGrid),f(nGrid),mu(nGrid)) + + call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO) + call density(nGrid,nBas,PHF(:,:,1),AO(:,:),rho(:)) + call MO_values_grid(nBas,nGrid,cHF(:,:,1),AO,dAO,MO,dMO) + call f_grid(nBas,nO(1),nGrid,MO,ERI_MO_basis,f) + call mu_grid(nGrid,rho,f,mu) + call ec_srlda(nGrid,weight,rho,mu) + !------------------------------------------------------------------------ ! End of QuAcK !------------------------------------------------------------------------ diff --git a/src/QuAcK/density.f90 b/src/QuAcK/density.f90 index 90a17c1..e6bfd22 100644 --- a/src/QuAcK/density.f90 +++ b/src/QuAcK/density.f90 @@ -1,51 +1,38 @@ -subroutine density(doDrift,nBas,nWalk,P,gAO,dgAO,g,dg) +subroutine density(nGrid,nBas,P,AO,rho) -! Calculate the Green functions +! Calculate one-electron density implicit none + include 'parameters.h' ! Input variables - logical,intent(in) :: doDrift - integer,intent(in) :: nBas,nWalk - double precision,intent(in) :: P(nBas,nBas),gAO(nWalk,2,nBas),dgAO(nWalk,2,3,nBas) + double precision,parameter :: thresh = 1d-15 + + integer,intent(in) :: nGrid + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: AO(nBas,nGrid) ! Local variables - integer :: iW,iEl,ixyz,mu,nu + integer :: iG,mu,nu ! Output variables - double precision,intent(out) :: g(nWalk,2),dg(nWalk,2,3) + double precision,intent(out) :: rho(nGrid) - g = 0d0 - do iW=1,nWalk - do iEl=1,2 - do mu=1,nBas - do nu=1,nBas - g(iW,iEl) = g(iW,iEl) + gAO(iW,iEl,mu)*P(mu,nu)*gAO(iW,iEl,nu) - enddo - enddo - enddo - enddo - - if(doDrift) then - - dg = 0d0 - do iW=1,nWalk - do iEl=1,2 - do ixyz=1,3 - do mu=1,nBas - do nu=1,nBas - dg(iW,iEl,ixyz) = dg(iW,iEl,ixyz) & - + P(mu,nu)*(dgAO(iW,iEl,ixyz,mu)*gAO(iW,iEl,nu) & - + gAO(iW,iEl,mu)*dgAO(iW,iEl,ixyz,nu)) - enddo - enddo - enddo + rho(:) = 0d0 + do iG=1,nGrid + do mu=1,nBas + do nu=1,nBas + rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG) enddo enddo + enddo - endif +! do iG=1,nGrid +! rho(iG) = max(rho(iG),thresh) +! enddo end subroutine density diff --git a/src/eDFT/dft_grid.f90 b/src/QuAcK/dft_grid.f similarity index 100% rename from src/eDFT/dft_grid.f90 rename to src/QuAcK/dft_grid.f diff --git a/src/QuAcK/ec_srlda.f90 b/src/QuAcK/ec_srlda.f90 new file mode 100644 index 0000000..85b653a --- /dev/null +++ b/src/QuAcK/ec_srlda.f90 @@ -0,0 +1,46 @@ +subroutine ec_srlda(nGrid,weight,rho,mu) + +! Compute sr-lda ec + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: mu(nGrid) + +! Local variables + + integer :: iG + double precision :: r + double precision :: rs + double precision :: ecsr + double precision :: ec + double precision,parameter :: thres = 1d-15 + +! Initialization + + ec = 0d0 + + do iG=1,ngrid + + r = max(0d0,rho(iG)) + + if(r > thres) then + + rs = (4d0*pi*r/3d0)**(-1d0/3d0) + + call srlda(rs,mu(iG),ecsr) + + ec = ec + weight(iG)*ecsr*rho(iG) + + end if + + end do + + print*, 'ec = ',ec + +end subroutine ec_srlda diff --git a/src/QuAcK/excitation_density_SOSEX_RI.f90 b/src/QuAcK/excitation_density_SOSEX_RI.f90 new file mode 100644 index 0000000..8f3f6b5 --- /dev/null +++ b/src/QuAcK/excitation_density_SOSEX_RI.f90 @@ -0,0 +1,65 @@ +subroutine excitation_density_SOSEX_RI(nBas,nC,nO,nR,nS,c,G,XpY,rho) + +! Compute excitation densities + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nR,nS + double precision,intent(in) :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS) + +! Local variables + + double precision,allocatable :: scr(:,:,:) + integer :: mu,nu,la,si,ia,jb,x,y,j,b + +! Output variables + + double precision,intent(out) :: rho(nBas,nBas,nS) + +! Memory allocation + allocate(scr(nBas,nBas,nS)) + + rho(:,:,:) = 0d0 + do nu=1,nBas + do si=1,nBas + do ia=1,nS + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b) + enddo + enddo + enddo + enddo + enddo + + scr(:,:,:) = 0d0 + do mu=1,nBas + do la=1,nBas + do ia=1,nS + do nu=1,nBas + do si=1,nBas + scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia) + enddo + enddo + enddo + enddo + enddo + + rho(:,:,:) = 0d0 + do ia=1,nS + do x=nC+1,nBas-nR + do y=nC+1,nBas-nR + do mu=1,nBas + do la=1,nBas + rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y) + enddo + enddo + enddo + enddo + enddo + +end subroutine excitation_density_SOSEX_RI diff --git a/src/QuAcK/f_grid.f90 b/src/QuAcK/f_grid.f90 new file mode 100644 index 0000000..63bd3f8 --- /dev/null +++ b/src/QuAcK/f_grid.f90 @@ -0,0 +1,44 @@ +subroutine f_grid(nBas,nO,nGrid,MO,ERI,f) + +! Compute f + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: nGrid + double precision,intent(in) :: MO(nBas,nGrid) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: p,q + integer :: i,j + integer :: iG + +! Output variables + + double precision,intent(out) :: f(nGrid) + +! Initialization + + f(:) = 0d0 + + do p=1,nBas + do q=1,nBas + do i=1,nO + do j=1,nO + do iG=1,ngrid + + f(iG) = f(iG) + MO(p,iG)*MO(q,iG)*ERI(p,q,i,j)*MO(i,iG)*MO(j,iG) + + end do + end do + end do + end do + end do + +end subroutine f_grid diff --git a/src/QuAcK/mu_grid.f90 b/src/QuAcK/mu_grid.f90 new file mode 100644 index 0000000..23a6ff2 --- /dev/null +++ b/src/QuAcK/mu_grid.f90 @@ -0,0 +1,46 @@ +subroutine mu_grid(nGrid,rho,f,mu) + +! Compute mu + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nGrid + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: f(nGrid) + +! Local variables + + integer :: iG + double precision,parameter :: thres = 1d-10 + double precision :: n2 + +! Output variables + + double precision,intent(out) :: mu(nGrid) + +! Initialization + + mu(:) = 0d0 + + do iG=1,ngrid + + n2 = rho(iG)**2 + + if(abs(n2) > thres) then + + mu(iG) = f(iG)/n2 + + else + + mu(iG) = 1d0/thres + + end if + + end do + + mu(:) = 0.5d0*sqrt(pi)*mu(:) + +end subroutine mu_grid diff --git a/src/QuAcK/quadrature_grid.f90 b/src/QuAcK/quadrature_grid.f90 new file mode 100644 index 0000000..420e80a --- /dev/null +++ b/src/QuAcK/quadrature_grid.f90 @@ -0,0 +1,77 @@ +subroutine quadrature_grid(nRad,nAng,nGrid,root,weight) + +! Build roots and weights of quadrature grid + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nRad,nAng,nGrid + +! Local variables + + integer :: i,j,k + double precision :: scale + double precision,allocatable :: Radius(:) + double precision,allocatable :: RadWeight(:) + double precision,allocatable :: XYZ(:,:) + double precision,allocatable :: XYZWeight(:) + +! Output variables + + double precision,intent(out) :: root(3,nGrid) + double precision,intent(out) :: weight(nGrid) + +! Memory allocation + + allocate(Radius(nRad),RadWeight(nRad),XYZ(3,nAng),XYZWeight(nAng)) + +! Findthe radial grid + + scale = 1d0 + call EulMac(Radius,RadWeight,nRad,scale) + + write(*,20) + write(*,30) + write(*,20) + do i = 1,nRad + write(*,40) i,Radius(i),RadWeight(i) + end do + write(*,20) + write(*,*) + +! Find the angular grid + + call Lebdev(XYZ,XYZWeight,nAng) + + write(*,20) + write(*,50) + write(*,20) + do j = 1,nAng + write(*,60) j,(XYZ(k,j),k=1,3), XYZWeight(j) + end do + write(*,20) + +! Form the roots and weights + + k = 0 + do i=1,nRad + do j=1,nAng + k = k + 1 + root(:,k) = Radius(i)*XYZ(:,j) + weight(k) = RadWeight(i)*XYZWeight(j) + enddo + enddo + +! Compute values of the basis functions (and the its gradient if required) at each grid point + +20 format(T2,58('-')) +30 format(T20,'Radial Quadrature',/, & + T6,'I',T26,'Radius',T50,'Weight') +40 format(T3,I4,T18,F17.10,T35,F25.10) +50 format(T20,'Angular Quadrature',/, & + T6,'I',T19,'X',T29,'Y',T39,'Z',T54,'Weight') +60 format(T3,I4,T13,3F10.5,T50,F10.5) + +end subroutine quadrature_grid diff --git a/src/QuAcK/read_grid.f90 b/src/QuAcK/read_grid.f90 new file mode 100644 index 0000000..5cf4391 --- /dev/null +++ b/src/QuAcK/read_grid.f90 @@ -0,0 +1,47 @@ +subroutine read_grid(SGn,nRad,nAng,nGrid) + +! Read grid type + + implicit none + +! Input variables + + integer,intent(in) :: SGn + +! Output variables + + integer,intent(out) :: nRad + integer,intent(out) :: nAng + integer,intent(out) :: nGrid + + write(*,*)'----------------------------------------------------------' + write(*,'(A22,I1)')' Quadrature grid: SG-',SGn + write(*,*)'----------------------------------------------------------' + + select case (SGn) + + case(0) + nRad = 23 + nAng = 170 + + case(1) + nRad = 50 + nAng = 194 + + case(2) + nRad = 75 + nAng = 302 + + case(3) + nRad = 99 + nAng = 590 + + case default + call print_warning('!!! Quadrature grid not available !!!') + stop + + end select + + nGrid = nRad*nAng + +end subroutine read_grid diff --git a/src/QuAcK/srlda.f90 b/src/QuAcK/srlda.f90 new file mode 100644 index 0000000..aa8609e --- /dev/null +++ b/src/QuAcK/srlda.f90 @@ -0,0 +1,125 @@ +subroutine srlda(rs,mu,ecsr) + +! Correlation energy of a spin unpolarized uniform electron gas +! with short-range interaction erfc(mu*r)/r +! See Zecca et al. PRB 70, 205127 (2004) + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: rs + double precision,intent(in) :: mu + +! Local variables + + double precision :: ec + double precision :: cf + double precision :: b1 + double precision :: b2 + double precision :: b3 + double precision :: b4 + double precision :: a0 + double precision :: bb + double precision :: m1 + +! Ouput variables + + double precision,intent(out) :: ecsr + +! Compute PW LDA correlation energy + + call ecPW(rs,0d0,ec) + +! Define various stuff + + cf = (9d0*pi/4d0)**(1d0/3d0) + bb = 1.27329d0 + m1 = 0.0357866d0 + a0 = ec + b3 = bb*rs**(7d0/2d0) + b2 = -3d0/2d0/pi/cf*rs/a0 + b1 = (b3-1d0/sqrt(3d0*pi)*rs**(3d0/2d0)/a0)/b2 + b4 = -a0*b1*rs**3/m1 + +! Compute short-range correlation energy + + ecsr = a0*(1d0 + b1*mu)/(1d0 + b1*mu+b2*mu**2 + b3*mu**3 + b4*mu**4) + +end subroutine srlda + +!========================================================================================== + +subroutine ecPW(x,y,ec) + +! Correlation energy of the 3D electron gas of density rs and spin polarization z +! Perdew & Wang, PRB 45, 13244 (1992) + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: x + double precision,intent(in) :: y + +! Local variables + + double precision :: f02 + double precision :: ff + double precision :: aaa + double precision :: G + double precision :: ec0 + double precision :: ec1 + double precision :: alfac + +! Ouput variables + + double precision,intent(out) :: ec + + f02 = 4d0/(9d0*(2d0**(1d0/3d0) - 1d0)) + + ff = ((1d0+y)**(4d0/3d0) + (1d0-y)**(4d0/3d0)-2d0)/(2d0**(4d0/3d0) - 2d0) + + aaa = (1d0 - log(2d0))/pi**2 + + call GPW(x,aaa,0.21370d0,7.5957d0,3.5876d0,1.6382d0,0.49294d0,G) + ec0 = G + + aaa=aaa/2d0 + call GPW(x,aaa,0.20548d0,14.1189d0,6.1977d0,3.3662d0,0.62517d0,G) + ec1 = G + + call GPW(x,0.016887d0,0.11125d0,10.357d0,3.6231d0,0.88026d0,0.49671d0,G) + alfac = -G + + ec = ec0 + alfac*ff/f02*(1d0 - y**4) + (ec1 - ec0)*ff*y**4 + +end subroutine ecPW + +!========================================================================================== + +subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G) + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: Ac + double precision,intent(in) :: alfa1 + double precision,intent(in) :: beta1 + double precision,intent(in) :: beta2 + double precision,intent(in) :: beta3 + double precision,intent(in) :: beta4 + double precision,intent(in) :: x + +! Ouput variables + + double precision,intent(out) :: G + + G = -2d0*Ac*(1d0 + alfa1*x)*log(1d0 & + + 1d0/(2d0*Ac*(beta1*sqrt(x) + beta2*x + beta3*x*sqrt(x) + beta4*x**2))) + +end subroutine GPW diff --git a/src/eDFT/DIIS_extrapolation.f90 b/src/eDFT/DIIS_extrapolation.f90 deleted file mode 100644 index e53db06..0000000 --- a/src/eDFT/DIIS_extrapolation.f90 +++ /dev/null @@ -1,54 +0,0 @@ -subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout) - -! Perform DIIS extrapolation - - implicit none - - include 'parameters.h' - -! Input variables - - integer,intent(in) :: n_err,n_e - double precision,intent(in) :: error_in(n_err),error(n_err,n_diis),e(n_e,n_diis) - -! Local variables - - double precision,allocatable :: A(:,:),b(:),w(:) - -! Output variables - - double precision,intent(out) :: rcond - integer,intent(inout) :: n_diis - double precision,intent(inout):: e_inout(n_e) - -! Memory allocaiton - - allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1)) - -! Update DIIS "history" - - call prepend(n_err,n_diis,error,error_in) - call prepend(n_e,n_diis,e,e_inout) - -! Build A matrix - - A(1:n_diis,1:n_diis) = matmul(transpose(error),error) - - A(1:n_diis,n_diis+1) = -1d0 - A(n_diis+1,1:n_diis) = -1d0 - A(n_diis+1,n_diis+1) = +0d0 - -! Build x matrix - - b(1:n_diis) = +0d0 - b(n_diis+1) = -1d0 - -! Solve linear system - - call linear_solve(n_diis+1,A,b,w,rcond) - -! Extrapolate - - e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis))) - -end subroutine DIIS_extrapolation diff --git a/src/eDFT/Makefile b/src/eDFT/Makefile index a72dc51..2c3316c 100644 --- a/src/eDFT/Makefile +++ b/src/eDFT/Makefile @@ -1,12 +1,17 @@ IDIR =../../include BDIR =../../bin ODIR = obj +OODIR = ../IntPak/obj SDIR =. -FC = gfortran -I$(IDIR) +FC = gfortran -I$(IDIR) -Wall -g -Wno-unused -Wno-unused-dummy-argument ifeq ($(DEBUG),1) FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant else -FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O2 +FFLAGS = -O3 +endif + +ifeq ($(PROFIL),1) + FC += -p -fno-inline endif LIBS = ~/Dropbox/quack/lib/*.a @@ -18,7 +23,7 @@ SRC = $(wildcard *.f) OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC)) -$(ODIR)/%.o: %.f90 +$(ODIR)/%.o: %.f90 $(FC) -c -o $@ $< $(FFLAGS) $(ODIR)/%.o: %.f @@ -29,7 +34,9 @@ $(BDIR)/eDFT: $(OBJ) debug: DEBUG=1 make $(BDIR)/eDFT -#DEBUG=1 make clean $(BDIR)/eDFT + +profil: + PROFIL=1 make $(BDIR)/eDFT clean: rm -f $(ODIR)/*.o $(BDIR)/eDFT $(BDIR)/debug diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 9db226b..9ebae91 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -2,9 +2,11 @@ program eDFT ! exchange-correlation density-functional theory calculations + implicit none include 'parameters.h' - integer :: nNuc,nBas,nEl(nspin),nO(nspin),nV(nspin) + integer :: nNuc,nBas + integer :: nEl(nspin),nC(nspin),nO(nspin),nV(nspin),nR(nspin) double precision :: ENuc,EKS double precision,allocatable :: ZNuc(:),rNuc(:,:) @@ -35,7 +37,9 @@ program eDFT integer :: maxSCF,max_diis double precision :: thresh - logical :: DIIS,guess_type,ortho_type + logical :: DIIS + integer :: guess_type + integer :: ortho_type ! Hello World @@ -50,12 +54,14 @@ program eDFT !------------------------------------------------------------------------ ! Read number of atoms, number of electrons of the system +! nC = number of core orbitals ! nO = number of occupied orbitals ! nV = number of virtual orbitals (see below) +! nR = number of Rydberg orbitals ! nBas = number of basis functions (see below) ! = nO + nV - call read_molecule(nNuc,nEl,nO) + call read_molecule(nNuc,nEl(:),nO(:),nC(:),nR(:)) allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) ! Read geometry @@ -71,23 +77,6 @@ program eDFT call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) -!------------------------------------------------------------------------ -! Read one- and two-electron integrals -!------------------------------------------------------------------------ - -! Memory allocation for one- and two-electron integrals - - allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), & - ERI(nBas,nBas,nBas,nBas)) - -! Read integrals - - call read_integrals(nBas,S,T,V,Hc,ERI) - -! Orthogonalization X = S^(-1/2) - - call orthogonalization_matrix(nBas,S,X) - !------------------------------------------------------------------------ ! DFT options !------------------------------------------------------------------------ @@ -95,9 +84,24 @@ program eDFT ! Allocate ensemble weights allocate(wEns(maxEns)) - call read_options(x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type) +!------------------------------------------------------------------------ +! Read one- and two-electron integrals +!------------------------------------------------------------------------ + +! Memory allocation for one- and two-electron integrals + + allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI(nBas,nBas,nBas,nBas)) + +! Read integrals + + call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI) + +! Orthogonalization X = S^(-1/2) + + call orthogonalization_matrix(ortho_type,nBas,S,X) + !------------------------------------------------------------------------ ! Construct quadrature grid !------------------------------------------------------------------------ @@ -111,8 +115,7 @@ program eDFT !------------------------------------------------------------------------ allocate(AO(nBas,nGrid),dAO(ncart,nBas,nGrid)) - call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - nGrid,root,AO,dAO) + call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell,nGrid,root,AO,dAO) !------------------------------------------------------------------------ ! Compute KS energy diff --git a/src/eDFT/orthogonalization_matrix.f90 b/src/eDFT/orthogonalization_matrix.f90 index ed837d4..15ea4ac 100644 --- a/src/eDFT/orthogonalization_matrix.f90 +++ b/src/eDFT/orthogonalization_matrix.f90 @@ -1,12 +1,12 @@ -subroutine orthogonalization_matrix(nBas,S,X) +subroutine orthogonalization_matrix(ortho_type,nBas,S,X) -! Compute the orthogonalization matrix X = S^(-1/2) +! Compute the orthogonalization matrix X implicit none ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas,ortho_type double precision,intent(in) :: S(nBas,nBas) ! Local variables @@ -23,30 +23,87 @@ subroutine orthogonalization_matrix(nBas,S,X) debug = .false. +! Type of orthogonalization ortho_type +! +! 1 = Lowdin +! 2 = Canonical +! 3 = SVD +! + allocate(Uvec(nBas,nBas),Uval(nBas)) - write(*,*) - write(*,*) ' *** Lowdin orthogonalization X = S^(-1/2) *** ' - write(*,*) + if(ortho_type == 1) then - Uvec = S - call diagonalize_matrix(nBas,Uvec,Uval) + write(*,*) + write(*,*) ' Lowdin orthogonalization' + write(*,*) - do i=1,nBas + Uvec = S + call diagonalize_matrix(nBas,Uvec,Uval) - if(Uval(i) > thresh) then + do i=1,nBas - Uval(i) = 1d0/sqrt(Uval(i)) + if(Uval(i) > thresh) then - else + Uval(i) = 1d0/sqrt(Uval(i)) - write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization' + else - endif + write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization' - enddo - - call ADAt(nBas,Uvec,Uval,X) + endif + + enddo + + call ADAt(nBas,Uvec,Uval,X) + + elseif(ortho_type == 2) then + + write(*,*) + write(*,*) 'Canonical orthogonalization' + write(*,*) + + Uvec = S + call diagonalize_matrix(nBas,Uvec,Uval) + + do i=1,nBas + + if(Uval(i) > thresh) then + + Uval(i) = 1d0/sqrt(Uval(i)) + + else + + write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' + + endif + + enddo + + call AD(nBas,Uvec,Uval) + X = Uvec + + elseif(ortho_type == 3) then + + write(*,*) + write(*,*) ' SVD-based orthogonalization NYI' + write(*,*) + +! Uvec = S +! call diagonalize_matrix(nBas,Uvec,Uval) + +! do i=1,nBas +! if(Uval(i) > thresh) then +! Uval(i) = 1d0/sqrt(Uval(i)) +! else +! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization' +! endif +! enddo +! +! call AD(nBas,Uvec,Uval) +! X = Uvec + + endif ! Print results diff --git a/src/utils/Makefile b/src/utils/Makefile index 09a8268..d30e9a9 100644 --- a/src/utils/Makefile +++ b/src/utils/Makefile @@ -5,7 +5,7 @@ ODIR = obj OODIR = ../IntPak/obj SDIR =. FC = gfortran -I$(IDIR) -AR = ar +AR = libtool ifeq ($(DEBUG),1) FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant else @@ -28,7 +28,7 @@ $(ODIR)/%.o: %.f $(FC) -c -o $@ $< $(FFLAGS) $(LDIR)/utils.a: $(OBJ) - $(AR) -cru $@ $^ + $(AR) -static -o $@ $^ debug: DEBUG=1 make $(LDIR)/utils.a diff --git a/src/utils/map_module.f90 b/src/utils/map_module.f90 deleted file mode 100644 index 98e7347..0000000 --- a/src/utils/map_module.f90 +++ /dev/null @@ -1,902 +0,0 @@ -module map_module - -! A map is an array of maps (cache_maps) -! A cache map is an array of keys and values sorted by keys -! A cache map has its own OpenMP lock -! To access a (key,value) pair in the map, the -! index of the cache_map in the map array is obtained -! by removing the first 15 bits of the key. -! The key in the cache_map is composed of the first -! 15 bits of the key. Therefore, it can be stored -! as integer*2 and is found by applying the map_mask -! to the initial key. The element are found in the -! cache_map using a binary search -! -! When using the map_update subroutine to build the map, -! the map_merge subroutine -! should be called before getting data from the map. - - use omp_lib - - integer, parameter :: integral_kind = 8 - - integer, parameter :: cache_key_kind = 2 - integer, parameter :: cache_map_size_kind = 4 - - integer, parameter :: key_kind = 8 - integer, parameter :: map_size_kind = 8 - - integer, parameter :: map_shift = 15 - integer*8, parameter :: map_mask = ibset(0_8,15)-1_8 - - type cache_map_type - real(integral_kind), pointer :: value(:) - integer(cache_key_kind), pointer :: key(:) - logical :: sorted - integer(cache_map_size_kind) :: map_size - integer(cache_map_size_kind) :: n_elements - integer(omp_lock_kind) :: lock - end type cache_map_type - - type map_type - type(cache_map_type), allocatable :: map(:) - real(integral_kind), pointer :: consolidated_value(:) - integer(cache_key_kind), pointer :: consolidated_key(:) - integer*8, pointer :: consolidated_idx(:) - logical :: sorted - logical :: consolidated - integer(map_size_kind) :: map_size - integer(map_size_kind) :: n_elements - integer(omp_lock_kind) :: lock - end type map_type - -end module map_module - - -double precision function map_mb(map) - use map_module - use omp_lib - implicit none - type (map_type), intent(in) :: map - integer(map_size_kind) :: i - - map_mb = dble(8+map_size_kind+map_size_kind+omp_lock_kind+4) - do i=0,map%map_size - map_mb = map_mb + dble(map%map(i)%map_size*(cache_key_kind+integral_kind) +& - 8+8+4+cache_map_size_kind+cache_map_size_kind+omp_lock_kind) - enddo - map_mb = map_mb / (1024.d0*1024.d0) -end - -subroutine cache_map_init(map,sze) - use map_module - implicit none - type (cache_map_type), intent(inout) :: map - integer(cache_map_size_kind) :: sze - call omp_set_lock(map%lock) - map%n_elements = 0_8 - map%map_size = 0_8 - map%sorted = .True. - NULLIFY(map%value, map%key) - call cache_map_reallocate(map,sze) - call omp_unset_lock(map%lock) -end - -subroutine map_init(map,keymax) - use map_module - implicit none - integer*8, intent(in) :: keymax - type (map_type), intent(inout) :: map - integer(map_size_kind) :: i - integer(cache_map_size_kind) :: sze - integer :: err - - call omp_init_lock(map%lock) - call omp_set_lock(map%lock) - - map%n_elements = 0_8 - map%map_size = shiftr(keymax,map_shift) - map%consolidated = .False. - - allocate(map%map(0_8:map%map_size),stat=err) - if (err /= 0) then - print *, 'Unable to allocate map' - stop 5 - endif - sze = 2 - do i=0_8,map%map_size - call omp_init_lock(map%map(i)%lock) - enddo - !$OMP PARALLEL DEFAULT(NONE) SHARED(map,sze) PRIVATE(i) - !$OMP DO SCHEDULE(STATIC,512) - do i=0_8,map%map_size - call cache_map_init(map%map(i),sze) - enddo - !$OMP ENDDO - !$OMP END PARALLEL - map%sorted = .True. - - call omp_unset_lock(map%lock) - -end - -subroutine cache_map_reallocate(map,sze) - use map_module - implicit none - integer(cache_map_size_kind), intent(in) :: sze - type (cache_map_type), intent(inout) :: map - - integer(cache_key_kind), pointer :: key_new(:) - real(integral_kind), pointer :: value_new(:) - integer(map_size_kind) :: i - integer :: err - !DIR$ ATTRIBUTES ALIGN : 64 :: key_new, value_new - - if (sze < map%n_elements) then - print *, 'Unable to resize map : map too large' - stop 3 - endif - - ! Resize keys - allocate( key_new(sze), stat=err ) - if (err /= 0) then - print *, 'Unable to allocate map', sze - stop 1 - endif - if (associated(map%key)) then - do i=1_8,min(size(map%key),map%n_elements) - key_new(i) = map%key(i) - enddo - deallocate(map%key) - endif - - ! Resize values - allocate( value_new(sze), stat=err ) - if (err /= 0) then - print *, 'Unable to allocate map', sze - stop 2 - endif - if (associated(map%value)) then - do i=1_8,min(size(map%key),map%n_elements) - value_new(i) = map%value(i) - enddo - deallocate(map%value) - endif - - ! Set new pointers - map%key => key_new - map%value => value_new - map%map_size = sze - -end - - -subroutine cache_map_deinit(map) - use map_module - implicit none - type (cache_map_type), intent(inout) :: map - - integer :: err - - if (associated( map % value )) then - deallocate( map % value, stat=err ) - if (err /= 0) then - print *, 'Unable to deallocate map' - stop 2 - endif - NULLIFY(map%value) - endif - - if (associated( map % key )) then - deallocate( map % key, stat=err ) - if (err /= 0) then - print *, 'Unable to deallocate map' - stop 4 - endif - NULLIFY(map%key) - endif - - map%n_elements = 0_8 - map%map_size = 0_8 - call omp_destroy_lock(map%lock) -end - -subroutine map_deinit(map) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer :: err - integer(map_size_kind) :: i - - if (allocated( map % map )) then - do i=0_8, map%map_size - call cache_map_deinit(map%map(i)) - enddo - deallocate( map % map, stat=err ) - if (err /= 0) then - print *, 'Unable to deallocate map' - stop 6 - endif - endif - - map%n_elements = 0_8 - map%map_size = 0_8 - call omp_destroy_lock(map%lock) -end - -subroutine cache_map_sort(map) - use map_module - implicit none - type (cache_map_type), intent(inout) :: map - integer(cache_map_size_kind), allocatable :: iorder(:) - integer(cache_map_size_kind) :: i - !DIR$ ATTRIBUTES ALIGN : 64 :: iorder - - if (.not.map%sorted) then - allocate(iorder(map%n_elements)) - do i=1,map%n_elements - iorder(i) = i - enddo - if (cache_key_kind == 2) then - call i2radix_sort(map%key,iorder,map%n_elements,-1) - else if (cache_key_kind == 4) then - call iradix_sort(map%key,iorder,map%n_elements,-1) - else if (cache_key_kind == 8) then - call i8radix_sort(map%key,iorder,map%n_elements,-1) - endif - if (integral_kind == 4) then - call set_order(map%value,iorder,map%n_elements) - else if (integral_kind == 8) then - call dset_order(map%value,iorder,map%n_elements) - endif - deallocate(iorder) - map%sorted = .True. - endif - -end - -subroutine map_sort(map) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer(map_size_kind) :: i - - if (.not.map%sorted) then - !$OMP PARALLEL DO SCHEDULE(static,1024) DEFAULT(SHARED) PRIVATE(i) - do i=0_8,map%map_size - call omp_set_lock(map%map(i)%lock) - call cache_map_sort(map%map(i)) - call omp_unset_lock(map%map(i)%lock) - enddo - !$OMP END PARALLEL DO - map%sorted = .True. - endif - -end - -subroutine cache_map_merge(map) - use map_module - implicit none - type (cache_map_type), intent(inout) :: map - integer(cache_key_kind) :: prev_key - integer(cache_map_size_kind) :: i, j - - call cache_map_sort(map) - prev_key = -1_8 - j=0 - do i=1,map%n_elements - if (map%key(i) /= prev_key) then - j = j+1 - map%value(j) = map%value(i) - map%key(j) = map%key(i) - prev_key = map%key(i) - else - map%value(j) = map%value(j)+map%value(i) - endif - enddo - map%n_elements = j - -end - -subroutine cache_map_unique(map) - use map_module - implicit none - type (cache_map_type), intent(inout) :: map - integer(cache_key_kind) :: prev_key - integer(cache_map_size_kind) :: i, j - - call cache_map_sort(map) - prev_key = -1_8 - j=0 - do i=1,map%n_elements - if (map%key(i) /= prev_key) then - j = j+1 - map%value(j) = map%value(i) - map%key(j) = map%key(i) - prev_key = map%key(i) - endif - enddo - map%n_elements = j - -end - -subroutine cache_map_shrink(map,thr) - use map_module - implicit none - type (cache_map_type), intent(inout) :: map - real(integral_kind) , intent(in) :: thr - integer(cache_map_size_kind) :: i,j - - j=0 - do i=1,map%n_elements - if (abs(map%value(i)) > thr) then - j = j+1 - map%value(j) = map%value(i) - map%key(j) = map%key(i) - endif - enddo - map%n_elements = j - -end - -subroutine map_unique(map) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer(map_size_kind) :: i - integer(map_size_kind) :: icount - - icount = 0_8 - !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)& - !$OMP REDUCTION(+:icount) - do i=0_8,map%map_size - call omp_set_lock(map%map(i)%lock) - call cache_map_unique(map%map(i)) - call omp_unset_lock(map%map(i)%lock) - icount = icount + map%map(i)%n_elements - enddo - !$OMP END PARALLEL DO - map%n_elements = icount - -end - -subroutine map_merge(map) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer(map_size_kind) :: i - integer(map_size_kind) :: icount - - icount = 0_8 - !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)& - !$OMP REDUCTION(+:icount) - do i=0_8,map%map_size - call omp_set_lock(map%map(i)%lock) - call cache_map_merge(map%map(i)) - call omp_unset_lock(map%map(i)%lock) - icount = icount + map%map(i)%n_elements - enddo - !$OMP END PARALLEL DO - map%n_elements = icount - -end - -subroutine map_shrink(map,thr) - use map_module - implicit none - type (map_type), intent(inout) :: map - real(integral_kind), intent(in) :: thr - integer(map_size_kind) :: i - integer(map_size_kind) :: icount - - icount = 0_8 - !$OMP PARALLEL DO SCHEDULE(dynamic,1000) DEFAULT(SHARED) PRIVATE(i)& - !$OMP REDUCTION(+:icount) - do i=0_8,map%map_size - call omp_set_lock(map%map(i)%lock) - call cache_map_shrink(map%map(i),thr) - call omp_unset_lock(map%map(i)%lock) - icount = icount + map%map(i)%n_elements - enddo - !$OMP END PARALLEL DO - map%n_elements = icount - -end - -subroutine map_update(map, key, value, sze, thr) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer, intent(in) :: sze - integer(key_kind), intent(inout) :: key(sze) - real(integral_kind), intent(inout) :: value(sze) - real(integral_kind), intent(in) :: thr - - integer :: i - integer(map_size_kind) :: idx_cache, idx_cache_new - integer(cache_map_size_kind) :: idx - integer :: sze2 - integer(cache_key_kind) :: cache_key - integer(map_size_kind) :: n_elements_temp - type (cache_map_type) :: local_map - logical :: map_sorted - - sze2 = sze - map_sorted = .True. - - n_elements_temp = 0_8 - n_elements_temp = n_elements_temp + 1_8 - do while (sze2>0) - i=1 - do while (i<=sze) - if (key(i) /= 0_8) then - idx_cache = shiftr(key(i),map_shift) - if (omp_test_lock(map%map(idx_cache)%lock)) then - local_map%key => map%map(idx_cache)%key - local_map%value => map%map(idx_cache)%value - local_map%sorted = map%map(idx_cache)%sorted - local_map%map_size = map%map(idx_cache)%map_size - local_map%n_elements = map%map(idx_cache)%n_elements - do - !DIR$ FORCEINLINE - call search_key_big_interval(key(i),local_map%key, local_map%n_elements, idx, 1, local_map%n_elements) - if (idx > 0_8) then - local_map%value(idx) = local_map%value(idx) + value(i) - else - ! Assert that the map has a proper size - if (local_map%n_elements == local_map%map_size) then - call cache_map_merge(local_map) - call cache_map_reallocate(local_map, local_map%n_elements + local_map%n_elements) - call cache_map_shrink(local_map,thr) - endif - cache_key = int(iand(key(i),map_mask),2) - local_map%n_elements = local_map%n_elements + 1 - local_map%value(local_map%n_elements) = value(i) - local_map%key(local_map%n_elements) = cache_key - local_map%sorted = .False. - n_elements_temp = n_elements_temp + 1_8 - endif ! idx > 0 - key(i) = 0_8 - i = i+1 - sze2 = sze2-1 - if (i>sze) then - i=1 - endif - if ( (shiftr(key(i),map_shift) /= idx_cache).or.(key(i)==0_8)) then - exit - endif - enddo - map%map(idx_cache)%key => local_map%key - map%map(idx_cache)%value => local_map%value - map%map(idx_cache)%sorted = local_map%sorted - map%map(idx_cache)%n_elements = local_map%n_elements - map%map(idx_cache)%map_size = local_map%map_size - map_sorted = map_sorted .and. local_map%sorted - call omp_unset_lock(map%map(idx_cache)%lock) - endif ! omp_test_lock - else - i=i+1 - endif ! key = 0 - enddo ! i -enddo ! sze2 > 0 -call omp_set_lock(map%lock) -map%n_elements = map%n_elements + n_elements_temp -map%sorted = map%sorted .and. map_sorted -call omp_unset_lock(map%lock) - -end - -subroutine map_append(map, key, value, sze) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer, intent(in) :: sze - integer(key_kind), intent(inout) :: key(sze) - real(integral_kind), intent(inout) :: value(sze) - - integer :: i - integer(cache_map_size_kind) :: n_elements - integer(map_size_kind) :: idx_cache - integer(cache_key_kind) :: cache_key - - do i=1,sze - idx_cache = shiftr(key(i),map_shift) - call omp_set_lock(map%map(idx_cache)%lock) - n_elements = map%map(idx_cache)%n_elements + 1 - ! Assert that the map has a proper size - if (n_elements == map%map(idx_cache)%map_size) then - call cache_map_reallocate(map%map(idx_cache), n_elements+ shiftr(n_elements,1)) - endif - cache_key = int(iand(key(i),map_mask),2) - map%map(idx_cache)%value(n_elements) = value(i) - map%map(idx_cache)%key(n_elements) = cache_key - map%map(idx_cache)%n_elements = n_elements - if (map%map(idx_cache)%sorted.and.n_elements > 1) then - map%map(idx_cache)%sorted = (map%map(idx_cache)%key(n_elements-1) <= cache_key) - map%sorted = map%sorted .and. map%map(idx_cache)%sorted - endif - call omp_unset_lock(map%map(idx_cache)%lock) - enddo - call omp_set_lock(map%lock) - map%n_elements = map%n_elements + sze - call omp_unset_lock(map%lock) - -end - -subroutine map_get(map, key, value) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer(key_kind), intent(in) :: key - real(integral_kind), intent(out) :: value - integer(map_size_kind) :: idx_cache - integer(cache_map_size_kind) :: idx - - ! index in tha pointers array - idx_cache = shiftr(key,map_shift) - !DIR$ FORCEINLINE - call cache_map_get_interval(map%map(idx_cache), key, value, 1, map%map(idx_cache)%n_elements,idx) -end - -subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx) - use map_module - implicit none - type (cache_map_type), intent(inout) :: map - integer(key_kind), intent(in) :: key - integer(cache_map_size_kind), intent(in) :: ibegin, iend - real(integral_kind), intent(out) :: value - integer(cache_map_size_kind), intent(inout) :: idx - double precision, pointer :: v(:) - integer :: i - - call search_key_big_interval(key,map%key, map%n_elements, idx, ibegin, iend) - if (idx > 0) then - value = map%value(idx) - else - value = 0._integral_kind - endif -! call search_key_value_big_interval(key, value, map%key, map%value, map%n_elements, idx, ibegin, iend) -end - - -subroutine map_get_many(map, key, value, sze) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer, intent(in) :: sze - integer(key_kind), intent(in) :: key(sze) - real(integral_kind), intent(out) :: value(sze) - integer :: i - integer(map_size_kind) :: idx_cache - integer(cache_map_size_kind) :: ibegin, iend - integer(cache_map_size_kind), allocatable :: idx(:) - !DIR$ ATTRIBUTES ALIGN : 64 :: idx - - allocate(idx(sze)) - do i=1,sze - idx_cache = shiftr(key(i),map_shift) - iend = map%map(idx_cache)%n_elements - !DIR$ FORCEINLINE - call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend) - enddo - do i=1,sze - idx_cache = shiftr(key(i),map_shift) - if (idx(i) > 0) then - value(i) = map%map(idx_cache)%value(idx(i)) - else - value(i) = 0. - endif - enddo - deallocate(idx) -end - -subroutine map_exists_many(map, key, sze) - use map_module - implicit none - type (map_type), intent(inout) :: map - integer, intent(in) :: sze - integer(key_kind), intent(inout) :: key(sze) - integer :: i - integer(map_size_kind) :: idx_cache, idx_cache_prev - integer(cache_map_size_kind) :: ibegin, iend - integer(cache_map_size_kind), allocatable :: idx(:) - !DIR$ ATTRIBUTES ALIGN : 64 :: idx - - idx_cache_prev = -1_map_size_kind - allocate(idx(sze)) - do i=1,sze - idx_cache = shiftr(key(i),map_shift) - iend = map%map(idx_cache)%n_elements - if (idx_cache == idx_cache_prev) then - if ((idx(i-1) > 0_cache_map_size_kind).and.(idx(i-1) < iend)) then - if ((key(i) == key(i-1)+1).and.(map%map(idx_cache)%key(idx(i-1))+1) == key(i)) then - idx(i) = idx(i-1)+1 - cycle - endif - endif - endif - !DIR$ FORCEINLINE - call search_key_big_interval(key(i),map%map(idx_cache)%key, iend, idx(i), 1, iend) - idx_cache_prev = idx_cache - enddo - do i=1,sze - idx_cache = shiftr(key(i),map_shift) - if (idx(i) <= 0) then - key(i) = 0_key_kind - endif - enddo - deallocate(idx) -end - -subroutine search_key_big(key,X,sze,idx) - use map_module - implicit none - integer(cache_map_size_kind), intent(in) :: sze - integer(key_kind) , intent(in) :: key - integer(cache_key_kind) , intent(in) :: X(sze) - integer(cache_map_size_kind), intent(out) :: idx - - call search_key_big_interval(key,X,sze,idx,1,sze) -end - - -subroutine search_key_big_interval(key,X,sze,idx,ibegin_in,iend_in) - use map_module - implicit none - integer(cache_map_size_kind), intent(in) :: sze - integer(key_kind) , intent(in) :: key - integer(cache_key_kind) , intent(in) :: X(sze) - integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in - integer(cache_map_size_kind), intent(out) :: idx - - integer(cache_map_size_kind) :: istep, ibegin, iend, i - integer(cache_key_kind) :: cache_key - - if (sze /= 0) then - continue - else - idx = -1 - return - endif - cache_key = int(iand(key,map_mask),2) - ibegin = min(ibegin_in,sze) - iend = min(iend_in,sze) - if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then - - istep = shiftr(iend-ibegin,1) - idx = ibegin + istep - do while (istep > 4) - idx = ibegin + istep - ! TODO : Cache misses - if (cache_key < X(idx)) then - iend = idx - istep = shiftr(idx-ibegin,1) - idx = ibegin + istep - if (cache_key < X(idx)) then - iend = idx - istep = shiftr(idx-ibegin,1) - cycle - else if (cache_key > X(idx)) then - ibegin = idx - istep = shiftr(iend-idx,1) - cycle - else - return - endif - else if (cache_key > X(idx)) then - ibegin = idx - istep = shiftr(iend-idx,1) - idx = idx + istep - if (cache_key < X(idx)) then - iend = idx - istep = shiftr(idx-ibegin,1) - cycle - else if (cache_key > X(idx)) then - ibegin = idx - istep = shiftr(iend-idx,1) - cycle - else - return - endif - else - return - endif - enddo - idx = ibegin - if (min(iend_in,sze) > ibegin+4) then - iend = ibegin+4 - !DIR$ LOOP COUNT MAX(4) - do while (cache_key > X(idx)) - idx = idx+1 - end do - else - !DIR$ LOOP COUNT MAX(4) - do while (cache_key > X(idx)) - idx = idx+1 - if (idx == iend) then - exit - endif - end do - endif - if (cache_key /= X(idx)) then - idx = 1-idx - endif - return - - else - - if (cache_key < X(ibegin)) then - idx = -ibegin - return - endif - if (cache_key > X(iend)) then - idx = -iend - return - endif - if (cache_key == X(ibegin)) then - idx = ibegin - return - endif - if (cache_key == X(iend)) then - idx = iend - return - endif - endif - -end - -subroutine search_key_value_big_interval(key,value,X,Y,sze,idx,ibegin_in,iend_in) - use map_module - implicit none - integer(cache_map_size_kind), intent(in) :: sze - integer(key_kind) , intent(in) :: key - real(integral_kind) , intent(out) :: value - integer(cache_key_kind) , intent(in) :: X(sze) - real(integral_kind) , intent(in) :: Y(sze) - integer(cache_map_size_kind), intent(in) :: ibegin_in, iend_in - integer(cache_map_size_kind), intent(out) :: idx - - integer(cache_map_size_kind) :: istep, ibegin, iend, i - integer(cache_key_kind) :: cache_key - - if (sze /= 0) then - continue - else - idx = -1 - value = 0.d0 - return - endif - cache_key = int(iand(key,map_mask),2) - ibegin = min(ibegin_in,sze) - iend = min(iend_in,sze) - if ((cache_key > X(ibegin)) .and. (cache_key < X(iend))) then - - istep = shiftr(iend+ibegin,1) - idx = ibegin + istep - do while (istep > 4) - idx = ibegin + istep - ! TODO : Cache misses - if (cache_key < X(idx)) then - iend = idx - istep = shiftr(idx-ibegin,1) - idx = ibegin + istep - if (cache_key < X(idx)) then - iend = idx - istep = shiftr(idx-ibegin,1) - cycle - else if (cache_key > X(idx)) then - ibegin = idx - istep = shiftr(iend-idx,1) - cycle - else - value = Y(idx) - return - endif - else if (cache_key > X(idx)) then - ibegin = idx - istep = shiftr(iend-idx,1) - idx = idx + istep - if (cache_key < X(idx)) then - iend = idx - istep = shiftr(idx-ibegin,1) - cycle - else if (cache_key > X(idx)) then - ibegin = idx - istep = shiftr(iend-idx,1) - cycle - else - value = Y(idx) - return - endif - else - value = Y(idx) - return - endif - enddo - idx = ibegin - if (min(iend_in,sze) > ibegin+4) then - iend = ibegin+4 - !DIR$ LOOP COUNT MAX(4) - do while (cache_key > X(idx)) - idx = idx+1 - end do - else - !DIR$ LOOP COUNT MAX(4) - do while (cache_key > X(idx)) - idx = idx+1 - if (idx == iend) then - exit - endif - end do - endif - if (cache_key /= X(idx)) then - idx = 1-idx - value = 0.d0 - else - value = Y(idx) - endif - return - - else - - if (cache_key < X(ibegin)) then - idx = -ibegin - value = 0.d0 - return - endif - if (cache_key > X(iend)) then - idx = -iend - value = 0.d0 - return - endif - if (cache_key == X(ibegin)) then - idx = ibegin - value = Y(idx) - return - endif - if (cache_key == X(iend)) then - idx = iend - value = Y(idx) - return - endif - endif - -end - - -subroutine get_cache_map_n_elements_max(map,n_elements_max) - use map_module - implicit none - ! Returns the size of the largest cache_map - type (map_type), intent(in) :: map - integer(cache_map_size_kind), intent(out) :: n_elements_max - integer(map_size_kind) :: i - n_elements_max = 0_cache_map_size_kind - do i=0_8,map%map_size - n_elements_max = max(n_elements_max, map%map(i)%n_elements) - enddo -end - - - -subroutine get_cache_map(map,map_idx,keys,values,n_elements) - use map_module - implicit none - type (map_type), intent(in) :: map - integer(map_size_kind), intent(in) :: map_idx - integer(cache_map_size_kind), intent(inout) :: n_elements - integer(key_kind), intent(out) :: keys(n_elements) - double precision, intent(out) :: values(n_elements) - integer(cache_map_size_kind) :: i - integer(key_kind) :: shift - - shift = shiftl(map_idx,map_shift) - - n_elements = map%map(map_idx)%n_elements - do i=1,n_elements - keys(i) = map%map(map_idx)%key(i) + shift - values(i) = map%map(map_idx)%value(i) - enddo - -end - diff --git a/src/utils/map_module.mod b/src/utils/map_module.mod deleted file mode 100644 index a1277b96060691ab0fdd3a4803ba17d802a4cbc5..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 4548 zcmV;#5j*Z5iwFP!000006Wv`~bJ{o?zSpnN9pqxa(ZR<{N@ypu3CWP8yZvrFl+eag z0vRCL&G+j+$?`!b8#|a_I@_HVZ0zTe-jk$PufBY~p4?3aWB9ascv;Nh>wLY5S4-H# zZV%4Z^XYc}0Dr7kPw;8_9EqF$^CLu{-#>%NACUcH&p&5?J&FC_@6D#O$9crxw|#q_ z_h7Ghbq^ADZ?EJUW!FG&_4FM5oNuG)ayH*=SL3SP)8TaUUF`do%?Hl=i5BmJo`R&&*f8dT zi*y%(n92{_&Zke&CjOs#mjkiOIwSdJ!|r;R|Cqik_}{PbOuf~EIOmew`nH@t#WU$1 zfCGp_`JSigztMEIjbG={V*WZ`C`1T?1j~0XS*?FRn)8+l%XQ2I(w+#9%la-QJ4#020F9ulXB;$M5MJME|8!U zC9IDURHCE^_BN|E7Z=#Sy#8`AWR2&fIJ3_J64=#rJ3*ZRD+Da$MCerBSdOi+a;w)m z2qwMDXug;~&6nHF@nBHcTdZckBmhALfW1V-hTsg|+<&^f80Nj4z8vL0V zxEp-B92IJ4Q$NtBF6_-#%gt&LKX7LsnHCN1wrF5KwXAI`JPoYXx`DdZ?E_3O^)xVmADesTB5?cnQ3_?F4t#b@>vBXT#Gd>P$|R@S^6l)aqND=HJWj7hb_ zJ}d{?jDyEHyl$u8Ke#6ZboiC|`BddHz>eRHL11HW&>pD39yY`73+xPN)8X@A z(qHes-Eal&yUagCkJHT~Gm3f8Q5CN56LK&E$lwm{MzBA;9^c*#7*u_X08|6g%R=Cu zu%LaH`H~&pti_>LikH8aIge=m?{>cAgZ+M`n6~0U@)$=3!L5B}YRtI&03Zr8Q-rk= z9&W{Ft7jg;7Qb>`!>P3c=Jx*6ZVmtz>D`KOCzXsoKqtV|&Ng;nkqysKQLHJ~F z{gnl%7w5?s<7U8|^zF@P$mznHn;};*lt>&mjz8k%bfJuBUEd7UTL|L{Z#x+c?kEVf7N7TXq|eT2pKsLf(~ z&TnP0J=bcnJ+E^XI}1Upr%Dlro(%%u*CPH7+YRh}-^g|Yd(YSFTG@7kQrm6Jaa7uF zP;R?D=FG$*-LLsu^ebLI*v<9LU^2KGG2pcHLsc^^h6}AvjUt` z2&IKbZgEeFm{Us3N3eQDwb2yE`Dl0vt@`_AG-l>_DBhjS7iUnj$>W@rIg0a{yNSPV z=G%tzFT1Dcb9{rzWbiFwS~Z$P=cAj^_t{$@QjX5L5&uuA*_*rf|W!Xdy`=?RF+g6G3@ zmad%1)&*dv{+1ePTrA(ttuTUB)PDKdSDpR zf<}2BUPrVhbP6K+N@FOgjhQ9bnXV}GKGgL-^a}=j8YyO>dUv}bD2Vtvk+2}r)Bwx6 z4fet>>9%Ff8lxgYWO!8q7)i=qnNHdfONCLP;tV{{=}?uy|HdGw@8jbGSDlF+aHuM}z zdR}Zmsi!E%aw9}icpPeO#JBOWZ>VX9k2L(PGt#OW24liAZ~HuDJmI5J^Qx&GVz4Ex z0KT6T6{H+1DNH1bkdx`6A`!_V5F^UBkT|%F-gjh=i7&cFjO}$*3Q9Rui`vmTC>hve zDn|B!Vd```wjFV;Pl@eq6R5BWlyAbL))`9{0F^L1L_f7K-G`|oIJxYfo^{d>Q`wMn zI^V?!O8%7^{#AV6!_{WwNIs*})n=XTBWIKuNdiRS@kfbCfLbkplNLSbB#c8$?9lHW;Z#uxDuPbcdh z7grjBy2cRHr3RCZ>WV9kL0w}E>PpJ;&o(~9npiCg_6?5Ml@#4I72Q|BKOJipxeUDI z4PC< zZB**@(eWlLd?9l_8=3f?J^b}k3mk8{v7T3H2z<#N`-VNva$N*u9sI{9ap3^vFp``n zX>6WlW7h%40mv2yz7kD(7d3x6SwkcX7X<*MeYo37&OoBuw z_)8gdmvQS%>REXbG7hNDZMj;`Ply|;;O1DkZJy_|7oNQGT7zY0q|>R8PFH7nniPWy zik?pKWqFD=23Jt@C5k6zUdl6Daga`y7lF=i#4U{y@@d#;wpeXWN#$s{JcdXo_hY&^ zDRQWS9M;L5m{e9c^pNP>o~E1ae0@SaR}qS337iGH%`=~IJ0W%zt~_*fcH8yKDFtJ= zLN3!VC0QDYrZ3yo3Gu4ZH&5ne24%F$;T^W)zO7xA~LTDqDC(X2PW}3)UJe9 zyT=i9TzmHf2A$6KizkTq)_x&KV$yKtyygPaZW0L~P|YG)GbR|>S>S?1qiqVumotn7 zF?#|dK^4S0n!ZRd6_o7hizLgu$cm14cfp2;lxBiUbcd=$4zX*YRZ$BBo6z1bL}$}g zizc#}39;|vC9-%OGSR%vt|7vC5^C=n$4aTQ{Ucl_v9CfCDPB8?1V19ums_D;ScR0tULQ!082 z6QOC4_Og6OnD&9PB~oFDQxc^)Sf7Hlp z%_^4d$ZUK4=rUP#ea2GkEy_wD<`j7{&E|EwPT;dY0g`HZv*)q4wl7^s$ieyj)s>cQ zUR43hY-pZ2v6NBJqo808F{`K`VE8JGofx3N0BZG(6$4)fRWXDuU})A3n2*|&y_&-J z5V_tv+}$TYJ|4x=J5?=%+0O^JxA#}0^XT^Czg4AkOJ{XWPlWlP$Ht7d%dD!s55+8y zm=(mDO0Z3fnn^I`2YFeLG?nB1^_g$hvml}P2yH@PFMimRV5n!_I;Ki|BxGq^=Lp^G zX^rc&+qg&2W`5#Ox=fgI3Y&AAg3VJzRY(t>m}MG@#I{9~X9VK3d;VMFmT zV!=xy%Mm@_rtzDjBPZ($dr43U~c)>Iq z6w(5eW?e&gv9oNcy;##p@zm)hWr113Q#H$M!vj)0egz({1w27BJiMlbY|;3V`4L9u zgcm#DOV#<}yFi1Rh4QUPFLvxu=*LSh&r*`j8(&&JO7Ad_ldn~$NP;~~TvO#qVeo_v z4F+ju-s~{y>P!ZMnKyNfIK|S2eQPMJeBIiLrI`n6MOyh9tsKGgCD0FdgM6CA4B(qW zX{B!{(%Ykrz_}?1>kzLglvd<}pOmfm8hVsjU_y1PN{?ez9o_fR`a@>bfzA;-bN z!syLvfgk7}4G zutNlS5uqRPGLB}es)lY8X%=)~H3->>l&Z?&HVZy>ie~+zkrUD6XPnB)qQhxHS(+m9 zbY*e(Rp6S+Vg=ejS$w;)sPMa>EJ2B~91iLGl;sS7eos%)K+LH1uB4B)GvX0+D5!Y^ zJG<8@yH|dXI0}8oG2%>I?8l8OPYUWBiOnBnKRa7_R%jmehoBb#`(b(PQfMAGhafMk zlb-L!K4739o5sE*^1);5GnG&7VxJu~>wxXy?PM7QtIwW^gg_NfQ4c;rP(Bnz-;~6z(&Cr-W*7wJ zlRNY}IS9(9Vq{E|6DBMh9}qN7aL|rCK*_lR>K)kt6hY;Lfo!~vpjxkAs&x+#H0t1I zRk0B?E6a;&zx;}XK`L^NpjrN0LWlk0SN2le{3dI?ZaK4nNs{e3YMp4NCem|?#Z<&g z8w2x7L9V>ev3PEv|GKEM%i-(4mMhir=9Ye46n#ZdDF<8RS@R3#_sO%AGNnabPXx_^ zpJD*YuGlWHhjBTO>rpC291*ltc+B4AN|KP48bnYpK=_otJ*z8~7sumYh9LX0s5gh8 zREI7utb|$EO)9s=Stupb6a!OIc}y0%Maz^>*qPkSJBa?j_qLIe=fJ`Y*-Cdwwa~mLct1=U*$^*$exwZwSwJ zP4k~9+o=oG%3kDb&n}WS^IF|Y`HtAsW~3dmQQePai)MPiu|<$YTqr0WV)ZF=a{HMxmoT5sE;lK#Roc19qm?w2Tt@F!)P i5ai5#{h*BP5avb