LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-25 17:54:44 +02:00
Code Releases Activity

Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

Browse Source
This commit is contained in:
Anthony Scemama 2025-03-14 18:11:46 +01:00
commit aa70f6fada
15 changed files with 369 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/irpf90 vendored

@ -1 +1 @@
Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b
Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6

4
plugins/local/extra_basis_int/ao_overlap.irp.f
View File

@ -4,7 +4,7 @@
BEGIN_PROVIDER [double precision, ao_extra_overlap , (ao_extra_num, ao_extra_num)]
BEGIN_DOC
! Overlap between atomic basis functions:
! Overlap between atomic basis functions belonging to the EXTRA BASIS
!
! :math:`\int \chi_i(r) \chi_j(r) dr`
END_DOC
@ -71,6 +71,8 @@ END_PROVIDER
BEGIN_DOC
! Overlap between atomic basis functions:
!
! first index belongs to the REGULAR AO basis, second to the EXTRA basis
!
! <AO_i|AO_j extra basis>
END_DOC

10
plugins/local/extra_basis_int/coul_1s.irp.f
View File

@ -11,6 +11,11 @@ double precision function coul_full_ao_pq_r_1s(p,q,R,R_p,R_q)
double precision, intent(in) :: R(3),R_p(3),R_q(3)
integer, intent(in) :: p,q
double precision :: coef,dist,P_pq(3),coefaos
if(.not.ao_extra_only_1s)then
print*,'You are using a function assuming that the extra basis is fitted on 1s functions'
print*,'But this is not the case apparently ... stopping'
stop
endif
coefaos= ao_extra_coef_normalized(p,1) * ao_extra_coef_normalized(q,1)
coef = inv_pi_gamma_pq_3_2_ao_extra(p,q) * E_pq_ao_extra(p,q)
P_pq = ao_extra_expo(p,1) * R_p + ao_extra_expo(q,1) * R_q
@ -40,6 +45,11 @@ double precision function coul_pq_r_1s(p,q,R,R_p,R_q)
double precision, intent(in) :: R(3),R_p(3),R_q(3)
integer, intent(in) :: p,q
double precision :: dist,P_pq(3)
if(.not.ao_extra_only_1s)then
print*,'You are using a function assuming that the extra basis is fitted on 1s functions'
print*,'But this is not the case apparently ... stopping'
stop
endif
P_pq = ao_extra_expo(p,1) * R_p + ao_extra_expo(q,1) * R_q
P_pq = P_pq * inv_gamma_pq_ao_extra(q,p)
dist = (P_pq(1)-R(1)) * (P_pq(1)-R(1))

56
plugins/local/extra_basis_int/extra_basis_int.irp.f
View File

@ -11,18 +11,30 @@ program extra_basis_int
! call routine_pot_ne
! call routine_test_pot_ne_extra_mixed
! call routine_test_coul_1s
call print_v_ne_extra_basis
call print_v_ne_basis
! call print_v_ne_extra_basis
! call print_v_ne_basis
! call test_v_ne_a_extra_basis
! call print_v_ee_mixed_direct
call print_v_ee_mixed_exchange
end
subroutine test_v_ne_a_extra_basis
implicit none
integer :: i,j
do i = 1, ao_extra_num
write(*,'(100(F16.10,X))')pot_vne_A_extra_basis(1:ao_extra_num,i)
enddo
end
subroutine test_overlap
implicit none
integer :: i,j
do i = 1, ao_extra_num
do j = 1, ao_extra_num
write(33,*)ao_extra_overlap(j,i)
enddo
do i = 1, ao_num
! do j = 1, ao_num
write(33,'(100(F16.10,X))')ao_extra_overlap_mixed(i,1:ao_extra_num)
! enddo
enddo
end
@ -189,3 +201,35 @@ subroutine print_v_ne_basis
print*,'accu = ',accu
end
subroutine print_v_ee_mixed_direct
implicit none
integer :: i,j,k,l
double precision :: ao_two_e_integral_mixed_direct
do i = 1, ao_num
do j = 1, ao_num
do k = 1, ao_extra_num
do l = 1, ao_extra_num
write(34,*)ao_two_e_integral_mixed_direct(i, j, k, l)
enddo
enddo
enddo
enddo
end
subroutine print_v_ee_mixed_exchange
implicit none
integer :: i,j,k,l
double precision :: ao_two_e_integral_mixed_exchange
do i = 1, ao_num
do j = 1, ao_extra_num
do k = 1, ao_num
do l = 1, ao_extra_num
write(34,*)ao_two_e_integral_mixed_exchange(i, j, k, l)
enddo
enddo
enddo
enddo
end

