From 1363b436e399d666cb00effc887c4e843178d19b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 16 Oct 2015 14:33:57 +0200 Subject: [PATCH] Accelerated pseudopotentials --- ocaml/Qputils.ml | 1 - plugins/Full_CI/README.rst | 4 ++ plugins/Full_CI/target_pt2.irp.f | 2 + plugins/Hartree_Fock/README.rst | 4 ++ plugins/Hartree_Fock/huckel.irp.f | 3 - src/Integrals_Monoelec/README.rst | 5 ++ src/Integrals_Monoelec/pseudopot.f90 | 83 ++++++++++++++-------------- src/Pseudo/EZFIO.cfg | 2 +- src/Pseudo/README.rst | 2 +- 9 files changed, 60 insertions(+), 46 deletions(-) diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index ed112de3..7f840fde 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -12,7 +12,6 @@ let rec transpose = function ;; *) - let input_to_sexp s = let result = String.split_lines s diff --git a/plugins/Full_CI/README.rst b/plugins/Full_CI/README.rst index bc2307cd..a2910125 100644 --- a/plugins/Full_CI/README.rst +++ b/plugins/Full_CI/README.rst @@ -23,6 +23,10 @@ Documentation .. by the `update_README.py` script. +`e_curve `_ + Undocumented + + `full_ci `_ Undocumented diff --git a/plugins/Full_CI/target_pt2.irp.f b/plugins/Full_CI/target_pt2.irp.f index f86a21a4..7e7c8fcf 100644 --- a/plugins/Full_CI/target_pt2.irp.f +++ b/plugins/Full_CI/target_pt2.irp.f @@ -27,6 +27,8 @@ program var_pt2_ratio_run call diagonalize_CI ratio = (CI_energy(1) - HF_energy) / (CI_energy(1)+pt2(1) - HF_energy) if (N_det > 20000) then + N_det = 20000 + TOUCH N_det exit endif enddo diff --git a/plugins/Hartree_Fock/README.rst b/plugins/Hartree_Fock/README.rst index ffe80c75..5ebd9bf7 100644 --- a/plugins/Hartree_Fock/README.rst +++ b/plugins/Hartree_Fock/README.rst @@ -161,6 +161,10 @@ Documentation optional: mo_basis.mo_coef +`simple_scf `_ + Undocumented + + `thresh_scf `_ Threshold on the convergence of the Hartree Fock energy diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f index a0987bc9..6139ac46 100644 --- a/plugins/Hartree_Fock/huckel.irp.f +++ b/plugins/Hartree_Fock/huckel.irp.f @@ -27,11 +27,8 @@ subroutine huckel_guess Fock_matrix_ao(j,j) = Fock_matrix_alpha_ao(j,j) enddo TOUCH Fock_matrix_ao - print *, "Huckel matrix computed" mo_coef = eigenvectors_fock_matrix_mo SOFT_TOUCH mo_coef - print *, "Saving MOs" call save_mos - print *, "Saving MOs saved" end diff --git a/src/Integrals_Monoelec/README.rst b/src/Integrals_Monoelec/README.rst index 13aceb0e..16230cbf 100644 --- a/src/Integrals_Monoelec/README.rst +++ b/src/Integrals_Monoelec/README.rst @@ -93,6 +93,11 @@ Documentation : sum of the kinetic and nuclear electronic potential +`ao_mono_elec_integral_diag `_ + array of the mono electronic hamiltonian on the AOs basis + : sum of the kinetic and nuclear electronic potential + + `ao_nucl_elec_integral `_ interaction nuclear electron diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90 index bd00dc51..1f2d4b2f 100644 --- a/src/Integrals_Monoelec/pseudopot.f90 +++ b/src/Integrals_Monoelec/pseudopot.f90 @@ -1829,11 +1829,12 @@ end double precision function binom_gen(alpha,n) implicit none integer :: n,k - double precision :: fact,alpha,prod + double precision :: fact,alpha,prod, factn_inv prod=1.d0 + factn_inv = 1.d0/fact(n) do k=0,n-1 prod=prod*(alpha-k) - binom_gen = prod/(fact(n)) + binom_gen = prod*factn_inv enddo end @@ -1881,6 +1882,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) double precision :: term_A, term_B, term_rap, expo double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square double precision :: int_prod_bessel_loc + double precision :: inverses(0:300) logical done @@ -1927,18 +1929,32 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) ! Initialise the first recurence terme for the q loop s_q_0 = s_0_0 + ! Loop over q for the convergence of the sequence do while (.not.done) ! Init - sum=0 s_q_k=s_q_0 + sum=s_q_0 - ! Iteration of k - do k=0,q + if (q>300) then + stop 'pseudopot.f90 : q > 300' + endif + + do k=0,q-1 + s_q_k = ( dble(2*(q-k+m)+1)*dble(q-k)*inverses(k) ) * s_q_k sum=sum+s_q_k - s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)/dble(2*(k+n)+3) ) * ( dble(q-k)/dble(k+1)) * s_q_k enddo + inverses(q) = a_over_b_square/(dble(2*(q+n)+3) * dble(q+1)) +! do k=0,q +! sum=sum+s_q_k +! s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k +! enddo + ! Iteration of k +! do k=0,q +! sum=sum+s_q_k +! s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k +! enddo int=int+sum @@ -2120,15 +2136,14 @@ parameter (ntot_max=14) integer l,m double precision a(3),g_a,c(3) double precision prod,binom_func,accu,bigI,ylm,bessel_mod -double precision theta_AC0,phi_AC0,ac,factor,fourpi,arg,r,areal +double precision theta_AC0,phi_AC0,ac,ac2,factor,fourpi,arg,r,areal integer ntotA,mu,k1,k2,k3,lambda integer n_a(3) -double precision & -array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max) double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max), y +ac2=(a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2 ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) -arg=g_a*(ac**2+r**2) +arg=g_a*(ac2+r*r) fourpi=4.d0*dacos(-1.d0) factor=fourpi*dexp(-arg) areal=2.d0*g_a*ac @@ -2154,39 +2169,27 @@ else enddo enddo - do lambda=0,l+ntotA - do mu=-lambda,lambda - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) - enddo - enddo - enddo - enddo - enddo - accu=0.d0 do lambda=0,l+ntotA - do mu=-lambda,lambda - y = ylm(lambda,mu,theta_AC0,phi_AC0) - if (y == 0.d0) then - cycle - endif - do k1=0,n_a(1) - do k2=0,n_a(2) - do k3=0,n_a(3) - prod=y*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3) - if (prod == 0.d0) then - cycle - endif - if (areal*r < 100.d0) then ! overflow! - accu=accu+prod*bessel_mod(areal*r,lambda) - endif + do mu=-lambda,lambda + y = ylm(lambda,mu,theta_AC0,phi_AC0) + if (y == 0.d0) then + cycle + endif + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + prod=y*array_coefs_A(k1,k2,k3)*bigI(lambda,mu,l,m,k1,k2,k3) + if (prod == 0.d0) then + cycle + endif + if (areal*r < 100.d0) then ! overflow! + accu=accu+prod*bessel_mod(areal*r,lambda) + endif + enddo + enddo + enddo enddo - enddo - enddo - enddo enddo ylm_orb=factor*accu return diff --git a/src/Pseudo/EZFIO.cfg b/src/Pseudo/EZFIO.cfg index 80d857e3..58df78a0 100644 --- a/src/Pseudo/EZFIO.cfg +++ b/src/Pseudo/EZFIO.cfg @@ -51,7 +51,7 @@ size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax) [do_pseudo] type: logical -doc: Using pseudo potential integral of not +doc: Using pseudo potential integral or not interface: ezfio,provider,ocaml default: False diff --git a/src/Pseudo/README.rst b/src/Pseudo/README.rst index 062a9465..cba187aa 100644 --- a/src/Pseudo/README.rst +++ b/src/Pseudo/README.rst @@ -29,7 +29,7 @@ Documentation `do_pseudo `_ - Using pseudo potential integral of not + Using pseudo potential integral or not `pseudo_dz_k `_