LCPQ/eplf
10
0
Fork 0
mirror of https://gitlab.com/scemama/eplf synced 2024-12-23 04:43:52 +01:00
Code Issues Releases Wiki Activity

Using MO basis => fast!!!

Browse Source
This commit is contained in:
Anthony Scemama 2009-05-15 17:18:39 +02:00
parent b61b8e51e3
commit 79949f7bca
5 changed files with 101 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
debug.irp.f
View File

@ -34,11 +34,11 @@ program debug
print *, 'Overlap integral :', ao_overlap(i,j) print *, 'Overlap integral :', ao_overlap(i,j)
print *, 'Overlap integral N :', ao_overlap_numeric(i,j) print *, 'Overlap integral N :', ao_overlap_numeric(i,j)
double precision :: eplf_integral, eplf_integral_numeric double precision :: ao_eplf_integral, ao_eplf_integral_numeric
print *, '' print *, ''
print *, 'EPLF gamma : ', eplf_gamma print *, 'EPLF gamma : ', eplf_gamma
print *, 'EPLF integral :', eplf_integral(i,j,eplf_gamma,point) print *, 'EPLF integral :', ao_eplf_integral(i,j,eplf_gamma,point)
print *, 'EPLF integral N :', eplf_integral_numeric(i,j,eplf_gamma,point) print *, 'EPLF integral N :', ao_eplf_integral_numeric(i,j,eplf_gamma,point)
print *, '' print *, ''
print *, 'EPLF grid Npoints :', grid_eplf_x_num, grid_eplf_y_num, grid_eplf_z_num print *, 'EPLF grid Npoints :', grid_eplf_x_num, grid_eplf_y_num, grid_eplf_z_num

116
eplf.irp.f
View File