50
plugins/local/extra_basis_int/pot_ao_ints.irp.f
View File

@ -1,3 +1,53 @@
BEGIN_PROVIDER [ double precision, pot_vne_A_extra_basis, (ao_extra_num,ao_extra_num)]
implicit none
BEGIN_DOC
!
! Computes the following integral :
! $\sum_{R in the USUAL nuclei} -Z <chi_i|1/|r-R||chi_j>$
!
! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis
END_DOC
integer :: mu,nu
double precision :: v_nucl_extra_ao
pot_vne_A_extra_basis = 0.d0
do mu = 1, ao_extra_num
do nu = 1, ao_extra_num
pot_vne_A_extra_basis(nu,mu)= v_nucl_extra_ao(mu,nu)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, pot_vne_extra_basis, (ao_num,ao_num)]
implicit none
BEGIN_DOC
!
! Computes the following integral :
! $\sum_{R in EXTRA nuclei} -Z <chi_i|1/|r-R||chi_j>$
!
!
! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the USUAL basis
END_DOC
integer :: mu,nu,k_nucl
double precision :: mu_in, R_nucl(3),charge_nucl, integral
double precision :: NAI_pol_mult_erf_ao
mu_in = 10.d0**10
pot_vne_extra_basis = 0.d0
do mu = 1, ao_num
do nu = 1, ao_num
do k_nucl = 1, extra_nucl_num
R_nucl(1:3) = extra_nucl_coord_transp(1:3,k_nucl)
charge_nucl = extra_nucl_charge(k_nucl)
integral = NAI_pol_mult_erf_ao(mu, nu, mu_in, R_nucl)
pot_vne_extra_basis(nu,mu) += -integral * charge_nucl
enddo
enddo
enddo
END_PROVIDER
double precision function NAI_pol_mult_erf_ao_extra(i_ao, j_ao, mu_in, C_center)

36
plugins/local/extra_basis_int/ref_extra_basis.irp.f
View File

@ -6,7 +6,9 @@ program pouet
! call routine_pot_ne_extra
! call ref_pot_ne_mixed
! call ref_pot_ne
call ref_pot_ne_extra_mixed
! call ref_pot_ne_extra_mixed
! call ref_v_ee_mixed_direct
call ref_v_ee_mixed_exchange
end
@ -113,3 +115,35 @@ subroutine ref_pot_ne_extra_mixed
enddo
enddo
end
subroutine ref_v_ee_mixed_direct
implicit none
integer :: i,j,k,l
double precision :: ao_two_e_integral
do i = 1, 15
do j = 1, 15
do k = 16, ao_num
do l = 16, ao_num
write(33,*)ao_two_e_integral(i, j, k, l)
enddo
enddo
enddo
enddo
end
subroutine ref_v_ee_mixed_exchange
implicit none
integer :: i,j,k,l
double precision :: ao_two_e_integral
do i = 1, 15
do j = 16, ao_num
do k = 1, 15
do l = 16, ao_num
write(33,*)ao_two_e_integral(i, j, k, l)
enddo
enddo
enddo
enddo
end

145
plugins/local/extra_basis_int/two_e_int.irp.f Normal file
View File

