10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
This commit is contained in:
Anthony Scemama 2017-09-14 11:36:27 +02:00
parent 091e9a0208
commit 1e40552708
8 changed files with 61 additions and 9 deletions

View File

@ -1,7 +1,7 @@
BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ] BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Alpha density matrix in the AO basis ! S^{-1}.P_alpha.S^{-1}
END_DOC END_DOC
call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
@ -14,7 +14,7 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num,ao_num) ] BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num,ao_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Beta density matrix in the AO basis ! S^{-1}.P_beta.S^{-1}
END_DOC END_DOC
call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
@ -27,7 +27,7 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ] BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Density matrix in the AO basis ! S^{-1}.P.S^{-1} where P = C.C^t
END_DOC END_DOC
ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1)) ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
if (elec_alpha_num== elec_beta_num) then if (elec_alpha_num== elec_beta_num) then

View File

@ -7,11 +7,11 @@ subroutine huckel_guess
double precision :: accu double precision :: accu
double precision :: c double precision :: c
character*(64) :: label character*(64) :: label
double precision, allocatable :: A(:,:), T(:,:) double precision, allocatable :: A(:,:)
label = "Guess" label = "Guess"
c = 0.25d0 * 1.75d0 c = 0.5d0 * 1.75d0
allocate (A(ao_num, ao_num),T(ao_num, ao_num)) allocate (A(ao_num, ao_num))
A = 0.d0 A = 0.d0
do j=1,ao_num do j=1,ao_num
do i=1,ao_num do i=1,ao_num
@ -29,7 +29,6 @@ subroutine huckel_guess
mo_coef = eigenvectors_fock_matrix_mo mo_coef = eigenvectors_fock_matrix_mo
SOFT_TOUCH mo_coef SOFT_TOUCH mo_coef
call save_mos call save_mos
print *, 'E=', HF_energy deallocate(A)
deallocate(A,T)
end end

View File

@ -34,6 +34,8 @@
power_A(1) = ao_power( j, 1 ) power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 ) power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 ) power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num do i= 1,ao_num
ao_overlap(i,j)= 0.d0 ao_overlap(i,j)= 0.d0
ao_overlap_x(i,j)= 0.d0 ao_overlap_x(i,j)= 0.d0
@ -47,6 +49,7 @@
power_B(3) = ao_power( i, 3 ) power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j) do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j) alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i) do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i) beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
@ -100,6 +103,8 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
power_A(1) = ao_power( j, 1 ) power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 ) power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 ) power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num do i= 1,ao_num
ao_overlap_abs(i,j)= 0.d0 ao_overlap_abs(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 ) B_center(1) = nucl_coord( ao_nucl(i), 1 )
@ -110,6 +115,7 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
power_B(3) = ao_power( i, 3 ) power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j) do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j) alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i) do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i) beta = ao_expo_ordered_transp(l,i)
call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1) call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)

View File

@ -50,7 +50,7 @@ default: 0.99
type: Threshold type: Threshold
doc: Thresholds on selectors (fraction of the norm) doc: Thresholds on selectors (fraction of the norm)
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 0.995 default: 0.999
[n_int] [n_int]
interface: ezfio interface: ezfio

View File

@ -7,6 +7,7 @@
! : sum of the kinetic and nuclear electronic potential ! : sum of the kinetic and nuclear electronic potential
END_DOC END_DOC
do j = 1, ao_num do j = 1, ao_num
!DIR$ VECTOR ALIGNED
do i = 1, ao_num do i = 1, ao_num
ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j) + ao_pseudo_integral(i,j) ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j) + ao_pseudo_integral(i,j)
enddo enddo

View File

@ -45,6 +45,8 @@
power_A(1) = ao_power( j, 1 ) power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 ) power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 ) power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num do i= 1,ao_num
ao_deriv2_x(i,j)= 0.d0 ao_deriv2_x(i,j)= 0.d0
ao_deriv2_y(i,j)= 0.d0 ao_deriv2_y(i,j)= 0.d0
@ -57,6 +59,7 @@
power_B(3) = ao_power( i, 3 ) power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j) do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j) alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i) do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i) beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
@ -136,6 +139,8 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
!$OMP PRIVATE(i,j) & !$OMP PRIVATE(i,j) &
!$OMP SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z) !$OMP SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
do j = 1, ao_num do j = 1, ao_num
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
do i = 1, ao_num do i = 1, ao_num
ao_kinetic_integral(i,j) = -0.5d0 * (ao_deriv2_x(i,j) + ao_deriv2_y(i,j) + ao_deriv2_z(i,j) ) ao_kinetic_integral(i,j) = -0.5d0 * (ao_deriv2_x(i,j) + ao_deriv2_y(i,j) + ao_deriv2_z(i,j) )
enddo enddo

View File

@ -35,6 +35,8 @@
power_A(1) = ao_power( j, 1 ) power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 ) power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 ) power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 ) B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 ) B_center(2) = nucl_coord( ao_nucl(i), 2 )
@ -47,6 +49,7 @@
accu_z = 0.d0 accu_z = 0.d0
do n = 1,ao_prim_num(j) do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j) alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i) do l = 1, ao_prim_num(i)
c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,i) c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,i)
beta = ao_expo_ordered_transp(l,i) beta = ao_expo_ordered_transp(l,i)
@ -106,6 +109,8 @@
power_A(1) = ao_power( j, 1 ) power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 ) power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 ) power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 ) B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 ) B_center(2) = nucl_coord( ao_nucl(i), 2 )
@ -118,6 +123,7 @@
accu_z = 0.d0 accu_z = 0.d0
do n = 1,ao_prim_num(j) do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j) alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i) do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i) beta = ao_expo_ordered_transp(l,i)
c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j) c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j)
@ -177,6 +183,8 @@
power_A(1) = ao_power( j, 1 ) power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 ) power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 ) power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 ) B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 ) B_center(2) = nucl_coord( ao_nucl(i), 2 )
@ -189,6 +197,7 @@
accu_z = 0.d0 accu_z = 0.d0
do n = 1,ao_prim_num(j) do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j) alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i) do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i) beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)

View File

@ -49,6 +49,9 @@ END_PROVIDER
do j=1,ao_num do j=1,ao_num
mo_coef(j,i) = buffer(j,i) mo_coef(j,i) = buffer(j,i)
enddo enddo
do j=ao_num+1,ao_num_align
mo_coef(j,i) = 0.d0
enddo
enddo enddo
deallocate(buffer) deallocate(buffer)
else else
@ -57,6 +60,9 @@ END_PROVIDER
do j=1,ao_num do j=1,ao_num
mo_coef(j,i) = ao_ortho_canonical_coef(j,i) mo_coef(j,i) = ao_ortho_canonical_coef(j,i)
enddo enddo
do j=ao_num+1,ao_num_align
mo_coef(j,i) = 0.d0
enddo
enddo enddo
endif endif
END_PROVIDER END_PROVIDER
@ -202,6 +208,32 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
deallocate(T) deallocate(T)
end end
subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the MO basis to the S^-1 AO basis
! Useful for density matrix
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo,mo_tot_num)
double precision, intent(out) :: A_ao(LDA_ao,ao_num)
double precision, allocatable :: T(:,:)
allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,size(A_mo,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, size(T,1))
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, size(T,1), &
0.d0, A_ao, size(A_ao,1))
deallocate(T)
end
subroutine mix_mo_jk(j,k) subroutine mix_mo_jk(j,k)
implicit none implicit none