From e5c2d5ee65bb6d716caf744402b5a1bce6ef1974 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 3 Apr 2015 13:34:27 +0200 Subject: [PATCH 01/27] add michel routine --- ocaml/Input.ml | 2 +- .../{convert_ezfio.sh => upgrade_1.0_2.0.sh} | 0 src/NEEDED_MODULES | 2 +- src/Pseudo/ASSUMPTIONS.rst | 0 src/Pseudo/Makefile | 6 + src/Pseudo/NEEDED_MODULES | 0 src/Pseudo/README.rst | 4 + src/Pseudo/int.f90 | 2022 +++++++++++++++++ 8 files changed, 2034 insertions(+), 2 deletions(-) rename scripts/ezfio_interface/{convert_ezfio.sh => upgrade_1.0_2.0.sh} (100%) create mode 100644 src/Pseudo/ASSUMPTIONS.rst create mode 100644 src/Pseudo/Makefile create mode 100644 src/Pseudo/NEEDED_MODULES create mode 100644 src/Pseudo/README.rst create mode 100644 src/Pseudo/int.f90 diff --git a/ocaml/Input.ml b/ocaml/Input.ml index 01bb54a0..fe155f80 100644 --- a/ocaml/Input.ml +++ b/ocaml/Input.ml @@ -8,4 +8,4 @@ include Input_determinants;; include Input_electrons;; include Input_mo_basis;; include Input_nuclei;; -include Input_auto_generated;; \ No newline at end of file +include Input_auto_generated;; diff --git a/scripts/ezfio_interface/convert_ezfio.sh b/scripts/ezfio_interface/upgrade_1.0_2.0.sh similarity index 100% rename from scripts/ezfio_interface/convert_ezfio.sh rename to scripts/ezfio_interface/upgrade_1.0_2.0.sh diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index a144c42c..35aa5ec3 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected MRCC +AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected MRCC Pseudo diff --git a/src/Pseudo/ASSUMPTIONS.rst b/src/Pseudo/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/Pseudo/Makefile b/src/Pseudo/Makefile new file mode 100644 index 00000000..5cf11b78 --- /dev/null +++ b/src/Pseudo/Makefile @@ -0,0 +1,6 @@ +# Define here all new external source files and objects.Don't forget to prefix the +# object files with IRPF90_temp/ +SRC=int.f90 +OBJ=IRPF90_temp/int.o + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/Pseudo/NEEDED_MODULES b/src/Pseudo/NEEDED_MODULES new file mode 100644 index 00000000..e69de29b diff --git a/src/Pseudo/README.rst b/src/Pseudo/README.rst new file mode 100644 index 00000000..ebb762c1 --- /dev/null +++ b/src/Pseudo/README.rst @@ -0,0 +1,4 @@ +======= + Module +======= + diff --git a/src/Pseudo/int.f90 b/src/Pseudo/int.f90 new file mode 100644 index 00000000..16aca8e0 --- /dev/null +++ b/src/Pseudo/int.f90 @@ -0,0 +1,2022 @@ +!! +!! Computation of Vps, matrix element of the +!! pseudo-potential centered at point C +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! Vps= < Phi_A | Vloc(C) + Vpp(C) | Phi_B> +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! Phi_M (M=A,B) Cartesian gaussian orbital centered at point M : +!! Phi_M = (x-M_x)**n^M_x *(y-M_y)**n^M_y *(z-M_z)**n^M_z exp(-g_M rM**2) +!! with rM**2=(x-M_x)**2 + (y-M_y)**2 + (z-M_z)**2 +!! +!!** Vloc(C)= \sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2) +!! +!!** Vpp(C)= \sum_{l=0}^lmax v_l(rC) \sum_{m=-l}^{m=l} |Y_lm> : +!! function Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +!! lmax of formula above +!! kmax of formula above +!! v_kl = array v_kl(kmax_max,0:lmax_max) +!! n_kl = array n_kl(kmax_max,0:lmax_max) +!! dz_kl = array dz_kl(kmax_max,0:lmax_max) +!! n_a(1),n_a(2),n_a(3) +!! a(1),a(2),a(3) +!! g_a +!! n_b(1),n_b(2),n_b(3) +!! b(1),b(2),b(3) +!! g_b +!! c(1),c(2),c(3) +!! +!! Routine computing : +!! function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) +!! klocmax of formula above +!! v_k = array v_k(klocmax_max) +!! n_k = array n_k(klocmax_max) +!! dz_k= array dz_k(klocmax_max) +!! Routine total matrix element : +!! function Vps(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) +!! +!! Routines Vps_num, Vpseudo_num, and Vloc_num = brute force numerical +!! estimations of the same integrals + + +!! Vps= +!! +!! with: Vloc(C)=\sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2) +!! Vpp(C)=\sum_{l=0}^lmax\sum_{k=1}^kmax v_kl rC**n_kl exp(-dz_kl rC**2)*|l> 0 factor* P_l^|m|(cos(theta)) cos (|m| phi) +! m = 0 1/sqrt(2) *factor* P_l^0(cos(theta)) +! m < 0 factor* P_l^|m|(cos(theta)) sin (|m| phi) +! +! x=cos(theta) + + double precision function ylm_real(l,m,x,phi) + implicit double precision (a-h,o-z) + DIMENSION PM(0:100,0:100) + MM=100 + pi=dacos(-1.d0) + iabs_m=iabs(m) + if(iabs_m.gt.l)stop 'm must be between -l and l' + factor= dsqrt( ((2*l+1)*fact(l-iabs_m))/(4.d0*pi*fact(l+iabs_m)) ) + if(dabs(x).gt.1.d0)then + print*,'pb. in ylm_no' + print*,'x=',x + stop + endif + call LPMN(MM,l,l,X,PM) + plm=PM(iabs_m,l) + coef=factor*plm + if(m.gt.0)ylm_real=dsqrt(2.d0)*coef*dcos(iabs_m*phi) + if(m.eq.0)ylm_real=coef + if(m.lt.0)ylm_real=dsqrt(2.d0)*coef*dsin(iabs_m*phi) + + fourpi=4.d0*dacos(-1.d0) + if(l.eq.0)ylm_real=dsqrt(1.d0/fourpi) + + xchap=dsqrt(1.d0-x**2)*dcos(phi) + ychap=dsqrt(1.d0-x**2)*dsin(phi) + zchap=x + if(l.eq.1.and.m.eq.1)ylm_real=dsqrt(3.d0/fourpi)*xchap + if(l.eq.1.and.m.eq.0)ylm_real=dsqrt(3.d0/fourpi)*zchap + if(l.eq.1.and.m.eq.-1)ylm_real=dsqrt(3.d0/fourpi)*ychap + + if(l.eq.2.and.m.eq.2)ylm_real=dsqrt(15.d0/16.d0/pi)*(xchap**2-ychap**2) + if(l.eq.2.and.m.eq.1)ylm_real=dsqrt(15.d0/fourpi)*xchap*zchap + if(l.eq.2.and.m.eq.0)ylm_real=dsqrt(5.d0/16.d0/pi)*(-xchap**2-ychap**2+2.d0*zchap**2) + if(l.eq.2.and.m.eq.-1)ylm_real=dsqrt(15.d0/fourpi)*ychap*zchap + if(l.eq.2.and.m.eq.-2)ylm_real=dsqrt(15.d0/fourpi)*xchap*ychap + + if(l.gt.2)stop 'l > 2 not coded!' + + end + +!! Routine Vpseudo is based on formumla (66) +!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976): +!! +!! Vpseudo= (4pi)**2* \sum_{l=0}^lmax \sum_{m=-l}^{l} +!! \sum{lambda=0}^{l+nA} \sum_{mu=-lambda}^{lambda} +!! \sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz} +!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3)* Y_{lambda mu}(AC_unit) +!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)* +!! bigI(lambda,mu,l,m,k1,k2,k3) +!! \sum{lambdap=0}^{l+nB} \sum_{mup=-lambdap}^{lambdap} +!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz} +!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p)* Y_{lambdap mup}(BC_unit) +!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)* +!! bigI(lambdap,mup,l,m,k1p,k2p,k3p)* +!! \sum_{k=1}^{kmax} v_kl(k,l)* +!! bigR(lambda,lambdap,k1+k2+k3+k1p+k2p+k3p+n_kl(k,l),g_a,g_b,AC,BC,dz_kl(k,l)) +!! +!! nA=nAx+nAy+nAz +!! nB=nBx+nBy+nBz +!! AC=|A-C| +!! AC_x= A_x-C_x, etc. +!! BC=|B-C| +!! AC_unit= vect(AC)/AC +!! BC_unit= vect(BC)/BC +!! bigI(lambda,mu,l,m,k1,k2,k3)= +!! \int dOmega Y_{lambda mu}(Omega) xchap^k1 ychap^k2 zchap^k3 Y_{l m}(Omega) +!! +!! bigR(lambda,lambdap,N,g_a,g_b,gamm_k,AC,BC) +!! = exp(-g_a* AC**2 -g_b* BC**2) * int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) +!! /int dx x^{ktot} exp(-g_k)*x**2) M_lambda(2 g_k D x) + +double precision function Vpseudo & +(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +implicit none +integer kmax_max,lmax_max,ntot_max,nkl_max +parameter (kmax_max=2,lmax_max=2,nkl_max=4) +parameter (ntot_max=10) +integer lmax,kmax,n_kl(kmax_max,0:lmax_max),l,k +double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) +double precision a(3),g_a,b(3),g_b,c(3) +double precision fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm +double precision theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big +double precision areal,freal,breal,t1,t2,int_prod_bessel +integer ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot +integer n_a(3),n_b(3) +double precision array_R(0:ntot_max+nkl_max,kmax_max,0:lmax_max,0:lmax_max+ntot_max,0:lmax_max+ntot_max) +double precision & +array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max) +double precision & +array_I_B(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max) + +double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max) +double precision array_coefs_B(0:ntot_max,0:ntot_max,0:ntot_max) +double precision arg + +fourpi=4.d0*dacos(-1.d0) +ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) +bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) +arg=g_a*ac**2+g_b*bc**2 +if(arg.gt.-dlog(10.d-20))then +Vpseudo=0.d0 +return +endif +freal=dexp(-arg) + +areal=2.d0*g_a*ac +breal=2.d0*g_b*bc +ntotA=n_a(1)+n_a(2)+n_a(3) +ntotB=n_b(1)+n_b(2)+n_b(3) +ntot=ntotA+ntotB + +if(ntot.gt.ntot_max)stop 'increase ntot_max' + +if(ac.eq.0.d0.and.bc.eq.0.d0)then + + accu=0.d0 + do k=1,kmax + do l=0,lmax + ktot=ntot+n_kl(k,l) + do m=-l,l + prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) + prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) + accu=accu+prod*prodp*v_kl(k,l)*freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal) + enddo + enddo + enddo + Vpseudo=accu*fourpi + return +endif + +if(ac.ne.0.d0.and.bc.ne.0.d0)then + + f=fourpi**2 + + theta_AC0=dacos( (a(3)-c(3))/ac ) + phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) + theta_BC0=dacos( (b(3)-c(3))/bc ) + phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc) + + do ktot=0,ntotA+ntotB+nkl_max + do lambda=0,lmax+ntotA + do lambdap=0,lmax+ntotB + do k=1,kmax + do l=0,lmax + array_R(ktot,k,l,lambda,lambdap)= & + freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal) + enddo + enddo + enddo + enddo + enddo + + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & + *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) + enddo + enddo + enddo + + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & + *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) + enddo + enddo + enddo + + accu=0.d0 + do l=0,lmax + do m=-l,l + + do lambda=0,l+ntotA + do mu=-lambda,lambda + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) + enddo + enddo + enddo + enddo + enddo + + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) + enddo + enddo + enddo + enddo + enddo + + do lambda=0,l+ntotA + do mu=-lambda,lambda + + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + + prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3) + + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + + prodp=ylm(lambdap,mup,theta_BC0,phi_BC0)*array_coefs_B(k1p,k2p,k3p)*array_I_B(lambdap,mup,k1p,k2p,k3p) + + do k=1,kmax + ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l) + accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,lambdap) + enddo + + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + Vpseudo=f*accu + return +endif + +if(ac.eq.0.d0.and.bc.ne.0.d0)then + + f=fourpi**1.5d0 + theta_BC0=dacos( (b(3)-c(3))/bc ) + phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc) + + areal=2.d0*g_a*ac + breal=2.d0*g_b*bc + freal=dexp(-g_a*ac**2-g_b*bc**2) + do ktot=0,ntotA+ntotB+nkl_max + do lambdap=0,lmax+ntotB + do k=1,kmax + do l=0,lmax + array_R(ktot,k,l,0,lambdap)= & + freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal) + enddo + enddo + enddo + enddo + + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & + *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) + enddo + enddo + enddo + + accu=0.d0 + do l=0,lmax + do m=-l,l + + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) + enddo + enddo + enddo + enddo + enddo + + prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) + + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + + prodp=array_coefs_B(k1p,k2p,k3p)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(lambdap,mup,k1p,k2p,k3p) + + do k=1,kmax + ktot=ntotA+k1p+k2p+k3p+n_kl(k,l) + accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,0,lambdap) + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + Vpseudo=f*accu + return +endif + +if(ac.ne.0.d0.and.bc.eq.0.d0)then + + f=fourpi**1.5d0 + theta_AC0=dacos( (a(3)-c(3))/ac ) + phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) + + areal=2.d0*g_a*ac + breal=2.d0*g_b*bc + freal=dexp(-g_a*ac**2-g_b*bc**2) + do ktot=0,ntotA+ntotB+nkl_max + do lambda=0,lmax+ntotA + do k=1,kmax + do l=0,lmax + array_R(ktot,k,l,lambda,0)= & + freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal) + enddo + enddo + enddo + enddo + + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & + *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) + enddo + enddo + enddo + + accu=0.d0 + do l=0,lmax + do m=-l,l + + do lambda=0,l+ntotA + do mu=-lambda,lambda + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) + enddo + enddo + enddo + enddo + enddo + + do lambda=0,l+ntotA + do mu=-lambda,lambda + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + prod=array_coefs_A(k1,k2,k3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(lambda,mu,k1,k2,k3) + prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) + do k=1,kmax + ktot=k1+k2+k3+ntotB+n_kl(k,l) + accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,0) + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + Vpseudo=f*accu + return +endif +end + +double precision function Vpseudo_num(npts,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +implicit none +integer kmax_max,lmax_max +parameter (kmax_max=2,lmax_max=2) +integer lmax,kmax, n_kl(kmax_max,0:lmax_max),l,m,k,kk +double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) +double precision a(3),g_a,b(3),g_b,c(3),ac(3),bc(3) +integer n_a(3),n_b(3),npts +double precision rmax,dr,sum,rC +double precision overlap_orb_ylm_brute + +do l=1,3 + ac(l)=a(l)-c(l) + bc(l)=b(l)-c(l) +enddo +dr=rmax/npts +sum=0.d0 +do l=0,lmax + do m=-l,l + do k=1,npts + rC=(k-1)*dr+dr/2.d0 + do kk=1,kmax + sum=sum+dr*v_kl(kk,l)*rC**(n_kl(kk,l)+2)*dexp(-dz_kl(kk,l)*rC**2) & + *overlap_orb_ylm_brute(npts,rC,n_a,ac,g_a,l,m) & + *overlap_orb_ylm_brute(npts,rC,n_b,bc,g_b,l,m) + enddo + enddo + enddo +enddo +Vpseudo_num=sum +return +end +!! Routine Vloc is a variation of formumla (66) +!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976) +!! without the projection operator +!! +!! Vloc= (4pi)**3/2* \sum_{k=1}^{klocmax} \sum_{l=0}^lmax \sum_{m=-l}^{l} +!!\sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz} +!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3) +!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)* +!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz} +!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p) +!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)* +!!\sum_{l=0}^lmax \sum_{m=-l}^{l} + +!! bigI(0,0,l,m,k1+k1p,k2+k2p,k3+k3p)*Y_{l m}(D_unit) +!! *v_k(k)* bigR(lambda,k1+k2+k3+k1p+k2p+k3p+n_k(k),g_a,g_b,AC,BC,dz_k(k)) +!! +!! nA=nAx+nAy+nAz +!! nB=nBx+nBy+nBz +!! D=(g_a AC+g_b BC) +!! D_unit= vect(D)/D +!! AC_x= A_x-C_x, etc. +!! BC=|B-C| +!! AC_unit= vect(AC)/AC +!! BC_unit= vect(BC)/BCA +!! +!! bigR(lambda,g_a,g_b,g_k,AC,BC) +!! = exp(-g_a* AC**2 -g_b* BC**2)* +!! I_loc= \int dx x**l *exp(-gam*x**2) M_n(ax) l=ktot+2 gam=g_a+g_b+dz_k(k) a=dreal n=l +!! M_n(x) modified spherical bessel function + + +double precision function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) +implicit none +integer klocmax_max,lmax_max,ntot_max +parameter (klocmax_max=10,lmax_max=2) +parameter (ntot_max=10) +integer klocmax +double precision v_k(klocmax_max),dz_k(klocmax_max),crochet,bigA +integer n_k(klocmax_max) +double precision a(3),g_a,b(3),g_b,c(3),d(3) +integer n_a(3),n_b(3),ntotA,ntotB,ntot,m +integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p +double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0 +double precision,allocatable :: array_R_loc(:,:,:) +double precision,allocatable :: array_coefs(:,:,:,:,:,:) +double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg + + + fourpi=4.d0*dacos(-1.d0) + f=fourpi**1.5d0 + ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) + bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) + arg=g_a*ac**2+g_b*bc**2 + if(arg.gt.-dlog(10.d-20))then + Vloc=0.d0 + return + endif + freal=dexp(-g_a*ac**2-g_b*bc**2) + + ntotA=n_a(1)+n_a(2)+n_a(3) + ntotB=n_b(1)+n_b(2)+n_b(3) + ntot=ntotA+ntotB + + d2=0.d0 + do i=1,3 + d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i)) + d2=d2+d(i)**2 + enddo + d2=dsqrt(d2) + dreal=2.d0*d2 + + theta_DC0=dacos(d(3)/d2) + phi_DC0=datan2(d(2)/d2,d(1)/d2) + +allocate (array_R_loc(-2:ntot_max+klocmax_max,klocmax_max,0:ntot_max)) +allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) + + do ktot=-2,ntotA+ntotB+klocmax + do l=0,ntot + do k=1,klocmax + array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) + enddo + enddo + enddo + + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & + *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) & + *binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & + *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) + enddo + enddo + enddo + enddo + enddo + enddo + + accu=0.d0 + do k=1,klocmax + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + + do l=0,ntot + do m=-l,l + prod=ylm(l,m,theta_DC0,phi_DC0)*array_coefs(k1,k2,k3,k1p,k2p,k3p) & + *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) + ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) + accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) + enddo + enddo + + enddo + enddo + enddo + enddo + enddo + enddo + enddo + Vloc=f*accu + + deallocate (array_R_loc) + deallocate (array_coefs) +end + +double precision function bigA(i,j,k) +implicit none +integer i,j,k +double precision fourpi,dblefact +fourpi=4.d0*dacos(-1.d0) +bigA=0.d0 +if(mod(i,2).eq.1)return +if(mod(j,2).eq.1)return +if(mod(k,2).eq.1)return +bigA=fourpi*dblefact(i-1)*dblefact(j-1)*dblefact(k-1)/dblefact(i+j+k+1) +end +!! +!! I_{lambda,mu,l,m}^{k1,k2,k3} = /int dOmega Y_{lambda mu} xchap^k1 ychap^k2 zchap^k3 Y_{lm} +!! + +double precision function bigI(lambda,mu,l,m,k1,k2,k3) +implicit none +integer lambda,mu,l,m,k1,k2,k3 +integer k,i,kp,ip +double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom,fact,coef_pm +pi=dacos(-1.d0) + +if(mu.gt.0.and.m.gt.0)then +sum=0.d0 +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu)))) +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +do k=0,mu/2 + do i=0,lambda-mu + do kp=0,m/2 + do ip=0,l-m + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3) + enddo + enddo + enddo +enddo +bigI=sum +return +endif + +if(mu.eq.0.and.m.eq.0)then +factor1=dsqrt((2*lambda+1)/(4.d0*pi)) +factor2=dsqrt((2*l+1)/(4.d0*pi)) +sum=0.d0 +do i=0,lambda + do ip=0,l + cylm=factor1*coef_pm(lambda,i) + cylmp=factor2*coef_pm(l,ip) + sum=sum+cylm*cylmp*bigA(k1,k2,i+ip+k3) + enddo +enddo +bigI=sum +return +endif + +if(mu.eq.0.and.m.gt.0)then +factor1=dsqrt((2*lambda+1)/(4.d0*pi)) +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +sum=0.d0 +do i=0,lambda + do kp=0,m/2 + do ip=0,l-m + cylm=factor1*coef_pm(lambda,i) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3) + enddo + enddo +enddo +bigI=sum +return +endif + +if(mu.gt.0.and.m.eq.0)then +sum=0.d0 +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu)))) +factor2=dsqrt((2*l+1)/(4.d0*pi)) +do k=0,mu/2 + do i=0,lambda-mu + do ip=0,l + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=factor2*coef_pm(l,ip) + sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3) + enddo + enddo +enddo +bigI=sum +return +endif + +if(mu.lt.0.and.m.lt.0)then +mu=-mu +m=-m +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu)))) +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +sum=0.d0 +do k=0,(mu-1)/2 + do i=0,lambda-mu + do kp=0,(m-1)/2 + do ip=0,l-m + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3) + enddo + enddo + enddo +enddo +mu=-mu +m=-m +bigI=sum +return +endif + +if(mu.eq.0.and.m.lt.0)then +m=-m +factor1=dsqrt((2*lambda+1)/(4.d0*pi)) +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +sum=0.d0 +do i=0,lambda + do kp=0,(m-1)/2 + do ip=0,l-m + cylm=factor1*coef_pm(lambda,i) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3) + enddo + enddo +enddo +m=-m +bigI=sum +return +endif + +if(mu.lt.0.and.m.eq.0)then +sum=0.d0 +mu=-mu +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu)))) +factor2=dsqrt((2*l+1)/(4.d0*pi)) +do k=0,(mu-1)/2 + do i=0,lambda-mu + do ip=0,l + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=factor2*coef_pm(l,ip) + sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3) + enddo + enddo +enddo +mu=-mu +bigI=sum +return +endif + +if(mu.gt.0.and.m.lt.0)then +sum=0.d0 +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu)))) +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +m=-m +do k=0,mu/2 + do i=0,lambda-mu + do kp=0,(m-1)/2 + do ip=0,l-m + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3) + enddo + enddo + enddo +enddo +m=-m +bigI=sum +return +endif + +if(mu.lt.0.and.m.gt.0)then +mu=-mu +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(4.d0*pi*fact(lambda+iabs(mu)))) +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +sum=0.d0 +do k=0,(mu-1)/2 + do i=0,lambda-mu + do kp=0,m/2 + do ip=0,l-m + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3) + enddo + enddo + enddo +enddo +bigI=sum +mu=-mu +return +endif + +stop 'pb in bigI!' +end + +double precision function crochet(n,g) +implicit none +integer n +double precision g,pi,dblefact,expo +pi=dacos(-1.d0) +expo=0.5d0*dfloat(n+1) +crochet=dblefact(n-1)/(2.d0*g)**expo +if(mod(n,2).eq.0)crochet=crochet*dsqrt(pi/2.d0) +end + +!! +!! overlap= = /int dOmega Ylm (x-center_x)**nx*(y-center_y)**nx*(z-center)**nx +!! *exp(-g*(r-center)**2) +!! +double precision function overlap_orb_ylm_brute(npts,r,npower_orb,center_orb,g_orb,l,m) +implicit none +integer npower_orb(3),l,m,i,j,npts +double precision u,g_orb,du,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3) +double precision x_orb,y_orb,z_orb,twopi,r +twopi=2.d0*dacos(-1.d0) +du=2.d0/npts +dphi=twopi/npts +overlap_orb_ylm_brute=0.d0 +do i=1,npts + u=-1.d0+du*(i-1)+du/2.d0 + sintheta=dsqrt(1.d0-u**2) + do j=1,npts + phi=dphi*(j-1)+dphi/2.d0 + x_orb=r*dcos(phi)*sintheta + y_orb=r*dsin(phi)*sintheta + z_orb=r*u + term=orb_phi(x_orb,y_orb,z_orb,npower_orb,center_orb,g_orb)*ylm_real(l,m,u,phi) + overlap_orb_ylm_brute= overlap_orb_ylm_brute+term*du*dphi + enddo +enddo +end + +double precision function overlap_orb_ylm_grid(nptsgrid,r_orb,npower_orb,center_orb,g_orb,l,m) +implicit none +!! PSEUDOS +integer nptsgridmax,nptsgrid +double precision coefs_pseudo,ptsgrid +parameter(nptsgridmax=50) +common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) +!!!!! +integer npower_orb(3),l,m,i +double precision x,g_orb,two_pi,dx,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3) +double precision x_orb,y_orb,z_orb,twopi,pi,cosphi,sinphi,xbid +pi=dacos(-1.d0) +twopi=2.d0*pi +overlap_orb_ylm_grid=0.d0 +do i=1,nptsgrid + x_orb=r_orb*ptsgrid(i,1) + y_orb=r_orb*ptsgrid(i,2) + z_orb=r_orb*ptsgrid(i,3) + x=ptsgrid(i,3) + phi=datan2(ptsgrid(i,2),ptsgrid(i,1)) + term=orb_phi(x_orb,y_orb,z_orb,npower_orb,center_orb,g_orb)*ylm_real(l,m,x,phi) + overlap_orb_ylm_grid= overlap_orb_ylm_grid+coefs_pseudo(i)*term +enddo +overlap_orb_ylm_grid=2.d0*twopi*overlap_orb_ylm_grid +end + +! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) exp(i m phi) +! l=0,1,2,.... +! m=0,1,...,l +! Here: +! n=l (n=0,1,...) +! m=0,1,...,n +! x=cos(theta) 0 < x < 1 +! +! +! This routine computes: PM(m,n) for n=0,...,N (number N in input) and m=0,..,n + +! Exemples (see 'Associated Legendre Polynomilas wikipedia') +! P_{0}^{0}(x)=1 +! P_{1}^{-1}(x)=-1/2 P_{1}^{1}(x) +! P_{1}^{0}(x)=x +! P_{1}^{1}(x)=-(1-x^2)^{1/2} +! P_{2}^{-2}(x)=1/24 P_{2}^{2}(x) +! P_{2}^{-1}(x)=-1/6 P_{2}^{1}(x) +! P_{2}^{0}(x)=1/2 (3x^{2}-1) +! P_{2}^{1}(x)=-3x(1-x^2)^{1/2} +! P_{2}^{2}(x)=3(1-x^2) + + + SUBROUTINE LPMN(MM,M,N,X,PM) +! +! Here N = LMAX +! Here M= MMAX (we take M=LMAX in input) +! +! ===================================================== +! Purpose: Compute the associated Legendre functions Pmn(x) +! Input : x --- Argument of Pmn(x) +! m --- Order of Pmn(x), m = 0,1,2,...,n +! n --- Degree of Pmn(x), n = 0,1,2,...,N +! mm --- Physical dimension of PM +! Output: PM(m,n) --- Pmn(x) +! ===================================================== +! + IMPLICIT DOUBLE PRECISION (P,X) + DIMENSION PM(0:MM,0:(N+1)) + DO 10 I=0,N + DO 10 J=0,M +10 PM(J,I)=0.0D0 + PM(0,0)=1.0D0 + IF (DABS(X).EQ.1.0D0) THEN + DO 15 I=1,N +15 PM(0,I)=X**I + RETURN + ENDIF + LS=1 + IF (DABS(X).GT.1.0D0) LS=-1 + XQ=DSQRT(LS*(1.0D0-X*X)) + XS=LS*(1.0D0-X*X) + DO 30 I=1,M +30 PM(I,I)=-LS*(2.0D0*I-1.0D0)*XQ*PM(I-1,I-1) + DO 35 I=0,M +35 PM(I,I+1)=(2.0D0*I+1.0D0)*X*PM(I,I) + + DO 40 I=0,M + DO 40 J=I+2,N + PM(I,J)=((2.0D0*J-1.0D0)*X*PM(I,J-1)- (I+J-1.0D0)*PM(I,J-2))/(J-I) +40 CONTINUE + END + + +! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 +! P_l^|m|(cos(theta)) exp(i m phi) + + double precision function fact(n) + implicit double precision(a-h,o-z) + fact=1.d0 + if(n.eq.0)return + do i=1,n + fact=fact*dfloat(i) + enddo + end + + subroutine erreur(x,n,rmoy,error) + implicit double precision(a-h,o-z) + dimension x(n) +! calcul de la moyenne + rmoy=0.d0 + do i=1,n + rmoy=rmoy+x(i) + enddo + rmoy=rmoy/dfloat(n) +! calcul de l'erreur + error=0.d0 + do i=1,n + error=error+(x(i)-rmoy)**2 + enddo + if(n.gt.1)then + rn=dfloat(n) + rn1=dfloat(n-1) + error=dsqrt(error)/dsqrt(rn*rn1) + else + write(2,*)'Seulement un block Erreur nondefinie' + error=0.d0 + endif + end + + subroutine initpseudos(nptsgrid) + implicit none + integer nptsgridmax,nptsgrid,ik + double precision coefs_pseudo,ptsgrid + double precision p,q,r,s + parameter(nptsgridmax=50) + common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) + + p=1.d0/dsqrt(2.d0) + q=1.d0/dsqrt(3.d0) + r=1.d0/dsqrt(11.d0) + s=3.d0/dsqrt(11.d0) + + if(nptsgrid.eq.4)then + + ptsgrid(1,1)=q + ptsgrid(1,2)=q + ptsgrid(1,3)=q + + ptsgrid(2,1)=q + ptsgrid(2,2)=-q + ptsgrid(2,3)=-q + + ptsgrid(3,1)=-q + ptsgrid(3,2)=q + ptsgrid(3,3)=-q + + ptsgrid(4,1)=-q + ptsgrid(4,2)=-q + ptsgrid(4,3)=q + + do ik=1,4 + coefs_pseudo(ik)=1.d0/4.d0 + enddo + return + endif + + ptsgrid(1,1)=1.d0 + ptsgrid(1,2)=0.d0 + ptsgrid(1,3)=0.d0 + + ptsgrid(2,1)=-1.d0 + ptsgrid(2,2)=0.d0 + ptsgrid(2,3)=0.d0 + + ptsgrid(3,1)=0.d0 + ptsgrid(3,2)=1.d0 + ptsgrid(3,3)=0.d0 + + ptsgrid(4,1)=0.d0 + ptsgrid(4,2)=-1.d0 + ptsgrid(4,3)=0.d0 + + ptsgrid(5,1)=0.d0 + ptsgrid(5,2)=0.d0 + ptsgrid(5,3)=1.d0 + + ptsgrid(6,1)=0.d0 + ptsgrid(6,2)=0.d0 + ptsgrid(6,3)=-1.d0 + + do ik=1,6 + coefs_pseudo(ik)=1.d0/6.d0 + enddo + + if(nptsgrid.eq.6)return + + ptsgrid(7,1)=p + ptsgrid(7,2)=p + ptsgrid(7,3)=0.d0 + + ptsgrid(8,1)=p + ptsgrid(8,2)=-p + ptsgrid(8,3)=0.d0 + + ptsgrid(9,1)=-p + ptsgrid(9,2)=p + ptsgrid(9,3)=0.d0 + + ptsgrid(10,1)=-p + ptsgrid(10,2)=-p + ptsgrid(10,3)=0.d0 + + ptsgrid(11,1)=p + ptsgrid(11,2)=0.d0 + ptsgrid(11,3)=p + + ptsgrid(12,1)=p + ptsgrid(12,2)=0.d0 + ptsgrid(12,3)=-p + + ptsgrid(13,1)=-p + ptsgrid(13,2)=0.d0 + ptsgrid(13,3)=p + + ptsgrid(14,1)=-p + ptsgrid(14,2)=0.d0 + ptsgrid(14,3)=-p + + ptsgrid(15,1)=0.d0 + ptsgrid(15,2)=p + ptsgrid(15,3)=p + + ptsgrid(16,1)=0.d0 + ptsgrid(16,2)=p + ptsgrid(16,3)=-p + + ptsgrid(17,1)=0.d0 + ptsgrid(17,2)=-p + ptsgrid(17,3)=p + + ptsgrid(18,1)=0.d0 + ptsgrid(18,2)=-p + ptsgrid(18,3)=-p + + do ik=1,6 + coefs_pseudo(ik)=1.d0/30.d0 + enddo + do ik=7,18 + coefs_pseudo(ik)=1.d0/15.d0 + enddo + + if(nptsgrid.eq.18)return + + ptsgrid(19,1)=q + ptsgrid(19,2)=q + ptsgrid(19,3)=q + + ptsgrid(20,1)=-q + ptsgrid(20,2)=q + ptsgrid(20,3)=q + + ptsgrid(21,1)=q + ptsgrid(21,2)=-q + ptsgrid(21,3)=q + + ptsgrid(22,1)=q + ptsgrid(22,2)=q + ptsgrid(22,3)=-q + + ptsgrid(23,1)=-q + ptsgrid(23,2)=-q + ptsgrid(23,3)=q + + ptsgrid(24,1)=-q + ptsgrid(24,2)=q + ptsgrid(24,3)=-q + + ptsgrid(25,1)=q + ptsgrid(25,2)=-q + ptsgrid(25,3)=-q + + ptsgrid(26,1)=-q + ptsgrid(26,2)=-q + ptsgrid(26,3)=-q + + do ik=1,6 + coefs_pseudo(ik)=1.d0/21.d0 + enddo + do ik=7,18 + coefs_pseudo(ik)=4.d0/105.d0 + enddo + do ik=19,26 + coefs_pseudo(ik)=27.d0/840.d0 + enddo + + if(nptsgrid.eq.26)return + + ptsgrid(27,1)=r + ptsgrid(27,2)=r + ptsgrid(27,3)=s + + ptsgrid(28,1)=r + ptsgrid(28,2)=-r + ptsgrid(28,3)=s + + ptsgrid(29,1)=-r + ptsgrid(29,2)=r + ptsgrid(29,3)=s + + ptsgrid(30,1)=-r + ptsgrid(30,2)=-r + ptsgrid(30,3)=s + + ptsgrid(31,1)=r + ptsgrid(31,2)=r + ptsgrid(31,3)=-s + + ptsgrid(32,1)=r + ptsgrid(32,2)=-r + ptsgrid(32,3)=-s + + ptsgrid(33,1)=-r + ptsgrid(33,2)=r + ptsgrid(33,3)=-s + + ptsgrid(34,1)=-r + ptsgrid(34,2)=-r + ptsgrid(34,3)=-s + + ptsgrid(35,1)=r + ptsgrid(35,2)=s + ptsgrid(35,3)=r + + ptsgrid(36,1)=-r + ptsgrid(36,2)=s + ptsgrid(36,3)=r + + ptsgrid(37,1)=r + ptsgrid(37,2)=s + ptsgrid(37,3)=-r + + ptsgrid(38,1)=-r + ptsgrid(38,2)=s + ptsgrid(38,3)=-r + + ptsgrid(39,1)=r + ptsgrid(39,2)=-s + ptsgrid(39,3)=r + + ptsgrid(40,1)=r + ptsgrid(40,2)=-s + ptsgrid(40,3)=-r + + ptsgrid(41,1)=-r + ptsgrid(41,2)=-s + ptsgrid(41,3)=r + + ptsgrid(42,1)=-r + ptsgrid(42,2)=-s + ptsgrid(42,3)=-r + + ptsgrid(43,1)=s + ptsgrid(43,2)=r + ptsgrid(43,3)=r + + ptsgrid(44,1)=s + ptsgrid(44,2)=r + ptsgrid(44,3)=-r + + ptsgrid(45,1)=s + ptsgrid(45,2)=-r + ptsgrid(45,3)=r + + ptsgrid(46,1)=s + ptsgrid(46,2)=-r + ptsgrid(46,3)=-r + + ptsgrid(47,1)=-s + ptsgrid(47,2)=r + ptsgrid(47,3)=r + + ptsgrid(48,1)=-s + ptsgrid(48,2)=r + ptsgrid(48,3)=-r + + ptsgrid(49,1)=-s + ptsgrid(49,2)=-r + ptsgrid(49,3)=r + + ptsgrid(50,1)=-s + ptsgrid(50,2)=-r + ptsgrid(50,3)=-r + + do ik=1,6 + coefs_pseudo(ik)=4.d0/315.d0 + enddo + do ik=7,18 + coefs_pseudo(ik)=64.d0/2835.d0 + enddo + do ik=19,26 + coefs_pseudo(ik)=27.d0/1280.d0 + enddo + do ik=27,50 + coefs_pseudo(ik)=14641.d0/725760.d0 + enddo + + if(nptsgrid.eq.50)return + + write(*,*)'Grid for pseudos not available!' + write(*,*)'N=4-6-18-26-50 only!' + stop + end + +double precision function dblefact(n) +implicit none +integer :: n,k +double precision prod +dblefact=1.d0 +if(n.lt.0)return +if(mod(n,2).eq.1)then + prod=1.d0 + do k=1,n,2 + prod=prod*dfloat(k) + enddo + dblefact=prod + return + endif + if(mod(n,2).eq.0)then + prod=1.d0 + do k=2,n,2 + prod=prod*dfloat(k) + enddo + dblefact=prod + return + endif +end +!! +!! R_{lambda,lamba',N}= exp(-ga_a AC**2 -g_b BC**2) /int_{0}{+infty} r**(2+n) exp(-(g_a+g_b+g_k)r**2) +!! * M_{lambda}( 2g_a ac r) M_{lambda'}(2g_b bc r) +!! + double precision function bigR(lambda,lambdap,n,g_a,g_b,ac,bc,g_k) + implicit none + integer lambda,lambdap,n,npts,i + double precision g_a,g_b,ac,bc,g_k,arg,factor,delta1,delta2,cc,rmax,dr,sum,x1,x2,r + double precision bessel_mod + arg=g_a*ac**2+g_b*bc**2 + factor=dexp(-arg) + delta1=2.d0*g_a*ac + delta2=2.d0*g_b*bc + cc=g_a+g_b+g_k + if(cc.eq.0.d0)stop 'pb. in bigR' + rmax=dsqrt(-dlog(10.d-20)/cc) + npts=500 + dr=rmax/npts + sum=0.d0 + do i=1,npts + r=(i-1)*dr + x1=delta1*r + x2=delta2*r + sum=sum+dr*r**(n+2)*dexp(-cc*r**2)*bessel_mod(x1,lambda)*bessel_mod(x2,lambdap) + enddo + bigR=sum*factor + end + + double precision function bessel_mod(x,n) + implicit none + integer n + double precision x,bessel_mod_exp,bessel_mod_recur + if(x.le.0.8d0)then + bessel_mod=bessel_mod_exp(n,x) + else + bessel_mod=bessel_mod_recur(n,x) + endif + end + + recursive function bessel_mod_recur(n,x) result(a) + implicit none + integer n + double precision x,a,bessel_mod_exp + if(x.le.0.8d0)then + a=bessel_mod_exp(n,x) + return + endif + if(n.eq.0)a=dsinh(x)/x + if(n.eq.1)a=(x*dcosh(x)-dsinh(x))/x**2 + if(n.ge.2)a=bessel_mod_recur(n-2,x)-(2*n-1)/x*bessel_mod_recur(n-1,x) + end + + double precision function bessel_mod_exp(n,x) + implicit none + integer n,k + double precision x,coef,accu,fact,dblefact + accu=0.d0 + do k=0,10 + coef=1.d0/fact(k)/dblefact(2*(n+k)+1) + accu=accu+(x**2/2.d0)**k*coef + enddo + bessel_mod_exp=x**n*accu + end + +! double precision function bessel_mod(x,n) +! IMPLICIT DOUBLE PRECISION (A-H,O-Z) +! parameter(NBESSMAX=100) +! dimension SI(0:NBESSMAX),DI(0:NBESSMAX) +! if(n.lt.0.or.n.gt.NBESSMAX)stop 'pb with argument of bessel_mod' +! CALL SPHI(N,X,NBESSMAX,SI,DI) +! bessel_mod=si(n) +! end + + SUBROUTINE SPHI(N,X,NMAX,SI,DI) +! +! ======================================================== +! Purpose: Compute modified spherical Bessel functions +! of the first kind, in(x) and in'(x) +! Input : x --- Argument of in(x) +! n --- Order of in(x) ( n = 0,1,2,... ) +! Output: SI(n) --- in(x) +! DI(n) --- in'(x) +! NM --- Highest order computed +! Routines called: +! MSTA1 and MSTA2 for computing the starting +! point for backward recurrence +! ======================================================== +! + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + DIMENSION SI(0:NMAX),DI(0:NMAX) + NM=N + IF (DABS(X).LT.1.0D-100) THEN + DO 10 K=0,N + SI(K)=0.0D0 +10 DI(K)=0.0D0 + SI(0)=1.0D0 + DI(1)=0.333333333333333D0 + RETURN + ENDIF + SI(0)=DSINH(X)/X + SI(1)=-(DSINH(X)/X-DCOSH(X))/X + SI0=SI(0) + IF (N.GE.2) THEN + M=MSTA1(X,200) + + write(34,*)'m=',m + + IF (M.LT.N) THEN + NM=M + ELSE + M=MSTA2(X,N,15) + write(34,*)'m=',m + ENDIF + F0=0.0D0 + F1=1.0D0-100 + DO 15 K=M,0,-1 + F=(2.0D0*K+3.0D0)*F1/X+F0 + IF (K.LE.NM) SI(K)=F + F0=F1 +15 F1=F + CS=SI0/F + write(34,*)'cs=',cs + DO 20 K=0,NM +20 SI(K)=CS*SI(K) + ENDIF + DI(0)=SI(1) + DO 25 K=1,NM +25 DI(K)=SI(K-1)-(K+1.0D0)/X*SI(K) + RETURN + END + + + INTEGER FUNCTION MSTA1(X,MP) +! +! =================================================== +! Purpose: Determine the starting point for backward +! recurrence such that the magnitude of +! Jn(x) at that point is about 10^(-MP) +! Input : x --- Argument of Jn(x) +! MP --- Value of magnitude +! Output: MSTA1 --- Starting point +! =================================================== +! + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + A0=DABS(X) + N0=INT(1.1*A0)+1 + F0=ENVJ(N0,A0)-MP + N1=N0+5 + F1=ENVJ(N1,A0)-MP + DO 10 IT=1,20 + NN=N1-(N1-N0)/(1.0D0-F0/F1) + F=ENVJ(NN,A0)-MP + IF(ABS(NN-N1).LT.1) GO TO 20 + N0=N1 + F0=F1 + N1=NN + 10 F1=F + 20 MSTA1=NN + RETURN + END + + + INTEGER FUNCTION MSTA2(X,N,MP) +! +! =================================================== +! Purpose: Determine the starting point for backward +! recurrence such that all Jn(x) has MP +! significant digits +! Input : x --- Argument of Jn(x) +! n --- Order of Jn(x) +! MP --- Significant digit +! Output: MSTA2 --- Starting point +! =================================================== +! + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + A0=DABS(X) + HMP=0.5D0*MP + EJN=ENVJ(N,A0) + IF (EJN.LE.HMP) THEN + OBJ=MP + N0=INT(1.1*A0) + ELSE + OBJ=HMP+EJN + N0=N + ENDIF + F0=ENVJ(N0,A0)-OBJ + N1=N0+5 + F1=ENVJ(N1,A0)-OBJ + DO 10 IT=1,20 + NN=N1-(N1-N0)/(1.0D0-F0/F1) + F=ENVJ(NN,A0)-OBJ + IF (iABS(NN-N1).LT.1) GO TO 20 + N0=N1 + F0=F1 + N1=NN +10 F1=F +20 MSTA2=NN+10 + RETURN + END + + double precision FUNCTION ENVJ(N,X) + DOUBLE PRECISION X + integer N + ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N) + RETURN + END + +!c Computation of real spherical harmonics Ylm(theta,phi) +!c +!c l=0,1,.... +!c m=-l,l +!c +!c m>0: Y_lm = sqrt(2) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) cos(m phi) +!c m=0: Y_l0 = ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^0 (cos(theta)) +!c m<0: Y_lm = sqrt(2) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) sin(|m|phi) + +!Examples(wikipedia http://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics) + +! l = 0 + +! Y_00 = \sqrt{1/(4pi)} + +! l = 1 + +! Y_1-1= \sqrt{3/(4pi)} y/r +! Y_10 = \sqrt{3/(4pi)} z/r +! Y_11 = \sqrt{3/(4pi)} x/r +! +! l = 2 +! +! Y_2,-2= 1/2 \sqrt{15/pi} xy/r^2 +! Y_2,-1= 1/2 \sqrt{15/pi} yz/r^2 +! Y_20 = 1/4 \sqrt{15/pi} (-x^2-y^2 +2z^2)/r^2 +! Y_21 = 1/2 \sqrt{15/pi} zx/r^2 +! Y_22 = 1/4 \sqrt{15/pi} (x^2-y^2)/r^2 +! +!c +double precision function ylm(l,m,theta,phi) +implicit none +integer l,m +double precision theta,phi,pm,factor,pi,x,fact,sign +DIMENSION PM(0:100,0:100) +pi=dacos(-1.d0) +x=dcos(theta) +sign=(-1.d0)**m +CALL LPMN(100,l,l,X,PM) +factor=dsqrt( (2*l+1)*fact(l-iabs(m)) /(4.d0*pi*fact(l+iabs(m))) ) +if(m.gt.0)ylm=sign*dsqrt(2.d0)*factor*pm(m,l)*dcos(dfloat(m)*phi) +if(m.eq.0)ylm=factor*pm(m,l) +if(m.lt.0)ylm=sign*dsqrt(2.d0)*factor*pm(iabs(m),l)*dsin(dfloat(iabs(m))*phi) +end + +!c Explicit representation of Legendre polynomials P_n(x) +!! +!! P_n0(x) = P_n(x)= \sum_{k=0}^n a_k x^k +!! +!! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n ) +!! coef_pm(n,k) is the k_th coefficient of P_n(x) +double precision function coef_pm(n,k) +implicit none +integer n,k +double precision arg,binom,binom_gen +if(n.eq.0.and.k.ne.0)stop 'coef_pm not defined' +if(n.eq.0.and.k.eq.0)then +coef_pm=1.d0 +return +endif +arg=0.5d0*dfloat(n+k-1) +coef_pm=2.d0**n*binom(n,k)*binom_gen(arg,n) +end + +!! Ylm_bis uses the series expansion of Ylm in xchap^i ychap^j zchap^k +!! xchap=x/r etc. +!c m>0: Y_lm = sqrt(2)*factor* P_l^|m|(cos(theta)) cos(m phi) +!c m=0: Y_l0 = factor* P_l^0 (cos(theta)) +!c m<0: Y_lm = sqrt(2) factor* P_l^|m|(cos(theta)) sin(|m|phi) +!c factor= ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 + +!! P_l^m (x) = (-1)**m (1-x**2)^m/2 d^m/dx^m P_l(x) m >0 or 0 +!! the series expansion of P_m (x) is known +!! +!! sin(theta)**m cos(mphi) = \sum_0^[m/2] binom(m,2k) x^(m-2k) y^2k (-1)**k (easy to proove with +!! Moivre formula) +!! (here x = xchap...) +!! +!! Ylm m> 0 = \sum_{k=0}^[m/2] \sum_{i=0}^(l-m) c_ki x^(m-2k) y^2k z^i +!! +!! c_ki= (-1)^k sqrt(2)*factor*binom(m,2k)*(m+i)!/i!*coef_pm(l,i+m) +!! +!! Ylm m< 0 = \sum_{k=0}^[(m-1)/2] \sum_{i=0}^(l-m) c_ki x^(m-(2k+1)) y^(2k+1) z^i +!! +!! c_ki= (-1)^k sqrt(2)*factor*binom(m,2k+1)*(m+i)!/i!*coef_pm(l,i+m) + + +double precision function ylm_bis(l,m,theta,phi) +implicit none +integer l,m,k,i +double precision x,y,z,theta,phi,sum,factor,pi,binom,fact,coef_pm,cylm +pi=dacos(-1.d0) +x=dsin(theta)*dcos(phi) +y=dsin(theta)*dsin(phi) +z=dcos(theta) +factor=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +if(m.gt.0)then +sum=0.d0 +do k=0,m/2 + do i=0,l-m + cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom(m,2*k)*fact(m+i)/fact(i)*coef_pm(l,i+m) + sum=sum+cylm*x**(m-2*k)*y**(2*k)*z**i + enddo +enddo +ylm_bis=sum +return +endif +if(m.eq.0)then +sum=0.d0 +do i=0,l + sum=sum+factor*coef_pm(l,i)*z**i +enddo +ylm_bis=sum +return +endif +if(m.lt.0)then +m=-m +sum=0.d0 +do k=0,(m-1)/2 + do i=0,l-m + cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom(m,2*k+1)*fact(m+i)/fact(i)*coef_pm(l,i+m) + sum=sum+cylm*x**(m-(2*k+1))*y**(2*k+1)*z**i + enddo +enddo +ylm_bis=sum +m=-m +return +endif +end + +!c +!c Computation of associated Legendre Polynomials PM(m,n) for n=0,...,N +!c Here: +!c n=l (n=0,1,...) +!c m=0,1,...,n +!c x=cos(theta) 0 < x < 1 +!c +!c This routine computes: PM(m,n) for n=0,...,N (number N in input) and m=0,..,n +!c Exemples (see 'Associated Legendre Polynomilas wikipedia') +!c P_{0}^{0}(x)=1 +!c P_{1}^{-1}(x)=-1/2 P_{1}^{1}(x) +!c P_{1}^{0}(x)=x +!c P_{1}^{1}(x)=-(1-x^2)^{1/2} +!c P_{2}^{-2}(x)=1/24 P_{2}^{2}(x) +!c P_{2}^{-1}(x)=-1/6 P_{2}^{1}(x) +!c P_{2}^{0}(x)=1/2 (3x^{2}-1) +!c P_{2}^{1}(x)=-3x(1-x^2)^{1/2} +!c P_{2}^{2}(x)=3(1-x^2) +!c +!c Explicit representation: +!! +!! P_n0(x) = P_n(x)= \sum_{k=0}^n a_k x^k +!! +!! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n ) + +double precision function binom(i,j) + implicit none + integer :: i,j + double precision :: fact + binom = fact(i)/(fact(j)*fact(i-j)) +end + +double precision function binom_gen(alpha,n) + implicit none + integer :: n,k + double precision :: fact,alpha,prod + prod=1.d0 + do k=0,n-1 + prod=prod*(alpha-k) + binom_gen = prod/(fact(n)) + enddo +end + +double precision function test_int(g_a,g_b,g_c,ac,bc) +implicit none +double precision factor,g_a,g_b,g_c,ac,bc,x,dx,sum,alpha,beta,pi +integer i,npts +pi=dacos(-1.d0) +factor=0.5d0*pi/(g_a*g_b*ac*bc*dsqrt(g_a+g_b+g_c))*dexp(-g_a*ac**2-g_b*bc**2) +npts=2000 +dx=20.d0/npts +sum=0.d0 +alpha=(2.d0*g_a*ac+2.d0*g_b*bc)/dsqrt(g_c+g_a+g_b) +beta=(2.d0*g_b*bc-2.d0*g_b*bc)/dsqrt(g_c+g_a+g_b) +do i=1,npts + x=(i-1)*dx+0.5d0*dx + sum=sum+dx*dexp(-x**2)*(dcosh(alpha*x)-dcosh(beta*x)) +enddo +test_int=factor*sum +end + +recursive function fact1(n,a) result(x) +implicit none +integer n +double precision a,x,erf +if(n.eq.0)then +x=dsqrt(dacos(-1.d0))/2.d0*erf(a) +return +endif +if(n.eq.1)then +x=1.d0-dexp(-a**2) +return +endif +if(mod(n,2).eq.0)x=0.5d0*dfloat(n-1)*fact1(n-2,a)+a**n*dexp(-a**2) +if(mod(n,2).eq.1)x=0.5d0*dfloat(n-1)*fact1(n-2,a)+0.5d0*a**(n-1)*dexp(-a**2) +end + + double precision FUNCTION ERF(X) + implicit double precision(a-h,o-z) + IF(X.LT.0.d0)THEN + ERF=-GAMMP(.5d0,X**2) + ELSE + ERF=GAMMP(.5d0,X**2) + ENDIF + RETURN + END + + double precision FUNCTION GAMMLN(XX) + implicit double precision(a-h,o-z) + REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER + DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, & + -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ + DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ + X=XX-ONE + TMP=X+FPF + TMP=(X+HALF)*DLOG(TMP)-TMP + SER=ONE + DO 11 J=1,6 + X=X+ONE + SER=SER+COF(J)/X +11 CONTINUE + GAMMLN=TMP+DLOG(STP*SER) + RETURN + END + FUNCTION GAMMP(A,X) + implicit double precision(a-h,o-z) + IF(X.LT.0.d0.OR.A.LE.0.d0)PAUSE + IF(X.LT.A+1.d0)THEN + CALL GSER(GAMMP,A,X,GLN) + ELSE + CALL GCF(GAMMCF,A,X,GLN) + GAMMP=1.d0-GAMMCF + ENDIF + RETURN + END + SUBROUTINE GCF(GAMMCF,A,X,GLN) + implicit double precision(a-h,o-z) + PARAMETER (ITMAX=100,EPS=3.D-7) + GLN=GAMMLN(A) + GOLD=0.d0 + A0=1.d0 + A1=X + B0=0.d0 + B1=1.d0 + FAC=1.d0 + DO 11 N=1,ITMAX + AN=DFLOAT(N) + ANA=AN-A + A0=(A1+A0*ANA)*FAC + B0=(B1+B0*ANA)*FAC + ANF=AN*FAC + A1=X*A0+ANF*A1 + B1=X*B0+ANF*B1 + IF(A1.NE.0.d0)THEN + FAC=1.d0/A1 + G=B1*FAC + IF(DABS((G-GOLD)/G).LT.EPS)GO TO 1 + GOLD=G + ENDIF +11 CONTINUE + PAUSE 'A TOO LARGE, ITMAX TOO SMALL' +1 GAMMCF=DEXP(-X+A*DLOG(X)-GLN)*G + RETURN + END + SUBROUTINE GSER(GAMSER,A,X,GLN) + implicit double precision(a-h,o-z) + PARAMETER (ITMAX=100,EPS=3.D-7) + GLN=GAMMLN(A) + IF(X.LE.0.d0)THEN + IF(X.LT.0.d0)PAUSE + GAMSER=0.d0 + RETURN + ENDIF + AP=A + SUM=1.d0/A + DEL=SUM + DO 11 N=1,ITMAX + AP=AP+1.d0 + DEL=DEL*X/AP + SUM=SUM+DEL + IF(DABS(DEL).LT.DABS(SUM)*EPS)GO TO 1 +11 CONTINUE + PAUSE 'A TOO LARGE, ITMAX TOO SMALL' +1 GAMSER=SUM*DEXP(-X+A*DLOG(X)-GLN) + RETURN + END + + + double precision function coef_nk(n,k) + implicit none + integer n,k + double precision gam,dblefact,fact + gam=dblefact(2*(n+k)+1) + coef_nk=1.d0/(2.d0**k*fact(k)*gam) + end + +!! Calculation of +!! +!! I= \int dx x**l *exp(-gam*x**2) M_n(ax) M_m(bx) +!! +!! M_n(x) modified spherical bessel function +!! + double precision function int_prod_bessel(l,gam,n,m,a,b) + implicit none + integer n,k,m,q,l,kcp + double precision gam,dblefact,fact,pi,a,b + double precision int,intold,sum,coef_nk,crochet + logical done + + if(a.ne.0.d0.and.b.ne.0.d0)then + q=0 + intold=-1.d0 + int=0.d0 + done=.false. + kcp=0 + do while (.not.done) + kcp=kcp+1 + sum=0.d0 + do k=0,q + sum=sum+coef_nk(n,k)*coef_nk(m,q-k)*a**(n+2*k)*b**(m-2*k+2*q) + enddo + int=int+sum*crochet(2*q+n+m+l,gam) + if(dabs(int-intold).lt.1d-15)then + done=.true. + else + q=q+1 + intold=int + endif + enddo + int_prod_bessel=int + if(kcp.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel' + return + endif + + if(a.eq.0.d0.and.b.ne.0.d0)then + if(n.ne.0)then + int_prod_bessel=0.d0 + return + endif + q=0 + intold=-1.d0 + int=0.d0 + done=.false. + kcp=0 + do while (.not.done) + kcp=kcp+1 + int=int+coef_nk(m,q)*b**(m+2*q)*crochet(2*q+m+l,gam) + if(dabs(int-intold).lt.1d-15)then + done=.true. + else + q=q+1 + intold=int + endif + enddo + int_prod_bessel=int + if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel' + return + endif + + if(a.ne.0.d0.and.b.eq.0.d0)then + if(m.ne.0)then + int_prod_bessel=0.d0 + return + endif + q=0 + intold=-1.d0 + int=0.d0 + done=.false. + kcp=0 + do while (.not.done) + kcp=kcp+1 + int=int+coef_nk(n,q)*a**(n+2*q)*crochet(2*q+n+l,gam) + if(dabs(int-intold).lt.1d-15)then + done=.true. + else + q=q+1 + intold=int + endif + enddo + int_prod_bessel=int + if(kcp.gt.100)stop '**WARNING** bad convergence in int_prod_bessel' + return + endif + + if(a.eq.0.d0.and.b.eq.0.d0)then + if(n.ne.0.or.m.ne.0)then + int_prod_bessel=0.d0 + return + endif + int_prod_bessel=crochet(l,gam) + return + endif + + stop 'pb in int_prod_bessel!!' + end + +!! Calculation of +!! +!! I= \int dx x**l *exp(-gam*x**2) M_n(ax) +!! +!! M_n(x) modified spherical bessel function +!! + double precision function int_prod_bessel_loc(l,gam,n,a) + implicit none + integer n,k,l,kcp + double precision gam,a + double precision int,intold,coef_nk,crochet + logical done + k=0 + intold=-1.d0 + int=0.d0 + done=.false. + kcp=0 + do while (.not.done) + kcp=kcp+1 + int=int+coef_nk(n,k)*a**(n+2*k)*crochet(2*k+n+l,gam) + if(dabs(int-intold).lt.1d-15)then + done=.true. + else + k=k+1 + intold=int + endif + enddo + int_prod_bessel_loc=int + if(kcp.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel' + end + + double precision function int_prod_bessel_num(l,gam,n,m,a,b) + implicit none + integer n,m,l,i,npoints + double precision gam,a,b + double precision sum,dx,x,bessel_mod + sum=0.d0 + npoints=20000 + dx=30.d0/npoints + do i=1,npoints + x=(i-1)*dx+0.5d0*dx + sum=sum+dx*x**l*dexp(-gam*x**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m) + enddo + int_prod_bessel_num=sum + end + + + From 37e049ea1a745123144586473c5779f5b1a70fff Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 3 Apr 2015 14:12:51 +0200 Subject: [PATCH 02/27] Add ezfio support --- scripts/pseudo/put_pseudo_in_ezfio.py | 183 +++++++++++++++++++++++++ setup_environment.sh | 4 +- src/Pseudo/NEEDED_MODULES | 1 + src/Pseudo/int.f90 | 66 --------- src/Pseudo/pseudo.ezfio_config | 10 ++ src/Pseudo/test_michel.irp.f | 188 ++++++++++++++++++++++++++ 6 files changed, 385 insertions(+), 67 deletions(-) create mode 100755 scripts/pseudo/put_pseudo_in_ezfio.py create mode 100644 src/Pseudo/pseudo.ezfio_config create mode 100644 src/Pseudo/test_michel.irp.f diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py new file mode 100755 index 00000000..792df9de --- /dev/null +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -0,0 +1,183 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Create the pseudo potential for a given atom + +Usage: + put_pseudo_in_ezfio.py --ezfio= --atom=... +""" + + +import os +import sys +from docopt import docopt + +from subprocess import Popen, PIPE + +qpackage_root = os.environ['QPACKAGE_ROOT'] + +EZFIO = "{0}/EZFIO".format(qpackage_root) +sys.path = [EZFIO + "/Python"] + sys.path + +from ezfio import ezfio + +import re +p = re.compile(ur'\|(\d+)><\d+\|') + + +def get_pseudo_str(l_atom): + """ + Run EMSL_local for geting the str of the speudo potential + """ + + EMSL_root = "{0}/EMSL_Basis/".format(qpackage_root) + EMSL_path = "{0}/EMSL_api.py".format(EMSL_root) + db_path = "{0}/db/Pseudo.db".format(EMSL_root) + + l_cmd_atom = [] + for a in l_atom: + l_cmd_atom += ["--atom", a] + + l_cmd_head = [EMSL_path, "get_basis_data", + "--db_path", db_path, + "--basis", "BFD-Pseudo"] + + process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE) + + stdout, _ = process.communicate() + return stdout.strip() + + +def get_v_n_dz_local(str_ele): + """ + From a str_ele of the pseudo (aka only one ele in the str) + get the list ussefull for the Local potential : v_k n_k and dz_k + """ + l_v_k = [] + l_n_k = [] + l_dz_k = [] + + for l in str_ele.splitlines(): + try: + v, n, dz = l.split() + v = float(v) + n = int(n) + dz = float(dz) + except ValueError: + pass + else: + l_v_k.append(v) + l_n_k.append(n) + l_dz_k.append(dz) + + return l_v_k, l_n_k, l_dz_k + + +def get_v_n_dz_l_nonlocal(str_ele): + """ + From a str_ele of the pseudo (aka only one ele in the str) + get the list ussefull for the non Local potential + v_kl (v, l) + n_k (v, l) + dz_k (dz ,l) + """ + l_v_kl = [] + l_n_kl = [] + l_dz_kl = [] + + for l in str_ele.splitlines(): + try: + v, n, dz, proj = l.split() + v = float(v) + n = int(n) + dz = float(dz) + l = int(p.match(proj).group(1)) + + except ValueError: + pass + else: + l_v_kl.append((v, l)) + l_n_kl.append((n, l)) + l_dz_kl.append((dz, l)) + + return l_v_kl, l_n_kl, l_dz_kl + + +if __name__ == "__main__": + arguments = docopt(__doc__) + + # ___ + # | ._ o _|_ + # _|_ | | | |_ + # + + # ~#~#~#~#~ # + # E Z F I O # + # ~#~#~#~#~ # + + ezfio_path = arguments["--ezfio"] + ezfio_path = os.path.expanduser(ezfio_path) + ezfio_path = os.path.expandvars(ezfio_path) + ezfio_path = os.path.abspath(ezfio_path) + + ezfio.set_file("{0}".format(ezfio_path)) + + # ~#~#~#~#~#~#~#~#~#~#~ # + # P s e u d o _ d a t a # + # ~#~#~#~#~#~#~#~#~#~#~ # + + l_ele = arguments["--atom"] + str_ = get_pseudo_str(l_ele) + + # _ + # |_) _. ._ _ _ + # | (_| | _> (/_ + # + + l_str_ele = [str_ele for str_ele in str_.split("Element Symbol: ") + if str_ele] + + for str_ele in l_str_ele: + + # ~#~#~#~#~ # + # S p l i t # + # ~#~#~#~#~ # + + l = str_ele.find("Local component:") + nl = str_ele.find("Non-local component") + + # ~#~#~#~#~ # + # L o c a l # + # ~#~#~#~#~ # + + print "local" + + l_v, l_n, l_dz = get_v_n_dz_local(str_ele[l:nl]) + + print l_v + print l_n + print l_dz + + ezfio.pseudo_klocmax = len(l_v) + ezfio.pseudo_v_k = l_v + ezfio.pseudo_n_k = l_n + ezfio.pseudo_dz_k = l_dz + + # ~#~#~#~#~#~#~#~#~ # + # N o n _ L o c a l # + # ~#~#~#~#~#~#~#~#~ # + + print "non local" + + l_v_kl, l_n_kl, l_dz_kl = get_v_n_dz_l_nonlocal(str_ele[nl:]) + + print l_v_kl + print l_n_kl + print l_dz_kl + + if l_v_kl: + ezfio.pseudo_lmax = max([i[1] for i in l_v_kl]) + 1 + ezfio.pseudo_kmax = len(l_v_kl) + ezfio.pseudo_v_kl = l_v_kl + ezfio.pseudo_n_kl = l_n_kl + ezfio.pseudo_dz_kl = l_dz_kl diff --git a/setup_environment.sh b/setup_environment.sh index c3dc4194..59f3af70 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -17,7 +17,9 @@ export QPACKAGE_ROOT=\$( cd \$(dirname "\${BASH_SOURCE}") ; pwd -P ) export LD_LIBRARY_PATH="\${QPACKAGE_ROOT}"/lib:\${LD_LIBRARY_PATH} export LIBRARY_PATH="\${QPACKAGE_ROOT}"/lib:\${LIBRARY_PATH} export C_INCLUDE_PATH="\${QPACKAGE_ROOT}"/include:\${C_INCLUDE_PATH} -export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts:"\${QPACKAGE_ROOT}"/scripts/ezfio_interface +export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts +export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts/ezfio_interface +export PYTHONPATH=\${PYTHONPATH}:"\${QPACKAGE_ROOT}"/scripts/pseudo export PATH=\${PATH}:"\${QPACKAGE_ROOT}"/scripts:"\${QPACKAGE_ROOT}"/scripts/ezfio_interface export PATH=\${PATH}:"\${QPACKAGE_ROOT}"/bin export PATH=\${PATH}:"\${QPACKAGE_ROOT}"/ocaml diff --git a/src/Pseudo/NEEDED_MODULES b/src/Pseudo/NEEDED_MODULES index e69de29b..542760a9 100644 --- a/src/Pseudo/NEEDED_MODULES +++ b/src/Pseudo/NEEDED_MODULES @@ -0,0 +1 @@ +Ezfio_files Nuclei Output Utils \ No newline at end of file diff --git a/src/Pseudo/int.f90 b/src/Pseudo/int.f90 index 16aca8e0..b67e6102 100644 --- a/src/Pseudo/int.f90 +++ b/src/Pseudo/int.f90 @@ -1,60 +1,3 @@ -!! -!! Computation of Vps, matrix element of the -!! pseudo-potential centered at point C -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!! Vps= < Phi_A | Vloc(C) + Vpp(C) | Phi_B> -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -!! Phi_M (M=A,B) Cartesian gaussian orbital centered at point M : -!! Phi_M = (x-M_x)**n^M_x *(y-M_y)**n^M_y *(z-M_z)**n^M_z exp(-g_M rM**2) -!! with rM**2=(x-M_x)**2 + (y-M_y)**2 + (z-M_z)**2 -!! -!!** Vloc(C)= \sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2) -!! -!!** Vpp(C)= \sum_{l=0}^lmax v_l(rC) \sum_{m=-l}^{m=l} |Y_lm> : -!! function Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) -!! lmax of formula above -!! kmax of formula above -!! v_kl = array v_kl(kmax_max,0:lmax_max) -!! n_kl = array n_kl(kmax_max,0:lmax_max) -!! dz_kl = array dz_kl(kmax_max,0:lmax_max) -!! n_a(1),n_a(2),n_a(3) -!! a(1),a(2),a(3) -!! g_a -!! n_b(1),n_b(2),n_b(3) -!! b(1),b(2),b(3) -!! g_b -!! c(1),c(2),c(3) -!! -!! Routine computing : -!! function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) -!! klocmax of formula above -!! v_k = array v_k(klocmax_max) -!! n_k = array n_k(klocmax_max) -!! dz_k= array dz_k(klocmax_max) -!! Routine total matrix element : -!! function Vps(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) -!! -!! Routines Vps_num, Vpseudo_num, and Vloc_num = brute force numerical -!! estimations of the same integrals - - !! Vps= !! !! with: Vloc(C)=\sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2) @@ -1019,15 +962,6 @@ end ! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 ! P_l^|m|(cos(theta)) exp(i m phi) - double precision function fact(n) - implicit double precision(a-h,o-z) - fact=1.d0 - if(n.eq.0)return - do i=1,n - fact=fact*dfloat(i) - enddo - end - subroutine erreur(x,n,rmoy,error) implicit double precision(a-h,o-z) dimension x(n) diff --git a/src/Pseudo/pseudo.ezfio_config b/src/Pseudo/pseudo.ezfio_config new file mode 100644 index 00000000..1ed75773 --- /dev/null +++ b/src/Pseudo/pseudo.ezfio_config @@ -0,0 +1,10 @@ +pseudo + klocmax integer + v_k double precision (pseudo_klocmax) + n_k integer (pseudo_klocmax) + dz_k double precision (pseudo_klocmax) + lmax integer + kmax integer + v_kl double precision (pseudo_kmax,pseudo_lmax) + n_kl integer (pseudo_kmax,pseudo_lmax) + dz_kl double precision (pseudo_kmax,pseudo_lmax) \ No newline at end of file diff --git a/src/Pseudo/test_michel.irp.f b/src/Pseudo/test_michel.irp.f new file mode 100644 index 00000000..e290ea4a --- /dev/null +++ b/src/Pseudo/test_michel.irp.f @@ -0,0 +1,188 @@ +!! +!! Computation of Vps, matrix element of the +!! pseudo-potential centered at point C +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! Vps= < Phi_A | Vloc(C) + Vpp(C) | Phi_B> +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! Phi_M (M=A,B) Cartesian gaussian orbital centered at point M : +!! Phi_M = (x-M_x)**n^M_x *(y-M_y)**n^M_y *(z-M_z)**n^M_z exp(-g_M rM**2) +!! with rM**2=(x-M_x)**2 + (y-M_y)**2 + (z-M_z)**2 +!! +!!** Vloc(C)= \sum_{k=1}^klocmax v_k rC**n_k exp(-dz_k rC**2) +!! +!!** Vpp(C)= \sum_{l=0}^lmax v_l(rC) \sum_{m=-l}^{m=l} |Y_lm> : +!! function Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +!! lmax of formula above +!! kmax of formula above +!! v_kl = array v_kl(kmax_max,0:lmax_max) +!! n_kl = array n_kl(kmax_max,0:lmax_max) +!! dz_kl = array dz_kl(kmax_max,0:lmax_max) +!! n_a(1),n_a(2),n_a(3) +!! a(1),a(2),a(3) +!! g_a +!! n_b(1),n_b(2),n_b(3) +!! b(1),b(2),b(3) +!! g_b +!! c(1),c(2),c(3) +!! +!! Routine computing : +!! function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) +!! klocmax of formula above +!! v_k = array v_k(klocmax_max) +!! n_k = array n_k(klocmax_max) +!! dz_k= array dz_k(klocmax_max) +!! Routine total matrix element : +!! function Vps(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) +!! +!! Routines Vps_num, Vpseudo_num, and Vloc_num = brute force numerical +!! estimations of the same integrals + + +program compute_integrals_pseudo + implicit none + integer n_a(3),n_b(3),npts + double precision g_a,g_b,a(3),b(3),c(3) + double precision Vpseudo,Vpseudo_num,Vloc,Vloc_num + double precision v3,v4 + + + double precision vps,vps_num + + ! PSEUDOS + integer nptsgridmax,nptsgrid + double precision coefs_pseudo,ptsgrid + + double precision rmax + double precision time_1,time_2,time_3,time_4,time_5 + integer kga,kgb,na1,na2,na3,nb1,nb2,nb3 + + CALL RANDOM_SEED() + + nptsgrid=50 + call initpseudos(nptsgrid) + + PROVIDE ezfio_filename + + ! + ! | _ _ _. | + ! |_ (_) (_ (_| | + ! + + integer klocmax + integer, allocatable :: n_k(:) + double precision, allocatable :: v_k(:), dz_k(:) + + call ezfio_get_pseudo_klocmax(klocmax) + + allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax)) + + call ezfio_get_pseudo_v_k(v_k) + call ezfio_get_pseudo_n_k(n_k) + call ezfio_get_pseudo_dz_k(dz_k) + + print*, "klocmax", klocmax + + print*, "n_k_ezfio", n_k + print*, "v_k_ezfio",v_k + print*, "dz_k_ezfio", dz_k + + ! + ! |\ | _ ._ | _ _ _. | + ! | \| (_) | | | (_) (_ (_| | + ! + + !! Parameters of non local part of pseudo: + + integer :: kmax,lmax + integer, allocatable :: n_kl(:,:) + double precision, allocatable :: v_kl(:,:), dz_kl(:,:) + + call ezfio_get_pseudo_lmax(lmax) + call ezfio_get_pseudo_kmax(kmax) + lmax = lmax - 1 + + allocate(n_kl(kmax,0:lmax), v_kl(kmax,0:lmax), dz_kl(kmax,0:lmax)) + + call ezfio_get_pseudo_n_kl(n_kl) + call ezfio_get_pseudo_v_kl(v_kl) + call ezfio_get_pseudo_dz_kl(dz_kl) + + + print*, "lmax",lmax + print*, "kmax", kmax + + print*,"n_kl_ezfio", n_kl + print*,"v_kl_ezfio", v_kl + print*,"dz_kl_ezfio", dz_kl + + ! _ + ! / _. | _ | + ! \_ (_| | (_ |_| | + ! + + write(*,*)'a?' + read*,a(1),a(2),a(3) + !write(*,*)'b?' + !read*,b(1),b(2),b(3) + b(1)=-0.1d0 + b(2)=-0.2d0 + b(3)=0.3d0 + !write(*,*)'a?' + !read*,c(1),c(2),c(3) + c(1)=0.1d0 + c(2)=0.2d0 + c(3)=0.3d0 + + print*,'ntps? rmax for brute force integration' + read*,npts,rmax + + do kga=0,5 + g_a=10.d0**kga + do kgb=0,5 + g_b=10.d0**kgb + + do na1=0,1 + do na2=0,1 + do na3=0,1 + do nb1=0,1 + do nb2=0,1 + do nb3=0,1 + n_a(1)=na1 + n_a(2)=na2 + n_a(3)=na3 + n_b(1)=nb1 + n_b(2)=nb2 + n_b(3)=nb3 + + v3=Vps(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) + v4=Vps_num(npts,rmax,a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) + print*,'Vps= ',v3,' Vps_num=',v4,' diff=',v4-v3 + write(33,'(3f10.6)')v3,v4,v4-v3 + + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + +end From 5240e59d40d6db7e6335bf5e0bcfc2ca7c6fab80 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 7 Apr 2015 08:38:55 +0200 Subject: [PATCH 03/27] Mv Pseudo in MonoInts --- src/MonoInts/Makefile | 6 +++--- src/{Pseudo => MonoInts}/int.f90 | 0 src/{Pseudo => MonoInts}/pseudo.ezfio_config | 0 src/{Pseudo => MonoInts}/test_michel.irp.f | 0 src/NEEDED_MODULES | 2 +- src/Pseudo/ASSUMPTIONS.rst | 0 src/Pseudo/Makefile | 6 ------ src/Pseudo/NEEDED_MODULES | 1 - src/Pseudo/README.rst | 4 ---- 9 files changed, 4 insertions(+), 15 deletions(-) rename src/{Pseudo => MonoInts}/int.f90 (100%) rename src/{Pseudo => MonoInts}/pseudo.ezfio_config (100%) rename src/{Pseudo => MonoInts}/test_michel.irp.f (100%) delete mode 100644 src/Pseudo/ASSUMPTIONS.rst delete mode 100644 src/Pseudo/Makefile delete mode 100644 src/Pseudo/NEEDED_MODULES delete mode 100644 src/Pseudo/README.rst diff --git a/src/MonoInts/Makefile b/src/MonoInts/Makefile index 06dc50ff..8ae5c9fb 100644 --- a/src/MonoInts/Makefile +++ b/src/MonoInts/Makefile @@ -1,6 +1,6 @@ # Define here all new external source files and objects.Don't forget to prefix the # object files with IRPF90_temp/ -SRC= -OBJ= +SRC=int.f90 +OBJ=IRPF90_temp/int.o -include $(QPACKAGE_ROOT)/src/Makefile.common +include $(QPACKAGE_ROOT)/src/Makefile.common \ No newline at end of file diff --git a/src/Pseudo/int.f90 b/src/MonoInts/int.f90 similarity index 100% rename from src/Pseudo/int.f90 rename to src/MonoInts/int.f90 diff --git a/src/Pseudo/pseudo.ezfio_config b/src/MonoInts/pseudo.ezfio_config similarity index 100% rename from src/Pseudo/pseudo.ezfio_config rename to src/MonoInts/pseudo.ezfio_config diff --git a/src/Pseudo/test_michel.irp.f b/src/MonoInts/test_michel.irp.f similarity index 100% rename from src/Pseudo/test_michel.irp.f rename to src/MonoInts/test_michel.irp.f diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 35aa5ec3..a144c42c 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected MRCC Pseudo +AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected MRCC diff --git a/src/Pseudo/ASSUMPTIONS.rst b/src/Pseudo/ASSUMPTIONS.rst deleted file mode 100644 index e69de29b..00000000 diff --git a/src/Pseudo/Makefile b/src/Pseudo/Makefile deleted file mode 100644 index 5cf11b78..00000000 --- a/src/Pseudo/Makefile +++ /dev/null @@ -1,6 +0,0 @@ -# Define here all new external source files and objects.Don't forget to prefix the -# object files with IRPF90_temp/ -SRC=int.f90 -OBJ=IRPF90_temp/int.o - -include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/Pseudo/NEEDED_MODULES b/src/Pseudo/NEEDED_MODULES deleted file mode 100644 index 542760a9..00000000 --- a/src/Pseudo/NEEDED_MODULES +++ /dev/null @@ -1 +0,0 @@ -Ezfio_files Nuclei Output Utils \ No newline at end of file diff --git a/src/Pseudo/README.rst b/src/Pseudo/README.rst deleted file mode 100644 index ebb762c1..00000000 --- a/src/Pseudo/README.rst +++ /dev/null @@ -1,4 +0,0 @@ -======= - Module -======= - From ef65e3e51362dc5454f3b82cee5e611b0de6f41e Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 7 Apr 2015 10:27:22 +0200 Subject: [PATCH 04/27] New int.f90 working with A=B=C --- src/MonoInts/int.f90 | 12 +++++++++- src/MonoInts/pot_ao_ints.irp.f | 42 ++++++++++++++++++++++++++++++---- src/MonoInts/test_michel.irp.f | 32 ++++++++++++++++++-------- 3 files changed, 70 insertions(+), 16 deletions(-) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index b67e6102..640688d0 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -560,12 +560,22 @@ double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg Vloc=0.d0 return endif - freal=dexp(-g_a*ac**2-g_b*bc**2) ntotA=n_a(1)+n_a(2)+n_a(3) ntotB=n_b(1)+n_b(2)+n_b(3) ntot=ntotA+ntotB + if(ac.eq.0.d0.and.bc.eq.0.d0)then + accu=0.d0 + do k=1,klocmax + accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k)) + enddo + Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3)) + return + endif + + freal=dexp(-g_a*ac**2-g_b*bc**2) + d2=0.d0 do i=1,3 d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i)) diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index f430ace9..6ad38d6a 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -8,19 +8,43 @@ double precision :: A_center(3),B_center(3),C_center(3) integer :: power_A(3),power_B(3) integer :: i,j,k,l,n_pt_in,m - double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc integer :: nucl_numC ! Important for OpenMP ao_nucl_elec_integral = 0.d0 + ! + ! | _ _ _. | + ! |_ (_) (_ (_| | + ! + + integer klocmax + integer, allocatable :: n_k(:) + double precision, allocatable :: v_k(:), dz_k(:) + + call ezfio_get_pseudo_klocmax(klocmax) + + allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax)) + + call ezfio_get_pseudo_v_k(v_k) + call ezfio_get_pseudo_n_k(n_k) + call ezfio_get_pseudo_dz_k(dz_k) + + print*, "klocmax", klocmax + + print*, "n_k_ezfio", n_k + print*, "v_k_ezfio",v_k + print*, "dz_k_ezfio", dz_k + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, & - !$OMP num_A,num_B,Z,c,n_pt_in) & + !$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, & + !$OMP num_A, num_B, Z, c, n_pt_in,& + !$OMP v_k, n_k, dz_k, klocmax) & !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & - !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge) + !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge) n_pt_in = n_pt_max_integrals !$OMP DO SCHEDULE (guided) do j = 1, ao_num @@ -50,7 +74,15 @@ C_center(1) = nucl_coord(k,1) C_center(2) = nucl_coord(k,2) C_center(3) = nucl_coord(k,3) - c = c+Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) + c = c + Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) + + print*, A_center + print*, B_center + print*, C_center + + print*, Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center) +! c = c + Z*Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center) + enddo ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - & ao_coef_transp(l,j)*ao_coef_transp(m,i)*c diff --git a/src/MonoInts/test_michel.irp.f b/src/MonoInts/test_michel.irp.f index e290ea4a..c7e3c685 100644 --- a/src/MonoInts/test_michel.irp.f +++ b/src/MonoInts/test_michel.irp.f @@ -137,19 +137,31 @@ program compute_integrals_pseudo ! \_ (_| | (_ |_| | ! - write(*,*)'a?' - read*,a(1),a(2),a(3) +! write(*,*)'a?' +! read*,a(1),a(2),a(3) !write(*,*)'b?' !read*,b(1),b(2),b(3) - b(1)=-0.1d0 - b(2)=-0.2d0 - b(3)=0.3d0 - !write(*,*)'a?' - !read*,c(1),c(2),c(3) - c(1)=0.1d0 - c(2)=0.2d0 - c(3)=0.3d0 +! b(1)=-0.1d0 +! b(2)=-0.2d0 +! b(3)=0.3d0 +! !write(*,*)'a?' +! !read*,c(1),c(2),c(3) +! c(1)=0.1d0 +! c(2)=0.2d0 +! c(3)=0.3d0 + a(1)= 0.d0 + a(2)= 0.d0 + a(3)= 0.d0 + + b(1)= 0.d0 + b(2)= 0.d0 + b(3)= 0.d0 + + c(1)= 0.d0 + c(2)= 0.d0 + c(3)= 0.d0 + print*,'ntps? rmax for brute force integration' read*,npts,rmax From cca5ebc404ab3c6689a677b895afbe0b0ee80329 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 7 Apr 2015 14:24:05 +0200 Subject: [PATCH 05/27] Move need.irp.f into MonoInts --- src/MonoInts/need.irp.f | 289 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 289 insertions(+) create mode 100644 src/MonoInts/need.irp.f diff --git a/src/MonoInts/need.irp.f b/src/MonoInts/need.irp.f new file mode 100644 index 00000000..22cb6a48 --- /dev/null +++ b/src/MonoInts/need.irp.f @@ -0,0 +1,289 @@ + + double precision function SABpartial(zA,zB,A,B,nA,nB,gamA,gamB) + implicit double precision(a-h,o-z) + dimension nA(3),nB(3) + dimension A(3),B(3) + gamtot=gamA+gamB + SABpartial=1.d0 + + l=3 + u=gamA/gamtot*A(l)+gamB/gamtot*B(l) + arg=gamtot*u**2-gamA*A(l)**2-gamB*B(l)**2 + alpha=dexp(arg) + &/gamtot**((1.d0+dfloat(nA(l))+dfloat(nB(l)))/2.d0) + wA=dsqrt(gamtot)*(u-A(l)) + wB=dsqrt(gamtot)*(u-B(l)) + boundA=dsqrt(gamtot)*(zA-u) + boundB=dsqrt(gamtot)*(zB-u) + + accu=0.d0 + do n=0,nA(l) + do m=0,nB(l) + integ=nA(l)+nB(l)-n-m + accu=accu + & +wA**n*wB**m*binom(nA(l),n)*binom(nB(l),m) + & *(rinteg(integ,boundB)-rinteg(integ,boundA)) + enddo + enddo + SABpartial=SABpartial*accu*alpha + end + + double precision function rintgauss(n) + implicit double precision(a-h,o-z) + rintgauss=dsqrt(dacos(-1.d0)) + if(n.eq.0)return + if(n.eq.1)then + rintgauss=0.d0 + return + endif + if(iand(n,1).eq.1)then + rintgauss=0.d0 + return + endif + rintgauss=rintgauss/2.d0**(n/2) + rintgauss=rintgauss*ddfact2(n-1) + end + + double precision function rinteg(n,u) + implicit double precision(a-h,o-z) + include 'constants.F' +! pi=dacos(-1.d0) + ichange=1 + factor=1.d0 + if(u.lt.0.d0)then + u=-u + factor=(-1.d0)**(n+1) + ichange=-1 + endif + if(iand(n,1).eq.0)then + rinteg=0.d0 + do l=0,n-2,2 + prod=b_coef(l,u) + do k=l+2,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + enddo + prod=dsqrt(pi)/2.d0*erf0(u) + do k=0,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + endif + + if(iand(n,1).eq.1)then + rinteg=0.d0 + do l=1,n-2,2 + prod=b_coef(l,u) + do k=l+2,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + enddo + prod=0.5d0*(1.d0-dexp(-u**2)) + do k=1,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + endif + + rinteg=rinteg*factor + + if(ichange.eq.-1)u=-u + + end + +! +! +! +! +! +! +! +! + double precision function erf0(x) + implicit double precision (a-h,o-z) + if(x.lt.0.d0)then + erf0=-gammp(0.5d0,x**2) + else + erf0=gammp(0.5d0,x**2) + endif + end + + +! +! +! +! +! +! +! +! +! +! +! +! gcf +! gser +! +! +! + double precision function gammp(a,x) + implicit double precision (a-h,o-z) + if(x.lt.0..or.a.le.0.)stop 'error in gammp' + if(x.lt.a+1.)then + call gser(gammp,a,x,gln) + else + call gcf(gammcf,a,x,gln) + gammp=1.-gammcf + endif + return + end +! +! + + +! +! +! +! +! +! +! +! +! +! +! +! gammp +! +! +! + subroutine gser(gamser,a,x,gln) + implicit double precision (a-h,o-z) + parameter (itmax=100,eps=3.e-7) + gln=gammln(a) + if(x.le.0.)then + if(x.lt.0.) stop 'error in gser' + gamser=0. + return + endif + ap=a + sum=1./a + del=sum + do 11 n=1,itmax + ap=ap+1. + del=del*x/ap + sum=sum+del + if(abs(del).lt.abs(sum)*eps)go to 1 +11 continue + stop 'a too large, itmax too small' +1 gamser=sum*exp(-x+a*log(x)-gln) + return + end +! + +! +! +! +! +! +! +! +! +! +! +! +! +! gammp +! +! +! + subroutine gcf(gammcf,a,x,gln) + implicit double precision (a-h,o-z) + parameter (itmax=100,eps=3.e-7) + gln=gammln(a) + gold=0. + a0=1. + a1=x + b0=0. + b1=1. + fac=1. + do 11 n=1,itmax + an=float(n) + ana=an-a + a0=(a1+a0*ana)*fac + b0=(b1+b0*ana)*fac + anf=an*fac + a1=x*a0+anf*a1 + b1=x*b0+anf*b1 + if(a1.ne.0.)then + fac=1./a1 + g=b1*fac + if(abs((g-gold)/g).lt.eps)go to 1 + gold=g + endif +11 continue + stop 'a too large, itmax too small' +1 gammcf=exp(-x+a*log(x)-gln)*g + return + end + +! +! + double precision function ddfact2(n) + implicit double precision(a-h,o-z) + if(iand(n,1).eq.0)stop 'error in ddfact2' + ddfact2=1.d0 + do i=1,n,2 + ddfact2=ddfact2*dfloat(i) + enddo + end + + double precision function a_coef(n) + implicit double precision(a-h,o-z) + a_coef=dfloat(n+1)/2.d0 + end + + double precision function b_coef(n,u) + implicit double precision(a-h,o-z) + b_coef=-0.5d0*u**(n+1)*dexp(-u**2) + end + +! +! +! +! +! +! +! +! + double precision function gammln(xx) + implicit double precision (a-h,o-z) + real*8 cof(6),stp,half,one,fpf,x,tmp,ser + data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0, + * -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/ + data half,one,fpf/0.5d0,1.0d0,5.5d0/ + x=xx-one + tmp=x+fpf + tmp=(x+half)*log(tmp)-tmp + ser=one + do 11 j=1,6 + x=x+one + ser=ser+cof(j)/x +11 continue + gammln=tmp+log(stp*ser) + return + end +! +! From 64ea60c7277af739abeed790ab9d23b032535bff Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 7 Apr 2015 16:59:42 +0200 Subject: [PATCH 06/27] Working Local --- scripts/get_basis.sh | 6 +- src/MonoInts/int.f90 | 83 +--------- src/MonoInts/pot_ao_ints.irp.f | 63 +++++-- src/Properties/need.irp.f | 289 --------------------------------- 4 files changed, 54 insertions(+), 387 deletions(-) delete mode 100644 src/Properties/need.irp.f diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index 51b0a4f0..b5be03ac 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -42,9 +42,9 @@ then echo "ERROR" exit 1 fi -${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" $atoms - - +#${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" $atoms +cp He.dz_filipi.basis ${tmpfile} +echo ${tmpfile} diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index 640688d0..dc07fa54 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -550,7 +550,6 @@ double precision,allocatable :: array_R_loc(:,:,:) double precision,allocatable :: array_coefs(:,:,:,:,:,:) double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg - fourpi=4.d0*dacos(-1.d0) f=fourpi**1.5d0 ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) @@ -567,6 +566,7 @@ double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg if(ac.eq.0.d0.and.bc.eq.0.d0)then accu=0.d0 + do k=1,klocmax accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k)) enddo @@ -1727,87 +1727,6 @@ end RETURN END - double precision FUNCTION GAMMLN(XX) - implicit double precision(a-h,o-z) - REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER - DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, & - -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ - DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ - X=XX-ONE - TMP=X+FPF - TMP=(X+HALF)*DLOG(TMP)-TMP - SER=ONE - DO 11 J=1,6 - X=X+ONE - SER=SER+COF(J)/X -11 CONTINUE - GAMMLN=TMP+DLOG(STP*SER) - RETURN - END - FUNCTION GAMMP(A,X) - implicit double precision(a-h,o-z) - IF(X.LT.0.d0.OR.A.LE.0.d0)PAUSE - IF(X.LT.A+1.d0)THEN - CALL GSER(GAMMP,A,X,GLN) - ELSE - CALL GCF(GAMMCF,A,X,GLN) - GAMMP=1.d0-GAMMCF - ENDIF - RETURN - END - SUBROUTINE GCF(GAMMCF,A,X,GLN) - implicit double precision(a-h,o-z) - PARAMETER (ITMAX=100,EPS=3.D-7) - GLN=GAMMLN(A) - GOLD=0.d0 - A0=1.d0 - A1=X - B0=0.d0 - B1=1.d0 - FAC=1.d0 - DO 11 N=1,ITMAX - AN=DFLOAT(N) - ANA=AN-A - A0=(A1+A0*ANA)*FAC - B0=(B1+B0*ANA)*FAC - ANF=AN*FAC - A1=X*A0+ANF*A1 - B1=X*B0+ANF*B1 - IF(A1.NE.0.d0)THEN - FAC=1.d0/A1 - G=B1*FAC - IF(DABS((G-GOLD)/G).LT.EPS)GO TO 1 - GOLD=G - ENDIF -11 CONTINUE - PAUSE 'A TOO LARGE, ITMAX TOO SMALL' -1 GAMMCF=DEXP(-X+A*DLOG(X)-GLN)*G - RETURN - END - SUBROUTINE GSER(GAMSER,A,X,GLN) - implicit double precision(a-h,o-z) - PARAMETER (ITMAX=100,EPS=3.D-7) - GLN=GAMMLN(A) - IF(X.LE.0.d0)THEN - IF(X.LT.0.d0)PAUSE - GAMSER=0.d0 - RETURN - ENDIF - AP=A - SUM=1.d0/A - DEL=SUM - DO 11 N=1,ITMAX - AP=AP+1.d0 - DEL=DEL*X/AP - SUM=SUM+DEL - IF(DABS(DEL).LT.DABS(SUM)*EPS)GO TO 1 -11 CONTINUE - PAUSE 'A TOO LARGE, ITMAX TOO SMALL' -1 GAMSER=SUM*DEXP(-X+A*DLOG(X)-GLN) - RETURN - END - - double precision function coef_nk(n,k) implicit none integer n,k diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index 6ad38d6a..7fd24174 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -8,7 +8,7 @@ double precision :: A_center(3),B_center(3),C_center(3) integer :: power_A(3),power_B(3) integer :: i,j,k,l,n_pt_in,m - double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc + double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc, Vpseudo integer :: nucl_numC ! Important for OpenMP @@ -38,13 +38,55 @@ print*, "dz_k_ezfio", dz_k + call ezfio_get_pseudo_v_k(v_k) + call ezfio_get_pseudo_n_k(n_k) + call ezfio_get_pseudo_dz_k(dz_k) + + print*, "klocmax", klocmax + + print*, "n_k_ezfio", n_k + print*, "v_k_ezfio",v_k + print*, "dz_k_ezfio", dz_k + + + ! + ! |\ | _ ._ | _ _ _. | + ! | \| (_) | | | (_) (_ (_| | + ! + + !! Parameters of non local part of pseudo: + + integer :: kmax,lmax + integer, allocatable :: n_kl(:,:) + double precision, allocatable :: v_kl(:,:), dz_kl(:,:) + + call ezfio_get_pseudo_lmax(lmax) + call ezfio_get_pseudo_kmax(kmax) + lmax = lmax - 1 + + allocate(n_kl(kmax,0:lmax), v_kl(kmax,0:lmax), dz_kl(kmax,0:lmax)) + + call ezfio_get_pseudo_n_kl(n_kl) + call ezfio_get_pseudo_v_kl(v_kl) + call ezfio_get_pseudo_dz_kl(dz_kl) + + + print*, "lmax",lmax + print*, "kmax", kmax + + print*,"n_kl_ezfio", n_kl + print*,"v_kl_ezfio", v_kl + print*,"dz_kl_ezfio", dz_kl + + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, & - !$OMP num_A, num_B, Z, c, n_pt_in,& - !$OMP v_k, n_k, dz_k, klocmax) & + !$OMP num_A, num_B, Z, c, n_pt_in) & !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & - !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge) + !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge, & + !$OMP v_k, n_k, dz_k, klocmax, & + !$OMP lmax,kmax,v_kl,n_kl,dz_kl) n_pt_in = n_pt_max_integrals !$OMP DO SCHEDULE (guided) do j = 1, ao_num @@ -74,18 +116,13 @@ C_center(1) = nucl_coord(k,1) C_center(2) = nucl_coord(k,2) C_center(3) = nucl_coord(k,3) + c = c + Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - - print*, A_center - print*, B_center - print*, C_center - - print*, Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center) -! c = c + Z*Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center) - + c = c - Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center) + c = c - Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,C_center) enddo ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - & - ao_coef_transp(l,j)*ao_coef_transp(m,i)*c + ao_coef_transp(l,j)*ao_coef_transp(m,i)*c enddo enddo enddo diff --git a/src/Properties/need.irp.f b/src/Properties/need.irp.f deleted file mode 100644 index 22cb6a48..00000000 --- a/src/Properties/need.irp.f +++ /dev/null @@ -1,289 +0,0 @@ - - double precision function SABpartial(zA,zB,A,B,nA,nB,gamA,gamB) - implicit double precision(a-h,o-z) - dimension nA(3),nB(3) - dimension A(3),B(3) - gamtot=gamA+gamB - SABpartial=1.d0 - - l=3 - u=gamA/gamtot*A(l)+gamB/gamtot*B(l) - arg=gamtot*u**2-gamA*A(l)**2-gamB*B(l)**2 - alpha=dexp(arg) - &/gamtot**((1.d0+dfloat(nA(l))+dfloat(nB(l)))/2.d0) - wA=dsqrt(gamtot)*(u-A(l)) - wB=dsqrt(gamtot)*(u-B(l)) - boundA=dsqrt(gamtot)*(zA-u) - boundB=dsqrt(gamtot)*(zB-u) - - accu=0.d0 - do n=0,nA(l) - do m=0,nB(l) - integ=nA(l)+nB(l)-n-m - accu=accu - & +wA**n*wB**m*binom(nA(l),n)*binom(nB(l),m) - & *(rinteg(integ,boundB)-rinteg(integ,boundA)) - enddo - enddo - SABpartial=SABpartial*accu*alpha - end - - double precision function rintgauss(n) - implicit double precision(a-h,o-z) - rintgauss=dsqrt(dacos(-1.d0)) - if(n.eq.0)return - if(n.eq.1)then - rintgauss=0.d0 - return - endif - if(iand(n,1).eq.1)then - rintgauss=0.d0 - return - endif - rintgauss=rintgauss/2.d0**(n/2) - rintgauss=rintgauss*ddfact2(n-1) - end - - double precision function rinteg(n,u) - implicit double precision(a-h,o-z) - include 'constants.F' -! pi=dacos(-1.d0) - ichange=1 - factor=1.d0 - if(u.lt.0.d0)then - u=-u - factor=(-1.d0)**(n+1) - ichange=-1 - endif - if(iand(n,1).eq.0)then - rinteg=0.d0 - do l=0,n-2,2 - prod=b_coef(l,u) - do k=l+2,n-2,2 - prod=prod*a_coef(k) - enddo - rinteg=rinteg+prod - enddo - prod=dsqrt(pi)/2.d0*erf0(u) - do k=0,n-2,2 - prod=prod*a_coef(k) - enddo - rinteg=rinteg+prod - endif - - if(iand(n,1).eq.1)then - rinteg=0.d0 - do l=1,n-2,2 - prod=b_coef(l,u) - do k=l+2,n-2,2 - prod=prod*a_coef(k) - enddo - rinteg=rinteg+prod - enddo - prod=0.5d0*(1.d0-dexp(-u**2)) - do k=1,n-2,2 - prod=prod*a_coef(k) - enddo - rinteg=rinteg+prod - endif - - rinteg=rinteg*factor - - if(ichange.eq.-1)u=-u - - end - -! -! -! -! -! -! -! -! - double precision function erf0(x) - implicit double precision (a-h,o-z) - if(x.lt.0.d0)then - erf0=-gammp(0.5d0,x**2) - else - erf0=gammp(0.5d0,x**2) - endif - end - - -! -! -! -! -! -! -! -! -! -! -! -! gcf -! gser -! -! -! - double precision function gammp(a,x) - implicit double precision (a-h,o-z) - if(x.lt.0..or.a.le.0.)stop 'error in gammp' - if(x.lt.a+1.)then - call gser(gammp,a,x,gln) - else - call gcf(gammcf,a,x,gln) - gammp=1.-gammcf - endif - return - end -! -! - - -! -! -! -! -! -! -! -! -! -! -! -! gammp -! -! -! - subroutine gser(gamser,a,x,gln) - implicit double precision (a-h,o-z) - parameter (itmax=100,eps=3.e-7) - gln=gammln(a) - if(x.le.0.)then - if(x.lt.0.) stop 'error in gser' - gamser=0. - return - endif - ap=a - sum=1./a - del=sum - do 11 n=1,itmax - ap=ap+1. - del=del*x/ap - sum=sum+del - if(abs(del).lt.abs(sum)*eps)go to 1 -11 continue - stop 'a too large, itmax too small' -1 gamser=sum*exp(-x+a*log(x)-gln) - return - end -! - -! -! -! -! -! -! -! -! -! -! -! -! -! gammp -! -! -! - subroutine gcf(gammcf,a,x,gln) - implicit double precision (a-h,o-z) - parameter (itmax=100,eps=3.e-7) - gln=gammln(a) - gold=0. - a0=1. - a1=x - b0=0. - b1=1. - fac=1. - do 11 n=1,itmax - an=float(n) - ana=an-a - a0=(a1+a0*ana)*fac - b0=(b1+b0*ana)*fac - anf=an*fac - a1=x*a0+anf*a1 - b1=x*b0+anf*b1 - if(a1.ne.0.)then - fac=1./a1 - g=b1*fac - if(abs((g-gold)/g).lt.eps)go to 1 - gold=g - endif -11 continue - stop 'a too large, itmax too small' -1 gammcf=exp(-x+a*log(x)-gln)*g - return - end - -! -! - double precision function ddfact2(n) - implicit double precision(a-h,o-z) - if(iand(n,1).eq.0)stop 'error in ddfact2' - ddfact2=1.d0 - do i=1,n,2 - ddfact2=ddfact2*dfloat(i) - enddo - end - - double precision function a_coef(n) - implicit double precision(a-h,o-z) - a_coef=dfloat(n+1)/2.d0 - end - - double precision function b_coef(n,u) - implicit double precision(a-h,o-z) - b_coef=-0.5d0*u**(n+1)*dexp(-u**2) - end - -! -! -! -! -! -! -! -! - double precision function gammln(xx) - implicit double precision (a-h,o-z) - real*8 cof(6),stp,half,one,fpf,x,tmp,ser - data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0, - * -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/ - data half,one,fpf/0.5d0,1.0d0,5.5d0/ - x=xx-one - tmp=x+fpf - tmp=(x+half)*log(tmp)-tmp - ser=one - do 11 j=1,6 - x=x+one - ser=ser+cof(j)/x -11 continue - gammln=tmp+log(stp*ser) - return - end -! -! From 1cb79a58a41b2abe46a9df4855bf1f2db308fade Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 10 Apr 2015 09:43:28 +0200 Subject: [PATCH 07/27] Realy ugly (need to use create_ez and set alpha and beta), but working --- scripts/get_basis.sh | 2 +- scripts/pseudo/create_ez.sh | 14 +++ scripts/pseudo/put_pseudo_in_ezfio.py | 27 +++--- src/MonoInts/int.f90 | 15 ++-- src/MonoInts/pot_ao_ints.irp.f | 122 ++++++++++++++++---------- src/MonoInts/pseudo.ezfio_config | 8 +- src/MonoInts/test_michel.irp.f | 2 +- 7 files changed, 124 insertions(+), 66 deletions(-) create mode 100755 scripts/pseudo/create_ez.sh diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index b5be03ac..1d4e7472 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -45,6 +45,6 @@ fi #${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" $atoms -cp He.dz_filipi.basis ${tmpfile} +cp /home/razoa/quantum_package/scripts/pseudo/burkatzki_dz.basis ${tmpfile} echo ${tmpfile} diff --git a/scripts/pseudo/create_ez.sh b/scripts/pseudo/create_ez.sh new file mode 100755 index 00000000..3821d6c2 --- /dev/null +++ b/scripts/pseudo/create_ez.sh @@ -0,0 +1,14 @@ +#!/bin/bash +# $1 name +# $2 mult + +echo "name" $1 +echo "mul" $2 +echo "Zeff" $3 +echo "\`get_basis.sh\` need to be changed" + +rm -R $1.ezfio +qp_create_ezfio_from_xyz $1.xyz -b "cc-pvdz" -m $2 + +~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py --ezfio $1.ezfio/ --atom $1 --zeff $3 +qp_edit -c $1.ezfio \ No newline at end of file diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 792df9de..7d15f090 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -4,7 +4,7 @@ Create the pseudo potential for a given atom Usage: - put_pseudo_in_ezfio.py --ezfio= --atom=... + put_pseudo_in_ezfio.py --ezfio= --atom=... --zeff=... """ @@ -96,9 +96,14 @@ def get_v_n_dz_l_nonlocal(str_ele): except ValueError: pass else: - l_v_kl.append((v, l)) - l_n_kl.append((n, l)) - l_dz_kl.append((dz, l)) + l_v_kl.append([v]) + l_n_kl.append([n]) + l_dz_kl.append([dz]) + + if not l_v_kl: + l_v_kl.append([0.]) + l_n_kl.append([0]) + l_dz_kl.append([0.]) return l_v_kl, l_n_kl, l_dz_kl @@ -175,9 +180,11 @@ if __name__ == "__main__": print l_n_kl print l_dz_kl - if l_v_kl: - ezfio.pseudo_lmax = max([i[1] for i in l_v_kl]) + 1 - ezfio.pseudo_kmax = len(l_v_kl) - ezfio.pseudo_v_kl = l_v_kl - ezfio.pseudo_n_kl = l_n_kl - ezfio.pseudo_dz_kl = l_dz_kl + ezfio.pseudo_lmaxpo = len(l_v_kl) + ezfio.pseudo_kmax = len(l_v_kl[0]) + ezfio.pseudo_v_kl = l_v_kl + ezfio.pseudo_n_kl = l_n_kl + ezfio.pseudo_dz_kl = l_dz_kl + + if arguments["--zeff"]: + ezfio.nuclei_nucl_charge = map(int, arguments["--zeff"]) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index dc07fa54..cb95c3c9 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -8,7 +8,7 @@ implicit none integer n_a(3),n_b(3) double precision g_a,g_b,a(3),b(3),c(3) integer kmax_max,lmax_max -parameter (kmax_max=2,lmax_max=2) +parameter (kmax_max=4,lmax_max=2) integer lmax,kmax,n_kl(kmax_max,0:lmax_max) double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) integer klocmax_max @@ -18,7 +18,7 @@ double precision v_k(klocmax_max),dz_k(klocmax_max) double precision Vloc,Vpseudo Vps=Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) & - +Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) + +Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) end !! !! Vps_num: brute force numerical evaluation of the same matrix element Vps @@ -29,7 +29,7 @@ implicit none integer n_a(3),n_b(3) double precision g_a,g_b,a(3),b(3),c(3),rmax integer kmax_max,lmax_max -parameter (kmax_max=2,lmax_max=2) +parameter (kmax_max=4,lmax_max=2) integer lmax,kmax,n_kl(kmax_max,0:lmax_max) double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) integer klocmax_max;parameter (klocmax_max=10) @@ -170,12 +170,13 @@ end double precision function Vpseudo & (lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) implicit none +double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) + integer kmax_max,lmax_max,ntot_max,nkl_max -parameter (kmax_max=2,lmax_max=2,nkl_max=4) +parameter (kmax_max=4,lmax_max=2,nkl_max=4) parameter (ntot_max=10) integer lmax,kmax,n_kl(kmax_max,0:lmax_max),l,k double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) -double precision a(3),g_a,b(3),g_b,c(3) double precision fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm double precision theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big double precision areal,freal,breal,t1,t2,int_prod_bessel @@ -474,7 +475,7 @@ end double precision function Vpseudo_num(npts,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) implicit none integer kmax_max,lmax_max -parameter (kmax_max=2,lmax_max=2) +parameter (kmax_max=4,lmax_max=2) integer lmax,kmax, n_kl(kmax_max,0:lmax_max),l,m,k,kk double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) double precision a(3),g_a,b(3),g_b,c(3),ac(3),bc(3) @@ -486,6 +487,7 @@ do l=1,3 ac(l)=a(l)-c(l) bc(l)=b(l)-c(l) enddo + dr=rmax/npts sum=0.d0 do l=0,lmax @@ -571,6 +573,7 @@ double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k)) enddo Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3)) + !bigI frequantly is null return endif diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index 7fd24174..820ae937 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -8,7 +8,8 @@ double precision :: A_center(3),B_center(3),C_center(3) integer :: power_A(3),power_B(3) integer :: i,j,k,l,n_pt_in,m - double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc, Vpseudo + double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc, Vpseudo, Vpseudo_num + double precision :: dump integer :: nucl_numC ! Important for OpenMP @@ -20,27 +21,34 @@ ! integer klocmax - integer, allocatable :: n_k(:) - double precision, allocatable :: v_k(:), dz_k(:) +! integer, allocatable :: n_k(:) +! double precision, allocatable :: v_k(:), dz_k(:) +! +! call ezfio_get_pseudo_klocmax(klocmax) +! +! allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax)) +! +! call ezfio_get_pseudo_v_k(v_k) +! call ezfio_get_pseudo_n_k(n_k) +! call ezfio_get_pseudo_dz_k(dz_k) - call ezfio_get_pseudo_klocmax(klocmax) + klocmax = 3 - allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax)) + integer :: n_k(3) + double precision :: v_k(3), dz_k(3) - call ezfio_get_pseudo_v_k(v_k) - call ezfio_get_pseudo_n_k(n_k) - call ezfio_get_pseudo_dz_k(dz_k) + v_k(1) = 1.00000000d0 + v_k(2) = 5.35838717 + v_k(3) = -2.07764789 - print*, "klocmax", klocmax + n_k(1) = -1 + n_k(2) = 1 + n_k(3) = 0 - print*, "n_k_ezfio", n_k - print*, "v_k_ezfio",v_k - print*, "dz_k_ezfio", dz_k - - - call ezfio_get_pseudo_v_k(v_k) - call ezfio_get_pseudo_n_k(n_k) - call ezfio_get_pseudo_dz_k(dz_k) + dz_k(1) = 5.35838717 + dz_k(2) = 3.67918975 + dz_k(3) = 1.60507673 + print*, "klocmax", klocmax @@ -56,72 +64,98 @@ !! Parameters of non local part of pseudo: - integer :: kmax,lmax - integer, allocatable :: n_kl(:,:) - double precision, allocatable :: v_kl(:,:), dz_kl(:,:) + integer :: kmax,lmax + integer, allocatable :: n_kl(:,:) + double precision, allocatable :: v_kl(:,:), dz_kl(:,:) - call ezfio_get_pseudo_lmax(lmax) + call ezfio_get_pseudo_lmaxpo(lmax) call ezfio_get_pseudo_kmax(kmax) + !lmax plus one -> lmax lmax = lmax - 1 - + allocate(n_kl(kmax,0:lmax), v_kl(kmax,0:lmax), dz_kl(kmax,0:lmax)) call ezfio_get_pseudo_n_kl(n_kl) call ezfio_get_pseudo_v_kl(v_kl) call ezfio_get_pseudo_dz_kl(dz_kl) + print*, "raw" - print*, "lmax",lmax print*, "kmax", kmax + print*, "lmax",lmax print*,"n_kl_ezfio", n_kl print*,"v_kl_ezfio", v_kl print*,"dz_kl_ezfio", dz_kl +! lmax = 1 +! kmax = 1 + +! integer :: n_kl(1,0:1) +! double precision :: v_kl(1,0:1), dz_kl(1,0:1) + +! v_kl(1,0) =10.69640234 +! n_kl(1,0) = 0 +! dz_kl(1,0) = 1.32389367 +! +! v_kl(1,1) = 10.11238853 +! n_kl(1,1) = 0 +! dz_kl(1,1) = 1.14052020 +! +! print*, "kmax", kmax +! print*, "lmax",lmax +! +! print*,"n_kl_ezfio", n_kl +! print*,"v_kl_ezfio", v_kl +! print*,"dz_kl_ezfio", dz_kl + + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, & - !$OMP num_A, num_B, Z, c, n_pt_in) & + !$OMP num_A, num_B, Z, c, n_pt_in, dump) & !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge, & !$OMP v_k, n_k, dz_k, klocmax, & !$OMP lmax,kmax,v_kl,n_kl,dz_kl) + n_pt_in = n_pt_max_integrals !$OMP DO SCHEDULE (guided) do j = 1, ao_num - power_A(1)= ao_power(j,1) - power_A(2)= ao_power(j,2) - power_A(3)= ao_power(j,3) + num_A = ao_nucl(j) - A_center(1) = nucl_coord(num_A,1) - A_center(2) = nucl_coord(num_A,2) - A_center(3) = nucl_coord(num_A,3) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + do i = 1, ao_num - power_B(1)= ao_power(i,1) - power_B(2)= ao_power(i,2) - power_B(3)= ao_power(i,3) + num_B = ao_nucl(i) - B_center(1) = nucl_coord(num_B,1) - B_center(2) = nucl_coord(num_B,2) - B_center(3) = nucl_coord(num_B,3) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + do l=1,ao_prim_num(j) alpha = ao_expo_transp(l,j) + do m=1,ao_prim_num(i) beta = ao_expo_transp(m,i) + + double precision :: c c = 0.d0 do k = 1, nucl_num - double precision :: Z,c + double precision :: Z Z = nucl_charge(k) - C_center(1) = nucl_coord(k,1) - C_center(2) = nucl_coord(k,2) - C_center(3) = nucl_coord(k,3) - c = c + Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - c = c - Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center) - c = c - Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,C_center) + C_center(1:3) = nucl_coord(k,1:3) + + c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) + + c = c + Vloc( klocmax ,v_k ,n_k ,dz_k, A_center,power_A,alpha,B_center,power_B,beta,C_center) + c = c + Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,beta,C_center) +! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) + enddo - ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - & + ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + & ao_coef_transp(l,j)*ao_coef_transp(m,i)*c enddo enddo diff --git a/src/MonoInts/pseudo.ezfio_config b/src/MonoInts/pseudo.ezfio_config index 1ed75773..941c8ee7 100644 --- a/src/MonoInts/pseudo.ezfio_config +++ b/src/MonoInts/pseudo.ezfio_config @@ -3,8 +3,8 @@ pseudo v_k double precision (pseudo_klocmax) n_k integer (pseudo_klocmax) dz_k double precision (pseudo_klocmax) - lmax integer + lmaxpo integer kmax integer - v_kl double precision (pseudo_kmax,pseudo_lmax) - n_kl integer (pseudo_kmax,pseudo_lmax) - dz_kl double precision (pseudo_kmax,pseudo_lmax) \ No newline at end of file + v_kl double precision (pseudo_kmax,pseudo_lmaxpo) + n_kl integer (pseudo_kmax,pseudo_lmaxpo) + dz_kl double precision (pseudo_kmax,pseudo_lmaxpo) \ No newline at end of file diff --git a/src/MonoInts/test_michel.irp.f b/src/MonoInts/test_michel.irp.f index c7e3c685..ef905479 100644 --- a/src/MonoInts/test_michel.irp.f +++ b/src/MonoInts/test_michel.irp.f @@ -114,7 +114,7 @@ program compute_integrals_pseudo integer, allocatable :: n_kl(:,:) double precision, allocatable :: v_kl(:,:), dz_kl(:,:) - call ezfio_get_pseudo_lmax(lmax) + call ezfio_get_pseudo_lmaxpo(lmax) call ezfio_get_pseudo_kmax(kmax) lmax = lmax - 1 From 0ebbe6f82070481ef47aa2229461c3573c6968ef Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 10 Apr 2015 10:52:55 +0200 Subject: [PATCH 08/27] Working put_pseudo_in_ezfio --- scripts/pseudo/create_ez.sh | 5 +- scripts/pseudo/elts_num_ele.py | 118 ++++++++++++++++++++++++++ scripts/pseudo/put_pseudo_in_ezfio.py | 110 ++++++++++++++++++++---- 3 files changed, 211 insertions(+), 22 deletions(-) create mode 100644 scripts/pseudo/elts_num_ele.py diff --git a/scripts/pseudo/create_ez.sh b/scripts/pseudo/create_ez.sh index 3821d6c2..ca94dc95 100755 --- a/scripts/pseudo/create_ez.sh +++ b/scripts/pseudo/create_ez.sh @@ -4,11 +4,8 @@ echo "name" $1 echo "mul" $2 -echo "Zeff" $3 echo "\`get_basis.sh\` need to be changed" rm -R $1.ezfio qp_create_ezfio_from_xyz $1.xyz -b "cc-pvdz" -m $2 - -~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py --ezfio $1.ezfio/ --atom $1 --zeff $3 -qp_edit -c $1.ezfio \ No newline at end of file +~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py --ezfio $1.ezfio/ --atom $1 \ No newline at end of file diff --git a/scripts/pseudo/elts_num_ele.py b/scripts/pseudo/elts_num_ele.py new file mode 100644 index 00000000..3c4ad09f --- /dev/null +++ b/scripts/pseudo/elts_num_ele.py @@ -0,0 +1,118 @@ +name_to_elec = {"H": 1, + "He": 2, + "Li": 3, + "Be": 4, + "B": 5, + "C": 6, + "N": 7, + "O": 8, + "F": 9, + "Ne": 10, + "Na": 11, + "Mg": 12, + "Al": 13, + "Si": 14, + "P": 15, + "S": 16, + "Cl": 17, + "Ar": 18, + "K": 19, + "Ca": 20, + "Sc": 21, + "Ti": 22, + "V": 23, + "Cr": 24, + "Mn": 25, + "Fe": 26, + "Co": 27, + "Ni": 28, + "Cu": 29, + "Zn": 30, + "Ga": 31, + "Ge": 32, + "As": 33, + "Se": 34, + "Br": 35, + "Kr": 36, + "Rb": 37, + "Sr": 38, + "Y": 39, + "Zr": 40, + "Nb": 41, + "Mo": 42, + "Tc": 43, + "Ru": 44, + "Rh": 45, + "Pd": 46, + "Ag": 47, + "Cd": 48, + "In": 49, + "Sn": 50, + "Sb": 51, + "Te": 52, + "I": 53, + "Xe": 54, + "Cs": 55, + "Ba": 56, + "La": 57, + "Ce": 58, + "Pr": 59, + "Nd": 60, + "Pm": 61, + "Sm": 62, + "Eu": 63, + "Gd": 64, + "Tb": 65, + "Dy": 66, + "Ho": 67, + "Er": 68, + "Tm": 69, + "Yb": 70, + "Lu": 71, + "Hf": 72, + "Ta": 73, + "W": 74, + "Re": 75, + "Os": 76, + "Ir": 77, + "Pt": 78, + "Au": 79, + "Hg": 80, + "Tl": 81, + "Pb": 82, + "Bi": 83, + "Po": 84, + "At": 85, + "Rn": 86, + "Fr": 87, + "Ra": 88, + "Ac": 89, + "Th": 90, + "Pa": 91, + "U": 92, + "Np": 93, + "Pu": 94, + "Am": 95, + "Cm": 96, + "Bk": 97, + "Cf": 98, + "Es": 99, + "Fm": 100, + "Md": 101, + "No": 102, + "Lr": 103, + "Rf": 104, + "Db": 105, + "Sg": 106, + "Bh": 107, + "Hs": 108, + "Mt": 109, + "Ds": 110, + "Rg": 111, + "Cn": 112, + "Uut": 113, + "Fl": 114, + "Uup": 115, + "Lv": 116, + "Uus": 117, + "Uuo": 118} diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 7d15f090..aa4d008b 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -4,7 +4,10 @@ Create the pseudo potential for a given atom Usage: - put_pseudo_in_ezfio.py --ezfio= --atom=... --zeff=... + put_pseudo_in_ezfio.py --ezfio= --atom=... + +Help: + atom is the Abreviation of the atom """ @@ -28,6 +31,24 @@ p = re.compile(ur'\|(\d+)><\d+\|') def get_pseudo_str(l_atom): """ Run EMSL_local for geting the str of the speudo potential + + str_ele : + Element Symbol: Na + Number of replaced protons: 10 + Number of projectors: 2 + + Pseudopotential data: + + Local component: + Coeff. r^n Exp. + 1.00000000 -1 5.35838717 + 5.35838717 1 3.67918975 + -2.07764789 0 1.60507673 + + Non-local component: + Coeff. r^n Exp. Proj. + 10.69640234 0 1.32389367 |0><0| + 10.11238853 0 1.14052020 |1><1| """ EMSL_root = "{0}/EMSL_Basis/".format(qpackage_root) @@ -101,13 +122,63 @@ def get_v_n_dz_l_nonlocal(str_ele): l_dz_kl.append([dz]) if not l_v_kl: - l_v_kl.append([0.]) - l_n_kl.append([0]) - l_dz_kl.append([0.]) + l_v_kl.append([0.]) + l_n_kl.append([0]) + l_dz_kl.append([0.]) return l_v_kl, l_n_kl, l_dz_kl +def get_zeff_alpha_beta(str_ele): + """ + Return the the zeff, alpha num elec and beta num elec + Assert ezfio_set_file alredy defined + """ + + import re + + # ___ + # | ._ o _|_ + # _|_ | | | |_ + # + + # ~#~#~#~#~#~#~ # + # s t r _ e l e # + # ~#~#~#~#~#~#~ # + + m = re.search('Element Symbol: ([a-zA-Z]+)', str_ele) + name = m.group(1).capitalize() + + m = re.search('Number of replaced protons: (\d+)', str_ele) + z_remove = int(m.group(1)) + + # ~#~#~#~#~#~#~#~#~#~ # + # F r o m _ e z f i o # + # ~#~#~#~#~#~#~#~#~#~ # + + alpha = ezfio.get_electrons_elec_alpha_num() + beta = ezfio.get_electrons_elec_beta_num() + + # _ + # |_) _. ._ _ _ + # | (_| | _> (/_ + # + + from elts_num_ele import name_to_elec + z = name_to_elec[name] + + z_eff = z - z_remove + + alpha = alpha - (z_remove / 2) + beta = beta - (z_remove / 2) + + # _ + # |_) _ _|_ ._ ._ + # | \ (/_ |_ |_| | | | + # + + return [z_eff, alpha, beta] + if __name__ == "__main__": arguments = docopt(__doc__) @@ -142,6 +213,10 @@ if __name__ == "__main__": l_str_ele = [str_ele for str_ele in str_.split("Element Symbol: ") if str_ele] + l_zeff = [] + alpha_tot = 0 + beta_tot = 0 + for str_ele in l_str_ele: # ~#~#~#~#~ # @@ -155,14 +230,8 @@ if __name__ == "__main__": # L o c a l # # ~#~#~#~#~ # - print "local" - l_v, l_n, l_dz = get_v_n_dz_local(str_ele[l:nl]) - print l_v - print l_n - print l_dz - ezfio.pseudo_klocmax = len(l_v) ezfio.pseudo_v_k = l_v ezfio.pseudo_n_k = l_n @@ -172,19 +241,24 @@ if __name__ == "__main__": # N o n _ L o c a l # # ~#~#~#~#~#~#~#~#~ # - print "non local" - l_v_kl, l_n_kl, l_dz_kl = get_v_n_dz_l_nonlocal(str_ele[nl:]) - print l_v_kl - print l_n_kl - print l_dz_kl - ezfio.pseudo_lmaxpo = len(l_v_kl) ezfio.pseudo_kmax = len(l_v_kl[0]) ezfio.pseudo_v_kl = l_v_kl ezfio.pseudo_n_kl = l_n_kl ezfio.pseudo_dz_kl = l_dz_kl - if arguments["--zeff"]: - ezfio.nuclei_nucl_charge = map(int, arguments["--zeff"]) + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + # Z _ e f f , a l p h a / b e t a _ e l e c # + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + + zeff, alpha, beta = get_zeff_alpha_beta(str_) + + alpha_tot += alpha + beta_tot += beta + l_zeff.append(zeff) + + ezfio.electrons_elec_alpha_num = alpha_tot + ezfio.electrons_elec_beta_num = beta_tot + ezfio.nuclei_nucl_charge = l_zeff From c1d69942c07f6d84e97e28d53ada6dcdd0e50685 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 10 Apr 2015 16:02:21 +0200 Subject: [PATCH 09/27] working bug ugly pot_ao_ints --- scripts/get_basis.sh | 6 ++--- src/MonoInts/pot_ao_ints.irp.f | 48 +++++++++++++++++----------------- 2 files changed, 27 insertions(+), 27 deletions(-) diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index 1d4e7472..9f959110 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -43,8 +43,8 @@ then exit 1 fi -#${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" $atoms -cp /home/razoa/quantum_package/scripts/pseudo/burkatzki_dz.basis ${tmpfile} -echo ${tmpfile} +#${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" + +${EMSL_API_ROOT}/EMSL_api.py get_basis_data --save --path="${tmpfile}" --basis="${basis}" --db_path="${EMSL_API_ROOT}/db/Pseudo.db" diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index 820ae937..6f35a117 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -21,33 +21,33 @@ ! integer klocmax -! integer, allocatable :: n_k(:) -! double precision, allocatable :: v_k(:), dz_k(:) + integer, allocatable :: n_k(:) + double precision, allocatable :: v_k(:), dz_k(:) + + call ezfio_get_pseudo_klocmax(klocmax) + + allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax)) + + call ezfio_get_pseudo_v_k(v_k) + call ezfio_get_pseudo_n_k(n_k) + call ezfio_get_pseudo_dz_k(dz_k) + +! klocmax = 3 ! -! call ezfio_get_pseudo_klocmax(klocmax) +! integer :: n_k(3) +! double precision :: v_k(3), dz_k(3) ! -! allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax)) +! v_k(1) = 1.00000000d0 +! v_k(2) = 5.35838717 +! v_k(3) = -2.07764789 ! -! call ezfio_get_pseudo_v_k(v_k) -! call ezfio_get_pseudo_n_k(n_k) -! call ezfio_get_pseudo_dz_k(dz_k) - - klocmax = 3 - - integer :: n_k(3) - double precision :: v_k(3), dz_k(3) - - v_k(1) = 1.00000000d0 - v_k(2) = 5.35838717 - v_k(3) = -2.07764789 - - n_k(1) = -1 - n_k(2) = 1 - n_k(3) = 0 - - dz_k(1) = 5.35838717 - dz_k(2) = 3.67918975 - dz_k(3) = 1.60507673 +! n_k(1) = -1 +! n_k(2) = 1 +! n_k(3) = 0 +! +! dz_k(1) = 5.35838717 +! dz_k(2) = 3.67918975 +! dz_k(3) = 1.60507673 print*, "klocmax", klocmax From 96fed344efeb57c3d26a9b487672742d99ef357a Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 10 Apr 2015 16:08:06 +0200 Subject: [PATCH 10/27] Beter automisation for put_pseudo_in_ezfio.py --- scripts/pseudo/put_pseudo_in_ezfio.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index aa4d008b..182e699b 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -4,7 +4,7 @@ Create the pseudo potential for a given atom Usage: - put_pseudo_in_ezfio.py --ezfio= --atom=... + put_pseudo_in_ezfio.py Help: atom is the Abreviation of the atom @@ -181,7 +181,6 @@ def get_zeff_alpha_beta(str_ele): if __name__ == "__main__": arguments = docopt(__doc__) - # ___ # | ._ o _|_ # _|_ | | | |_ @@ -191,7 +190,7 @@ if __name__ == "__main__": # E Z F I O # # ~#~#~#~#~ # - ezfio_path = arguments["--ezfio"] + ezfio_path = arguments[""] ezfio_path = os.path.expanduser(ezfio_path) ezfio_path = os.path.expandvars(ezfio_path) ezfio_path = os.path.abspath(ezfio_path) @@ -202,7 +201,7 @@ if __name__ == "__main__": # P s e u d o _ d a t a # # ~#~#~#~#~#~#~#~#~#~#~ # - l_ele = arguments["--atom"] + l_ele = ezfio.get_nuclei_nucl_label() str_ = get_pseudo_str(l_ele) # _ From 409660b0c6fba64a858e5fd8791f56b77ae87108 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 10 Apr 2015 16:20:28 +0200 Subject: [PATCH 11/27] mend Beter automisation for put_pseudo_in_ezfio.py --- scripts/pseudo/create_ez.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/pseudo/create_ez.sh b/scripts/pseudo/create_ez.sh index ca94dc95..8de384b6 100755 --- a/scripts/pseudo/create_ez.sh +++ b/scripts/pseudo/create_ez.sh @@ -8,4 +8,4 @@ echo "\`get_basis.sh\` need to be changed" rm -R $1.ezfio qp_create_ezfio_from_xyz $1.xyz -b "cc-pvdz" -m $2 -~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py --ezfio $1.ezfio/ --atom $1 \ No newline at end of file +~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py --ezfio $1.ezfio \ No newline at end of file From 5962aff7eeee0b08d9813b1817d733fb28d2c20a Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Fri, 10 Apr 2015 16:51:46 +0200 Subject: [PATCH 12/27] mend Beter automisation for put_pseudo_in_ezfio.py --- scripts/pseudo/create_ez.sh | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/scripts/pseudo/create_ez.sh b/scripts/pseudo/create_ez.sh index 8de384b6..65c12c5e 100755 --- a/scripts/pseudo/create_ez.sh +++ b/scripts/pseudo/create_ez.sh @@ -3,9 +3,10 @@ # $2 mult echo "name" $1 -echo "mul" $2 +echo "basis" $2 +echo "mul" $3 echo "\`get_basis.sh\` need to be changed" rm -R $1.ezfio -qp_create_ezfio_from_xyz $1.xyz -b "cc-pvdz" -m $2 -~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py --ezfio $1.ezfio \ No newline at end of file +qp_create_ezfio_from_xyz $1.xyz -b $2 -m $3 +~/quantum_package/scripts/pseudo/put_pseudo_in_ezfio.py $1.ezfio \ No newline at end of file From e07351c729d3d396eeebac40987e64f922136176 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 13 Apr 2015 13:36:01 +0200 Subject: [PATCH 13/27] Pseduo readme --- src/MonoInts/README.rst | 58 +++++++++++++++++++++++++++++++++-------- 1 file changed, 47 insertions(+), 11 deletions(-) diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index fdbb086b..688bc647 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -63,44 +63,77 @@ Documentation array of the mono electronic hamiltonian on the MOs basis : sum of the kinetic and nuclear electronic potential +`a_coef `_ + Undocumented + +`b_coef `_ + Undocumented + +`ddfact2 `_ + Undocumented + +`erf0 `_ + Undocumented + +`gammln `_ + Undocumented + +`gammp `_ + Undocumented + +`gcf `_ + Undocumented + +`gser `_ + Undocumented + +`rinteg `_ + Undocumented + +`rintgauss `_ + Undocumented + +`sabpartial `_ + Undocumented + `orthonormalize_mos `_ Undocumented `ao_nucl_elec_integral `_ interaction nuclear electron -`ao_nucl_elec_integral_per_atom `_ +`ao_nucl_elec_integral_per_atom `_ ao_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom -`give_polynom_mult_center_mono_elec `_ +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented -`nai_pol_mult `_ +`nai_pol_mult `_ Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented `mo_nucl_elec_integral `_ @@ -221,5 +254,8 @@ Documentation array of the integrals of MO_i * y^2 MO_j array of the integrals of MO_i * z^2 MO_j +`compute_integrals_pseudo `_ + Undocumented + From 25f3b2ee01edc7ad1db610dd719207baef90ea89 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 13 Apr 2015 15:04:27 +0200 Subject: [PATCH 14/27] Working for Li2 but with **WARNING** bad convergence in int_prod_bessel --- scripts/pseudo/put_pseudo_in_ezfio.py | 43 ++++++++++++++------------- src/MonoInts/pot_ao_ints.irp.f | 4 +-- 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 182e699b..47e3d0d8 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -55,18 +55,21 @@ def get_pseudo_str(l_atom): EMSL_path = "{0}/EMSL_api.py".format(EMSL_root) db_path = "{0}/db/Pseudo.db".format(EMSL_root) - l_cmd_atom = [] + str_ = "" + for a in l_atom: - l_cmd_atom += ["--atom", a] + l_cmd_atom = ["--atom", a] - l_cmd_head = [EMSL_path, "get_basis_data", - "--db_path", db_path, - "--basis", "BFD-Pseudo"] + l_cmd_head = [EMSL_path, "get_basis_data", + "--db_path", db_path, + "--basis", "BFD-Pseudo"] - process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE) + process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE) - stdout, _ = process.communicate() - return stdout.strip() + stdout, _ = process.communicate() + str_ += stdout.strip() + "\n" + + return str_ def get_v_n_dz_local(str_ele): @@ -146,19 +149,13 @@ def get_zeff_alpha_beta(str_ele): # s t r _ e l e # # ~#~#~#~#~#~#~ # - m = re.search('Element Symbol: ([a-zA-Z]+)', str_ele) - name = m.group(1).capitalize() +# m = re.search('Element Symbol: ([a-zA-Z]+)', str_ele) +# name = m.group(1).capitalize() + name = str_ele.split("\n")[0].strip().capitalize() m = re.search('Number of replaced protons: (\d+)', str_ele) z_remove = int(m.group(1)) - # ~#~#~#~#~#~#~#~#~#~ # - # F r o m _ e z f i o # - # ~#~#~#~#~#~#~#~#~#~ # - - alpha = ezfio.get_electrons_elec_alpha_num() - beta = ezfio.get_electrons_elec_beta_num() - # _ # |_) _. ._ _ _ # | (_| | _> (/_ @@ -169,8 +166,8 @@ def get_zeff_alpha_beta(str_ele): z_eff = z - z_remove - alpha = alpha - (z_remove / 2) - beta = beta - (z_remove / 2) + alpha = (z_remove / 2) + beta = (z_remove / 2) # _ # |_) _ _|_ ._ ._ @@ -252,12 +249,16 @@ if __name__ == "__main__": # Z _ e f f , a l p h a / b e t a _ e l e c # # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # - zeff, alpha, beta = get_zeff_alpha_beta(str_) + zeff, alpha, beta = get_zeff_alpha_beta(str_ele) alpha_tot += alpha beta_tot += beta l_zeff.append(zeff) + ezfio.nuclei_nucl_charge = l_zeff + + alpha_tot = ezfio.get_electrons_elec_alpha_num() - alpha_tot + beta_tot = ezfio.get_electrons_elec_beta_num() - beta_tot + ezfio.electrons_elec_alpha_num = alpha_tot ezfio.electrons_elec_beta_num = beta_tot - ezfio.nuclei_nucl_charge = l_zeff diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index 6f35a117..bf01649d 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -150,8 +150,8 @@ c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - c = c + Vloc( klocmax ,v_k ,n_k ,dz_k, A_center,power_A,alpha,B_center,power_B,beta,C_center) - c = c + Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,beta,C_center) + c = c + Vloc( klocmax ,v_k ,n_k ,dz_k, A_center,power_A,alpha,B_center,power_B,beta,C_center) + c = c + Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,beta,C_center) ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) enddo From 268b252e11d1bbafc60c49092537197897edf074 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 13 Apr 2015 18:38:46 +0200 Subject: [PATCH 15/27] Working for Li2 but with **WARNING** bad convergence in int_prod_bessel --- scripts/pseudo/put_pseudo_in_ezfio.py | 32 ++++++++++++++++++++++----- src/MonoInts/README.rst | 22 +++++++++--------- src/MonoInts/pot_ao_ints.irp.f | 17 ++++++-------- src/MonoInts/pseudo.ezfio_config | 6 ++--- 4 files changed, 48 insertions(+), 29 deletions(-) diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 47e3d0d8..2ddde4b8 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -209,7 +209,9 @@ if __name__ == "__main__": l_str_ele = [str_ele for str_ele in str_.split("Element Symbol: ") if str_ele] - l_zeff = [] + for i in "l_zeff v_k n_k dz_k v_kl n_kl dz_kl".split(): + exec("{0} = []".format(i)) + alpha_tot = 0 beta_tot = 0 @@ -228,10 +230,9 @@ if __name__ == "__main__": l_v, l_n, l_dz = get_v_n_dz_local(str_ele[l:nl]) - ezfio.pseudo_klocmax = len(l_v) - ezfio.pseudo_v_k = l_v - ezfio.pseudo_n_k = l_n - ezfio.pseudo_dz_k = l_dz + v_k.append(l_v) + n_k.append(l_n) + dz_k.append(l_dz) # ~#~#~#~#~#~#~#~#~ # # N o n _ L o c a l # @@ -255,10 +256,31 @@ if __name__ == "__main__": beta_tot += beta l_zeff.append(zeff) + # _ + # /\ _| _| _|_ _ _ _ _|_ o _ + # /--\ (_| (_| |_ (_) (/_ /_ | | (_) + # + + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + # Z _ e f f , a l p h a / b e t a _ e l e c # + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + ezfio.nuclei_nucl_charge = l_zeff + print "alpha tot", alpha_tot + print "beta tot", beta_tot + alpha_tot = ezfio.get_electrons_elec_alpha_num() - alpha_tot beta_tot = ezfio.get_electrons_elec_beta_num() - beta_tot ezfio.electrons_elec_alpha_num = alpha_tot ezfio.electrons_elec_beta_num = beta_tot + + # ~#~#~#~#~ # + # L o c a l # + # ~#~#~#~#~ # + + ezfio.pseudo_klocmax = len(v_k[0]) + ezfio.pseudo_v_k = zip(*v_k) + ezfio.pseudo_n_k = zip(*n_k) + ezfio.pseudo_dz_k = zip(*dz_k) diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index 688bc647..751407c7 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -102,38 +102,38 @@ Documentation `ao_nucl_elec_integral `_ interaction nuclear electron -`ao_nucl_elec_integral_per_atom `_ +`ao_nucl_elec_integral_per_atom `_ ao_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom -`give_polynom_mult_center_mono_elec `_ +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented -`nai_pol_mult `_ +`nai_pol_mult `_ Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented `mo_nucl_elec_integral `_ diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index bf01649d..e57ffaab 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -21,12 +21,12 @@ ! integer klocmax - integer, allocatable :: n_k(:) - double precision, allocatable :: v_k(:), dz_k(:) + integer, allocatable :: n_k(:,:) + double precision, allocatable :: v_k(:,:), dz_k(:,:) call ezfio_get_pseudo_klocmax(klocmax) - allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax)) + allocate(n_k(nucl_num,klocmax),v_k(nucl_num,klocmax), dz_k(nucl_num,klocmax)) call ezfio_get_pseudo_v_k(v_k) call ezfio_get_pseudo_n_k(n_k) @@ -49,14 +49,13 @@ ! dz_k(2) = 3.67918975 ! dz_k(3) = 1.60507673 - + print*, "nucl_num", nucl_num print*, "klocmax", klocmax print*, "n_k_ezfio", n_k print*, "v_k_ezfio",v_k print*, "dz_k_ezfio", dz_k - ! ! |\ | _ ._ | _ _ _. | ! | \| (_) | | | (_) (_ (_| | @@ -79,8 +78,6 @@ call ezfio_get_pseudo_v_kl(v_kl) call ezfio_get_pseudo_dz_kl(dz_kl) - print*, "raw" - print*, "kmax", kmax print*, "lmax",lmax @@ -116,7 +113,7 @@ !$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, & !$OMP num_A, num_B, Z, c, n_pt_in, dump) & !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & - !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge, & + !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge,nucl_label, & !$OMP v_k, n_k, dz_k, klocmax, & !$OMP lmax,kmax,v_kl,n_kl,dz_kl) @@ -149,8 +146,8 @@ C_center(1:3) = nucl_coord(k,1:3) c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - - c = c + Vloc( klocmax ,v_k ,n_k ,dz_k, A_center,power_A,alpha,B_center,power_B,beta,C_center) + + c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) c = c + Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,beta,C_center) ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) diff --git a/src/MonoInts/pseudo.ezfio_config b/src/MonoInts/pseudo.ezfio_config index 941c8ee7..75e620ce 100644 --- a/src/MonoInts/pseudo.ezfio_config +++ b/src/MonoInts/pseudo.ezfio_config @@ -1,8 +1,8 @@ pseudo klocmax integer - v_k double precision (pseudo_klocmax) - n_k integer (pseudo_klocmax) - dz_k double precision (pseudo_klocmax) + v_k double precision (nuclei_nucl_num,pseudo_klocmax) + n_k integer (nuclei_nucl_num,pseudo_klocmax) + dz_k double precision (nuclei_nucl_num,pseudo_klocmax) lmaxpo integer kmax integer v_kl double precision (pseudo_kmax,pseudo_lmaxpo) From 79ee1f5a087728acb50f6e98b38b4ff985d8c567 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 13 Apr 2015 18:50:37 +0200 Subject: [PATCH 16/27] Working for the LOCAL part of NaCl --- scripts/pseudo/put_pseudo_in_ezfio.py | 22 ++++++++++++++-------- src/MonoInts/pot_ao_ints.irp.f | 8 ++++---- src/MonoInts/pseudo.ezfio_config | 6 +++--- 3 files changed, 21 insertions(+), 15 deletions(-) diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 2ddde4b8..53170271 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -240,11 +240,9 @@ if __name__ == "__main__": l_v_kl, l_n_kl, l_dz_kl = get_v_n_dz_l_nonlocal(str_ele[nl:]) - ezfio.pseudo_lmaxpo = len(l_v_kl) - ezfio.pseudo_kmax = len(l_v_kl[0]) - ezfio.pseudo_v_kl = l_v_kl - ezfio.pseudo_n_kl = l_n_kl - ezfio.pseudo_dz_kl = l_dz_kl + v_kl.append(l_v_kl) + n_kl.append(l_n_kl) + dz_kl.append(l_dz_kl) # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # # Z _ e f f , a l p h a / b e t a _ e l e c # @@ -267,9 +265,6 @@ if __name__ == "__main__": ezfio.nuclei_nucl_charge = l_zeff - print "alpha tot", alpha_tot - print "beta tot", beta_tot - alpha_tot = ezfio.get_electrons_elec_alpha_num() - alpha_tot beta_tot = ezfio.get_electrons_elec_beta_num() - beta_tot @@ -284,3 +279,14 @@ if __name__ == "__main__": ezfio.pseudo_v_k = zip(*v_k) ezfio.pseudo_n_k = zip(*n_k) ezfio.pseudo_dz_k = zip(*dz_k) + + # ~#~#~#~#~#~#~#~#~ # + # N o n _ L o c a l # + # ~#~#~#~#~#~#~#~#~ # + + ezfio.pseudo_lmaxpo = len(v_kl[0]) + ezfio.pseudo_kmax = len(v_kl[0][0]) + + ezfio.pseudo_v_kl = zip(*v_kl) + ezfio.pseudo_n_kl = zip(*n_kl) + ezfio.pseudo_dz_kl = zip(*dz_kl) diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index e57ffaab..4517bb51 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -64,15 +64,15 @@ !! Parameters of non local part of pseudo: integer :: kmax,lmax - integer, allocatable :: n_kl(:,:) - double precision, allocatable :: v_kl(:,:), dz_kl(:,:) + integer, allocatable :: n_kl(:,:,:) + double precision, allocatable :: v_kl(:,:,:), dz_kl(:,:,:) call ezfio_get_pseudo_lmaxpo(lmax) call ezfio_get_pseudo_kmax(kmax) !lmax plus one -> lmax lmax = lmax - 1 - allocate(n_kl(kmax,0:lmax), v_kl(kmax,0:lmax), dz_kl(kmax,0:lmax)) + allocate(n_kl(nucl_num,kmax,0:lmax), v_kl(nucl_num,kmax,0:lmax), dz_kl(nucl_num,kmax,0:lmax)) call ezfio_get_pseudo_n_kl(n_kl) call ezfio_get_pseudo_v_kl(v_kl) @@ -148,7 +148,7 @@ c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) - c = c + Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,beta,C_center) + c = c + Vpseudo(lmax,kmax,v_kl(k,:,:),n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,beta,C_center) ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) enddo diff --git a/src/MonoInts/pseudo.ezfio_config b/src/MonoInts/pseudo.ezfio_config index 75e620ce..97f9e1be 100644 --- a/src/MonoInts/pseudo.ezfio_config +++ b/src/MonoInts/pseudo.ezfio_config @@ -5,6 +5,6 @@ pseudo dz_k double precision (nuclei_nucl_num,pseudo_klocmax) lmaxpo integer kmax integer - v_kl double precision (pseudo_kmax,pseudo_lmaxpo) - n_kl integer (pseudo_kmax,pseudo_lmaxpo) - dz_kl double precision (pseudo_kmax,pseudo_lmaxpo) \ No newline at end of file + v_kl double precision (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo) + n_kl integer (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo) + dz_kl double precision (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo) \ No newline at end of file From 4a823050429887ebedd5bfd935f670d309b5f39b Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 13 Apr 2015 19:01:05 +0200 Subject: [PATCH 17/27] Print debug for n !! --- src/MonoInts/int.f90 | 5 +++++ src/MonoInts/pot_ao_ints.irp.f | 21 ++++++++++++++++++++- 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index cb95c3c9..2c944991 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -1288,6 +1288,11 @@ implicit none integer :: n,k double precision prod dblefact=1.d0 + +if (n.ge.3000) then + print*, n +endif + if(n.lt.0)return if(mod(n,2).eq.1)then prod=1.d0 diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index 4517bb51..cc699574 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -145,10 +145,29 @@ C_center(1:3) = nucl_coord(k,1:3) + print*, j, "j /", ao_num + print*, l, "l /", ao_prim_num(j) + print*, i, "i /", ao_num + print*, m, "m /", ao_prim_num(i) + c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) - c = c + Vpseudo(lmax,kmax,v_kl(k,:,:),n_kl,dz_kl,A_center,power_A,alpha,B_center,power_B,beta,C_center) + + print*, "lmax",lmax + print*, "kmax",kmax + print*, "v_kl",v_kl(k,:,:) + print*, "n_kl",n_kl(k,:,:) + print*, "dz_kl",dz_kl(k,:,:) + print*, "A_center", A_center + print*, "power_A",power_A + print*, "Alpha_B", alpha + print*, "B_center", B_center + print*, "power_B", power_B + print*, "beta", beta + print*, "C_center",C_center + + c = c + Vpseudo(lmax,kmax,v_kl(k,:,:),n_kl(k,:,:),dz_kl(k,:,:),A_center,power_A,alpha,B_center,power_B,beta,C_center) ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) enddo From b964301fd64e0d99dc6d8b502167e4a08d7c8018 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 14 Apr 2015 16:14:58 +0200 Subject: [PATCH 18/27] Put all the same size and allocate in int.f90 --- src/MonoInts/int.f90 | 475 +++++++++++++++++++++++++++---------------- 1 file changed, 303 insertions(+), 172 deletions(-) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index 2c944991..c7d2ac84 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -8,7 +8,7 @@ implicit none integer n_a(3),n_b(3) double precision g_a,g_b,a(3),b(3),c(3) integer kmax_max,lmax_max -parameter (kmax_max=4,lmax_max=2) +parameter (kmax_max=2,lmax_max=2) integer lmax,kmax,n_kl(kmax_max,0:lmax_max) double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) integer klocmax_max @@ -29,7 +29,7 @@ implicit none integer n_a(3),n_b(3) double precision g_a,g_b,a(3),b(3),c(3),rmax integer kmax_max,lmax_max -parameter (kmax_max=4,lmax_max=2) +parameter (kmax_max=2,lmax_max=2) integer lmax,kmax,n_kl(kmax_max,0:lmax_max) double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) integer klocmax_max;parameter (klocmax_max=10) @@ -135,6 +135,14 @@ end if(l.gt.2)stop 'l > 2 not coded!' end +! _ +! | | +! __ __ _ __ ___ ___ _ _ __| | ___ +! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ +! \ V / | |_) \__ \ __/ |_| | (_| | (_) | +! \_/ | .__/|___/\___|\__,_|\__,_|\___/ +! | | +! |_| !! Routine Vpseudo is based on formumla (66) !! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976): @@ -170,27 +178,60 @@ end double precision function Vpseudo & (lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) implicit none + +! ___ +! | ._ ._ _|_ +! _|_ | | |_) |_| |_ +! | double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) integer kmax_max,lmax_max,ntot_max,nkl_max -parameter (kmax_max=4,lmax_max=2,nkl_max=4) +parameter (kmax_max=2,lmax_max=2,nkl_max=4) parameter (ntot_max=10) -integer lmax,kmax,n_kl(kmax_max,0:lmax_max),l,k -double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) -double precision fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm -double precision theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big -double precision areal,freal,breal,t1,t2,int_prod_bessel -integer ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot -integer n_a(3),n_b(3) -double precision array_R(0:ntot_max+nkl_max,kmax_max,0:lmax_max,0:lmax_max+ntot_max,0:lmax_max+ntot_max) -double precision & -array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max) -double precision & -array_I_B(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max) +integer, intent(in) :: lmax,kmax,n_kl(kmax_max,0:lmax_max) +integer, intent(in) :: n_a(3),n_b(3) +double precision, intent(in) :: v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) + + +! +! | _ _ _. | _ +! |_ (_) (_ (_| | (/_ +! + +double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm +double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big +double precision :: areal,freal,breal,t1,t2,int_prod_bessel +double precision :: arg + +integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot +integer :: l,k + +! _ +! |_) o _ _. ._ ._ _. +! |_) | (_| (_| | | (_| \/ +! _| / double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max) double precision array_coefs_B(0:ntot_max,0:ntot_max,0:ntot_max) -double precision arg + +double precision, allocatable :: array_R(:,:,:,:,:) +double precision, allocatable :: array_I_A(:,:,:,:,:) +double precision, allocatable :: array_I_B(:,:,:,:,:) + +!=!=!=!=!=!=!=!=!=! +! A l l o c a t e ! +!=!=!=!=!=!=!=!=!=! + +allocate (array_R(0:ntot_max+nkl_max,kmax_max,0:lmax_max,0:lmax_max+ntot_max,0:lmax_max+ntot_max)) + +allocate (array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) + +allocate (array_I_B(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) + +! _ +! / _. | _ | +! \_ (_| | (_ |_| | +! fourpi=4.d0*dacos(-1.d0) ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) @@ -212,7 +253,17 @@ if(ntot.gt.ntot_max)stop 'increase ntot_max' if(ac.eq.0.d0.and.bc.eq.0.d0)then + + !=!=!=!=!=! + ! I n i t ! + !=!=!=!=!=! + accu=0.d0 + + !=!=!=!=!=!=!=! + ! c a l c u l ! + !=!=!=!=!=!=!=! + do k=1,kmax do l=0,lmax ktot=ntot+n_kl(k,l) @@ -223,11 +274,18 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then enddo enddo enddo - Vpseudo=accu*fourpi - return -endif -if(ac.ne.0.d0.and.bc.ne.0.d0)then + !=!=!=!=! + ! E n d ! + !=!=!=!=! + + Vpseudo=accu*fourpi + +else if(ac.ne.0.d0.and.bc.ne.0.d0)then + + !=!=!=!=!=! + ! I n i t ! + !=!=!=!=!=! f=fourpi**2 @@ -236,24 +294,27 @@ if(ac.ne.0.d0.and.bc.ne.0.d0)then theta_BC0=dacos( (b(3)-c(3))/bc ) phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc) + + + do ktot=0,ntotA+ntotB+nkl_max - do lambda=0,lmax+ntotA - do lambdap=0,lmax+ntotB - do k=1,kmax - do l=0,lmax - array_R(ktot,k,l,lambda,lambdap)= & - freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal) - enddo - enddo - enddo - enddo + do lambda=0,lmax+ntotA + do lambdap=0,lmax+ntotB + do k=1,kmax + do l=0,lmax + array_R(ktot,k,l,lambda,lambdap)= freal & + *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal) + enddo + enddo + enddo + enddo enddo do k1=0,n_a(1) do k2=0,n_a(2) do k3=0,n_a(3) array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & - *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) + *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) enddo enddo enddo @@ -262,79 +323,94 @@ if(ac.ne.0.d0.and.bc.ne.0.d0)then do k2p=0,n_b(2) do k3p=0,n_b(3) array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & - *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) + *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) enddo enddo enddo + !=!=!=!=!=!=!=! + ! c a l c u l ! + !=!=!=!=!=!=!=! + accu=0.d0 do l=0,lmax do m=-l,l - do lambda=0,l+ntotA - do mu=-lambda,lambda - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) - enddo - enddo - enddo - enddo - enddo - - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) - enddo - enddo - enddo - enddo - enddo - - do lambda=0,l+ntotA - do mu=-lambda,lambda - - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - - prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3) - - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - - prodp=ylm(lambdap,mup,theta_BC0,phi_BC0)*array_coefs_B(k1p,k2p,k3p)*array_I_B(lambdap,mup,k1p,k2p,k3p) - - do k=1,kmax - ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l) - accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,lambdap) - enddo - - enddo - enddo + do lambda=0,l+ntotA + do mu=-lambda,lambda + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) + enddo + enddo + enddo enddo enddo - enddo - enddo - enddo - enddo - enddo - enddo + + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) + enddo + enddo + enddo + enddo + enddo + + do lambda=0,l+ntotA + do mu=-lambda,lambda + + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + + prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3) + + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + + prodp=ylm(lambdap,mup,theta_BC0,phi_BC0)*array_coefs_B(k1p,k2p,k3p)*array_I_B(lambdap,mup,k1p,k2p,k3p) + + do k=1,kmax + ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l) + accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,lambdap) + enddo + + enddo + enddo + enddo + + enddo + enddo + + enddo + enddo + enddo + + enddo + enddo + enddo enddo - Vpseudo=f*accu - return -endif -if(ac.eq.0.d0.and.bc.ne.0.d0)then + !=!=!=!=! + ! E n d ! + !=!=!=!=! + + Vpseudo=f*accu + +else if(ac.eq.0.d0.and.bc.ne.0.d0)then + + !=!=!=!=!=! + ! I n i t ! + !=!=!=!=!=! f=fourpi**1.5d0 theta_BC0=dacos( (b(3)-c(3))/bc ) @@ -343,68 +419,85 @@ if(ac.eq.0.d0.and.bc.ne.0.d0)then areal=2.d0*g_a*ac breal=2.d0*g_b*bc freal=dexp(-g_a*ac**2-g_b*bc**2) + do ktot=0,ntotA+ntotB+nkl_max - do lambdap=0,lmax+ntotB - do k=1,kmax - do l=0,lmax - array_R(ktot,k,l,0,lambdap)= & - freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal) - enddo - enddo - enddo + do lambdap=0,lmax+ntotB + do k=1,kmax + do l=0,lmax + + array_R(ktot,k,l,0,lambdap)= freal & + *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal) + enddo + enddo + enddo enddo do k1p=0,n_b(1) do k2p=0,n_b(2) do k3p=0,n_b(3) + array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p) & - *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) + *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) enddo enddo enddo + !=!=!=!=!=!=!=! + ! c a l c u l ! + !=!=!=!=!=!=!=! + accu=0.d0 do l=0,lmax do m=-l,l - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) - enddo - enddo - enddo - enddo - enddo - - prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) - - do lambdap=0,l+ntotB - do mup=-lambdap,lambdap - do k1p=0,n_b(1) - do k2p=0,n_b(2) - do k3p=0,n_b(3) - - prodp=array_coefs_B(k1p,k2p,k3p)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(lambdap,mup,k1p,k2p,k3p) - - do k=1,kmax - ktot=ntotA+k1p+k2p+k3p+n_kl(k,l) - accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,0,lambdap) + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_I_B(lambdap,mup,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) + enddo + enddo + enddo enddo enddo - enddo + + prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) + + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + + prodp=array_coefs_B(k1p,k2p,k3p)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(lambdap,mup,k1p,k2p,k3p) + + do k=1,kmax + + ktot=ntotA+k1p+k2p+k3p+n_kl(k,l) + accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,0,lambdap) + + enddo + + enddo + enddo + enddo enddo enddo - enddo enddo enddo - Vpseudo=f*accu - return -endif -if(ac.ne.0.d0.and.bc.eq.0.d0)then + !=!=!=!=! + ! E n d ! + !=!=!=!=! + + Vpseudo=f*accu + +else if(ac.ne.0.d0.and.bc.eq.0.d0)then + + !=!=!=!=!=! + ! I n i t ! + !=!=!=!=!=! f=fourpi**1.5d0 theta_AC0=dacos( (a(3)-c(3))/ac ) @@ -413,69 +506,102 @@ if(ac.ne.0.d0.and.bc.eq.0.d0)then areal=2.d0*g_a*ac breal=2.d0*g_b*bc freal=dexp(-g_a*ac**2-g_b*bc**2) + do ktot=0,ntotA+ntotB+nkl_max - do lambda=0,lmax+ntotA - do k=1,kmax - do l=0,lmax - array_R(ktot,k,l,lambda,0)= & - freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal) - enddo - enddo - enddo + do lambda=0,lmax+ntotA + do k=1,kmax + do l=0,lmax + + array_R(ktot,k,l,lambda,0)= freal & + *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal) + + enddo + enddo + enddo enddo do k1=0,n_a(1) do k2=0,n_a(2) do k3=0,n_a(3) + array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3) & - *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) + *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) + enddo enddo enddo + !=!=!=!=!=!=!=! + ! c a l c u l ! + !=!=!=!=!=!=!=! + accu=0.d0 do l=0,lmax do m=-l,l - do lambda=0,l+ntotA - do mu=-lambda,lambda - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) - enddo - enddo - enddo - enddo - enddo - - do lambda=0,l+ntotA - do mu=-lambda,lambda - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - prod=array_coefs_A(k1,k2,k3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(lambda,mu,k1,k2,k3) - prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) - do k=1,kmax - ktot=k1+k2+k3+ntotB+n_kl(k,l) - accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,0) + do lambda=0,l+ntotA + do mu=-lambda,lambda + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) + enddo + enddo + enddo enddo - enddo - enddo + enddo + + do lambda=0,l+ntotA + do mu=-lambda,lambda + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + + prod=array_coefs_A(k1,k2,k3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(lambda,mu,k1,k2,k3) + prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) + + do k=1,kmax + ktot=k1+k2+k3+ntotB+n_kl(k,l) + accu=accu+prod*prodp*v_kl(k,l)*array_R(ktot,k,l,lambda,0) + enddo + + enddo + enddo + enddo enddo enddo - enddo + enddo enddo + + !=!=!=!=! + ! E n d ! + !=!=!=!=! + Vpseudo=f*accu - return endif + +! _ +! |_ o ._ _. | o _ _ +! | | | | (_| | | _> (/_ +! + deallocate (array_R, array_I_A, array_I_B) + return end +! _ +! | | +!__ __ _ __ ___ ___ _ _ __| | ___ _ __ _ _ _ __ ___ +!\ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ | '_ \| | | | '_ ` _ \ +! \ V / | |_) \__ \ __/ |_| | (_| | (_) | | | | | |_| | | | | | | +! \_/ | .__/|___/\___|\__,_|\__,_|\___/ |_| |_|\__,_|_| |_| |_| +! | | +! |_| + double precision function Vpseudo_num(npts,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) implicit none integer kmax_max,lmax_max -parameter (kmax_max=4,lmax_max=2) +parameter (kmax_max=2,lmax_max=2) integer lmax,kmax, n_kl(kmax_max,0:lmax_max),l,m,k,kk double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) double precision a(3),g_a,b(3),g_b,c(3),ac(3),bc(3) @@ -1289,10 +1415,6 @@ integer :: n,k double precision prod dblefact=1.d0 -if (n.ge.3000) then - print*, n -endif - if(n.lt.0)return if(mod(n,2).eq.1)then prod=1.d0 @@ -1777,7 +1899,16 @@ end endif enddo int_prod_bessel=int - if(kcp.gt.100)print*,'**WARNING** bad convergence in int_prod_bessel' + if(kcp.gt.100) then + print*,"l",l + print*, "gam", gam + print*, "n", n + print*, "m", m + print*, "a", a + print*, "b", b + print*, "kcp", kcp + print*,'**WARNING** bad convergence in int_prod_bessel' + endif return endif From 1fc1124d959b89048255a9abc692890d72e39c1a Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Thu, 16 Apr 2015 09:41:16 +0200 Subject: [PATCH 19/27] Remove open mp and add tmp array --- src/MonoInts/README.rst | 22 +++++----- src/MonoInts/pot_ao_ints.irp.f | 77 ++++++++++++++++++++-------------- 2 files changed, 56 insertions(+), 43 deletions(-) diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index 751407c7..ec92eada 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -102,38 +102,38 @@ Documentation `ao_nucl_elec_integral `_ interaction nuclear electron -`ao_nucl_elec_integral_per_atom `_ +`ao_nucl_elec_integral_per_atom `_ ao_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom -`give_polynom_mult_center_mono_elec `_ +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented -`nai_pol_mult `_ +`nai_pol_mult `_ Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented `mo_nucl_elec_integral `_ diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index cc699574..ef2e8c8f 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -48,7 +48,10 @@ ! dz_k(1) = 5.35838717 ! dz_k(2) = 3.67918975 ! dz_k(3) = 1.60507673 - + print*, "=======================" + print*, "=======================" + print*, "=======================" + print*, "nucl_num", nucl_num print*, "klocmax", klocmax @@ -85,6 +88,10 @@ print*,"v_kl_ezfio", v_kl print*,"dz_kl_ezfio", dz_kl + print*, "=======================" + print*, "=======================" + print*, "=======================" + ! lmax = 1 ! kmax = 1 @@ -108,29 +115,31 @@ ! print*,"dz_kl_ezfio", dz_kl - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, & - !$OMP num_A, num_B, Z, c, n_pt_in, dump) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & - !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge,nucl_label, & - !$OMP v_k, n_k, dz_k, klocmax, & - !$OMP lmax,kmax,v_kl,n_kl,dz_kl) + integer, allocatable :: n_kl_dump(:,:) + double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:) + + allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax)) + n_pt_in = n_pt_max_integrals - !$OMP DO SCHEDULE (guided) do j = 1, ao_num num_A = ao_nucl(j) power_A(1:3)= ao_power(j,1:3) A_center(1:3) = nucl_coord(num_A,1:3) + print*, "J", j, "/", ao_num + print*,"===================" + do i = 1, ao_num num_B = ao_nucl(i) power_B(1:3)= ao_power(i,1:3) B_center(1:3) = nucl_coord(num_B,1:3) + + print*, "i", i, "/", ao_num + do l=1,ao_prim_num(j) alpha = ao_expo_transp(l,j) @@ -145,31 +154,37 @@ C_center(1:3) = nucl_coord(k,1:3) - print*, j, "j /", ao_num - print*, l, "l /", ao_prim_num(j) - print*, i, "i /", ao_num - print*, m, "m /", ao_prim_num(i) - c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) - - print*, "lmax",lmax - print*, "kmax",kmax - print*, "v_kl",v_kl(k,:,:) - print*, "n_kl",n_kl(k,:,:) - print*, "dz_kl",dz_kl(k,:,:) - print*, "A_center", A_center - print*, "power_A",power_A - print*, "Alpha_B", alpha - print*, "B_center", B_center - print*, "power_B", power_B - print*, "beta", beta - print*, "C_center",C_center - c = c + Vpseudo(lmax,kmax,v_kl(k,:,:),n_kl(k,:,:),dz_kl(k,:,:),A_center,power_A,alpha,B_center,power_B,beta,C_center) -! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) + n_kl_dump = n_kl(k,1:kmax,0:lmax) + v_kl_dump = v_kl(k,1:kmax,0:lmax) + dz_kl_dump = dz_kl(k,1:kmax,0:lmax) + +! print*, "lmax",lmax +! print*, "kmax",kmax +! print*, "v_kl",v_kl_dump +! print*, "n_kl",n_kl_dump +! print*, n_kl_dump(1,0) +! print*, n_kl_dump(1,1) +! print*, "dz_kl",dz_kl_dump +! print*, dz_kl_dump(1,0) +! print*, dz_kl_dump(1,1) +! print*, "A_center", A_center +! print*, "power_A",power_A +! print*, "alpha", alpha +! print*, "B_center", B_center +! print*, "power_B", power_B +! print*, "beta", beta +! print*, "C_center",C_center + + ! c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) + +! print*, "#################" +! print*, "#################" enddo ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + & ao_coef_transp(l,j)*ao_coef_transp(m,i)*c @@ -177,8 +192,6 @@ enddo enddo enddo - !$OMP END DO - !$OMP END PARALLEL END_PROVIDER From 8054b1e58ace6286a628fad6983c7084ee7a946e Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 20 Apr 2015 17:17:36 +0200 Subject: [PATCH 20/27] New version of int.f90 for big alpha but not to much --- src/MonoInts/int.f90 | 70 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 62 insertions(+), 8 deletions(-) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index c7d2ac84..85b5c71e 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -140,7 +140,7 @@ end ! __ __ _ __ ___ ___ _ _ __| | ___ ! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ ! \ V / | |_) \__ \ __/ |_| | (_| | (_) | -! \_/ | .__/|___/\___|\__,_|\__,_|\___/ +! \_/ | .__/|___/\___|\__,_|\____|\___/ ! | | ! |_| @@ -200,7 +200,7 @@ double precision, intent(in) :: v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_ double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big -double precision :: areal,freal,breal,t1,t2,int_prod_bessel +double precision :: areal,freal,breal,t1,t2,int_prod_bessel, int_prod_bessel_num_soph_p double precision :: arg integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot @@ -270,7 +270,9 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then do m=-l,l prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) - accu=accu+prod*prodp*v_kl(k,l)*freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal) + + accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg) + enddo enddo enddo @@ -303,7 +305,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then do k=1,kmax do l=0,lmax array_R(ktot,k,l,lambda,lambdap)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal) + *int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg) enddo enddo enddo @@ -426,8 +428,8 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then do l=0,lmax array_R(ktot,k,l,0,lambdap)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal) - enddo + * int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg) + enddo enddo enddo enddo @@ -513,8 +515,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then do l=0,lmax array_R(ktot,k,l,lambda,0)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal) - + * int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg) enddo enddo enddo @@ -1974,6 +1975,59 @@ end stop 'pb in int_prod_bessel!!' end + + double precision function int_prod_bessel_num_soph_p(l,gam,n,m,a,b,arg) + implicit none + integer n,m,l + double precision gam,a,b,arg,arg_new + double precision bessel_mod,factor +logical not_done +double precision bigA,xold,x,dx,accu,intnew,intold,intold2,u,v,freal +integer iter +double precision n0,nI,i,eps + u=(a+b)/(2.d0*dsqrt(gam)) + arg_new=u**2-arg + freal=dexp(arg_new) + v=u/dsqrt(gam) + +bigA=v+dsqrt(-dlog(1.d-15)/gam) +n0=5 +accu=0.d0 +dx=bigA/(n0-1.d0) +iter=0 +do i=1.d0,n0 + x=(i-1.d0)*dx + accu=accu+x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) +enddo +accu=accu*freal +intold=accu*dx + +eps=1.d-08 +nI=n0-1.d0 +dx=dx/2.d0 +not_done=.true. + +do while(not_done) + iter=iter+1 + accu=0.d0 + do i=1.d0,nI + x=dx+(i-1.d0)*2.d0*dx + accu=accu+dx*x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) + enddo + accu=accu*freal + intnew=intold/2.d0+accu + if(iter.gt.1.and.dabs(intnew-intold).lt.eps.and.dabs(intnew-intold2).lt.eps)then + not_done=.false. + else + intold2=intold + intold=intnew + dx=dx/2.d0 + nI=2.d0*nI + endif +enddo +int_prod_bessel_num_soph_p=intold +end + !! Calculation of !! !! I= \int dx x**l *exp(-gam*x**2) M_n(ax) From fbab2d613eee25e6f611072b4fa4c978673f3ff1 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 20 Apr 2015 17:17:36 +0200 Subject: [PATCH 21/27] New version of int.f90 for big alpha but not to much --- src/MonoInts/int.f90 | 77 ++++++++++++++++++++++++++++------ src/MonoInts/pot_ao_ints.irp.f | 8 ++-- 2 files changed, 70 insertions(+), 15 deletions(-) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index c7d2ac84..750dbaa7 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -140,7 +140,7 @@ end ! __ __ _ __ ___ ___ _ _ __| | ___ ! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ ! \ V / | |_) \__ \ __/ |_| | (_| | (_) | -! \_/ | .__/|___/\___|\__,_|\__,_|\___/ +! \_/ | .__/|___/\___|\__,_|\____|\___/ ! | | ! |_| @@ -200,7 +200,7 @@ double precision, intent(in) :: v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_ double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big -double precision :: areal,freal,breal,t1,t2,int_prod_bessel +double precision :: areal,freal,breal,t1,t2,int_prod_bessel, int_prod_bessel_num_soph_p double precision :: arg integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot @@ -237,7 +237,7 @@ fourpi=4.d0*dacos(-1.d0) ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) arg=g_a*ac**2+g_b*bc**2 -if(arg.gt.-dlog(10.d-20))then +if(arg.gt.-dlog(1.d-20))then Vpseudo=0.d0 return endif @@ -270,7 +270,9 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then do m=-l,l prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) - accu=accu+prod*prodp*v_kl(k,l)*freal*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal) + + accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg) + enddo enddo enddo @@ -302,8 +304,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then do lambdap=0,lmax+ntotB do k=1,kmax do l=0,lmax - array_R(ktot,k,l,lambda,lambdap)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal) + array_R(ktot,k,l,lambda,lambdap)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg) enddo enddo enddo @@ -425,9 +426,8 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then do k=1,kmax do l=0,lmax - array_R(ktot,k,l,0,lambdap)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal) - enddo + array_R(ktot,k,l,0,lambdap)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg) + enddo enddo enddo enddo @@ -512,9 +512,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then do k=1,kmax do l=0,lmax - array_R(ktot,k,l,lambda,0)= freal & - *int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal) - + array_R(ktot,k,l,lambda,0)= int_prod_bessel_num_soph_p(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg) enddo enddo enddo @@ -1974,6 +1972,61 @@ end stop 'pb in int_prod_bessel!!' end + +double precision function int_prod_bessel_num_soph_p(l,gam,n,m,a,b,arg) + implicit none + integer :: n,m,l + double precision :: gam,a,b,arg,arg_new + double precision :: bessel_mod,factor + logical :: not_done + double precision :: bigA,xold,x,dx,accu,intnew,intold,intold2,u,v,freal + integer :: iter, i, nI, n0 + double precision :: eps + + u=(a+b)/(2.d0*dsqrt(gam)) + arg_new=u**2-arg + freal=dexp(arg_new) + v=u/dsqrt(gam) + + bigA=v+dsqrt(-dlog(1.d-15)/gam) + n0=5 + accu=0.d0 + dx=bigA/(float(n0)-1.d0) + iter=0 + do i=1,n0 + x=(float(i)-1.d0)*dx + accu=accu+x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) + enddo + +accu=accu*freal +intold=accu*dx + +eps=1.d-08 +nI=n0-1 +dx=dx/2.d0 +not_done=.true. + +do while(not_done) + iter=iter+1 + accu=0.d0 + do i=1,nI + x=dx+(float(i)-1.d0)*2.d0*dx + accu=accu+dx*x**l*dexp(-gam*(x-v)**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) + enddo + accu=accu*freal + intnew=intold/2.d0+accu + if(iter.gt.1.and.dabs(intnew-intold).lt.eps.and.dabs(intnew-intold2).lt.eps)then + not_done=.false. + else + intold2=intold + intold=intnew + dx=dx/2.d0 + nI=2*nI + endif +enddo +int_prod_bessel_num_soph_p=intold +end + !! Calculation of !! !! I= \int dx x**l *exp(-gam*x**2) M_n(ax) diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index ef2e8c8f..badf4afd 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -156,7 +156,7 @@ c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) +! c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) n_kl_dump = n_kl(k,1:kmax,0:lmax) @@ -167,7 +167,7 @@ ! print*, "kmax",kmax ! print*, "v_kl",v_kl_dump ! print*, "n_kl",n_kl_dump -! print*, n_kl_dump(1,0) +! print*, n_kl_ump(1,0) ! print*, n_kl_dump(1,1) ! print*, "dz_kl",dz_kl_dump ! print*, dz_kl_dump(1,0) @@ -180,7 +180,9 @@ ! print*, "beta", beta ! print*, "C_center",C_center - ! c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + dump = Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + print*, dump ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) ! print*, "#################" From cda3a5ee9c969e5d735ce0f3a6123543a47107f7 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 21 Apr 2015 15:05:39 +0200 Subject: [PATCH 22/27] Working ! --- scripts/get_basis.sh | 2 +- scripts/pseudo/put_pseudo_in_ezfio.py | 1 - src/MonoInts/README.rst | 22 +- src/MonoInts/int.f90 | 4 +- src/MonoInts/pot_ao_ints.irp.f | 303 +++++++++++++++----------- src/MonoInts/pseudo.ezfio_config | 3 +- 6 files changed, 188 insertions(+), 147 deletions(-) diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index 9f959110..7cfe8305 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -47,4 +47,4 @@ fi #${EMSL_API_ROOT}/EMSL_api.py get_basis_data --treat_l --save --path="${tmpfile}" --basis="${basis}" -${EMSL_API_ROOT}/EMSL_api.py get_basis_data --save --path="${tmpfile}" --basis="${basis}" --db_path="${EMSL_API_ROOT}/db/Pseudo.db" +${EMSL_API_ROOT}/EMSL_api.py get_basis_data --save --path="${tmpfile}" --basis="${basis}" --db_path="${EMSL_API_ROOT}/db/Pseudo.db" \ No newline at end of file diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 53170271..0ba71b0c 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -253,7 +253,6 @@ if __name__ == "__main__": alpha_tot += alpha beta_tot += beta l_zeff.append(zeff) - # _ # /\ _| _| _|_ _ _ _ _|_ o _ # /--\ (_| (_| |_ (_) (/_ /_ | | (_) diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index ec92eada..69da98ed 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -102,38 +102,38 @@ Documentation `ao_nucl_elec_integral `_ interaction nuclear electron -`ao_nucl_elec_integral_per_atom `_ +`ao_nucl_elec_integral_per_atom `_ ao_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom -`give_polynom_mult_center_mono_elec `_ +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented -`nai_pol_mult `_ +`nai_pol_mult `_ Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented `mo_nucl_elec_integral `_ diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index 750dbaa7..31b51c51 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -188,9 +188,9 @@ double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) integer kmax_max,lmax_max,ntot_max,nkl_max parameter (kmax_max=2,lmax_max=2,nkl_max=4) parameter (ntot_max=10) -integer, intent(in) :: lmax,kmax,n_kl(kmax_max,0:lmax_max) +integer, intent(in) :: lmax,kmax,n_kl(kmax,0:lmax) integer, intent(in) :: n_a(3),n_b(3) -double precision, intent(in) :: v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) +double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) ! diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index badf4afd..065785c2 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -4,21 +4,93 @@ END_DOC implicit none double precision :: alpha, beta, gama, delta - integer :: i_c,num_A,num_B - double precision :: A_center(3),B_center(3),C_center(3) - integer :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt_in,m - double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc, Vpseudo, Vpseudo_num - double precision :: dump - integer :: nucl_numC - ! Important for OpenMP + integer :: num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt_in,m + double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - ao_nucl_elec_integral = 0.d0 + ao_nucl_elec_integral = ao_nucl_elec_integral_pseudo ! 0.d0 + + ! _ + ! /| / |_) + ! | / | \ + ! + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, & + !$OMP num_A,num_B,Z,c,n_pt_in) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & + !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge) + + n_pt_in = n_pt_max_integrals + + !$OMP DO SCHEDULE (guided) + + do j = 1, ao_num + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_transp(m,i) + + double precision :: c + c = 0.d0 + + do k = 1, nucl_num + double precision :: Z + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) + + enddo + ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + & + ao_coef_transp(l,j)*ao_coef_transp(m,i)*c + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + + END_PROVIDER + + BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_pseudo, (ao_num_align,ao_num)] + BEGIN_DOC +! interaction nuclear electron + END_DOC + implicit none + double precision :: alpha, beta, gama, delta + integer :: num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt_in,m + double precision :: Vloc, Vpseudo + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + integer :: thread_num + + ao_nucl_elec_integral_pseudo = 0.d0 ! ! | _ _ _. | ! |_ (_) (_ (_| | ! + !! Parameters of the local part of pseudo: integer klocmax integer, allocatable :: n_k(:,:) @@ -32,43 +104,22 @@ call ezfio_get_pseudo_n_k(n_k) call ezfio_get_pseudo_dz_k(dz_k) -! klocmax = 3 -! -! integer :: n_k(3) -! double precision :: v_k(3), dz_k(3) -! -! v_k(1) = 1.00000000d0 -! v_k(2) = 5.35838717 -! v_k(3) = -2.07764789 -! -! n_k(1) = -1 -! n_k(2) = 1 -! n_k(3) = 0 -! -! dz_k(1) = 5.35838717 -! dz_k(2) = 3.67918975 -! dz_k(3) = 1.60507673 - print*, "=======================" - print*, "=======================" - print*, "=======================" + !! Dump array + integer, allocatable :: n_k_dump(:) + double precision, allocatable :: v_k_dump(:), dz_k_dump(:) + + allocate(n_k_dump(1:klocmax), v_k_dump(1:klocmax), dz_k_dump(1:klocmax)) - print*, "nucl_num", nucl_num - print*, "klocmax", klocmax - - print*, "n_k_ezfio", n_k - print*, "v_k_ezfio",v_k - print*, "dz_k_ezfio", dz_k ! ! |\ | _ ._ | _ _ _. | ! | \| (_) | | | (_) (_ (_| | ! - !! Parameters of non local part of pseudo: - integer :: kmax,lmax - integer, allocatable :: n_kl(:,:,:) - double precision, allocatable :: v_kl(:,:,:), dz_kl(:,:,:) + integer :: kmax,lmax + integer, allocatable :: n_kl(:,:,:) + double precision, allocatable :: v_kl(:,:,:), dz_kl(:,:,:) call ezfio_get_pseudo_lmaxpo(lmax) call ezfio_get_pseudo_kmax(kmax) @@ -81,55 +132,43 @@ call ezfio_get_pseudo_v_kl(v_kl) call ezfio_get_pseudo_dz_kl(dz_kl) - print*, "kmax", kmax - print*, "lmax",lmax - - print*,"n_kl_ezfio", n_kl - print*,"v_kl_ezfio", v_kl - print*,"dz_kl_ezfio", dz_kl - - print*, "=======================" - print*, "=======================" - print*, "=======================" - - -! lmax = 1 -! kmax = 1 - -! integer :: n_kl(1,0:1) -! double precision :: v_kl(1,0:1), dz_kl(1,0:1) - -! v_kl(1,0) =10.69640234 -! n_kl(1,0) = 0 -! dz_kl(1,0) = 1.32389367 -! -! v_kl(1,1) = 10.11238853 -! n_kl(1,1) = 0 -! dz_kl(1,1) = 1.14052020 -! -! print*, "kmax", kmax -! print*, "lmax",lmax -! -! print*,"n_kl_ezfio", n_kl -! print*,"v_kl_ezfio", v_kl -! print*,"dz_kl_ezfio", dz_kl - + !! Dump array integer, allocatable :: n_kl_dump(:,:) double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:) allocate(n_kl_dump(kmax,0:lmax), v_kl_dump(kmax,0:lmax), dz_kl_dump(kmax,0:lmax)) + ! _ + ! / _. | _ | + ! \_ (_| | (_ |_| | + ! - n_pt_in = n_pt_max_integrals - do j = 1, ao_num + write(output_monoints,*) 'Providing the nuclear electron pseudo integrals ' - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) + call wall_time(wall_1) + call cpu_time(cpu_1) - print*, "J", j, "/", ao_num - print*,"===================" + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, & + !$OMP num_A,num_B,Z,c,n_pt_in, & + !$OMP v_k_dump,n_k_dump, dz_k_dump, n_kl_dump, v_kl_dump, dz_kl_dump, & + !$OMP wall_0,wall_2,thread_num, output_monoints) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & + !$OMP n_pt_max_integrals,ao_nucl_elec_integral_pseudo,nucl_num,nucl_charge, & + !$OMP klocmax,lmax,kmax,v_k,n_k, dz_k, n_kl, v_kl, dz_kl, & + !$OMP wall_1) + + n_pt_in = n_pt_max_integrals + + !$OMP DO SCHEDULE (guided) + + do j = 1, ao_num + + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) do i = 1, ao_num @@ -137,63 +176,66 @@ power_B(1:3)= ao_power(i,1:3) B_center(1:3) = nucl_coord(num_B,1:3) - - print*, "i", i, "/", ao_num - - do l=1,ao_prim_num(j) - alpha = ao_expo_transp(l,j) - + do l=1,ao_prim_num(j) + alpha = ao_expo_transp(l,j) + do m=1,ao_prim_num(i) - beta = ao_expo_transp(m,i) + beta = ao_expo_transp(m,i) + double precision :: c + c = 0.d0 + + do k = 1, nucl_num + double precision :: Z + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + v_k_dump = v_k(k,1:klocmax) + n_k_dump = n_k(k,1:klocmax) + dz_k_dump = dz_k(k,1:klocmax) - double precision :: c - c = 0.d0 - do k = 1, nucl_num - double precision :: Z - Z = nucl_charge(k) + c = c + Vloc(klocmax, v_k_dump,n_k_dump, dz_k_dump, & + A_center,power_A,alpha,B_center,power_B,beta,C_center) + - C_center(1:3) = nucl_coord(k,1:3) - - c = c - Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) - -! c = c + Vloc( klocmax ,v_k(k,:) ,n_k(k,:) ,dz_k(k,:), A_center,power_A,alpha,B_center,power_B,beta,C_center) - - - n_kl_dump = n_kl(k,1:kmax,0:lmax) - v_kl_dump = v_kl(k,1:kmax,0:lmax) - dz_kl_dump = dz_kl(k,1:kmax,0:lmax) - -! print*, "lmax",lmax -! print*, "kmax",kmax -! print*, "v_kl",v_kl_dump -! print*, "n_kl",n_kl_dump -! print*, n_kl_ump(1,0) -! print*, n_kl_dump(1,1) -! print*, "dz_kl",dz_kl_dump -! print*, dz_kl_dump(1,0) -! print*, dz_kl_dump(1,1) -! print*, "A_center", A_center -! print*, "power_A",power_A -! print*, "alpha", alpha -! print*, "B_center", B_center -! print*, "power_B", power_B -! print*, "beta", beta -! print*, "C_center",C_center - - c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) - dump = Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) - print*, dump - ! c = c - Vps(A_center,power_A,alpha,B_center,power_B,beta,C_center,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) - -! print*, "#################" -! print*, "#################" + n_kl_dump = n_kl(k,1:kmax,0:lmax) + v_kl_dump = v_kl(k,1:kmax,0:lmax) + dz_kl_dump = dz_kl(k,1:kmax,0:lmax) + + c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + + enddo + ao_nucl_elec_integral_pseudo(i,j) = ao_nucl_elec_integral_pseudo(i,j) + & + ao_coef_transp(l,j)*ao_coef_transp(m,i)*c + enddo enddo - ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + & - ao_coef_transp(l,j)*ao_coef_transp(m,i)*c - enddo - enddo enddo - enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + write(output_monoints,*) 100.*float(j)/float(ao_num), '% in ', & + wall_2-wall_1, 's' + endif + endif + enddo + + !$OMP END DO + !$OMP END PARALLEL + + +! _ +! | \ _ _. | | _ _ _. _|_ _ +! |_/ (/_ (_| | | (_) (_ (_| |_ (/_ +! + + deallocate(n_k,v_k, dz_k) + deallocate(n_k_dump,v_k_dump, dz_k_dump) + + deallocate(n_kl,v_kl, dz_kl) + deallocate(n_kl_dump,v_kl_dump, dz_kl_dump) + END_PROVIDER @@ -210,7 +252,6 @@ END_PROVIDER integer :: power_A(3),power_B(3) integer :: i,j,k,l,n_pt_in,m double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - integer :: nucl_numC ! Important for OpenMP ao_nucl_elec_integral_per_atom = 0.d0 diff --git a/src/MonoInts/pseudo.ezfio_config b/src/MonoInts/pseudo.ezfio_config index 97f9e1be..db0da938 100644 --- a/src/MonoInts/pseudo.ezfio_config +++ b/src/MonoInts/pseudo.ezfio_config @@ -1,4 +1,5 @@ pseudo + do_pseudo logical klocmax integer v_k double precision (nuclei_nucl_num,pseudo_klocmax) n_k integer (nuclei_nucl_num,pseudo_klocmax) @@ -7,4 +8,4 @@ pseudo kmax integer v_kl double precision (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo) n_kl integer (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo) - dz_kl double precision (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo) \ No newline at end of file + dz_kl double precision (nuclei_nucl_num,pseudo_kmax,pseudo_lmaxpo) From 185616f7e304ba0b493a54f973d2a37d3e363841 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 28 Apr 2015 14:03:44 +0200 Subject: [PATCH 23/27] Add capitzalise to EZFIO_convert_output for consitenciy --- .../qp_convert_output_to_ezfio.py | 2 +- src/MonoInts/README.rst | 25 +++++++++++-------- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py index 6b5c5fcd..9d67611e 100755 --- a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py +++ b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py @@ -105,7 +105,7 @@ def write_ezfio(res, filename): # Transformt H1 into H import re p = re.compile(ur'(\d*)$') - label = [p.sub("", x.name) for x in res.geometry] + label = [p.sub("", x.name).capitalize() for x in res.geometry] ezfio.set_nuclei_nucl_label(label) ezfio.set_nuclei_nucl_coord(coord_x + coord_y + coord_z) diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index 69da98ed..ffdcdc94 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -102,38 +102,41 @@ Documentation `ao_nucl_elec_integral `_ interaction nuclear electron -`ao_nucl_elec_integral_per_atom `_ +`ao_nucl_elec_integral_per_atom `_ ao_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom -`give_polynom_mult_center_mono_elec `_ +`ao_nucl_elec_integral_pseudo `_ + interaction nuclear electron + +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented -`nai_pol_mult `_ +`nai_pol_mult `_ Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented `mo_nucl_elec_integral `_ From be4db7c56d28505d3417fd7a1526a9cf3fb7247e Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 28 Apr 2015 17:08:59 +0200 Subject: [PATCH 24/27] delete the **** parameters --- src/MonoInts/int.f90 | 78 +++++++++++++++++++++++----------- src/MonoInts/pot_ao_ints.irp.f | 2 +- 2 files changed, 54 insertions(+), 26 deletions(-) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index 31b51c51..98247553 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -184,15 +184,10 @@ implicit none ! _|_ | | |_) |_| |_ ! | double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) - -integer kmax_max,lmax_max,ntot_max,nkl_max -parameter (kmax_max=2,lmax_max=2,nkl_max=4) -parameter (ntot_max=10) integer, intent(in) :: lmax,kmax,n_kl(kmax,0:lmax) integer, intent(in) :: n_a(3),n_b(3) double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) - ! ! | _ _ _. | _ ! |_ (_) (_ (_| | (/_ @@ -204,35 +199,36 @@ double precision :: areal,freal,breal,t1,t2,int_prod_bessel, int_prod_bessel_num double precision :: arg integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot -integer :: l,k +integer :: l,k, nkl_max ! _ ! |_) o _ _. ._ ._ _. ! |_) | (_| (_| | | (_| \/ ! _| / -double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max) -double precision array_coefs_B(0:ntot_max,0:ntot_max,0:ntot_max) +double precision, allocatable :: array_coefs_A(:,:,:) +double precision, allocatable :: array_coefs_B(:,:,:) double precision, allocatable :: array_R(:,:,:,:,:) double precision, allocatable :: array_I_A(:,:,:,:,:) double precision, allocatable :: array_I_B(:,:,:,:,:) -!=!=!=!=!=!=!=!=!=! -! A l l o c a t e ! -!=!=!=!=!=!=!=!=!=! - -allocate (array_R(0:ntot_max+nkl_max,kmax_max,0:lmax_max,0:lmax_max+ntot_max,0:lmax_max+ntot_max)) - -allocate (array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) - -allocate (array_I_B(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)) ! _ ! / _. | _ | ! \_ (_| | (_ |_| | ! + +print*, "lmax",lmax +print*, "kmax",kmax +print*, "n_kl",n_kl +print*, "n_a",n_a +print*, "n_b",n_b +print*, "v_kl",v_kl +print*, "dz_kl",dz_kl + + fourpi=4.d0*dacos(-1.d0) ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) @@ -249,8 +245,21 @@ ntotA=n_a(1)+n_a(2)+n_a(3) ntotB=n_b(1)+n_b(2)+n_b(3) ntot=ntotA+ntotB -if(ntot.gt.ntot_max)stop 'increase ntot_max' +nkl_max=4 +!=!=!=!=!=!=!=!=!=! +! A l l o c a t e ! +!=!=!=!=!=!=!=!=!=! +allocate (array_coefs_A(0:ntot,0:ntot,0:ntot)) +allocate (array_coefs_B(0:ntot,0:ntot,0:ntot)) + +allocate (array_R(0:ntot+nkl_max,kmax,0:lmax,0:lmax+ntot,0:lmax+ntot)) + +allocate (array_I_A(0:lmax+ntot,-(lmax+ntot):lmax+ntot,0:ntot,0:ntot,0:ntot)) + +allocate (array_I_B(0:lmax+ntot,-(lmax+ntot):lmax+ntot,0:ntot,0:ntot,0:ntot)) + +print*, ac,bc if(ac.eq.0.d0.and.bc.eq.0.d0)then @@ -584,6 +593,7 @@ endif ! | | | | (_| | | _> (/_ ! deallocate (array_R, array_I_A, array_I_B) + deallocate (array_coefs_A, array_coefs_B) return end @@ -598,15 +608,33 @@ end double precision function Vpseudo_num(npts,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) implicit none -integer kmax_max,lmax_max -parameter (kmax_max=2,lmax_max=2) -integer lmax,kmax, n_kl(kmax_max,0:lmax_max),l,m,k,kk -double precision v_kl(kmax_max,0:lmax_max),dz_kl(kmax_max,0:lmax_max) -double precision a(3),g_a,b(3),g_b,c(3),ac(3),bc(3) -integer n_a(3),n_b(3),npts -double precision rmax,dr,sum,rC + + +! ___ +! | ._ ._ _|_ +! _|_ | | |_) |_| |_ +! | +double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) +integer, intent(in) :: lmax,kmax,npts +integer, intent(in) :: n_a(3),n_b(3), n_kl(kmax,0:lmax) +double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) +double precision, intent(in) :: rmax + +! +! | _ _ _. | _ +! |_ (_) (_ (_| | (/_ +! + +integer :: l,m,k,kk +double precision ac(3),bc(3) +double precision dr,sum,rC double precision overlap_orb_ylm_brute +! _ +! / _. | _ | +! \_ (_| | (_ |_| | +! + do l=1,3 ac(l)=a(l)-c(l) bc(l)=b(l)-c(l) diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index 065785c2..da9f1d68 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -201,7 +201,7 @@ n_kl_dump = n_kl(k,1:kmax,0:lmax) v_kl_dump = v_kl(k,1:kmax,0:lmax) dz_kl_dump = dz_kl(k,1:kmax,0:lmax) - + c = c + Vpseudo(lmax,kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) enddo From 1938425cfed0ef69498bb07a2906def6aee70687 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Tue, 28 Apr 2015 17:17:34 +0200 Subject: [PATCH 25/27] Add a return condition for kmax=1, lmax=1 and v_kl=0.d0 in vpseudo --- scripts/pseudo/put_pseudo_in_ezfio.py | 7 +++++++ src/MonoInts/int.f90 | 20 ++++++++------------ 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 0ba71b0c..68c16729 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -275,6 +275,13 @@ if __name__ == "__main__": # ~#~#~#~#~ # ezfio.pseudo_klocmax = len(v_k[0]) + + # Change all the array 'cause EZFIO + # v_kl (v, l) => v_kl(l,v) + # v_kl => zip(*_v_kl) + # [[7.0, 79.74474797, -49.45159098], [1.0, 5.41040609, -4.60151975]] + # [(7.0, 1.0), (79.74474797, 5.41040609), (-49.45159098, -4.60151975)] + ezfio.pseudo_v_k = zip(*v_k) ezfio.pseudo_n_k = zip(*n_k) ezfio.pseudo_dz_k = zip(*dz_k) diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index 98247553..be806b3f 100644 --- a/src/MonoInts/int.f90 +++ b/src/MonoInts/int.f90 @@ -219,24 +219,21 @@ double precision, allocatable :: array_I_B(:,:,:,:,:) ! \_ (_| | (_ |_| | ! - -print*, "lmax",lmax -print*, "kmax",kmax -print*, "n_kl",n_kl -print*, "n_a",n_a -print*, "n_b",n_b -print*, "v_kl",v_kl -print*, "dz_kl",dz_kl - +if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then + Vpseudo=0.d0 + return +end if fourpi=4.d0*dacos(-1.d0) ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) arg=g_a*ac**2+g_b*bc**2 + if(arg.gt.-dlog(1.d-20))then -Vpseudo=0.d0 -return + Vpseudo=0.d0 + return endif + freal=dexp(-arg) areal=2.d0*g_a*ac @@ -259,7 +256,6 @@ allocate (array_I_A(0:lmax+ntot,-(lmax+ntot):lmax+ntot,0:ntot,0:ntot,0:ntot)) allocate (array_I_B(0:lmax+ntot,-(lmax+ntot):lmax+ntot,0:ntot,0:ntot,0:ntot)) -print*, ac,bc if(ac.eq.0.d0.and.bc.eq.0.d0)then From 3badfe5a7fb0677135d203eb09969e6cd47e616b Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Wed, 29 Apr 2015 14:34:14 +0200 Subject: [PATCH 26/27] Solve put_pseudo_in_ezfio when len(v_kl) is not equal for all the atom --- scripts/pseudo/put_pseudo_in_ezfio.py | 54 +++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 8 deletions(-) diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index 68c16729..6a7aaef7 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -176,6 +176,36 @@ def get_zeff_alpha_beta(str_ele): return [z_eff, alpha, beta] + +def add_zero(array, size, type): + for add in xrange(len(array), size): + array.append([type(0)]) + + return array + + +def make_it_square(matrix, dim, type=float): + """ + matix the matrix to squate + dim array [lmax, kmax] + type the null value you want + [[[28.59107316], [19.37583724]], [[50.25646328]]] + => + [[[28.59107316], [19.37583724]], [[50.25646328], [0.0]]] + """ + + lmax = dim[0] + kmax = dim[1] + + for l_list in matrix: + + l_list = add_zero(l_list, lmax, type) + + for k_list in list_: + k_list = add_zero(k_list, kmax, type) + + return matrix + if __name__ == "__main__": arguments = docopt(__doc__) # ___ @@ -270,18 +300,19 @@ if __name__ == "__main__": ezfio.electrons_elec_alpha_num = alpha_tot ezfio.electrons_elec_beta_num = beta_tot - # ~#~#~#~#~ # - # L o c a l # - # ~#~#~#~#~ # - - ezfio.pseudo_klocmax = len(v_k[0]) - # Change all the array 'cause EZFIO # v_kl (v, l) => v_kl(l,v) # v_kl => zip(*_v_kl) # [[7.0, 79.74474797, -49.45159098], [1.0, 5.41040609, -4.60151975]] # [(7.0, 1.0), (79.74474797, 5.41040609), (-49.45159098, -4.60151975)] + # ~#~#~#~#~ # + # L o c a l # + # ~#~#~#~#~ # + + klocmax = max([len(i) for i in v_k]) + ezfio.pseudo_klocmax = klocmax + ezfio.pseudo_v_k = zip(*v_k) ezfio.pseudo_n_k = zip(*n_k) ezfio.pseudo_dz_k = zip(*dz_k) @@ -290,8 +321,15 @@ if __name__ == "__main__": # N o n _ L o c a l # # ~#~#~#~#~#~#~#~#~ # - ezfio.pseudo_lmaxpo = len(v_kl[0]) - ezfio.pseudo_kmax = len(v_kl[0][0]) + lmax = max([len(i) for i in v_kl]) + kmax = max([len(sublist) for list_ in v_kl for sublist in list_]) + + ezfio.pseudo_lmaxpo = lmax + ezfio.pseudo_kmax = kmax + + v_kl = make_it_square(v_kl, [lmax, kmax]) + n_kl = make_it_square(n_kl, [lmax, kmax], int) + dz_kl = make_it_square(dz_kl, [lmax, kmax]) ezfio.pseudo_v_kl = zip(*v_kl) ezfio.pseudo_n_kl = zip(*n_kl) From 24b4e1b7ca8ab87456c415cac6ef13ba051b5cdb Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Thu, 30 Apr 2015 14:00:43 +0200 Subject: [PATCH 27/27] Move install script into install --- scripts/{ => install}/install_curl.sh | 0 scripts/{ => install}/install_emsl.sh | 0 scripts/{ => install}/install_ezfio.sh | 0 scripts/{ => install}/install_irpf90.sh | 0 scripts/{ => install}/install_m4.sh | 0 scripts/{ => install}/install_ocaml.sh | 0 scripts/{ => install}/install_resultsFile.sh | 0 scripts/{ => install}/install_zlib.sh | 0 scripts/upgrade_ezfio.sh | 2 +- scripts/upgrade_irpf90.sh | 2 +- setup_environment.sh | 16 ++++++++-------- 11 files changed, 10 insertions(+), 10 deletions(-) rename scripts/{ => install}/install_curl.sh (100%) rename scripts/{ => install}/install_emsl.sh (100%) rename scripts/{ => install}/install_ezfio.sh (100%) rename scripts/{ => install}/install_irpf90.sh (100%) rename scripts/{ => install}/install_m4.sh (100%) rename scripts/{ => install}/install_ocaml.sh (100%) rename scripts/{ => install}/install_resultsFile.sh (100%) rename scripts/{ => install}/install_zlib.sh (100%) diff --git a/scripts/install_curl.sh b/scripts/install/install_curl.sh similarity index 100% rename from scripts/install_curl.sh rename to scripts/install/install_curl.sh diff --git a/scripts/install_emsl.sh b/scripts/install/install_emsl.sh similarity index 100% rename from scripts/install_emsl.sh rename to scripts/install/install_emsl.sh diff --git a/scripts/install_ezfio.sh b/scripts/install/install_ezfio.sh similarity index 100% rename from scripts/install_ezfio.sh rename to scripts/install/install_ezfio.sh diff --git a/scripts/install_irpf90.sh b/scripts/install/install_irpf90.sh similarity index 100% rename from scripts/install_irpf90.sh rename to scripts/install/install_irpf90.sh diff --git a/scripts/install_m4.sh b/scripts/install/install_m4.sh similarity index 100% rename from scripts/install_m4.sh rename to scripts/install/install_m4.sh diff --git a/scripts/install_ocaml.sh b/scripts/install/install_ocaml.sh similarity index 100% rename from scripts/install_ocaml.sh rename to scripts/install/install_ocaml.sh diff --git a/scripts/install_resultsFile.sh b/scripts/install/install_resultsFile.sh similarity index 100% rename from scripts/install_resultsFile.sh rename to scripts/install/install_resultsFile.sh diff --git a/scripts/install_zlib.sh b/scripts/install/install_zlib.sh similarity index 100% rename from scripts/install_zlib.sh rename to scripts/install/install_zlib.sh diff --git a/scripts/upgrade_ezfio.sh b/scripts/upgrade_ezfio.sh index 4a9af403..c35a2dbd 100755 --- a/scripts/upgrade_ezfio.sh +++ b/scripts/upgrade_ezfio.sh @@ -12,7 +12,7 @@ fi cd -- ${QPACKAGE_ROOT} mv -- ${QPACKAGE_ROOT}/EZFIO ${QPACKAGE_ROOT}/EZFIO.old -${QPACKAGE_ROOT}/scripts/install_ezfio.sh +${QPACKAGE_ROOT}/scripts/install/install_ezfio.sh if [[ $? -eq 0 ]] then diff --git a/scripts/upgrade_irpf90.sh b/scripts/upgrade_irpf90.sh index 5735754f..dea48014 100755 --- a/scripts/upgrade_irpf90.sh +++ b/scripts/upgrade_irpf90.sh @@ -12,7 +12,7 @@ fi cd -- ${QPACKAGE_ROOT} mv -f -- ${QPACKAGE_ROOT}/irpf90 ${QPACKAGE_ROOT}/irpf90.old -${QPACKAGE_ROOT}/scripts/install_irpf90.sh +${QPACKAGE_ROOT}/scripts/install/install_irpf90.sh if [[ $? -eq 0 ]] then diff --git a/setup_environment.sh b/setup_environment.sh index 12fcf4a5..72d2834c 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -30,7 +30,7 @@ EOF source quantum_package.rc echo "${BLUE}===== Installing IRPF90 ===== ${BLACK}" -${QPACKAGE_ROOT}/scripts/install_irpf90.sh | tee install_irpf90.log +${QPACKAGE_ROOT}/scripts/install/install_irpf90.sh | tee install_irpf90.log if [[ ! -d ${QPACKAGE_ROOT}/irpf90 ]] then echo $RED "Error in IRPF90 installation" $BLACK @@ -51,16 +51,16 @@ fi echo "${BLUE}===== Installing Zlib ===== ${BLACK}" -${QPACKAGE_ROOT}/scripts/install_zlib.sh | tee install_zlib.log +${QPACKAGE_ROOT}/scripts/install/install_zlib.sh | tee install_zlib.log echo "${BLUE}===== Installing Curl ===== ${BLACK}" -${QPACKAGE_ROOT}/scripts/install_curl.sh | tee install_curl.log +${QPACKAGE_ROOT}/scripts/install/install_curl.sh | tee install_curl.log echo "${BLUE}===== Installing M4 ===== ${BLACK}" -${QPACKAGE_ROOT}/scripts/install_m4.sh | tee install_m4.log +${QPACKAGE_ROOT}/scripts/install/install_m4.sh | tee install_m4.log echo "${BLUE}===== Installing EMSL Basis set library ===== ${BLACK}" -${QPACKAGE_ROOT}/scripts/install_emsl.sh | tee install_emsl.log +${QPACKAGE_ROOT}/scripts/install/install_emsl.sh | tee install_emsl.log if [[ ! -d ${QPACKAGE_ROOT}/EMSL_Basis ]] then @@ -70,7 +70,7 @@ fi echo "${BLUE}===== Installing EZFIO ===== ${BLACK}" -${QPACKAGE_ROOT}/scripts/install_ezfio.sh | tee install_ezfio.log +${QPACKAGE_ROOT}/scripts/install/install_ezfio.sh | tee install_ezfio.log if [[ ! -d ${QPACKAGE_ROOT}/EZFIO ]] then echo $RED "Error in EZFIO installation" $BLACK @@ -80,7 +80,7 @@ fi echo "${BLUE}===== Installing Ocaml compiler and libraries ===== ${BLACK}" rm -f -- ocaml/Qptypes.ml -${QPACKAGE_ROOT}/scripts/install_ocaml.sh | tee install_ocaml.log +${QPACKAGE_ROOT}/scripts/install/install_ocaml.sh | tee install_ocaml.log if [[ ! -f ${QPACKAGE_ROOT}/ocaml/Qptypes.ml ]] then @@ -89,7 +89,7 @@ then fi echo "${BLUE}===== Installing resultsFile Python library ===== ${BLACK}" -${QPACKAGE_ROOT}/scripts/install_resultsFile.sh +${QPACKAGE_ROOT}/scripts/install/install_resultsFile.sh if [[ ! -d ${QPACKAGE_ROOT}/resultsFile ]] then echo $RED "Error in resultsFile installation" $BLACK