@ -6,17 +6,38 @@ BEGIN_PROVIDER [ real, eplf_gamma ]
eplf_gamma = 1000. eplf_gamma = 1000.
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, eplf_integral_matrix, (ao_num,ao_num) ] BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF ! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
END_DOC END_DOC
integer :: i, j integer :: i, j
double precision :: eplf_integral double precision :: ao_eplf_integral
do i=1,ao_num do i=1,ao_num
do j=i,ao_num do j=i,ao_num
eplf_integral_matrix(j,i) = eplf_integral(j,i,eplf_gamma,point) ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,eplf_gamma,point)
eplf_integral_matrix(i,j) = eplf_integral_matrix(j,i) ao_eplf_integral_matrix(i,j) = ao_eplf_integral_matrix(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_eplf_integral_matrix, (mo_occ_num,mo_occ_num) ]
implicit none
BEGIN_DOC
! Array of all the <chi_i chi_j | exp(-gamma r^2)> for EPLF
END_DOC
integer :: i, j, k, l
double precision :: ao_eplf_integral
do i=1,mo_occ_num
do j=i,mo_occ_num
mo_eplf_integral_matrix(j,i) = 0.
do k=1,ao_num
do l=1,ao_num
mo_eplf_integral_matrix(j,i) = mo_eplf_integral_matrix(j,i) + &
mo_coef(k,i)*mo_coef(l,j)*ao_eplf_integral_matrix(k,l)
enddo
enddo
mo_eplf_integral_matrix(i,j) = mo_eplf_integral_matrix(j,i)
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
@ -37,62 +58,29 @@ END_PROVIDER
PROVIDE mo_coef_transp PROVIDE mo_coef_transp
do m=1,ao_num
do n=1,ao_num
double precision :: ao_mn
ao_mn = eplf_integral_matrix(m,n)
if (abs(ao_mn) > thr) then
do k=1,ao_num
do l=1,ao_num
double precision :: ao_klmn
ao_klmn = ao_mn*ao_value_p(k)*ao_value_p(l)
if (abs(ao_klmn) > thr) then
do j=1,elec_alpha_num do j=1,elec_alpha_num
if (mo_coef_transp(j,n) /= 0.d0) then
do i=1,elec_alpha_num do i=1,elec_alpha_num
eplf_up_up = eplf_up_up + ao_klmn* & eplf_up_up = eplf_up_up + mo_value_p(i)* ( &
mo_coef_transp(i,k)*mo_coef_transp(j,n) * & mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
(mo_coef_transp(i,l)*mo_coef_transp(j,m) - & mo_value_p(j)*mo_eplf_integral_matrix(i,j) )
mo_coef_transp(i,m)*mo_coef_transp(j,l) )
enddo enddo
endif
enddo enddo
do j=1,elec_beta_num do j=1,elec_beta_num
if (mo_coef_transp(j,n) /= 0.d0) then
do i=1,elec_beta_num do i=1,elec_beta_num
eplf_up_up = eplf_up_up + ao_klmn* & eplf_up_up = eplf_up_up + mo_value_p(i)* ( &
mo_coef_transp(i,k)*mo_coef_transp(j,n) * & mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
(mo_coef_transp(i,l)*mo_coef_transp(j,m) - & mo_value_p(j)*mo_eplf_integral_matrix(i,j) )
mo_coef_transp(i,m)*mo_coef_transp(j,l) )
enddo enddo
endif
enddo enddo
do j=1,elec_beta_num do j=1,elec_beta_num
if ( (mo_coef_transp(j,n) /= 0.d0).and. &
(mo_coef_transp(j,m) /= 0.d0) ) then
do i=1,elec_alpha_num do i=1,elec_alpha_num
eplf_up_dn = eplf_up_dn + 2.d0*ao_klmn* & eplf_up_dn = eplf_up_dn + mo_value_p(i)**2 * &
mo_coef_transp(i,k)*mo_coef_transp(j,n) * & mo_eplf_integral_matrix(j,j)
mo_coef_transp(i,l)*mo_coef_transp(j,m)
enddo
endif
enddo
endif
enddo
enddo
endif
enddo enddo
enddo enddo
eplf_up_dn = 2.d0*eplf_up_dn
END_PROVIDER END_PROVIDER
@ -120,7 +108,7 @@ BEGIN_PROVIDER [ real, eplf_value ]
END_PROVIDER END_PROVIDER
double precision function eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,gmma,xr) double precision function ao_eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,gmma,xr)
implicit none implicit none
include 'constants.F' include 'constants.F'
@ -150,22 +138,22 @@ double precision function eplf_integral_primitive_oneD_numeric(a,xa,i,b,xb,j,gmm
(x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2)) (x-xa)**i * (x-xb)**j * exp(-(a*(x-xa)**2+b*(x-xb)**2+gmma*(x-xr)**2))
x = x+dx x = x+dx
enddo enddo
eplf_integral_primitive_oneD_numeric = dtemp*dx ao_eplf_integral_primitive_oneD_numeric = dtemp*dx
end function end function
double precision function eplf_integral_numeric(i,j,gmma,center) double precision function ao_eplf_integral_numeric(i,j,gmma,center)
implicit none implicit none
integer, intent(in) :: i, j integer, intent(in) :: i, j
integer :: p,q,k integer :: p,q,k
double precision :: integral(ao_prim_num_max,ao_prim_num_max) double precision :: integral(ao_prim_num_max,ao_prim_num_max)
double precision :: eplf_integral_primitive_oneD_numeric double precision :: ao_eplf_integral_primitive_oneD_numeric
real :: gmma, center(3) real :: gmma, center(3)
do q=1,ao_prim_num(j) do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i) do p=1,ao_prim_num(i)
integral(p,q) = & integral(p,q) = &
eplf_integral_primitive_oneD_numeric( & ao_eplf_integral_primitive_oneD_numeric( &
ao_expo(p,i), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),1), & nucl_coord(ao_nucl(i),1), &
ao_power(i,1), & ao_power(i,1), &
@ -174,7 +162,7 @@ double precision function eplf_integral_numeric(i,j,gmma,center)
ao_power(j,1), & ao_power(j,1), &
gmma, & gmma, &
center(1)) * & center(1)) * &
eplf_integral_primitive_oneD_numeric( & ao_eplf_integral_primitive_oneD_numeric( &
ao_expo(p,i), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),2), & nucl_coord(ao_nucl(i),2), &
ao_power(i,2), & ao_power(i,2), &
@ -183,7 +171,7 @@ double precision function eplf_integral_numeric(i,j,gmma,center)
ao_power(j,2), & ao_power(j,2), &
gmma, & gmma, &
center(2)) * & center(2)) * &
eplf_integral_primitive_oneD_numeric( & ao_eplf_integral_primitive_oneD_numeric( &
ao_expo(p,i), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),3), & nucl_coord(ao_nucl(i),3), &
ao_power(i,3), & ao_power(i,3), &
@ -201,16 +189,16 @@ double precision function eplf_integral_numeric(i,j,gmma,center)
enddo enddo
enddo enddo
eplf_integral_numeric = 0. ao_eplf_integral_numeric = 0.
do q=1,ao_prim_num(j) do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i) do p=1,ao_prim_num(i)
eplf_integral_numeric = eplf_integral_numeric + integral(p,q) ao_eplf_integral_numeric = ao_eplf_integral_numeric + integral(p,q)
enddo enddo
enddo enddo
end function end function
double precision function eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) double precision function ao_eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
implicit none implicit none
include 'constants.F' include 'constants.F'
@ -266,16 +254,16 @@ double precision function eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr)
enddo enddo
enddo enddo
eplf_integral_primitive_oneD = S(i,j) ao_eplf_integral_primitive_oneD = S(i,j)
end function end function
double precision function eplf_integral(i,j,gmma,center) double precision function ao_eplf_integral(i,j,gmma,center)
implicit none implicit none
integer, intent(in) :: i, j integer, intent(in) :: i, j
integer :: p,q,k integer :: p,q,k
double precision :: integral(ao_prim_num_max,ao_prim_num_max) double precision :: integral(ao_prim_num_max,ao_prim_num_max)
double precision :: eplf_integral_primitive_oneD double precision :: ao_eplf_integral_primitive_oneD
real :: gmma, center(3) real :: gmma, center(3)
ASSERT(i>0) ASSERT(i>0)
@ -286,7 +274,7 @@ double precision function eplf_integral(i,j,gmma,center)
do q=1,ao_prim_num(j) do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i) do p=1,ao_prim_num(i)
integral(p,q) = & integral(p,q) = &
eplf_integral_primitive_oneD( & ao_eplf_integral_primitive_oneD( &
ao_expo(p,i), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),1), & nucl_coord(ao_nucl(i),1), &
ao_power(i,1), & ao_power(i,1), &
@ -295,7 +283,7 @@ double precision function eplf_integral(i,j,gmma,center)
ao_power(j,1), & ao_power(j,1), &
gmma, & gmma, &
center(1)) * & center(1)) * &
eplf_integral_primitive_oneD( & ao_eplf_integral_primitive_oneD( &
ao_expo(p,i), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),2), & nucl_coord(ao_nucl(i),2), &
ao_power(i,2), & ao_power(i,2), &
@ -304,7 +292,7 @@ double precision function eplf_integral(i,j,gmma,center)
ao_power(j,2), & ao_power(j,2), &
gmma, & gmma, &
center(2)) * & center(2)) * &
eplf_integral_primitive_oneD( & ao_eplf_integral_primitive_oneD( &
ao_expo(p,i), & ao_expo(p,i), &
nucl_coord(ao_nucl(i),3), & nucl_coord(ao_nucl(i),3), &
ao_power(i,3), & ao_power(i,3), &
@ -322,10 +310,10 @@ double precision function eplf_integral(i,j,gmma,center)
enddo enddo
enddo enddo
eplf_integral = 0. ao_eplf_integral = 0.
do q=1,ao_prim_num(j) do q=1,ao_prim_num(j)
do p=1,ao_prim_num(i) do p=1,ao_prim_num(i)
eplf_integral = eplf_integral + integral(p,q) ao_eplf_integral = ao_eplf_integral + integral(p,q)
enddo enddo
enddo enddo

8
eplf_hf.irp.f
View File

@ -1,5 +1,11 @@
program debug program eplf_hf
PROVIDE ao_prim_num_max PROVIDE ao_prim_num_max
call write_grid_eplf() call write_grid_eplf()
IRP_IF MPI
integer :: ierr
call MPI_FINALIZE(ierr)
IRP_ENDIF
end end

8
mo.irp.f
View File

@ -40,6 +40,14 @@ BEGIN_PROVIDER [ integer, mo_num ]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ integer, mo_occ_num ]
implicit none
BEGIN_DOC
! Number of occupied molecular orbitals
END_DOC
mo_occ_num = mo_closed_num + mo_active_num
END_PROVIDER
BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ] BEGIN_PROVIDER [ real, mo_coef, (ao_num,mo_num) ]
implicit none implicit none

19
mo_point.irp.f Normal file
View File

@ -0,0 +1,19 @@
BEGIN_PROVIDER [ real, mo_value_p, (mo_num) ]
implicit none
BEGIN_DOC
! Values of the molecular orbitals
END_DOC
integer :: i, j, k
do j=1,mo_num
mo_value_p(j) = 0.
do k=1,ao_num
mo_value_p(j) = mo_value_p(j)+mo_coef(k,j)*ao_value_p(k)
enddo
enddo
END_PROVIDER