From 79949f7bca242fc9f4853134300775b2f55f3f88 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 15 May 2009 17:18:39 +0200 Subject: [PATCH] Using MO basis => fast!!! --- debug.irp.f | 6 +-- eplf.irp.f | 140 ++++++++++++++++++++++--------------------------- eplf_hf.irp.f | 8 ++- mo.irp.f | 8 +++ mo_point.irp.f | 19 +++++++ 5 files changed, 101 insertions(+), 80 deletions(-) create mode 100644 mo_point.irp.f diff --git a/debug.irp.f b/debug.irp.f index ee3edb2..a030024 100644 --- a/debug.irp.f +++ b/debug.irp.f @@ -34,11 +34,11 @@ program debug print *, 'Overlap integral :', ao_overlap(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 *, 'EPLF gamma : ', eplf_gamma - print *, 'EPLF integral :', eplf_integral(i,j,eplf_gamma,point) - print *, 'EPLF integral N :', eplf_integral_numeric(i,j,eplf_gamma,point) + print *, 'EPLF integral :', ao_eplf_integral(i,j,eplf_gamma,point) + print *, 'EPLF integral N :', ao_eplf_integral_numeric(i,j,eplf_gamma,point) print *, '' print *, 'EPLF grid Npoints :', grid_eplf_x_num, grid_eplf_y_num, grid_eplf_z_num diff --git a/eplf.irp.f b/eplf.irp.f index ae0d8fe..956d24b 100644 --- a/eplf.irp.f +++ b/eplf.irp.f @@ -6,17 +6,38 @@ BEGIN_PROVIDER [ real, eplf_gamma ] eplf_gamma = 1000. 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 BEGIN_DOC ! Array of all the for EPLF END_DOC integer :: i, j - double precision :: eplf_integral + double precision :: ao_eplf_integral do i=1,ao_num do j=i,ao_num - eplf_integral_matrix(j,i) = eplf_integral(j,i,eplf_gamma,point) - eplf_integral_matrix(i,j) = eplf_integral_matrix(j,i) + ao_eplf_integral_matrix(j,i) = ao_eplf_integral(j,i,eplf_gamma,point) + 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 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 END_PROVIDER @@ -37,63 +58,30 @@ END_PROVIDER 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 - if (mo_coef_transp(j,n) /= 0.d0) then - do i=1,elec_alpha_num - eplf_up_up = eplf_up_up + ao_klmn* & - mo_coef_transp(i,k)*mo_coef_transp(j,n) * & - (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & - mo_coef_transp(i,m)*mo_coef_transp(j,l) ) - enddo - endif - enddo - - do j=1,elec_beta_num - if (mo_coef_transp(j,n) /= 0.d0) then - do i=1,elec_beta_num - eplf_up_up = eplf_up_up + ao_klmn* & - mo_coef_transp(i,k)*mo_coef_transp(j,n) * & - (mo_coef_transp(i,l)*mo_coef_transp(j,m) - & - mo_coef_transp(i,m)*mo_coef_transp(j,l) ) - enddo - endif - enddo - - 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 - eplf_up_dn = eplf_up_dn + 2.d0*ao_klmn* & - mo_coef_transp(i,k)*mo_coef_transp(j,n) * & - mo_coef_transp(i,l)*mo_coef_transp(j,m) - enddo - endif - enddo - - endif - - enddo - enddo - - endif - + do j=1,elec_alpha_num + do i=1,elec_alpha_num + eplf_up_up = eplf_up_up + mo_value_p(i)* ( & + mo_value_p(i)*mo_eplf_integral_matrix(j,j) - & + mo_value_p(j)*mo_eplf_integral_matrix(i,j) ) enddo enddo + do j=1,elec_beta_num + do i=1,elec_beta_num + eplf_up_up = eplf_up_up + mo_value_p(i)* ( & + mo_value_p(i)*mo_eplf_integral_matrix(j,j) - & + mo_value_p(j)*mo_eplf_integral_matrix(i,j) ) + enddo + enddo + + do j=1,elec_beta_num + do i=1,elec_alpha_num + eplf_up_dn = eplf_up_dn + mo_value_p(i)**2 * & + mo_eplf_integral_matrix(j,j) + enddo + enddo + eplf_up_dn = 2.d0*eplf_up_dn + END_PROVIDER @@ -120,7 +108,7 @@ BEGIN_PROVIDER [ real, eplf_value ] 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 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 = x+dx enddo - eplf_integral_primitive_oneD_numeric = dtemp*dx + ao_eplf_integral_primitive_oneD_numeric = dtemp*dx 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 integer, intent(in) :: i, j integer :: p,q,k 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) do q=1,ao_prim_num(j) do p=1,ao_prim_num(i) integral(p,q) = & - eplf_integral_primitive_oneD_numeric( & + ao_eplf_integral_primitive_oneD_numeric( & ao_expo(p,i), & nucl_coord(ao_nucl(i),1), & ao_power(i,1), & @@ -174,7 +162,7 @@ double precision function eplf_integral_numeric(i,j,gmma,center) ao_power(j,1), & gmma, & center(1)) * & - eplf_integral_primitive_oneD_numeric( & + ao_eplf_integral_primitive_oneD_numeric( & ao_expo(p,i), & nucl_coord(ao_nucl(i),2), & ao_power(i,2), & @@ -183,7 +171,7 @@ double precision function eplf_integral_numeric(i,j,gmma,center) ao_power(j,2), & gmma, & center(2)) * & - eplf_integral_primitive_oneD_numeric( & + ao_eplf_integral_primitive_oneD_numeric( & ao_expo(p,i), & nucl_coord(ao_nucl(i),3), & ao_power(i,3), & @@ -201,16 +189,16 @@ double precision function eplf_integral_numeric(i,j,gmma,center) enddo enddo - eplf_integral_numeric = 0. + ao_eplf_integral_numeric = 0. do q=1,ao_prim_num(j) 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 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 include 'constants.F' @@ -266,16 +254,16 @@ double precision function eplf_integral_primitive_oneD(a,xa,i,b,xb,j,gmma,xr) enddo enddo - eplf_integral_primitive_oneD = S(i,j) + ao_eplf_integral_primitive_oneD = S(i,j) end function -double precision function eplf_integral(i,j,gmma,center) +double precision function ao_eplf_integral(i,j,gmma,center) implicit none integer, intent(in) :: i, j integer :: p,q,k 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) ASSERT(i>0) @@ -286,7 +274,7 @@ double precision function eplf_integral(i,j,gmma,center) do q=1,ao_prim_num(j) do p=1,ao_prim_num(i) integral(p,q) = & - eplf_integral_primitive_oneD( & + ao_eplf_integral_primitive_oneD( & ao_expo(p,i), & nucl_coord(ao_nucl(i),1), & ao_power(i,1), & @@ -295,7 +283,7 @@ double precision function eplf_integral(i,j,gmma,center) ao_power(j,1), & gmma, & center(1)) * & - eplf_integral_primitive_oneD( & + ao_eplf_integral_primitive_oneD( & ao_expo(p,i), & nucl_coord(ao_nucl(i),2), & ao_power(i,2), & @@ -304,7 +292,7 @@ double precision function eplf_integral(i,j,gmma,center) ao_power(j,2), & gmma, & center(2)) * & - eplf_integral_primitive_oneD( & + ao_eplf_integral_primitive_oneD( & ao_expo(p,i), & nucl_coord(ao_nucl(i),3), & ao_power(i,3), & @@ -322,10 +310,10 @@ double precision function eplf_integral(i,j,gmma,center) enddo enddo - eplf_integral = 0. + ao_eplf_integral = 0. do q=1,ao_prim_num(j) 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 diff --git a/eplf_hf.irp.f b/eplf_hf.irp.f index b67765f..eb7b012 100644 --- a/eplf_hf.irp.f +++ b/eplf_hf.irp.f @@ -1,5 +1,11 @@ -program debug +program eplf_hf + PROVIDE ao_prim_num_max call write_grid_eplf() + + IRP_IF MPI + integer :: ierr + call MPI_FINALIZE(ierr) + IRP_ENDIF end diff --git a/mo.irp.f b/mo.irp.f index 72dbc3a..8632030 100644 --- a/mo.irp.f +++ b/mo.irp.f @@ -40,6 +40,14 @@ BEGIN_PROVIDER [ integer, mo_num ] 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) ] implicit none diff --git a/mo_point.irp.f b/mo_point.irp.f new file mode 100644 index 0000000..f9a3124 --- /dev/null +++ b/mo_point.irp.f @@ -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 + +