diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py index 9eb781cb..4cba032f 100755 --- a/scripts/ezfio_interface/ei_handler.py +++ b/scripts/ezfio_interface/ei_handler.py @@ -368,7 +368,10 @@ def create_ezfio_stuff(dict_ezfio_cfg, config_or_default="config"): except ValueError: a_size_raw.append(dim) else: - a_size_raw.append("{0}+{1}+1".format(begin, end)) + if begin[0] == '-': + a_size_raw.append("{0}+{1}+1".format(end, begin[1:])) + else: + a_size_raw.append("{0}-{1}+1".format(end, begin)) size_raw = ",".join(a_size_raw) diff --git a/src/AOs/tree_dependancy.png b/src/AOs/tree_dependancy.png index 41a281a5..35cbcbc3 100644 Binary files a/src/AOs/tree_dependancy.png and b/src/AOs/tree_dependancy.png differ diff --git a/src/Bitmask/tree_dependancy.png b/src/Bitmask/tree_dependancy.png index 0eb69ca7..dd1c0f54 100644 Binary files a/src/Bitmask/tree_dependancy.png and b/src/Bitmask/tree_dependancy.png differ diff --git a/src/CAS_SD/tree_dependancy.png b/src/CAS_SD/tree_dependancy.png index bd68d477..aae463ee 100644 Binary files a/src/CAS_SD/tree_dependancy.png and b/src/CAS_SD/tree_dependancy.png differ diff --git a/src/CID/tree_dependancy.png b/src/CID/tree_dependancy.png index 0e95f33d..aee0843d 100644 Binary files a/src/CID/tree_dependancy.png and b/src/CID/tree_dependancy.png differ diff --git a/src/CID_SC2_selected/tree_dependancy.png b/src/CID_SC2_selected/tree_dependancy.png index 688bd413..1b24cd17 100644 Binary files a/src/CID_SC2_selected/tree_dependancy.png and b/src/CID_SC2_selected/tree_dependancy.png differ diff --git a/src/CID_selected/tree_dependancy.png b/src/CID_selected/tree_dependancy.png index 01673a75..f6c4cc3e 100644 Binary files a/src/CID_selected/tree_dependancy.png and b/src/CID_selected/tree_dependancy.png differ diff --git a/src/CIS/tree_dependancy.png b/src/CIS/tree_dependancy.png index 5ca9a22a..87b2e1cd 100644 Binary files a/src/CIS/tree_dependancy.png and b/src/CIS/tree_dependancy.png differ diff --git a/src/CISD/tree_dependancy.png b/src/CISD/tree_dependancy.png index 31a6d3a6..2be99e68 100644 Binary files a/src/CISD/tree_dependancy.png and b/src/CISD/tree_dependancy.png differ diff --git a/src/CISD_SC2_selected/tree_dependancy.png b/src/CISD_SC2_selected/tree_dependancy.png index 3f67037a..8804d41f 100644 Binary files a/src/CISD_SC2_selected/tree_dependancy.png and b/src/CISD_SC2_selected/tree_dependancy.png differ diff --git a/src/CISD_selected/tree_dependancy.png b/src/CISD_selected/tree_dependancy.png index a6c1a604..282999d6 100644 Binary files a/src/CISD_selected/tree_dependancy.png and b/src/CISD_selected/tree_dependancy.png differ diff --git a/src/DDCI_selected/tree_dependancy.png b/src/DDCI_selected/tree_dependancy.png index 5cba3510..4ca5ff35 100644 Binary files a/src/DDCI_selected/tree_dependancy.png and b/src/DDCI_selected/tree_dependancy.png differ diff --git a/src/DensityFit/ASSUMPTIONS.rst b/src/DensityFit/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/DensityFit/Makefile b/src/DensityFit/Makefile new file mode 100644 index 00000000..06dc50ff --- /dev/null +++ b/src/DensityFit/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= +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/DensityFit/NEEDED_CHILDREN_MODULES b/src/DensityFit/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..13e425cc --- /dev/null +++ b/src/DensityFit/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +AOs diff --git a/src/DensityFit/README.rst b/src/DensityFit/README.rst new file mode 100644 index 00000000..23d6c06f --- /dev/null +++ b/src/DensityFit/README.rst @@ -0,0 +1,24 @@ +================= +DensityFit Module +================= + +In this module, the basis of all the products of atomic orbitals is built. + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +.. image:: tree_dependancy.png + +* `AOs `_ + diff --git a/src/DensityFit/aux_basis.irp.f b/src/DensityFit/aux_basis.irp.f new file mode 100644 index 00000000..1f28d985 --- /dev/null +++ b/src/DensityFit/aux_basis.irp.f @@ -0,0 +1,54 @@ + BEGIN_PROVIDER [ integer, aux_basis_num ] +&BEGIN_PROVIDER [ integer, aux_basis_num_8 ] + implicit none + BEGIN_DOC + ! Number of auxiliary basis functions + END_DOC + integer :: align_double + aux_basis_num = ao_num * (ao_num+1)/2 + aux_basis_num_8 = align_double(aux_basis_num) +END_PROVIDER + +BEGIN_PROVIDER [ integer, aux_basis_idx, (2,aux_basis_num) ] + implicit none + BEGIN_DOC +! aux_basis_idx(k) -> i,j + END_DOC + integer :: i,j,k + k=0 + do j=1,ao_num + do i=1,j + k = k+1 + aux_basis_idx(1,k) = i + aux_basis_idx(2,k) = j + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, aux_basis_overlap_matrix, (aux_basis_num_8,aux_basis_num) ] + implicit none + BEGIN_DOC +! Auxiliary basis set + END_DOC + integer :: m,n,i,j,k,l + double precision :: ao_four_overlap + + aux_basis_overlap_matrix(1,1) = ao_four_overlap(1,1,1,1) + !$OMP PARALLEL DO PRIVATE(i,j,k,l,m,n) SCHEDULE(GUIDED) + do m=1,aux_basis_num + i = aux_basis_idx(1,m) + j = aux_basis_idx(2,m) + do n=1,m + k = aux_basis_idx(1,n) + l = aux_basis_idx(2,n) + aux_basis_overlap_matrix(m,n) = ao_four_overlap(i,j,k,l) + aux_basis_overlap_matrix(n,m) = aux_basis_overlap_matrix(m,n) + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + + diff --git a/src/DensityFit/overlap.irp.f b/src/DensityFit/overlap.irp.f new file mode 100644 index 00000000..0394dbf9 --- /dev/null +++ b/src/DensityFit/overlap.irp.f @@ -0,0 +1,66 @@ +double precision function ao_four_overlap(i,j,k,l) + implicit none + BEGIN_DOC +! \int \chi_i(r) \chi_j(r) \chi_k(r) \chi_l(r) dr + END_DOC + integer,intent(in) :: i,j,k,l + integer :: p,q,r,s + double precision :: I_center(3),J_center(3),K_center(3),L_center(3) + integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3) + double precision :: overlap_x,overlap_y,overlap_z, overlap + include 'include/constants.F' + double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp + double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq + integer :: iorder_p(3), iorder_q(3) + + dim1 = n_pt_max_integrals + + num_i = ao_nucl(i) + num_j = ao_nucl(j) + num_k = ao_nucl(k) + num_l = ao_nucl(l) + ao_four_overlap = 0.d0 + + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + I_center(p) = nucl_coord(num_i,p) + J_center(p) = nucl_coord(num_j,p) + K_center(p) = nucl_coord(num_k,p) + L_center(p) = nucl_coord(num_l,p) + enddo + + do p = 1, ao_prim_num(i) + double precision :: coef1 + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & + I_power,J_power,I_center,J_center,dim1) + double precision :: coef2 + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)*fact_p + do r = 1, ao_prim_num(k) + double precision :: coef3 + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + double precision :: general_primitive_integral + call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& + ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & + K_power,L_power,K_center,L_center,dim1) + double precision :: coef4 + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)*fact_q + call overlap_gaussian_xyz(P_center,Q_center,pp,qq,iorder_p,iorder_q,overlap_x,overlap_y,overlap_z,overlap,dim1) + ao_four_overlap += coef4 * overlap + enddo ! s + enddo ! r + enddo ! q + enddo ! p + +end + +! TODO : Schwartz acceleration + + + diff --git a/src/Determinants/tree_dependancy.png b/src/Determinants/tree_dependancy.png index 1f56d1ee..fb18c1cd 100644 Binary files a/src/Determinants/tree_dependancy.png and b/src/Determinants/tree_dependancy.png differ diff --git a/src/Electrons/tree_dependancy.png b/src/Electrons/tree_dependancy.png index 8c7bb855..b90a1f83 100644 Binary files a/src/Electrons/tree_dependancy.png and b/src/Electrons/tree_dependancy.png differ diff --git a/src/Ezfio_files/tree_dependancy.png b/src/Ezfio_files/tree_dependancy.png index 4b0b9fd0..48f53991 100644 Binary files a/src/Ezfio_files/tree_dependancy.png and b/src/Ezfio_files/tree_dependancy.png differ diff --git a/src/FCIdump/README.rst b/src/FCIdump/README.rst index 308982d5..1f25ca6e 100644 --- a/src/FCIdump/README.rst +++ b/src/FCIdump/README.rst @@ -10,9 +10,6 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fcidump `_ - Undocumented - Needed Modules diff --git a/src/FCIdump/tree_dependancy.png b/src/FCIdump/tree_dependancy.png index e12a69a9..80c25a57 100644 Binary files a/src/FCIdump/tree_dependancy.png and b/src/FCIdump/tree_dependancy.png differ diff --git a/src/Full_CI/tree_dependancy.png b/src/Full_CI/tree_dependancy.png index 237d40ab..1297595b 100644 Binary files a/src/Full_CI/tree_dependancy.png and b/src/Full_CI/tree_dependancy.png differ diff --git a/src/Generators_CAS/tree_dependancy.png b/src/Generators_CAS/tree_dependancy.png index 348e50ca..2e2571d9 100644 Binary files a/src/Generators_CAS/tree_dependancy.png and b/src/Generators_CAS/tree_dependancy.png differ diff --git a/src/Generators_full/tree_dependancy.png b/src/Generators_full/tree_dependancy.png index 08fe461a..659db178 100644 Binary files a/src/Generators_full/tree_dependancy.png and b/src/Generators_full/tree_dependancy.png differ diff --git a/src/Hartree_Fock/tree_dependancy.png b/src/Hartree_Fock/tree_dependancy.png index 3ba1436b..aa20a902 100644 Binary files a/src/Hartree_Fock/tree_dependancy.png and b/src/Hartree_Fock/tree_dependancy.png differ diff --git a/src/Integrals_Bielec/tree_dependancy.png b/src/Integrals_Bielec/tree_dependancy.png index 8ecebb3d..e9430bf0 100644 Binary files a/src/Integrals_Bielec/tree_dependancy.png and b/src/Integrals_Bielec/tree_dependancy.png differ diff --git a/src/Integrals_Monoelec/README.rst b/src/Integrals_Monoelec/README.rst index 7ba2970f..775f45d2 100644 --- a/src/Integrals_Monoelec/README.rst +++ b/src/Integrals_Monoelec/README.rst @@ -103,6 +103,12 @@ Documentation `ao_pseudo_integral_non_local `_ Local pseudo-potential +`pseudo_grid `_ + Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C + .br + + `mo_nucl_elec_integral `_ interaction nuclear electron on the MO basis diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f index 093b9fc8..63828c8c 100644 --- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)] ! Pseudo-potential END_DOC if (do_pseudo) then - ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local + ao_pseudo_integral = ao_pseudo_integral_local !+ ao_pseudo_integral_non_local else ao_pseudo_integral = 0.d0 endif @@ -220,4 +220,67 @@ END_PROVIDER END_PROVIDER +BEGIN_PROVIDER [ double precision, pseudo_grid, (pseudo_grid_size,ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num) ] + implicit none + BEGIN_DOC +! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C +! +! + END_DOC + ! l,m : Y(l,m) parameters + ! c(3) : pseudopotential center + ! a(3) : Atomic Orbital center + ! n_a(3) : Powers of x,y,z in the Atomic Orbital + ! g_a : Atomic Orbital exponent + ! r : Distance between the Atomic Orbital center and the considered point + double precision, external :: ylm_orb + integer :: n_a(3) + double precision :: a(3), c(3), g_a + integer :: i,j,k,l,m,n,p + double precision :: r(pseudo_grid_size), dr, Ulc + double precision, parameter :: rmax= 10.d0 + + dr = rmax/dble(pseudo_grid_size) + r(1) = 0.d0 + do j=2,pseudo_grid_size + r(j) = r(j-1) + dr + enddo + + pseudo_grid = 0.d0 + do k=1,nucl_num + c(1:3) = nucl_coord(k,1:3) + do l=0,pseudo_lmax + do i=1,ao_num + a(1:3) = nucl_coord(ao_nucl(i),1:3) + n_a(1:3) = ao_power(i,1:3) + do j=1,pseudo_grid_size + do p=1,ao_prim_num(i) + g_a = ao_expo_ordered_transp(p,i) + do m=-l,l + double precision :: y + ! y = ylm_orb(l,m,c,a,n_a,g_a,r(j)) + ! if (y /= 0.d0) then + ! print *, 'y = ', y + ! print *, 'l = ', l + ! print *, 'm = ', m + ! print *, 'c = ', c + ! print *, 'a = ', a + ! print *, 'n = ', n_a + ! print *, 'g = ', g_a + ! print *, 'r = ', r(j) + ! print *, '' + ! endif + y = ylm_orb(l,m,c,a,n_a,g_a,r(j)) + pseudo_grid(j,i,m,l,k) = pseudo_grid(j,i,m,l,k) + & + ao_coef_normalized_ordered_transp(p,i)*y + enddo + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f index b6851f21..d34a05ff 100644 --- a/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f @@ -28,5 +28,7 @@ BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_n enddo enddo !$OMP END PARALLEL DO + call ezfio_set_pseudo_pseudo_matrix(mo_pseudo_integral(1:mo_tot_num,1:mo_tot_num)) END_PROVIDER + diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90 index 0cea5f66..d8ff7893 100644 --- a/src/Integrals_Monoelec/pseudopot.f90 +++ b/src/Integrals_Monoelec/pseudopot.f90 @@ -201,7 +201,7 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) ! |_ (_) (_ (_| | (/_ ! -double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm +double precision :: fourpi,f,prod,prodp,binom_func,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 @@ -327,7 +327,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then 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) & + array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(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 @@ -336,7 +336,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then 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) & + array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(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 @@ -448,7 +448,7 @@ else if(ac.eq.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) & + array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(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 @@ -534,7 +534,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then 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) & + array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(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 @@ -705,7 +705,7 @@ 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 +double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg fourpi=4.d0*dacos(-1.d0) f=fourpi**1.5d0 @@ -762,9 +762,9 @@ allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:n 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) & + array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(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) & + *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(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 @@ -823,7 +823,7 @@ 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 +double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm pi=dacos(-1.d0) if(mu.gt.0.and.m.gt.0)then @@ -834,8 +834,8 @@ 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) + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(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 @@ -868,7 +868,7 @@ 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) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(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 @@ -884,7 +884,7 @@ 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) + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(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 @@ -904,8 +904,8 @@ 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) + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(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 @@ -926,7 +926,7 @@ 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) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(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 @@ -944,7 +944,7 @@ 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) + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(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 @@ -964,8 +964,8 @@ 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) + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(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 @@ -985,8 +985,8 @@ 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) + cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(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 @@ -1485,6 +1485,7 @@ end a=bessel_mod_exp(n,x) return endif + print *, n,x 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) @@ -1699,14 +1700,14 @@ end double precision function coef_pm(n,k) implicit none integer n,k -double precision arg,binom,binom_gen +double precision arg,binom_func,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) +coef_pm=2.d0**n*binom_func(n,k)*binom_gen(arg,n) end !! Ylm_bis uses the series expansion of Ylm in xchap^i ychap^j zchap^k @@ -1735,7 +1736,7 @@ end 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 +double precision x,y,z,theta,phi,sum,factor,pi,binom_func,fact,coef_pm,cylm pi=dacos(-1.d0) x=dsin(theta)*dcos(phi) y=dsin(theta)*dsin(phi) @@ -1745,7 +1746,7 @@ 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) + cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(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 @@ -1765,7 +1766,7 @@ 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) + cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(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 @@ -1800,13 +1801,6 @@ end !! !! 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 @@ -2087,4 +2081,90 @@ end end - +! l,m : Y(l,m) parameters +! c(3) : pseudopotential center +! a(3) : Atomic Orbital center +! n_a(3) : Powers of x,y,z in the Atomic Orbital +! g_a : Atomic Orbital exponent +! r : Distance between the Atomic Orbital center and the considered point +double precision function ylm_orb(l,m,c,a,n_a,g_a,r) +implicit none +integer lmax_max,ntot_max +parameter (lmax_max=2) +parameter (ntot_max=14) +integer l,m +double precision a(3),g_a,c(3) +double precision prod,binom_func,accu,bigI,ylm,bessel_mod +double precision theta_AC0,phi_AC0,ac,factor,fourpi,arg,r,areal +integer ntotA,mu,k1,k2,k3,lambda +integer n_a(3) +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_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max), y + +ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) +arg=g_a*(ac**2+r**2) +fourpi=4.d0*dacos(-1.d0) +factor=fourpi*dexp(-arg) +areal=2.d0*g_a*ac +ntotA=n_a(1)+n_a(2)+n_a(3) + +if(ntotA.gt.ntot_max)stop 'increase ntot_max' + +if(ac.eq.0.d0)then + ylm_orb=dsqrt(fourpi)*r**ntotA*dexp(-g_a*r**2)*bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) + return +else + + theta_AC0=dacos( (a(3)-c(3))/ac ) + phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) + + 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_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(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) & + *r**(k1+k2+k3) + 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 + + accu=0.d0 + do lambda=0,l+ntotA + do mu=-lambda,lambda + y = ylm(lambda,mu,theta_AC0,phi_AC0) + if (y == 0.d0) then + cycle + endif + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + prod=y*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3) + if (prod == 0.d0) then + cycle + endif + if (areal*r < 100.d0) then ! overflow! + accu=accu+prod*bessel_mod(areal*r,lambda) + endif + enddo + enddo + enddo + enddo + enddo + ylm_orb=factor*accu + return +endif + +end diff --git a/src/Integrals_Monoelec/tree_dependancy.png b/src/Integrals_Monoelec/tree_dependancy.png index 628d4ed4..490d33b5 100644 Binary files a/src/Integrals_Monoelec/tree_dependancy.png and b/src/Integrals_Monoelec/tree_dependancy.png differ diff --git a/src/MOGuess/tree_dependancy.png b/src/MOGuess/tree_dependancy.png index 40c35021..14da5742 100644 Binary files a/src/MOGuess/tree_dependancy.png and b/src/MOGuess/tree_dependancy.png differ diff --git a/src/MOs/tree_dependancy.png b/src/MOs/tree_dependancy.png index 276be274..7250aa3f 100644 Binary files a/src/MOs/tree_dependancy.png and b/src/MOs/tree_dependancy.png differ diff --git a/src/MP2/tree_dependancy.png b/src/MP2/tree_dependancy.png index 9fb6376d..f0573cdf 100644 Binary files a/src/MP2/tree_dependancy.png and b/src/MP2/tree_dependancy.png differ diff --git a/src/MRCC/tree_dependancy.png b/src/MRCC/tree_dependancy.png index d2a82da4..9b5920af 100644 Binary files a/src/MRCC/tree_dependancy.png and b/src/MRCC/tree_dependancy.png differ diff --git a/src/Molden/tree_dependancy.png b/src/Molden/tree_dependancy.png index fd419b86..887372a5 100644 Binary files a/src/Molden/tree_dependancy.png and b/src/Molden/tree_dependancy.png differ diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 05dd72a9..b8ea4e38 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Determinants Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils \ No newline at end of file +AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Determinants Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils DensityFit diff --git a/src/Nuclei/tree_dependancy.png b/src/Nuclei/tree_dependancy.png index 91922cea..dfd6f57a 100644 Binary files a/src/Nuclei/tree_dependancy.png and b/src/Nuclei/tree_dependancy.png differ diff --git a/src/Pseudo/EZFIO.cfg b/src/Pseudo/EZFIO.cfg index 8a6afa2c..5a401d18 100644 --- a/src/Pseudo/EZFIO.cfg +++ b/src/Pseudo/EZFIO.cfg @@ -55,3 +55,22 @@ doc: Using pseudo potential integral of not interface: input default: False +[pseudo_grid_size] +type: integer +doc: Size of the QMC grid +interface: input +default: 100 + +[pseudo_grid] +type: double precision +doc: QMC grid +interface: output +size: (pseudo.pseudo_grid_size,ao_basis.ao_num,-pseudo.pseudo_lmax:pseudo.pseudo_lmax,0:pseudo.pseudo_lmax,nuclei.nucl_num) + +[pseudo_matrix] +type: double precision +doc: QMC grid +interface: output +size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num) + + diff --git a/src/Pseudo/tree_dependancy.png b/src/Pseudo/tree_dependancy.png index ce0dbd5b..389cb115 100644 Binary files a/src/Pseudo/tree_dependancy.png and b/src/Pseudo/tree_dependancy.png differ diff --git a/src/Selectors_full/tree_dependancy.png b/src/Selectors_full/tree_dependancy.png index cff76ed3..f88049bd 100644 Binary files a/src/Selectors_full/tree_dependancy.png and b/src/Selectors_full/tree_dependancy.png differ diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index ec1e4200..d356454b 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -131,13 +131,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,& integer :: iorder_p(3) call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) - if(fact_p.lt.1d-10)then - overlap_x = 0.d0 - overlap_y = 0.d0 - overlap_z = 0.d0 - overlap = 0.d0 - return - endif +! if(fact_p.lt.1d-20)then +! overlap_x = 0.d0 +! overlap_y = 0.d0 +! overlap_z = 0.d0 +! overlap = 0.d0 +! return +! endif integer :: nmax double precision :: F_integral nmax = maxval(iorder_p) diff --git a/src/Utils/tree_dependancy.png b/src/Utils/tree_dependancy.png index 68f0c2ca..38b04785 100644 Binary files a/src/Utils/tree_dependancy.png and b/src/Utils/tree_dependancy.png differ