@ -0,0 +1,145 @@
double precision function ao_two_e_integral_mixed_direct(i, j, k, l)
BEGIN_DOC
! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2)
! A A B B
!
! where i,j belong to the REGULAR AO basis (system A) and k,l to the EXTRA basis (system B)
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: i, j, k, l
integer :: p, q, r, s
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
integer :: iorder_p(3), iorder_q(3)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: integral
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
double precision :: general_primitive_integral
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_extra_nucl(k)
num_l = ao_extra_nucl(l)
ao_two_e_integral_mixed_direct = 0.d0
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_extra_power(k,p)
L_power(p) = ao_extra_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = nucl_coord(num_j,p)
K_center(p) = extra_nucl_coord(num_k,p)
L_center(p) = extra_nucl_coord(num_l,p)
enddo
double precision :: coef1, coef2, coef3, coef4
double precision :: p_inv,q_inv
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_prim_num(j)
coef2 = coef1*ao_coef_normalized_ordered_transp(q,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)
p_inv = 1.d0/pp
do r = 1, ao_extra_prim_num(k)
coef3 = coef2*ao_extra_coef_normalized_ordered_transp(r,k)
do s = 1, ao_extra_prim_num(l)
coef4 = coef3*ao_extra_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_extra_expo_ordered_transp(r,k),ao_extra_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = general_primitive_integral(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_two_e_integral_mixed_direct = ao_two_e_integral_mixed_direct + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
end
double precision function ao_two_e_integral_mixed_exchange(i, j, k, l)
BEGIN_DOC
! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2)
! A B A B
!
! where i,k belong to the REGULAR AO basis (system A) and j,l to the EXTRA basis (system B)
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: i, j, k, l
integer :: p, q, r, s
integer :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
integer :: iorder_p(3), iorder_q(3)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: integral
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
double precision :: general_primitive_integral
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_extra_nucl(j)
num_k = ao_nucl(k)
num_l = ao_extra_nucl(l)
ao_two_e_integral_mixed_exchange = 0.d0
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_extra_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_extra_power(l,p)
I_center(p) = nucl_coord(num_i,p)
J_center(p) = extra_nucl_coord(num_j,p)
K_center(p) = nucl_coord(num_k,p)
L_center(p) = extra_nucl_coord(num_l,p)
enddo
double precision :: coef1, coef2, coef3, coef4
double precision :: p_inv,q_inv
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p,i)
do q = 1, ao_extra_prim_num(j)
coef2 = coef1*ao_extra_coef_normalized_ordered_transp(q,j)
call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,&
ao_expo_ordered_transp(p,i),ao_extra_expo_ordered_transp(q,j), &
I_power,J_power,I_center,J_center,dim1)
p_inv = 1.d0/pp
do r = 1, ao_prim_num(k)
coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
do s = 1, ao_extra_prim_num(l)
coef4 = coef3*ao_extra_coef_normalized_ordered_transp(s,l)
call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,&
ao_expo_ordered_transp(r,k),ao_extra_expo_ordered_transp(s,l), &
K_power,L_power,K_center,L_center,dim1)
q_inv = 1.d0/qq
integral = general_primitive_integral(dim1, &
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
ao_two_e_integral_mixed_exchange = ao_two_e_integral_mixed_exchange + coef4 * integral
enddo ! s
enddo ! r
enddo ! q
enddo ! p
end

35
plugins/local/extra_basis_int/v_mixed_extra.irp.f
View File

@ -1,3 +1,6 @@
!!! TODO:: optimize when "ao_extra_only_1s" is True
double precision function v_extra_nucl_extra_ao(i_ao,j_ao)
implicit none
BEGIN_DOC
@ -6,9 +9,9 @@ double precision function v_extra_nucl_extra_ao(i_ao,j_ao)
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) v_ne^{extra}(r)$.
!
!
! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis
! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis (system B)
!
! and v_ne^{extra}(r) is the Coulomb potential coming from the EXTRA nuclei
! and v_ne^{extra}(r) is the Coulomb potential coming from the EXTRA nuclei (system B)
END_DOC
integer, intent(in) ::i_ao,j_ao
double precision :: mu_in,charge,coord(3)
@ -23,6 +26,30 @@ double precision function v_extra_nucl_extra_ao(i_ao,j_ao)
enddo
end
double precision function v_extra_nucl_ao(i_ao,j_ao)
implicit none
BEGIN_DOC
!
! Computes the following integral :
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) v_ne(r)$.
!
!
! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the REGULAR basis (system A)
!
! and v_ne(r) is the Coulomb potential coming from the EXTRA nuclei (system B)
END_DOC
integer, intent(in) ::i_ao,j_ao
integer :: i
double precision :: mu_in, coord(3),charge, integral
double precision :: NAI_pol_mult_erf_ao
mu_in = 1.d+10
do i = 1, extra_nucl_num
coord(1:3) = extra_nucl_coord_transp(1:3,i)
charge = extra_nucl_charge(i)
v_extra_nucl_ao += -NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, coord) * charge
enddo
end
double precision function v_nucl_extra_ao(i_ao,j_ao)
implicit none
@ -32,9 +59,9 @@ double precision function v_nucl_extra_ao(i_ao,j_ao)
! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) v_ne(r)$.
!
!
! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis
! where BOTH $\chi_i(r)$ AND $\chi_j(r)$ belongs to the EXTRA basis (system B)
!
! and v_ne(r) is the Coulomb potential coming from the REGULAR nuclei
! and v_ne(r) is the Coulomb potential coming from the REGULAR nuclei (system A)
END_DOC
integer, intent(in) ::i_ao,j_ao
double precision :: mu_in,charge,coord(3)

