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