From e9dccd2364f282397df9f3b5bc4e3373fe3bd7e6 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 25 Apr 2024 19:46:26 +0200 Subject: [PATCH 1/8] added spherical harmonics --- plugins/local/spher_harm/.gitignore | 59 +++++ plugins/local/spher_harm/NEED | 1 + plugins/local/spher_harm/README.rst | 4 + plugins/local/spher_harm/assoc_gaus_pol.irp.f | 50 ++++ plugins/local/spher_harm/spher_harm.irp.f | 217 ++++++++++++++++++ .../local/spher_harm/spher_harm_func.irp.f | 151 ++++++++++++ 6 files changed, 482 insertions(+) create mode 100644 plugins/local/spher_harm/.gitignore create mode 100644 plugins/local/spher_harm/NEED create mode 100644 plugins/local/spher_harm/README.rst create mode 100644 plugins/local/spher_harm/assoc_gaus_pol.irp.f create mode 100644 plugins/local/spher_harm/spher_harm.irp.f create mode 100644 plugins/local/spher_harm/spher_harm_func.irp.f diff --git a/plugins/local/spher_harm/.gitignore b/plugins/local/spher_harm/.gitignore new file mode 100644 index 00000000..1561915b --- /dev/null +++ b/plugins/local/spher_harm/.gitignore @@ -0,0 +1,59 @@ +IRPF90_temp/ +IRPF90_man/ +build.ninja +irpf90.make +ezfio_interface.irp.f +irpf90_entities +tags +Makefile +ao_basis +ao_one_e_ints +ao_two_e_erf_ints +ao_two_e_ints +aux_quantities +becke_numerical_grid +bitmask +cis +cisd +cipsi +davidson +davidson_dressed +davidson_undressed +density_for_dft +determinants +dft_keywords +dft_utils_in_r +dft_utils_one_e +dft_utils_two_body +dressing +dummy +electrons +ezfio_files +fci +generators_cas +generators_full +hartree_fock +iterations +kohn_sham +kohn_sham_rs +mo_basis +mo_guess +mo_one_e_ints +mo_two_e_erf_ints +mo_two_e_ints +mpi +mrpt_utils +nuclei +perturbation +pseudo +psiref_cas +psiref_utils +scf_utils +selectors_cassd +selectors_full +selectors_utils +single_ref_method +slave +tools +utils +zmq diff --git a/plugins/local/spher_harm/NEED b/plugins/local/spher_harm/NEED new file mode 100644 index 00000000..92df7f12 --- /dev/null +++ b/plugins/local/spher_harm/NEED @@ -0,0 +1 @@ +dft_utils_in_r diff --git a/plugins/local/spher_harm/README.rst b/plugins/local/spher_harm/README.rst new file mode 100644 index 00000000..bf897f73 --- /dev/null +++ b/plugins/local/spher_harm/README.rst @@ -0,0 +1,4 @@ +========== +spher_harm +========== + diff --git a/plugins/local/spher_harm/assoc_gaus_pol.irp.f b/plugins/local/spher_harm/assoc_gaus_pol.irp.f new file mode 100644 index 00000000..fa790307 --- /dev/null +++ b/plugins/local/spher_harm/assoc_gaus_pol.irp.f @@ -0,0 +1,50 @@ +double precision function plgndr(l,m,x) + integer, intent(in) :: l,m + double precision, intent(in) :: x + BEGIN_DOC + ! associated Legenre polynom P_l,m(x). Used for the Y_lm(theta,phi) + ! Taken from https://iate.oac.uncor.edu/~mario/materia/nr/numrec/f6-8.pdf + END_DOC + integer :: i,ll + double precision :: fact,pll,pmm,pmmp1,somx2 + if(m.lt.0.or.m.gt.l.or.dabs(x).gt.1.d0)then + print*,'bad arguments in plgndr' + pause + endif + pmm=1.d0 + if(m.gt.0) then + somx2=dsqrt((1.d0-x)*(1.d0+x)) + fact=1.d0 + do i=1,m + pmm=-pmm*fact*somx2 + fact=fact+2.d0 + enddo + endif ! m > 0 + if(l.eq.m) then + plgndr=pmm + else + pmmp1=x*(2*m+1)*pmm ! Compute P_m+1^m + if(l.eq.m+1) then + plgndr=pmmp1 + else ! Compute P_l^m, l> m+1 + do ll=m+2,l + pll=(x*dble(2*ll-1)*pmmp1-dble(ll+m-1)*pmm)/(ll-m) + pmm=pmmp1 + pmmp1=pll + enddo + plgndr=pll + endif ! l.eq.m+1 + endif ! l.eq.m + return +end + +double precision function ortho_assoc_gaus_pol(l1,m1,l2) + implicit none + integer, intent(in) :: l1,m1,l2 + double precision :: fact + if(l1.ne.l2)then + ortho_assoc_gaus_pol= 0.d0 + else + ortho_assoc_gaus_pol = 2.d0*fact(l1+m1) / (dble(2*l1+1)*fact(l1-m1)) + endif +end diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f new file mode 100644 index 00000000..40661db1 --- /dev/null +++ b/plugins/local/spher_harm/spher_harm.irp.f @@ -0,0 +1,217 @@ +program spher_harm + implicit none + call test_spher_harm +! call test_cart +! call test_brutal_spheric +end + +subroutine test_cart + implicit none + include 'constants.include.F' + double precision :: r(3),theta,phi,r_abs + print*,'' + r = 0.d0 + r(1) = 1.d0 + r(2) = 1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) =-1.d0 + r(2) = 1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) =-1.d0 + r(2) =-1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) = 1.d0 + r(2) =-1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi +end + +subroutine test_spher_harm + implicit none + include 'constants.include.F' + integer :: l1,m1,i,l2,m2,lmax + double precision :: r(3),weight,accu_re, accu_im,accu + double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 + l1 = 0 + m1 = 0 + l2 = 0 + m2 = 0 + lmax = 5 + do l1 = 0,lmax + do m1 = -l1 ,l1 + do l2 = 0,lmax + do m2 = -l2 ,l2 + accu_re = 0.d0 + accu_im = 0.d0 + ! = \int dOmega Y_l1,m1^* Y_l2,m2 + ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) + ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu = 0.d0 + do i = 1, n_points_integration_angular + double precision :: theta,phi,r_abs + r(1:3) = angular_quadrature_points(i,1:3) + weight = weights_angular_points(i) + call cartesian_to_spherical(r,theta,phi,r_abs) + if(theta.gt.pi.or.theta.lt.0.d0)then + print*,'pb with theta',theta + print*,r + endif + if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then + print*,'pb with phi',phi/pi + print*,r + endif + call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) + call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) +! call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) +! call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) +! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) +! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) +! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) +! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + write(33,'(10(F16.10,X))')phi/pi + enddo + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + enddo + enddo + enddo + enddo + double precision :: x,dx,xmax,xmin + integer:: nx + nx = 10000 + xmin = -5.d0 + xmax = 5.d0 + dx = (xmax - xmin)/dble(nx) + x = xmin + do i = 1, nx + write(34,*)x,datan(x),dacos(x) + x += dx + enddo +end + +subroutine test_brutal_spheric + implicit none + include 'constants.include.F' + integer :: itheta, iphi,ntheta,nphi + double precision :: theta_min, theta_max, dtheta,theta + double precision :: phi_min, phi_max, dphi,phi + double precision :: accu_re, accu_im,weight + double precision :: re_ylm_1, im_ylm_1 ,re_ylm_2, im_ylm_2,accu + integer :: l1,m1,i,l2,m2,lmax + phi_min = 0.d0 + phi_max = 2.D0 * pi + theta_min = 0.d0 + theta_max = 1.D0 * pi + ntheta = 1000 + nphi = 1000 + dphi = (phi_max - phi_min)/dble(nphi) + dtheta = (theta_max - theta_min)/dble(ntheta) + + lmax = 3 + do l1 = 0,lmax + do m1 = 0 ,l1 + do l2 = 0,lmax + do m2 = 0 ,l2 + accu_re = 0.d0 + accu_im = 0.d0 + accu = 0.d0 + theta = theta_min + do itheta = 1, ntheta + phi = phi_min + do iphi = 1, nphi +! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) +! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) + call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) + call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) + weight = dtheta * dphi * dsin(theta) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + phi += dphi + enddo + theta += dtheta + enddo + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + print*,'accu = ',accu + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + endif + endif + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + endif + endif + enddo + enddo + enddo + enddo + + +end + +subroutine test_assoc_leg_pol + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + integer :: l1,m1,ngrid,i,l2,m2 + l1 = 0 + m1 = 0 + l2 = 2 + m2 = 0 + double precision :: x, dx,xmax,accu,xmin + double precision :: plgndr,func_1,func_2,ortho_assoc_gaus_pol + ngrid = 100000 + xmax = 1.d0 + xmin = -1.d0 + dx = (xmax-xmin)/dble(ngrid) + do l2 = 0,10 + x = xmin + accu = 0.d0 + do i = 1, ngrid + func_1 = plgndr(l1,m1,x) + func_2 = plgndr(l2,m2,x) + write(33,*)x, func_1,func_2 + accu += func_1 * func_2 * dx + x += dx + enddo + print*,'l2 = ',l2 + print*,'accu = ',accu + print*,ortho_assoc_gaus_pol(l1,m1,l2) + enddo +end diff --git a/plugins/local/spher_harm/spher_harm_func.irp.f b/plugins/local/spher_harm/spher_harm_func.irp.f new file mode 100644 index 00000000..825bd8ac --- /dev/null +++ b/plugins/local/spher_harm/spher_harm_func.irp.f @@ -0,0 +1,151 @@ +subroutine spher_harm_func_r3(r,l,m,re_ylm, im_ylm) + implicit none + integer, intent(in) :: l,m + double precision, intent(in) :: r(3) + double precision, intent(out) :: re_ylm, im_ylm + + double precision :: theta, phi,r_abs + call cartesian_to_spherical(r,theta,phi,r_abs) + call spher_harm_func(l,m,theta,phi,re_ylm, im_ylm) +end + + +subroutine spher_harm_func_m_pos(l,m,theta,phi,re_ylm, im_ylm) + include 'constants.include.F' + implicit none + BEGIN_DOC +! Y_lm(theta,phi) with m >0 +! + END_DOC + double precision, intent(in) :: theta, phi + integer, intent(in) :: l,m + double precision, intent(out):: re_ylm,im_ylm + double precision :: prefact,fact,cos_theta,plgndr,p_lm + double precision :: tmp + prefact = dble(2*l+1)*fact(l-m)/(dfour_pi * fact(l+m)) + prefact = dsqrt(prefact) + cos_theta = dcos(theta) + p_lm = plgndr(l,m,cos_theta) + tmp = prefact * p_lm + re_ylm = dcos(dble(m)*phi) * tmp + im_ylm = dsin(dble(m)*phi) * tmp +end + +subroutine spher_harm_func(l,m,theta,phi,re_ylm, im_ylm) + implicit none + BEGIN_DOC + ! Y_lm(theta,phi) with -l l in spher_harm_func !! stopping ...' + stop + endif + if(m.ge.0)then + call spher_harm_func_m_pos(l,m,theta,phi,re_ylm_pos, im_ylm_pos) + re_ylm = re_ylm_pos + im_ylm = im_ylm_pos + else + minus_m = -m !> 0 + call spher_harm_func_m_pos(l,minus_m,theta,phi,re_ylm_pos, im_ylm_pos) + tmp = (-1)**minus_m + re_ylm = tmp * re_ylm_pos + im_ylm = -tmp * im_ylm_pos ! complex conjugate + endif +end + +subroutine cartesian_to_spherical(r,theta,phi,r_abs) + implicit none + double precision, intent(in) :: r(3) + double precision, intent(out):: theta, phi,r_abs + double precision :: r_2,x_2_y_2,tmp + include 'constants.include.F' + x_2_y_2 = r(1)*r(1) + r(2)*r(2) + r_2 = x_2_y_2 + r(3)*r(3) + r_abs = dsqrt(r_2) + + if(r_abs.gt.1.d-20)then + theta = dacos(r(3)/r_abs) + else + theta = 0.d0 + endif + + if(.true.)then + if(dabs(r(1)).gt.0.d0)then + tmp = datan(r(2)/r(1)) +! phi = datan2(r(2),r(1)) + endif + ! From Wikipedia on Spherical Harmonics + if(r(1).gt.0.d0)then + phi = tmp + else if(r(1).lt.0.d0.and.r(2).ge.0.d0)then + phi = tmp + pi + else if(r(1).lt.0.d0.and.r(2).lt.0.d0)then + phi = tmp - pi + else if(r(1)==0.d0.and.r(2).gt.0.d0)then + phi = 0.5d0*pi + else if(r(1)==0.d0.and.r(2).lt.0.d0)then + phi =-0.5d0*pi + else if(r(1)==0.d0.and.r(2)==0.d0)then + phi = 0.d0 + endif + if(r(2).lt.0.d0.and.r(1).le.0.d0)then + tmp = pi - dabs(phi) + phi = pi + tmp + else if(r(2).lt.0.d0.and.r(1).gt.0.d0)then + phi = dtwo_pi + phi + endif + endif + + if(.false.)then + x_2_y_2 = dsqrt(x_2_y_2) + if(dabs(x_2_y_2).gt.1.d-20.and.dabs(r(2)).gt.1.d-20)then + phi = dabs(r(2))/r(2) * dacos(r(1)/x_2_y_2) + else + phi = 0.d0 + endif + endif +end + + +subroutine spher_harm_func_expl(l,m,theta,phi,re_ylm, im_ylm) + implicit none + BEGIN_DOC + ! Y_lm(theta,phi) with -l Date: Thu, 25 Apr 2024 20:00:42 +0200 Subject: [PATCH 2/8] Added properly the routines for the test of the Spherical Harmonics --- plugins/local/spher_harm/routines_test.irp.f | 227 +++++++++++++++++++ plugins/local/spher_harm/spher_harm.irp.f | 210 ----------------- 2 files changed, 227 insertions(+), 210 deletions(-) create mode 100644 plugins/local/spher_harm/routines_test.irp.f diff --git a/plugins/local/spher_harm/routines_test.irp.f b/plugins/local/spher_harm/routines_test.irp.f new file mode 100644 index 00000000..6f7cbc1c --- /dev/null +++ b/plugins/local/spher_harm/routines_test.irp.f @@ -0,0 +1,227 @@ + +subroutine test_cart + implicit none + BEGIN_DOC + ! test for the cartesian --> spherical change of coordinates + ! + ! simple test such that the polar angle theta ranges in [0,pi] + ! + ! and the asymuthal angle phi ranges in [0,2pi] + END_DOC + include 'constants.include.F' + double precision :: r(3),theta,phi,r_abs + print*,'' + r = 0.d0 + r(1) = 1.d0 + r(2) = 1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) =-1.d0 + r(2) = 1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) =-1.d0 + r(2) =-1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) = 1.d0 + r(2) =-1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi +end + +subroutine test_spher_harm + implicit none + BEGIN_DOC + ! routine to test the spherical harmonics integration on a sphere with the grid. + ! + ! We test = delta_m1,m2 delta_l1,l2 + END_DOC + include 'constants.include.F' + integer :: l1,m1,i,l2,m2,lmax + double precision :: r(3),weight,accu_re, accu_im,accu + double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 + l1 = 0 + m1 = 0 + l2 = 0 + m2 = 0 + lmax = 5 + do l1 = 0,lmax + do m1 = -l1 ,l1 + do l2 = 0,lmax + do m2 = -l2 ,l2 + accu_re = 0.d0 + accu_im = 0.d0 + ! = \int dOmega Y_l1,m1^* Y_l2,m2 + ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) + ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu = 0.d0 + do i = 1, n_points_integration_angular + double precision :: theta,phi,r_abs + r(1:3) = angular_quadrature_points(i,1:3) + weight = weights_angular_points(i) + call cartesian_to_spherical(r,theta,phi,r_abs) + if(theta.gt.pi.or.theta.lt.0.d0)then + print*,'pb with theta',theta + print*,r + endif + if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then + print*,'pb with phi',phi/pi + print*,r + endif + call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) + call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + write(33,'(10(F16.10,X))')phi/pi + enddo + ! Test for the delta l1,l2 and delta m1,m2 + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + enddo + enddo + enddo + enddo + double precision :: x,dx,xmax,xmin + integer:: nx + nx = 10000 + xmin = -5.d0 + xmax = 5.d0 + dx = (xmax - xmin)/dble(nx) + x = xmin + do i = 1, nx + write(34,*)x,datan(x),dacos(x) + x += dx + enddo +end + +subroutine test_brutal_spheric + implicit none + include 'constants.include.F' + BEGIN_DOC + ! test for the = delta_m1,m2 delta_l1,l2 using a two dimentional integration + ! + ! \int_0^2pi d Phi \int_-1^+1 d(cos(Theta)) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) + ! + != \int_0^2pi d Phi \int_0^pi dTheta sin(Theta) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) + ! + ! Allows to test for the general functions spher_harm_func_m_pos with spher_harm_func_expl + END_DOC + integer :: itheta, iphi,ntheta,nphi + double precision :: theta_min, theta_max, dtheta,theta + double precision :: phi_min, phi_max, dphi,phi + double precision :: accu_re, accu_im,weight + double precision :: re_ylm_1, im_ylm_1 ,re_ylm_2, im_ylm_2,accu + integer :: l1,m1,i,l2,m2,lmax + phi_min = 0.d0 + phi_max = 2.D0 * pi + theta_min = 0.d0 + theta_max = 1.D0 * pi + ntheta = 1000 + nphi = 1000 + dphi = (phi_max - phi_min)/dble(nphi) + dtheta = (theta_max - theta_min)/dble(ntheta) + + lmax = 3 + do l1 = 0,lmax + do m1 = 0 ,l1 + do l2 = 0,lmax + do m2 = 0 ,l2 + accu_re = 0.d0 + accu_im = 0.d0 + accu = 0.d0 + theta = theta_min + do itheta = 1, ntheta + phi = phi_min + do iphi = 1, nphi +! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) +! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) + call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) + call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) + weight = dtheta * dphi * dsin(theta) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + phi += dphi + enddo + theta += dtheta + enddo + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + print*,'accu = ',accu + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + endif + endif + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + endif + endif + enddo + enddo + enddo + enddo + + +end + +subroutine test_assoc_leg_pol + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + integer :: l1,m1,ngrid,i,l2,m2 + l1 = 0 + m1 = 0 + l2 = 2 + m2 = 0 + double precision :: x, dx,xmax,accu,xmin + double precision :: plgndr,func_1,func_2,ortho_assoc_gaus_pol + ngrid = 100000 + xmax = 1.d0 + xmin = -1.d0 + dx = (xmax-xmin)/dble(ngrid) + do l2 = 0,10 + x = xmin + accu = 0.d0 + do i = 1, ngrid + func_1 = plgndr(l1,m1,x) + func_2 = plgndr(l2,m2,x) + write(33,*)x, func_1,func_2 + accu += func_1 * func_2 * dx + x += dx + enddo + print*,'l2 = ',l2 + print*,'accu = ',accu + print*,ortho_assoc_gaus_pol(l1,m1,l2) + enddo +end diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f index 40661db1..e8deafb9 100644 --- a/plugins/local/spher_harm/spher_harm.irp.f +++ b/plugins/local/spher_harm/spher_harm.irp.f @@ -5,213 +5,3 @@ program spher_harm ! call test_brutal_spheric end -subroutine test_cart - implicit none - include 'constants.include.F' - double precision :: r(3),theta,phi,r_abs - print*,'' - r = 0.d0 - r(1) = 1.d0 - r(2) = 1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi - print*,'' - r = 0.d0 - r(1) =-1.d0 - r(2) = 1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi - print*,'' - r = 0.d0 - r(1) =-1.d0 - r(2) =-1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi - print*,'' - r = 0.d0 - r(1) = 1.d0 - r(2) =-1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi -end - -subroutine test_spher_harm - implicit none - include 'constants.include.F' - integer :: l1,m1,i,l2,m2,lmax - double precision :: r(3),weight,accu_re, accu_im,accu - double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 - l1 = 0 - m1 = 0 - l2 = 0 - m2 = 0 - lmax = 5 - do l1 = 0,lmax - do m1 = -l1 ,l1 - do l2 = 0,lmax - do m2 = -l2 ,l2 - accu_re = 0.d0 - accu_im = 0.d0 - ! = \int dOmega Y_l1,m1^* Y_l2,m2 - ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) - ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu = 0.d0 - do i = 1, n_points_integration_angular - double precision :: theta,phi,r_abs - r(1:3) = angular_quadrature_points(i,1:3) - weight = weights_angular_points(i) - call cartesian_to_spherical(r,theta,phi,r_abs) - if(theta.gt.pi.or.theta.lt.0.d0)then - print*,'pb with theta',theta - print*,r - endif - if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then - print*,'pb with phi',phi/pi - print*,r - endif - call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) - call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) -! call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) -! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) -! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) - accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) - accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu += weight - write(33,'(10(F16.10,X))')phi/pi - enddo - if(l1.ne.l2.or.m1.ne.m2)then - if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then - print*,'pb OFF DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - if(l1==l2.and.m1==m2)then - if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then - print*,'pb DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - enddo - enddo - enddo - enddo - double precision :: x,dx,xmax,xmin - integer:: nx - nx = 10000 - xmin = -5.d0 - xmax = 5.d0 - dx = (xmax - xmin)/dble(nx) - x = xmin - do i = 1, nx - write(34,*)x,datan(x),dacos(x) - x += dx - enddo -end - -subroutine test_brutal_spheric - implicit none - include 'constants.include.F' - integer :: itheta, iphi,ntheta,nphi - double precision :: theta_min, theta_max, dtheta,theta - double precision :: phi_min, phi_max, dphi,phi - double precision :: accu_re, accu_im,weight - double precision :: re_ylm_1, im_ylm_1 ,re_ylm_2, im_ylm_2,accu - integer :: l1,m1,i,l2,m2,lmax - phi_min = 0.d0 - phi_max = 2.D0 * pi - theta_min = 0.d0 - theta_max = 1.D0 * pi - ntheta = 1000 - nphi = 1000 - dphi = (phi_max - phi_min)/dble(nphi) - dtheta = (theta_max - theta_min)/dble(ntheta) - - lmax = 3 - do l1 = 0,lmax - do m1 = 0 ,l1 - do l2 = 0,lmax - do m2 = 0 ,l2 - accu_re = 0.d0 - accu_im = 0.d0 - accu = 0.d0 - theta = theta_min - do itheta = 1, ntheta - phi = phi_min - do iphi = 1, nphi -! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) - call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) - call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) - weight = dtheta * dphi * dsin(theta) - accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) - accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu += weight - phi += dphi - enddo - theta += dtheta - enddo - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - print*,'accu = ',accu - if(l1.ne.l2.or.m1.ne.m2)then - if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then - print*,'pb OFF DIAG !!!!! ' - endif - endif - if(l1==l2.and.m1==m2)then - if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then - print*,'pb DIAG !!!!! ' - endif - endif - enddo - enddo - enddo - enddo - - -end - -subroutine test_assoc_leg_pol - implicit none - BEGIN_DOC -! TODO : Put the documentation of the program here - END_DOC - print *, 'Hello world' - integer :: l1,m1,ngrid,i,l2,m2 - l1 = 0 - m1 = 0 - l2 = 2 - m2 = 0 - double precision :: x, dx,xmax,accu,xmin - double precision :: plgndr,func_1,func_2,ortho_assoc_gaus_pol - ngrid = 100000 - xmax = 1.d0 - xmin = -1.d0 - dx = (xmax-xmin)/dble(ngrid) - do l2 = 0,10 - x = xmin - accu = 0.d0 - do i = 1, ngrid - func_1 = plgndr(l1,m1,x) - func_2 = plgndr(l2,m2,x) - write(33,*)x, func_1,func_2 - accu += func_1 * func_2 * dx - x += dx - enddo - print*,'l2 = ',l2 - print*,'accu = ',accu - print*,ortho_assoc_gaus_pol(l1,m1,l2) - enddo -end From c3483df9a16003065a41bfa92d37274a3eb466ee Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 25 Apr 2024 20:00:42 +0200 Subject: [PATCH 3/8] Added properly the routines for the test of the Spherical Harmonics --- plugins/local/spher_harm/routines_test.irp.f | 227 +++++++++++++++++++ plugins/local/spher_harm/spher_harm.irp.f | 210 ----------------- 2 files changed, 227 insertions(+), 210 deletions(-) create mode 100644 plugins/local/spher_harm/routines_test.irp.f diff --git a/plugins/local/spher_harm/routines_test.irp.f b/plugins/local/spher_harm/routines_test.irp.f new file mode 100644 index 00000000..6f7cbc1c --- /dev/null +++ b/plugins/local/spher_harm/routines_test.irp.f @@ -0,0 +1,227 @@ + +subroutine test_cart + implicit none + BEGIN_DOC + ! test for the cartesian --> spherical change of coordinates + ! + ! simple test such that the polar angle theta ranges in [0,pi] + ! + ! and the asymuthal angle phi ranges in [0,2pi] + END_DOC + include 'constants.include.F' + double precision :: r(3),theta,phi,r_abs + print*,'' + r = 0.d0 + r(1) = 1.d0 + r(2) = 1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) =-1.d0 + r(2) = 1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) =-1.d0 + r(2) =-1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi + print*,'' + r = 0.d0 + r(1) = 1.d0 + r(2) =-1.d0 + call cartesian_to_spherical(r,theta,phi,r_abs) + print*,r + print*,phi/pi +end + +subroutine test_spher_harm + implicit none + BEGIN_DOC + ! routine to test the spherical harmonics integration on a sphere with the grid. + ! + ! We test = delta_m1,m2 delta_l1,l2 + END_DOC + include 'constants.include.F' + integer :: l1,m1,i,l2,m2,lmax + double precision :: r(3),weight,accu_re, accu_im,accu + double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 + l1 = 0 + m1 = 0 + l2 = 0 + m2 = 0 + lmax = 5 + do l1 = 0,lmax + do m1 = -l1 ,l1 + do l2 = 0,lmax + do m2 = -l2 ,l2 + accu_re = 0.d0 + accu_im = 0.d0 + ! = \int dOmega Y_l1,m1^* Y_l2,m2 + ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) + ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu = 0.d0 + do i = 1, n_points_integration_angular + double precision :: theta,phi,r_abs + r(1:3) = angular_quadrature_points(i,1:3) + weight = weights_angular_points(i) + call cartesian_to_spherical(r,theta,phi,r_abs) + if(theta.gt.pi.or.theta.lt.0.d0)then + print*,'pb with theta',theta + print*,r + endif + if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then + print*,'pb with phi',phi/pi + print*,r + endif + call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) + call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + write(33,'(10(F16.10,X))')phi/pi + enddo + ! Test for the delta l1,l2 and delta m1,m2 + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + enddo + enddo + enddo + enddo + double precision :: x,dx,xmax,xmin + integer:: nx + nx = 10000 + xmin = -5.d0 + xmax = 5.d0 + dx = (xmax - xmin)/dble(nx) + x = xmin + do i = 1, nx + write(34,*)x,datan(x),dacos(x) + x += dx + enddo +end + +subroutine test_brutal_spheric + implicit none + include 'constants.include.F' + BEGIN_DOC + ! test for the = delta_m1,m2 delta_l1,l2 using a two dimentional integration + ! + ! \int_0^2pi d Phi \int_-1^+1 d(cos(Theta)) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) + ! + != \int_0^2pi d Phi \int_0^pi dTheta sin(Theta) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) + ! + ! Allows to test for the general functions spher_harm_func_m_pos with spher_harm_func_expl + END_DOC + integer :: itheta, iphi,ntheta,nphi + double precision :: theta_min, theta_max, dtheta,theta + double precision :: phi_min, phi_max, dphi,phi + double precision :: accu_re, accu_im,weight + double precision :: re_ylm_1, im_ylm_1 ,re_ylm_2, im_ylm_2,accu + integer :: l1,m1,i,l2,m2,lmax + phi_min = 0.d0 + phi_max = 2.D0 * pi + theta_min = 0.d0 + theta_max = 1.D0 * pi + ntheta = 1000 + nphi = 1000 + dphi = (phi_max - phi_min)/dble(nphi) + dtheta = (theta_max - theta_min)/dble(ntheta) + + lmax = 3 + do l1 = 0,lmax + do m1 = 0 ,l1 + do l2 = 0,lmax + do m2 = 0 ,l2 + accu_re = 0.d0 + accu_im = 0.d0 + accu = 0.d0 + theta = theta_min + do itheta = 1, ntheta + phi = phi_min + do iphi = 1, nphi +! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) +! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) + call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) + call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) + weight = dtheta * dphi * dsin(theta) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + phi += dphi + enddo + theta += dtheta + enddo + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + print*,'accu = ',accu + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + endif + endif + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + endif + endif + enddo + enddo + enddo + enddo + + +end + +subroutine test_assoc_leg_pol + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + integer :: l1,m1,ngrid,i,l2,m2 + l1 = 0 + m1 = 0 + l2 = 2 + m2 = 0 + double precision :: x, dx,xmax,accu,xmin + double precision :: plgndr,func_1,func_2,ortho_assoc_gaus_pol + ngrid = 100000 + xmax = 1.d0 + xmin = -1.d0 + dx = (xmax-xmin)/dble(ngrid) + do l2 = 0,10 + x = xmin + accu = 0.d0 + do i = 1, ngrid + func_1 = plgndr(l1,m1,x) + func_2 = plgndr(l2,m2,x) + write(33,*)x, func_1,func_2 + accu += func_1 * func_2 * dx + x += dx + enddo + print*,'l2 = ',l2 + print*,'accu = ',accu + print*,ortho_assoc_gaus_pol(l1,m1,l2) + enddo +end diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f index 40661db1..e8deafb9 100644 --- a/plugins/local/spher_harm/spher_harm.irp.f +++ b/plugins/local/spher_harm/spher_harm.irp.f @@ -5,213 +5,3 @@ program spher_harm ! call test_brutal_spheric end -subroutine test_cart - implicit none - include 'constants.include.F' - double precision :: r(3),theta,phi,r_abs - print*,'' - r = 0.d0 - r(1) = 1.d0 - r(2) = 1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi - print*,'' - r = 0.d0 - r(1) =-1.d0 - r(2) = 1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi - print*,'' - r = 0.d0 - r(1) =-1.d0 - r(2) =-1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi - print*,'' - r = 0.d0 - r(1) = 1.d0 - r(2) =-1.d0 - call cartesian_to_spherical(r,theta,phi,r_abs) - print*,r - print*,phi/pi -end - -subroutine test_spher_harm - implicit none - include 'constants.include.F' - integer :: l1,m1,i,l2,m2,lmax - double precision :: r(3),weight,accu_re, accu_im,accu - double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 - l1 = 0 - m1 = 0 - l2 = 0 - m2 = 0 - lmax = 5 - do l1 = 0,lmax - do m1 = -l1 ,l1 - do l2 = 0,lmax - do m2 = -l2 ,l2 - accu_re = 0.d0 - accu_im = 0.d0 - ! = \int dOmega Y_l1,m1^* Y_l2,m2 - ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) - ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu = 0.d0 - do i = 1, n_points_integration_angular - double precision :: theta,phi,r_abs - r(1:3) = angular_quadrature_points(i,1:3) - weight = weights_angular_points(i) - call cartesian_to_spherical(r,theta,phi,r_abs) - if(theta.gt.pi.or.theta.lt.0.d0)then - print*,'pb with theta',theta - print*,r - endif - if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then - print*,'pb with phi',phi/pi - print*,r - endif - call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) - call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) -! call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) -! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) -! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) - accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) - accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu += weight - write(33,'(10(F16.10,X))')phi/pi - enddo - if(l1.ne.l2.or.m1.ne.m2)then - if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then - print*,'pb OFF DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - if(l1==l2.and.m1==m2)then - if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then - print*,'pb DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - enddo - enddo - enddo - enddo - double precision :: x,dx,xmax,xmin - integer:: nx - nx = 10000 - xmin = -5.d0 - xmax = 5.d0 - dx = (xmax - xmin)/dble(nx) - x = xmin - do i = 1, nx - write(34,*)x,datan(x),dacos(x) - x += dx - enddo -end - -subroutine test_brutal_spheric - implicit none - include 'constants.include.F' - integer :: itheta, iphi,ntheta,nphi - double precision :: theta_min, theta_max, dtheta,theta - double precision :: phi_min, phi_max, dphi,phi - double precision :: accu_re, accu_im,weight - double precision :: re_ylm_1, im_ylm_1 ,re_ylm_2, im_ylm_2,accu - integer :: l1,m1,i,l2,m2,lmax - phi_min = 0.d0 - phi_max = 2.D0 * pi - theta_min = 0.d0 - theta_max = 1.D0 * pi - ntheta = 1000 - nphi = 1000 - dphi = (phi_max - phi_min)/dble(nphi) - dtheta = (theta_max - theta_min)/dble(ntheta) - - lmax = 3 - do l1 = 0,lmax - do m1 = 0 ,l1 - do l2 = 0,lmax - do m2 = 0 ,l2 - accu_re = 0.d0 - accu_im = 0.d0 - accu = 0.d0 - theta = theta_min - do itheta = 1, ntheta - phi = phi_min - do iphi = 1, nphi -! call spher_harm_func_expl(l1,m1,theta,phi,re_ylm_1, im_ylm_1) -! call spher_harm_func_expl(l2,m2,theta,phi,re_ylm_2, im_ylm_2) - call spher_harm_func_m_pos(l1,m1,theta,phi,re_ylm_1, im_ylm_1) - call spher_harm_func_m_pos(l2,m2,theta,phi,re_ylm_2, im_ylm_2) - weight = dtheta * dphi * dsin(theta) - accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) - accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu += weight - phi += dphi - enddo - theta += dtheta - enddo - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - print*,'accu = ',accu - if(l1.ne.l2.or.m1.ne.m2)then - if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then - print*,'pb OFF DIAG !!!!! ' - endif - endif - if(l1==l2.and.m1==m2)then - if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then - print*,'pb DIAG !!!!! ' - endif - endif - enddo - enddo - enddo - enddo - - -end - -subroutine test_assoc_leg_pol - implicit none - BEGIN_DOC -! TODO : Put the documentation of the program here - END_DOC - print *, 'Hello world' - integer :: l1,m1,ngrid,i,l2,m2 - l1 = 0 - m1 = 0 - l2 = 2 - m2 = 0 - double precision :: x, dx,xmax,accu,xmin - double precision :: plgndr,func_1,func_2,ortho_assoc_gaus_pol - ngrid = 100000 - xmax = 1.d0 - xmin = -1.d0 - dx = (xmax-xmin)/dble(ngrid) - do l2 = 0,10 - x = xmin - accu = 0.d0 - do i = 1, ngrid - func_1 = plgndr(l1,m1,x) - func_2 = plgndr(l2,m2,x) - write(33,*)x, func_1,func_2 - accu += func_1 * func_2 * dx - x += dx - enddo - print*,'l2 = ',l2 - print*,'accu = ',accu - print*,ortho_assoc_gaus_pol(l1,m1,l2) - enddo -end From 5c69a7c005ecabe8428c386bf17bad3327891578 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 26 Apr 2024 10:57:57 +0200 Subject: [PATCH 4/8] removed stupid stuffs in spher_harm --- plugins/local/spher_harm/README.rst | 3 + plugins/local/spher_harm/routines_test.irp.f | 172 ++++++++++--------- plugins/local/spher_harm/spher_harm.irp.f | 4 +- 3 files changed, 93 insertions(+), 86 deletions(-) diff --git a/plugins/local/spher_harm/README.rst b/plugins/local/spher_harm/README.rst index bf897f73..9c9b12a6 100644 --- a/plugins/local/spher_harm/README.rst +++ b/plugins/local/spher_harm/README.rst @@ -2,3 +2,6 @@ spher_harm ========== +Routines for spherical Harmonics evaluation in real space. +The main routine is "spher_harm_func_r3(r,l,m,re_ylm, im_ylm)". +The test routine is "test_spher_harm" where everything is explained in details. diff --git a/plugins/local/spher_harm/routines_test.irp.f b/plugins/local/spher_harm/routines_test.irp.f index 6f7cbc1c..fe8fc422 100644 --- a/plugins/local/spher_harm/routines_test.irp.f +++ b/plugins/local/spher_harm/routines_test.irp.f @@ -1,10 +1,93 @@ +subroutine test_spher_harm + implicit none + BEGIN_DOC + ! routine to test the generic spherical harmonics routine "spher_harm_func_r3" from R^3 --> C + ! + ! We test = delta_m1,m2 delta_l1,l2 + ! + ! The test is done through the integration on a sphere with the Lebedev grid. + END_DOC + include 'constants.include.F' + integer :: l1,m1,i,l2,m2,lmax + double precision :: r(3),weight,accu_re, accu_im,accu + double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 + double precision :: theta,phi,r_abs + lmax = 5 ! Maximum angular momentum until which we are going to test orthogonality conditions + do l1 = 0,lmax + do m1 = -l1 ,l1 + do l2 = 0,lmax + do m2 = -l2 ,l2 + accu_re = 0.d0 ! accumulator for the REAL part of + accu_im = 0.d0 ! accumulator for the IMAGINARY part of + accu = 0.d0 ! accumulator for the weights ==> should be \int dOmega == 4 pi + ! = \int dOmega Y_l1,m1^* Y_l2,m2 + ! \approx \sum_i W_i Y_l1,m1^*(r_i) Y_l2,m2(r_i) WITH r_i being on the spher of radius 1 + do i = 1, n_points_integration_angular + r(1:3) = angular_quadrature_points(i,1:3) ! ith Lebedev point (x,y,z) on the sphere of radius 1 + weight = weights_angular_points(i) ! associated Lebdev weight not necessarily positive + +!!!!!!!!!!! Test of the Cartesian --> Spherical coordinates + ! theta MUST belong to [0,pi] and phi to [0,2pi] + ! gets the cartesian to spherical change of coordinates + call cartesian_to_spherical(r,theta,phi,r_abs) + if(theta.gt.pi.or.theta.lt.0.d0)then + print*,'pb with theta, it should be in [0,pi]',theta + print*,r + endif + if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then + print*,'pb with phi, it should be in [0,2 pi]',phi/pi + print*,r + endif + +!!!!!!!!!!! Routines returning the Spherical harmonics on the grid point + call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) + call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) + +!!!!!!!!!!! Integration of Y_l1,m1^*(r) Y_l2,m2(r) + ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) + ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) + accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) + accu += weight + enddo + ! Test that the sum of the weights is 4 pi + if(dabs(accu - dfour_pi).gt.1.d-6)then + print*,'Problem !! The sum of the Lebedev weight is not 4 pi ..' + print*,accu + stop + endif + ! Test for the delta l1,l2 and delta m1,m2 + ! + ! Test for the off-diagonal part of the Kronecker delta + if(l1.ne.l2.or.m1.ne.m2)then + if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then + print*,'pb OFF DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + ! Test for the diagonal part of the Kronecker delta + if(l1==l2.and.m1==m2)then + if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then + print*,'pb DIAG !!!!! ' + print*,'l1,m1,l2,m2',l1,m1,l2,m2 + print*,'accu_re = ',accu_re + print*,'accu_im = ',accu_im + endif + endif + enddo + enddo + enddo + enddo +end subroutine test_cart implicit none BEGIN_DOC ! test for the cartesian --> spherical change of coordinates ! - ! simple test such that the polar angle theta ranges in [0,pi] + ! test the routine "cartesian_to_spherical" such that the polar angle theta ranges in [0,pi] ! ! and the asymuthal angle phi ranges in [0,2pi] END_DOC @@ -40,97 +123,18 @@ subroutine test_cart print*,phi/pi end -subroutine test_spher_harm - implicit none - BEGIN_DOC - ! routine to test the spherical harmonics integration on a sphere with the grid. - ! - ! We test = delta_m1,m2 delta_l1,l2 - END_DOC - include 'constants.include.F' - integer :: l1,m1,i,l2,m2,lmax - double precision :: r(3),weight,accu_re, accu_im,accu - double precision :: re_ylm_1, im_ylm_1,re_ylm_2, im_ylm_2 - l1 = 0 - m1 = 0 - l2 = 0 - m2 = 0 - lmax = 5 - do l1 = 0,lmax - do m1 = -l1 ,l1 - do l2 = 0,lmax - do m2 = -l2 ,l2 - accu_re = 0.d0 - accu_im = 0.d0 - ! = \int dOmega Y_l1,m1^* Y_l2,m2 - ! = \int dOmega (re_ylm_1 -i im_ylm_1) * (re_ylm_2 +i im_ylm_2) - ! = \int dOmega (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) +i (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu = 0.d0 - do i = 1, n_points_integration_angular - double precision :: theta,phi,r_abs - r(1:3) = angular_quadrature_points(i,1:3) - weight = weights_angular_points(i) - call cartesian_to_spherical(r,theta,phi,r_abs) - if(theta.gt.pi.or.theta.lt.0.d0)then - print*,'pb with theta',theta - print*,r - endif - if(phi.gt.2.d0*pi.or.phi.lt.0.d0)then - print*,'pb with phi',phi/pi - print*,r - endif - call spher_harm_func_r3(r,l1,m1,re_ylm_1, im_ylm_1) - call spher_harm_func_r3(r,l2,m2,re_ylm_2, im_ylm_2) - accu_re += weight * (re_ylm_1*re_ylm_2 + im_ylm_1*im_ylm_2) - accu_im += weight * (im_ylm_2*re_ylm_1 - im_ylm_1*re_ylm_2) - accu += weight - write(33,'(10(F16.10,X))')phi/pi - enddo - ! Test for the delta l1,l2 and delta m1,m2 - if(l1.ne.l2.or.m1.ne.m2)then - if(dabs(accu_re).gt.1.d-6.or.dabs(accu_im).gt.1.d-6)then - print*,'pb OFF DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - if(l1==l2.and.m1==m2)then - if(dabs(accu_re-1.d0).gt.1.d-5.or.dabs(accu_im).gt.1.d-6)then - print*,'pb DIAG !!!!! ' - print*,'l1,m1,l2,m2',l1,m1,l2,m2 - print*,'accu_re = ',accu_re - print*,'accu_im = ',accu_im - endif - endif - enddo - enddo - enddo - enddo - double precision :: x,dx,xmax,xmin - integer:: nx - nx = 10000 - xmin = -5.d0 - xmax = 5.d0 - dx = (xmax - xmin)/dble(nx) - x = xmin - do i = 1, nx - write(34,*)x,datan(x),dacos(x) - x += dx - enddo -end subroutine test_brutal_spheric implicit none include 'constants.include.F' BEGIN_DOC - ! test for the = delta_m1,m2 delta_l1,l2 using a two dimentional integration + ! Test for the = delta_m1,m2 delta_l1,l2 using the following two dimentional integration ! ! \int_0^2pi d Phi \int_-1^+1 d(cos(Theta)) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) ! != \int_0^2pi d Phi \int_0^pi dTheta sin(Theta) Y_l1,m1^*(Theta,Phi) Y_l2,m2(Theta,Phi) ! - ! Allows to test for the general functions spher_harm_func_m_pos with spher_harm_func_expl + ! Allows to test for the general functions "spher_harm_func_m_pos" with "spher_harm_func_expl" END_DOC integer :: itheta, iphi,ntheta,nphi double precision :: theta_min, theta_max, dtheta,theta @@ -147,7 +151,7 @@ subroutine test_brutal_spheric dphi = (phi_max - phi_min)/dble(nphi) dtheta = (theta_max - theta_min)/dble(ntheta) - lmax = 3 + lmax = 2 do l1 = 0,lmax do m1 = 0 ,l1 do l2 = 0,lmax @@ -196,7 +200,7 @@ end subroutine test_assoc_leg_pol implicit none BEGIN_DOC -! TODO : Put the documentation of the program here +! Test for the associated Legendre Polynoms. The test is done through the orthogonality condition. END_DOC print *, 'Hello world' integer :: l1,m1,ngrid,i,l2,m2 diff --git a/plugins/local/spher_harm/spher_harm.irp.f b/plugins/local/spher_harm/spher_harm.irp.f index e8deafb9..7a2eea06 100644 --- a/plugins/local/spher_harm/spher_harm.irp.f +++ b/plugins/local/spher_harm/spher_harm.irp.f @@ -1,7 +1,7 @@ program spher_harm implicit none - call test_spher_harm +! call test_spher_harm ! call test_cart -! call test_brutal_spheric + call test_brutal_spheric end From 8eea5d7f7f142103998d8bfa1b3bcc630935f69b Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 15 May 2024 15:41:35 +0200 Subject: [PATCH 5/8] fixed a bug in cholesk_ao_transp --- .../tuto_plugins/tuto_I/test_cholesky.irp.f | 53 +++++++++++++++++++ src/ao_two_e_ints/cholesky.irp.f | 2 +- 2 files changed, 54 insertions(+), 1 deletion(-) create mode 100644 plugins/local/tuto_plugins/tuto_I/test_cholesky.irp.f diff --git a/plugins/local/tuto_plugins/tuto_I/test_cholesky.irp.f b/plugins/local/tuto_plugins/tuto_I/test_cholesky.irp.f new file mode 100644 index 00000000..d09d100a --- /dev/null +++ b/plugins/local/tuto_plugins/tuto_I/test_cholesky.irp.f @@ -0,0 +1,53 @@ +program my_program_to_print_stuffs + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + integer :: i,j,k,l,m + double precision :: integral, accu, accu_tot, integral_cholesky + double precision :: get_ao_two_e_integral, get_two_e_integral ! declaration of the functions + print*,'AO integrals, physicist notations : ' + accu_tot = 0.D0 + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + integral = get_ao_two_e_integral(i, j, k, l, ao_integrals_map) + integral_cholesky = 0.D0 + do m = 1, cholesky_ao_num + integral_cholesky += cholesky_ao_transp(m,i,k) * cholesky_ao_transp(m,j,l) + enddo + accu = dabs(integral_cholesky-integral) + accu_tot += accu + if(accu.gt.1.d-10)then + print*,i,j,k,l + print*,accu, integral, integral_cholesky + endif + enddo + enddo + enddo + enddo + print*,'accu_tot',accu_tot + + print*,'MO integrals, physicist notations : ' + do i = 1, mo_num + do j = 1, mo_num + do k = 1, mo_num + do l = 1, mo_num + integral = get_two_e_integral(i, j, k, l, mo_integrals_map) + accu = 0.D0 + integral_cholesky = 0.D0 + do m = 1, cholesky_mo_num + integral_cholesky += cholesky_mo_transp(m,i,k) * cholesky_mo_transp(m,j,l) + enddo + accu = dabs(integral_cholesky-integral) + accu_tot += accu + if(accu.gt.1.d-10)then + print*,i,j,k,l + print*,accu, integral, integral_cholesky + endif + enddo + enddo + enddo + enddo +end diff --git a/src/ao_two_e_ints/cholesky.irp.f b/src/ao_two_e_ints/cholesky.irp.f index 33304026..5fbd166c 100644 --- a/src/ao_two_e_ints/cholesky.irp.f +++ b/src/ao_two_e_ints/cholesky.irp.f @@ -6,7 +6,7 @@ BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num, integer :: i,j,k do j=1,ao_num do i=1,ao_num - do k=1,ao_num + do k=1,cholesky_ao_num cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k) enddo enddo From c6a61639445229eca3ecb2e32556ddef646064d6 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 16 May 2024 17:57:00 +0200 Subject: [PATCH 6/8] added f_hf with cholesky by default --- src/dft_utils_in_r/mo_in_r.irp.f | 2 +- src/mu_of_r/f_cholesky.irp.f | 221 +++++++++++++++++++++++++++ src/mu_of_r/mu_of_r_conditions.irp.f | 46 +++++- 3 files changed, 264 insertions(+), 5 deletions(-) create mode 100644 src/mu_of_r/f_cholesky.irp.f diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 192cb25a..ad931402 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -48,7 +48,7 @@ integer :: i,j do i = 1, n_points_final_grid do j = 1, mo_num - mos_in_r_array_transp(i,j) = mos_in_r_array(j,i) + mos_in_r_array_transp(i,j) = mos_in_r_array_omp(j,i) enddo enddo END_PROVIDER diff --git a/src/mu_of_r/f_cholesky.irp.f b/src/mu_of_r/f_cholesky.irp.f new file mode 100644 index 00000000..1ad4ce36 --- /dev/null +++ b/src/mu_of_r/f_cholesky.irp.f @@ -0,0 +1,221 @@ +BEGIN_PROVIDER [integer, list_couple_orb_r1, (2,n_couple_orb_r1)] + implicit none + integer :: ii,i,mm,m,itmp + itmp = 0 + do ii = 1, n_occ_val_orb_for_hf(1) + i = list_valence_orb_for_hf(ii,1) + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + itmp += 1 + list_couple_orb_r1(1,itmp) = i + list_couple_orb_r1(2,itmp) = m + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [integer, list_couple_orb_r2, (2,n_couple_orb_r2)] + implicit none + integer :: ii,i,mm,m,itmp + itmp = 0 + do ii = 1, n_occ_val_orb_for_hf(2) + i = list_valence_orb_for_hf(ii,2) + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + itmp += 1 + list_couple_orb_r2(1,itmp) = i + list_couple_orb_r2(2,itmp) = m + enddo + enddo +END_PROVIDER + + +BEGIN_PROVIDER [integer, n_couple_orb_r1] + implicit none + BEGIN_DOC + ! number of couples of alpha occupied times any basis orbital + END_DOC + n_couple_orb_r1 = n_occ_val_orb_for_hf(1) * n_basis_orb +END_PROVIDER + +BEGIN_PROVIDER [integer, n_couple_orb_r2] + implicit none + BEGIN_DOC + ! number of couples of beta occupied times any basis orbital + END_DOC + n_couple_orb_r2 = n_occ_val_orb_for_hf(2) * n_basis_orb +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mos_times_cholesky_r1, (cholesky_mo_num,n_points_final_grid)] + implicit none + BEGIN_DOC + ! V1_AR = \sum_{I}V_AI Phi_IR where "R" specifies the index of the grid point and A the number of cholesky point + ! + ! here Phi_IR is phi_i(R)xphi_b(R) for r1 and V_AI = (ib|A) chollesky vector + END_DOC + double precision, allocatable :: mos_ib_r1(:,:),mo_chol_r1(:,:) + double precision, allocatable :: test(:,:) + double precision :: mo_i_r1,mo_b_r1 + integer :: ii,i,mm,m,itmp,ipoint,ll + allocate(mos_ib_r1(n_couple_orb_r1,n_points_final_grid)) + allocate(mo_chol_r1(cholesky_mo_num,n_couple_orb_r1)) + + do ipoint = 1, n_points_final_grid + itmp = 0 + do ii = 1, n_occ_val_orb_for_hf(1) + i = list_valence_orb_for_hf(ii,1) + mo_i_r1 = mos_in_r_array_omp(i,ipoint) + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + mo_b_r1 = mos_in_r_array_omp(m,ipoint) + itmp += 1 + mos_ib_r1(itmp,ipoint) = mo_i_r1 * mo_b_r1 + enddo + enddo + enddo + + itmp = 0 + do ii = 1, n_occ_val_orb_for_hf(1) + i = list_valence_orb_for_hf(ii,1) + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + itmp += 1 + do ll = 1, cholesky_mo_num + mo_chol_r1(ll,itmp) = cholesky_mo_transp(ll,m,i) + enddo + enddo + enddo + + call get_AB_prod(mo_chol_r1,cholesky_mo_num,n_couple_orb_r1,mos_ib_r1,n_points_final_grid,mos_times_cholesky_r1) + allocate(test(cholesky_mo_num,n_points_final_grid)) + test = 0.d0 + do ipoint = 1, n_points_final_grid + do itmp = 1, n_couple_orb_r1 + i = list_couple_orb_r1(1,itmp) + m = list_couple_orb_r1(2,itmp) + mo_i_r1 = mos_in_r_array_omp(i,ipoint) + mo_b_r1 = mos_in_r_array_omp(m,ipoint) + do mm = 1, cholesky_mo_num + test(mm,ipoint) += mo_i_r1 * mo_b_r1 * mo_chol_r1(mm,itmp) + enddo + enddo + enddo + double precision :: accu + accu = 0.d0 + do ipoint = 1, n_points_final_grid + do mm = 1, cholesky_mo_num + accu += dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint) ) + if(dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint)).gt.1.d-10)then + print*,'problem ! ',dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint)) & + , mos_times_cholesky_r1(mm,ipoint) , test(mm,ipoint) + endif + enddo + enddo + print*,'accu = ',accu + + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mos_times_cholesky_r2, (cholesky_mo_num,n_points_final_grid)] + implicit none + BEGIN_DOC + ! V1_AR = \sum_{I}V_AI Phi_IR where "R" specifies the index of the grid point and A the number of cholesky point + ! + ! here Phi_IR is phi_i(R)xphi_b(R) for r2 and V_AI = (ib|A) chollesky vector + END_DOC + double precision, allocatable :: mos_ib_r2(:,:),mo_chol_r2(:,:) + double precision, allocatable :: test(:,:) + double precision :: mo_i_r2,mo_b_r2 + integer :: ii,i,mm,m,itmp,ipoint,ll + allocate(mos_ib_r2(n_couple_orb_r2,n_points_final_grid)) + allocate(mo_chol_r2(cholesky_mo_num,n_couple_orb_r2)) + + do ipoint = 1, n_points_final_grid + itmp = 0 + do ii = 1, n_occ_val_orb_for_hf(2) + i = list_valence_orb_for_hf(ii,2) + mo_i_r2 = mos_in_r_array_omp(i,ipoint) + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + mo_b_r2 = mos_in_r_array_omp(m,ipoint) + itmp += 1 + mos_ib_r2(itmp,ipoint) = mo_i_r2 * mo_b_r2 + enddo + enddo + enddo + + itmp = 0 + do ii = 1, n_occ_val_orb_for_hf(2) + i = list_valence_orb_for_hf(ii,2) + do mm = 1, n_basis_orb ! electron 1 + m = list_basis(mm) + itmp += 1 + do ll = 1, cholesky_mo_num + mo_chol_r2(ll,itmp) = cholesky_mo_transp(ll,m,i) + enddo + enddo + enddo + + call get_AB_prod(mo_chol_r2,cholesky_mo_num,n_couple_orb_r2,mos_ib_r2,n_points_final_grid,mos_times_cholesky_r2) + allocate(test(cholesky_mo_num,n_points_final_grid)) + test = 0.d0 + do ipoint = 1, n_points_final_grid + do itmp = 1, n_couple_orb_r2 + i = list_couple_orb_r2(1,itmp) + m = list_couple_orb_r2(2,itmp) + mo_i_r2 = mos_in_r_array_omp(i,ipoint) + mo_b_r2 = mos_in_r_array_omp(m,ipoint) + do mm = 1, cholesky_mo_num + test(mm,ipoint) += mo_i_r2 * mo_b_r2 * mo_chol_r2(mm,itmp) + enddo + enddo + enddo + double precision :: accu + accu = 0.d0 + do ipoint = 1, n_points_final_grid + do mm = 1, cholesky_mo_num + accu += dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint) ) + if(dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint)).gt.1.d-10)then + print*,'problem ! ',dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint)) & + , mos_times_cholesky_r2(mm,ipoint) , test(mm,ipoint) + endif + enddo + enddo + print*,'accu = ',accu + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, f_hf_cholesky, (n_points_final_grid)] + implicit none + integer :: ipoint + !!f(R) = \sum_{I} \sum_{J} Phi_I(R) Phi_J(R) V_IJ + !! = \sum_{I}\sum_{J}\sum_A Phi_I(R) Phi_J(R) V_AI V_AJ + !! = \sum_A \sum_{I}Phi_I(R)V_AI \sum_{J}V_AJ Phi_J(R) + !! = \sum_A V_AR G_AR + !! V_AR = \sum_{I}Phi_IR V_AI = \sum_{I}Phi^t_RI V_AI + double precision :: u_dot_v + do ipoint = 1, n_points_final_grid + f_hf_cholesky(ipoint) = 2.D0 * u_dot_v(mos_times_cholesky_r2(1,ipoint),mos_times_cholesky_r1(1,ipoint),cholesky_mo_num) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)] + implicit none + integer :: ipoint,i,ii + double precision :: dm_a, dm_b + do ipoint = 1, n_points_final_grid + dm_a = 0.d0 + do ii = 1, n_occ_val_orb_for_hf(1) + i = list_valence_orb_for_hf(ii,1) + dm_a += mos_in_r_array_omp(i,ipoint)*mos_in_r_array_omp(i,ipoint) + enddo + dm_b = 0.d0 + do ii = 1, n_occ_val_orb_for_hf(2) + i = list_valence_orb_for_hf(ii,2) + dm_b += mos_in_r_array_omp(i,ipoint)*mos_in_r_array_omp(i,ipoint) + enddo + on_top_hf_grid(ipoint) = 2.D0 * dm_a*dm_b + enddo +END_PROVIDER + diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index 6b49b9df..5b4d4b83 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -61,7 +61,7 @@ END_DOC integer :: ipoint double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + PROVIDE f_hf_cholesky on_top_hf_grid print*,'providing mu_of_r_hf ...' call wall_time(wall0) sqpi = dsqrt(dacos(-1.d0)) @@ -69,10 +69,10 @@ !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & - !$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_psi_hf_ab,on_top_hf_mu_r,sqpi) + !$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_hf_cholesky,on_top_hf_grid,sqpi) do ipoint = 1, n_points_final_grid - f_hf = f_psi_hf_ab(ipoint) - on_top = on_top_hf_mu_r(ipoint) + f_hf = f_hf_cholesky(ipoint) + on_top = on_top_hf_grid(ipoint) if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then w_hf = 1.d+10 else @@ -85,6 +85,44 @@ print*,'Time to provide mu_of_r_hf = ',wall1-wall0 END_PROVIDER + BEGIN_PROVIDER [double precision, mu_of_r_hf_old, (n_points_final_grid) ] + implicit none + BEGIN_DOC + ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO) + ! + ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the two-body density matrix are excluded + END_DOC + integer :: ipoint + double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + print*,'providing mu_of_r_hf_old ...' + call wall_time(wall0) + sqpi = dsqrt(dacos(-1.d0)) + provide f_psi_hf_ab + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & + !$OMP ShARED (n_points_final_grid,mu_of_r_hf_old,f_psi_hf_ab,on_top_hf_mu_r,sqpi) + do ipoint = 1, n_points_final_grid + f_hf = f_psi_hf_ab(ipoint) + on_top = on_top_hf_mu_r(ipoint) + if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then + w_hf = 1.d+10 + else + w_hf = f_hf / on_top + endif + mu_of_r_hf_old(ipoint) = w_hf * sqpi * 0.5d0 + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide mu_of_r_hf_old = ',wall1-wall0 + END_PROVIDER + + BEGIN_PROVIDER [double precision, mu_of_r_psi_cas, (n_points_final_grid,N_states) ] implicit none BEGIN_DOC From ce042fbd787a21a600830596fa3caa5f7aa2cdb1 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 21 May 2024 12:01:28 +0200 Subject: [PATCH 7/8] basis set correction with cholesky works for hf --- .../local/basis_correction/51.basis_c.bats | 8 -- .../{01.convert.bats => convert_bats_old} | 0 src/hartree_fock/10.hf.bats | 13 -- src/mu_of_r/basis_def.irp.f | 45 +++++++ .../{f_cholesky.irp.f => f_hf_cholesky.irp.f} | 121 +++++++++--------- 5 files changed, 104 insertions(+), 83 deletions(-) rename src/ezfio_files/{01.convert.bats => convert_bats_old} (100%) rename src/mu_of_r/{f_cholesky.irp.f => f_hf_cholesky.irp.f} (67%) diff --git a/plugins/local/basis_correction/51.basis_c.bats b/plugins/local/basis_correction/51.basis_c.bats index 914b482b..1e20bae3 100644 --- a/plugins/local/basis_correction/51.basis_c.bats +++ b/plugins/local/basis_correction/51.basis_c.bats @@ -37,14 +37,6 @@ function run_sd() { eq $energy1 $1 $thresh } -@test "O2 CAS" { - qp set_file o2_cas.gms.ezfio - qp set_mo_class -c "[1-2]" -a "[3-10]" -d "[11-46]" - run -149.72435425 3.e-4 10000 - qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]" - run_md -0.1160222327 1.e-6 -} - @test "LiF RHF" { qp set_file lif.ezfio diff --git a/src/ezfio_files/01.convert.bats b/src/ezfio_files/convert_bats_old similarity index 100% rename from src/ezfio_files/01.convert.bats rename to src/ezfio_files/convert_bats_old diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index b496a089..214dfa86 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -115,9 +115,6 @@ rm -rf $EZFIO run hco.ezfio -113.1841002944744 } -@test "HBO" { # 0.805600 1.4543s - run hbo.ezfio -100.018582259096 -} @test "H2S" { # 1.655600 4.21402s run h2s.ezfio -398.6944130421982 @@ -127,9 +124,6 @@ rm -rf $EZFIO run h3coh.ezfio -114.9865030596373 } -@test "H2O" { # 1.811100 1.84387s - run h2o.ezfio -0.760270218692179E+02 -} @test "H2O2" { # 2.217000 8.50267s run h2o2.ezfio -150.7806608469964 @@ -187,13 +181,6 @@ rm -rf $EZFIO run oh.ezfio -75.42025413469165 } -@test "[Cu(NH3)4]2+" { # 59.610100 4.18766m - [[ -n $TRAVIS ]] && skip - qp set_file cu_nh3_4_2plus.ezfio - qp set scf_utils thresh_scf 1.e-10 - run cu_nh3_4_2plus.ezfio -1862.97590358903 -} - @test "SO2" { # 71.894900 3.22567m [[ -n $TRAVIS ]] && skip run so2.ezfio -41.55800401346361 diff --git a/src/mu_of_r/basis_def.irp.f b/src/mu_of_r/basis_def.irp.f index fff9f581..e433f4d8 100644 --- a/src/mu_of_r/basis_def.irp.f +++ b/src/mu_of_r/basis_def.irp.f @@ -114,3 +114,48 @@ BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_fi enddo enddo END_PROVIDER + +! BEGIN_PROVIDER [integer, n_docc_val_orb_for_cas] +!&BEGIN_PROVIDER [integer, n_max_docc_val_orb_for_cas] +! implicit none +! BEGIN_DOC +! ! Number of DOUBLY OCCUPIED VALENCE ORBITALS for the CAS wave function +! ! +! ! This determines the size of the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 +! END_DOC +! integer :: i +! n_docc_val_orb_for_cas = 0 +! ! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL +! do i = 1, elec_beta_num +! if( trim(mo_class(i))=="Inactive" & +! .or. trim(mo_class(i))=="Active" & +! .or. trim(mo_class(i))=="Virtual" )then +! n_docc_val_orb_for_cas +=1 +! endif +! enddo +! n_max_docc_val_orb_for_cas = maxval(n_docc_val_orb_for_cas) +! +!END_PROVIDER +! +!BEGIN_PROVIDER [integer, list_doc_valence_orb_for_cas, (n_max_docc_val_orb_for_cas)] +! implicit none +! BEGIN_DOC +! ! List of OCCUPIED valence orbitals for each spin to build the f_{HF}(r_1,r_2) function +! ! +! ! This corresponds to ALL OCCUPIED orbitals in the HF wave function, except those defined as "core" +! ! +! ! This determines the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937 +! END_DOC +! j = 0 +! ! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL +! do i = 1, elec_beta_num +! if( trim(mo_class(i))=="Inactive" & +! .or. trim(mo_class(i))=="Active" & +! .or. trim(mo_class(i))=="Virtual" )then +! j +=1 +! list_doc_valence_orb_for_cas(j) = i +! endif +! enddo +! +!END_PROVIDER + diff --git a/src/mu_of_r/f_cholesky.irp.f b/src/mu_of_r/f_hf_cholesky.irp.f similarity index 67% rename from src/mu_of_r/f_cholesky.irp.f rename to src/mu_of_r/f_hf_cholesky.irp.f index 1ad4ce36..84097f09 100644 --- a/src/mu_of_r/f_cholesky.irp.f +++ b/src/mu_of_r/f_hf_cholesky.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [integer, list_couple_orb_r1, (2,n_couple_orb_r1)] +BEGIN_PROVIDER [integer, list_couple_hf_orb_r1, (2,n_couple_orb_r1)] implicit none integer :: ii,i,mm,m,itmp itmp = 0 @@ -7,14 +7,14 @@ BEGIN_PROVIDER [integer, list_couple_orb_r1, (2,n_couple_orb_r1)] do mm = 1, n_basis_orb ! electron 1 m = list_basis(mm) itmp += 1 - list_couple_orb_r1(1,itmp) = i - list_couple_orb_r1(2,itmp) = m + list_couple_hf_orb_r1(1,itmp) = i + list_couple_hf_orb_r1(2,itmp) = m enddo enddo END_PROVIDER -BEGIN_PROVIDER [integer, list_couple_orb_r2, (2,n_couple_orb_r2)] +BEGIN_PROVIDER [integer, list_couple_hf_orb_r2, (2,n_couple_orb_r2)] implicit none integer :: ii,i,mm,m,itmp itmp = 0 @@ -23,8 +23,8 @@ BEGIN_PROVIDER [integer, list_couple_orb_r2, (2,n_couple_orb_r2)] do mm = 1, n_basis_orb ! electron 1 m = list_basis(mm) itmp += 1 - list_couple_orb_r2(1,itmp) = i - list_couple_orb_r2(2,itmp) = m + list_couple_hf_orb_r2(1,itmp) = i + list_couple_hf_orb_r2(2,itmp) = m enddo enddo END_PROVIDER @@ -87,31 +87,6 @@ BEGIN_PROVIDER [ double precision, mos_times_cholesky_r1, (cholesky_mo_num,n_poi enddo call get_AB_prod(mo_chol_r1,cholesky_mo_num,n_couple_orb_r1,mos_ib_r1,n_points_final_grid,mos_times_cholesky_r1) - allocate(test(cholesky_mo_num,n_points_final_grid)) - test = 0.d0 - do ipoint = 1, n_points_final_grid - do itmp = 1, n_couple_orb_r1 - i = list_couple_orb_r1(1,itmp) - m = list_couple_orb_r1(2,itmp) - mo_i_r1 = mos_in_r_array_omp(i,ipoint) - mo_b_r1 = mos_in_r_array_omp(m,ipoint) - do mm = 1, cholesky_mo_num - test(mm,ipoint) += mo_i_r1 * mo_b_r1 * mo_chol_r1(mm,itmp) - enddo - enddo - enddo - double precision :: accu - accu = 0.d0 - do ipoint = 1, n_points_final_grid - do mm = 1, cholesky_mo_num - accu += dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint) ) - if(dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint)).gt.1.d-10)then - print*,'problem ! ',dabs(mos_times_cholesky_r1(mm,ipoint) - test(mm,ipoint)) & - , mos_times_cholesky_r1(mm,ipoint) , test(mm,ipoint) - endif - enddo - enddo - print*,'accu = ',accu END_PROVIDER @@ -157,53 +132,72 @@ BEGIN_PROVIDER [ double precision, mos_times_cholesky_r2, (cholesky_mo_num,n_poi enddo call get_AB_prod(mo_chol_r2,cholesky_mo_num,n_couple_orb_r2,mos_ib_r2,n_points_final_grid,mos_times_cholesky_r2) - allocate(test(cholesky_mo_num,n_points_final_grid)) - test = 0.d0 - do ipoint = 1, n_points_final_grid - do itmp = 1, n_couple_orb_r2 - i = list_couple_orb_r2(1,itmp) - m = list_couple_orb_r2(2,itmp) - mo_i_r2 = mos_in_r_array_omp(i,ipoint) - mo_b_r2 = mos_in_r_array_omp(m,ipoint) - do mm = 1, cholesky_mo_num - test(mm,ipoint) += mo_i_r2 * mo_b_r2 * mo_chol_r2(mm,itmp) - enddo - enddo - enddo - double precision :: accu - accu = 0.d0 - do ipoint = 1, n_points_final_grid - do mm = 1, cholesky_mo_num - accu += dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint) ) - if(dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint)).gt.1.d-10)then - print*,'problem ! ',dabs(mos_times_cholesky_r2(mm,ipoint) - test(mm,ipoint)) & - , mos_times_cholesky_r2(mm,ipoint) , test(mm,ipoint) - endif - enddo - enddo - print*,'accu = ',accu END_PROVIDER BEGIN_PROVIDER [ double precision, f_hf_cholesky, (n_points_final_grid)] implicit none - integer :: ipoint + integer :: ipoint,m,k !!f(R) = \sum_{I} \sum_{J} Phi_I(R) Phi_J(R) V_IJ !! = \sum_{I}\sum_{J}\sum_A Phi_I(R) Phi_J(R) V_AI V_AJ !! = \sum_A \sum_{I}Phi_I(R)V_AI \sum_{J}V_AJ Phi_J(R) !! = \sum_A V_AR G_AR !! V_AR = \sum_{I}Phi_IR V_AI = \sum_{I}Phi^t_RI V_AI - double precision :: u_dot_v - do ipoint = 1, n_points_final_grid - f_hf_cholesky(ipoint) = 2.D0 * u_dot_v(mos_times_cholesky_r2(1,ipoint),mos_times_cholesky_r1(1,ipoint),cholesky_mo_num) - enddo + double precision :: u_dot_v,wall0,wall1 + if(elec_alpha_num == elec_beta_num)then + provide mos_times_cholesky_r1 + print*,'providing f_hf_cholesky ...' + call wall_time(wall0) + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,m) & + !$OMP ShARED (mos_times_cholesky_r1,cholesky_mo_num,f_hf_cholesky,n_points_final_grid) + do ipoint = 1, n_points_final_grid + f_hf_cholesky(ipoint) = 0.d0 + do m = 1, cholesky_mo_num + f_hf_cholesky(ipoint) = f_hf_cholesky(ipoint) + & + mos_times_cholesky_r1(m,ipoint) * mos_times_cholesky_r1(m,ipoint) + enddo + f_hf_cholesky(ipoint) *= 2.D0 + enddo + !$OMP END PARALLEL DO + + call wall_time(wall1) + print*,'Time to provide f_hf_cholesky = ',wall1-wall0 + free mos_times_cholesky_r1 + else + provide mos_times_cholesky_r2 mos_times_cholesky_r1 + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,m) & + !$OMP ShARED (mos_times_cholesky_r2,mos_times_cholesky_r1,cholesky_mo_num,f_hf_cholesky,n_points_final_grid) + do ipoint = 1, n_points_final_grid + f_hf_cholesky(ipoint) = 0.D0 + do m = 1, cholesky_mo_num + f_hf_cholesky(ipoint) = f_hf_cholesky(ipoint) + & + mos_times_cholesky_r2(m,ipoint)*mos_times_cholesky_r1(m,ipoint) + enddo + f_hf_cholesky(ipoint) *= 2.D0 + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide f_hf_cholesky = ',wall1-wall0 + free mos_times_cholesky_r2 mos_times_cholesky_r1 + endif END_PROVIDER BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)] implicit none integer :: ipoint,i,ii - double precision :: dm_a, dm_b + double precision :: dm_a, dm_b,wall0,wall1 + print*,'providing on_top_hf_grid ...' + provide mos_in_r_array_omp + call wall_time(wall0) + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,dm_a,dm_b,ii,i) & + !$OMP ShARED (n_points_final_grid,n_occ_val_orb_for_hf,mos_in_r_array_omp,list_valence_orb_for_hf,on_top_hf_grid) do ipoint = 1, n_points_final_grid dm_a = 0.d0 do ii = 1, n_occ_val_orb_for_hf(1) @@ -217,5 +211,8 @@ BEGIN_PROVIDER [ double precision, on_top_hf_grid, (n_points_final_grid)] enddo on_top_hf_grid(ipoint) = 2.D0 * dm_a*dm_b enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide on_top_hf_grid = ',wall1-wall0 END_PROVIDER From 112f113ccb3f363262930b53e21aed010a29f746 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 21 May 2024 12:26:30 +0200 Subject: [PATCH 8/8] fixed forgotten stuffs in normal_order_old/NEED --- plugins/local/normal_order_old/NEED | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/local/normal_order_old/NEED b/plugins/local/normal_order_old/NEED index 8b137891..e8c8c478 100644 --- a/plugins/local/normal_order_old/NEED +++ b/plugins/local/normal_order_old/NEED @@ -1 +1 @@ - +tc_scf