1
plugins/local/tuto_plugins/tuto_I/traces_one_e.irp.f
View File

@ -30,6 +30,7 @@ BEGIN_PROVIDER [ double precision, trace_ao_one_e_ints]
! have the same number of functions
END_DOC
integer :: i,j
double precision :: accu
double precision, allocatable :: inv_overlap_times_integrals(:,:) ! = h S^{-1}
allocate(inv_overlap_times_integrals(ao_num,ao_num))
! routine that computes the product of two matrices, you can check it with

5
src/ao_extra_basis/README.rst
View File

@ -5,7 +5,9 @@ extra_basis
Plugin to handle an extra basis, which is attached to the extra_nuclei.
It is essentially a duplication of all important quantities (coefficients, exponents and so on) of the usual |AO| basis.
An interesting feature is the possibility to fit any basis made at most with "p" functions onto a purely "s" basis.
Check in the directory "tuto" for a simple example of how to create a fictious system "B" attached independently to a system "A"
Another interesting feature is the possibility to fit any basis made at most with "p" functions onto a purely "s" basis.
This is done with the various scripts here:
- qp_fit_1s_basis : script that creates an |EZFIO| folder corresponding to an .xyz file and a basis fitted with only "s" functions
@ -13,3 +15,4 @@ This is done with the various scripts here:
Ex:
qp_add_extra_fit_system LiH.ezfio/ h2o.xyz # takes the EZFIO folder "LiH.ezfio" and creates all necessary additional basis and nuclei based on h2o.xyz, but only with 1s functions.

1
src/ao_extra_basis/fit_1s_basis.irp.f
View File

@ -31,6 +31,7 @@ program fit_1s_basis
call ezfio_set_extra_nuclei_extra_nucl_label(new_nucl_label_1s)
!
call ezfio_set_ao_extra_basis_ao_extra_num(n_func_tot)
call ezfio_set_ao_extra_basis_ao_extra_only_1s(.True.)
call ezfio_set_ao_extra_basis_ao_extra_center(ao_extra_center)
call ezfio_set_ao_extra_basis_ao_extra_nucl(new_ao_nucl_1s)
call ezfio_set_ao_extra_basis_ao_extra_prim_num(new_ao_prim_num_1s)

2
src/ao_extra_basis/qp_copy_extra_basis
View File

@ -58,7 +58,7 @@ do
done
i=primitives_normalized
newfile=primitives_normalized_extra
cp ${EZFIO_extra}/ao_basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
cp ${EZFIO_extra}/basis/$i ${EZFIO_target}/ao_extra_basis/$newfile
echo "COPYING ALL DATA FROM "$EZFIO_extra"/aux_quantities/ to "${EZFIO_target}"/ao_extra_basis/"
i=data_one_e_dm_tot_ao.gz

3
src/ao_extra_basis/tuto/He_A.xyz Normal file
View File

@ -0,0 +1,3 @@
1
He atom "A"
He 0. 0. 0.

26
src/ao_extra_basis/tuto/example_copy.sh Executable file
View File

@ -0,0 +1,26 @@
source ~/qp2/quantum_package.rc
## Example of how to generate an additional h2o molecule, stored as a extra basis/nuclei etc .. to an He
sys_B=h2o.xyz
basis_B=sto-3g
output_B=${sys_B%.xyz}_${basis_B}
sys_A=He_A.xyz
basis_A=cc-pvtz
output_A=${sys_A%.xyz}_${basis_A}_extra_${output_B}
# we create the system "B" that will be attached as an "extra system" to the syste "A"
qp create_ezfio -b $basis_B $sys_B -o ${output_B}
# we perform an HF calculation to obtain the AO density matrix
qp run scf
# we save the density matrix in the EZFIO
qp run save_one_e_dm
# we create the system "A"
qp create_ezfio -b $basis_A $sys_A -o ${output_A}
# We perform an SCF calculation
qp run scf
# we copy the system "B" information as extra nuclei/basis etc in the EZFIO of system "A"
qp_copy_extra_basis ${output_B} ${output_A}
# we execute an example of progra that prints a lot of useful integrals/information on the A-B interaction
qp run test_extra_basis | tee ${output_A}.test_extra_basis

7
src/ao_extra_basis/tuto/h2o.xyz Normal file
View File

@ -0,0 +1,7 @@
3
O 0.000000 -0.399441 3.000000
H 0.761232 0.199721 3.000000
H -0.761232 0.199721 3.000000