mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-07 22:53:52 +01:00
Merge pull request #337 from QuantumPackage/dev-spher_harm
Dev spher harm
This commit is contained in:
Emmanuel Giner
GitHub
commit
8c4183cf6e
@ -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
|
||||
|
||||
@ -1 +1 @@
|
||||
|
||||
tc_scf
|
||||
|
||||
59
plugins/local/spher_harm/.gitignore
vendored
Normal file
59
plugins/local/spher_harm/.gitignore
vendored
Normal file
@ -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
|
||||
1
plugins/local/spher_harm/NEED
Normal file
1
plugins/local/spher_harm/NEED
Normal file
@ -0,0 +1 @@
|
||||
dft_utils_in_r
|
||||
7
plugins/local/spher_harm/README.rst
Normal file
7
plugins/local/spher_harm/README.rst
Normal file
@ -0,0 +1,7 @@
|
||||
==========
|
||||
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.
|
||||
50
plugins/local/spher_harm/assoc_gaus_pol.irp.f
Normal file
50
plugins/local/spher_harm/assoc_gaus_pol.irp.f
Normal file
@ -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
|
||||
231
plugins/local/spher_harm/routines_test.irp.f
Normal file
231
plugins/local/spher_harm/routines_test.irp.f
Normal file
@ -0,0 +1,231 @@
|
||||
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 <Y_l1,m1|Y_l2,m2> = 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 <Y_l1,m1|Y_l2,m2>
|
||||
accu_im = 0.d0 ! accumulator for the IMAGINARY part of <Y_l1,m1|Y_l2,m2>
|
||||
accu = 0.d0 ! accumulator for the weights ==> should be \int dOmega == 4 pi
|
||||
! <l1,m1|l2,m2> = \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
|
||||
!
|
||||
! 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
|
||||
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_brutal_spheric
|
||||
implicit none
|
||||
include 'constants.include.F'
|
||||
BEGIN_DOC
|
||||
! Test for the <Y_l1,m1|Y_l2,m2> = 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"
|
||||
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 = 2
|
||||
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
|
||||
! 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
|
||||
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
|
||||
7
plugins/local/spher_harm/spher_harm.irp.f
Normal file
7
plugins/local/spher_harm/spher_harm.irp.f
Normal file
@ -0,0 +1,7 @@
|
||||
program spher_harm
|
||||
implicit none
|
||||
! call test_spher_harm
|
||||
! call test_cart
|
||||
call test_brutal_spheric
|
||||
end
|
||||
|
||||
151
plugins/local/spher_harm/spher_harm_func.irp.f
Normal file
151
plugins/local/spher_harm/spher_harm_func.irp.f
Normal file
@ -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<m<+l
|
||||
!
|
||||
END_DOC
|
||||
double precision, intent(in) :: theta, phi
|
||||
integer, intent(in) :: l,m
|
||||
double precision, intent(out):: re_ylm,im_ylm
|
||||
double precision :: re_ylm_pos,im_ylm_pos,tmp
|
||||
integer :: minus_m
|
||||
if(abs(m).gt.l)then
|
||||
print*,'|m| > 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<m<+l and 0<= l <=2
|
||||
!
|
||||
END_DOC
|
||||
double precision, intent(in) :: theta, phi
|
||||
integer, intent(in) :: l,m
|
||||
double precision, intent(out):: re_ylm,im_ylm
|
||||
double precision :: tmp
|
||||
include 'constants.include.F'
|
||||
if(l==0.and.m==0)then
|
||||
re_ylm = 0.5d0 * inv_sq_pi
|
||||
im_ylm = 0.d0
|
||||
else if(l==1.and.m==1)then
|
||||
tmp = - inv_sq_pi * dsqrt(3.d0/8.d0) * dsin(theta)
|
||||
re_ylm = tmp * dcos(phi)
|
||||
im_ylm = tmp * dsin(phi)
|
||||
else if(l==1.and.m==0)then
|
||||
tmp = inv_sq_pi * dsqrt(3.d0/4.d0) * dcos(theta)
|
||||
re_ylm = tmp
|
||||
im_ylm = 0.d0
|
||||
else if(l==2.and.m==2)then
|
||||
tmp = 0.25d0 * inv_sq_pi * dsqrt(0.5d0*15.d0) * dsin(theta)*dsin(theta)
|
||||
re_ylm = tmp * dcos(2.d0*phi)
|
||||
im_ylm = tmp * dsin(2.d0*phi)
|
||||
else if(l==2.and.m==1)then
|
||||
tmp = - inv_sq_pi * dsqrt(15.d0/8.d0) * dsin(theta) * dcos(theta)
|
||||
re_ylm = tmp * dcos(phi)
|
||||
im_ylm = tmp * dsin(phi)
|
||||
else if(l==2.and.m==0)then
|
||||
tmp = dsqrt(5.d0/4.d0) * inv_sq_pi* (1.5d0*dcos(theta)*dcos(theta)-0.5d0)
|
||||
re_ylm = tmp
|
||||
im_ylm = 0.d0
|
||||
endif
|
||||
end
|
||||
53
plugins/local/tuto_plugins/tuto_I/test_cholesky.irp.f
Normal file
53
plugins/local/tuto_plugins/tuto_I/test_cholesky.irp.f
Normal file
@ -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 : <i j|k l>'
|
||||
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 : <i j|k l>'
|
||||
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
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
218
src/mu_of_r/f_hf_cholesky.irp.f
Normal file
218
src/mu_of_r/f_hf_cholesky.irp.f
Normal file
@ -0,0 +1,218 @@
|
||||
BEGIN_PROVIDER [integer, list_couple_hf_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_hf_orb_r1(1,itmp) = i
|
||||
list_couple_hf_orb_r1(2,itmp) = m
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [integer, list_couple_hf_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_hf_orb_r2(1,itmp) = i
|
||||
list_couple_hf_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)
|
||||
|
||||
|
||||
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)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, f_hf_cholesky, (n_points_final_grid)]
|
||||
implicit none
|
||||
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,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,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)
|
||||
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
|
||||
!$OMP END PARALLEL DO
|
||||
call wall_time(wall1)
|
||||
print*,'Time to provide on_top_hf_grid = ',wall1-wall0
|
||||
END_PROVIDER
|
||||
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user