10
1
mirror of https://gitlab.com/scemama/qmcchem.git synced 2024-11-08 07:03:39 +01:00
qmcchem/src/pseudo.irp.f

384 lines
9.0 KiB
FortranFixed
Raw Normal View History

2015-12-19 03:29:52 +01:00
BEGIN_PROVIDER [ double precision, pseudo_v_k , (nucl_num,pseudo_klocmax) ]
implicit none
BEGIN_DOC
! V_k
END_DOC
call get_pseudo_pseudo_v_k(pseudo_v_k)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_v_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ]
implicit none
BEGIN_DOC
! V_Kl
END_DOC
call get_pseudo_pseudo_v_kl(pseudo_v_kl)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_kmax ]
implicit none
BEGIN_DOC
! kmax
END_DOC
call get_pseudo_pseudo_kmax(pseudo_kmax)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_dz_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ]
implicit none
BEGIN_DOC
! Exponents in the non-local part of the pseudo-potential
END_DOC
call get_pseudo_pseudo_dz_kl(pseudo_dz_kl)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_klocmax ]
implicit none
BEGIN_DOC
! klocmax
END_DOC
call get_pseudo_pseudo_klocmax(pseudo_klocmax)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_lmax ]
implicit none
BEGIN_DOC
! Max value of l in the pseudo-potential
END_DOC
call get_pseudo_pseudo_lmax(pseudo_lmax)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_grid_size ]
implicit none
BEGIN_DOC
! Size of the QMC grid (number of points)
END_DOC
call get_pseudo_pseudo_grid_size(pseudo_grid_size)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_grid_rmax ]
implicit none
BEGIN_DOC
! Size of the QMC grid (max distance)
END_DOC
call get_pseudo_pseudo_grid_rmax(pseudo_grid_rmax)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_non_loc_dim ]
&BEGIN_PROVIDER [ integer, pseudo_non_loc_dim_8 ]
&BEGIN_PROVIDER [ integer, pseudo_non_loc_dim_count, (nucl_num) ]
implicit none
BEGIN_DOC
! Dimension of of pseudo_non_local arrays
END_DOC
pseudo_non_loc_dim = 0
pseudo_non_loc_dim_count = 0
integer :: k,l,m
do k=1,nucl_num
do l=0,pseudo_lmax
pseudo_non_loc_dim_count(k) += 2*l+1
enddo
enddo
pseudo_non_loc_dim = sum(pseudo_non_loc_dim_count)
integer, external :: mod_align
pseudo_non_loc_dim_8 = mod_align(pseudo_non_loc_dim)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_n_kl , (nucl_num,pseudo_kmax,0:pseudo_lmax) ]
implicit none
BEGIN_DOC
! n_kl
END_DOC
call get_pseudo_pseudo_n_kl(pseudo_n_kl)
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_dz_k , (nucl_num,pseudo_klocmax) ]
implicit none
BEGIN_DOC
! dz_k
END_DOC
call get_pseudo_pseudo_dz_k(pseudo_dz_k)
END_PROVIDER
BEGIN_PROVIDER [ integer, pseudo_n_k , (nucl_num,pseudo_klocmax) ]
implicit none
BEGIN_DOC
! n_k
END_DOC
call get_pseudo_pseudo_n_k(pseudo_n_k)
END_PROVIDER
BEGIN_PROVIDER [ logical, do_pseudo ]
implicit none
BEGIN_DOC
! Using pseudo potential integral of not
END_DOC
call get_pseudo_do_pseudo(do_pseudo)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_grid, (ao_num, -pseudo_lmax:pseudo_lmax, 0:pseudo_lmax, nucl_num, pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Pseudopotential grid points
END_DOC
call get_pseudo_ao_pseudo_grid(ao_pseudo_grid)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_pseudo_grid, (ao_num, -pseudo_lmax:pseudo_lmax, 0:pseudo_lmax, nucl_num, pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Pseudopotential grid points
END_DOC
call get_pseudo_mo_pseudo_grid(mo_pseudo_grid)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_pseudo_grid_scaled, (pseudo_non_loc_dim_8,ao_num,pseudo_grid_size) ]
implicit none
BEGIN_DOC
! Pseudopotential grid points
END_DOC
integer :: i,k,l,m,kk,n
double precision :: c
c = 1.d0/mo_scale
do n=1,pseudo_grid_size
do i=1,ao_num
kk = 0
do k=1,nucl_num
do l=0,pseudo_lmax
do m=-l,l
kk = kk+1
mo_pseudo_grid_scaled(kk,i,n) = c * mo_pseudo_grid(i,m,l,k,n)
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, v_pseudo_local, (elec_num) ]
implicit none
BEGIN_DOC
! Local component of the pseudo-potential
END_DOC
integer :: i,k,l
double precision :: r, rn, alpha
do l=1,elec_num
v_pseudo_local(l) = 0.d0
do k=1,pseudo_klocmax
do i=1,nucl_num
r = nucl_elec_dist(i,l)
alpha = pseudo_dz_k(i,k)*r*r
if (alpha > 20.d0) then
cycle
endif
select case (pseudo_n_k(i,k))
case (-2)
rn = nucl_elec_dist_inv(i,l)*nucl_elec_dist_inv(i,l)
case (-1)
rn = nucl_elec_dist_inv(i,l)
case (0)
rn = 1.d0
case (1)
rn = r
case (2)
rn = r*r
case default
rn = r**pseudo_n_k(i,k)
end select
v_pseudo_local(l) += pseudo_v_k(i,k) * exp(-alpha) * rn
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, v_pseudo_non_local, (pseudo_non_loc_dim_8,elec_num) ]
implicit none
BEGIN_DOC
! Non-Local component of the pseudo-potential
END_DOC
integer :: i,k,j,l
integer :: kk,m
double precision :: r, rn, r2, alpha
double precision :: tmp(0:pseudo_lmax)
v_pseudo_non_local = 0.d0
do j=1,elec_num
kk = 0
do i=1,nucl_num
r = nucl_elec_dist(i,j)
r2 = r*r
tmp(0:pseudo_lmax) = 0.d0
do k=1,pseudo_kmax
do l=0,pseudo_lmax
alpha = pseudo_dz_kl(i,k,l)*r2
if (alpha > 20.d0) then
cycle
endif
select case (pseudo_n_kl(i,k,l))
case (0)
rn = 1.d0
case (-1)
rn = nucl_elec_dist_inv(i,j)
case (1)
rn = r
case (2)
rn = r*r
case (-2)
rn = nucl_elec_dist_inv(i,j)*nucl_elec_dist_inv(i,j)
case default
rn = r**pseudo_n_kl(i,k,l)
end select
tmp(l) += pseudo_v_kl(i,k,l) * exp(-alpha) * rn
enddo
enddo
do l=0,pseudo_lmax
do m=-l,l
kk += 1
v_pseudo_non_local(kk,j) = tmp(l)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, pseudo_mo_term, (mo_num,elec_num) ]
implicit none
BEGIN_DOC
! \sum < Ylm | MO > x Ylm(r) x V_nl(r)
END_DOC
integer :: ii,i,j,l,m,k,n,kk
double precision :: r, dr_inv
double precision :: tmp(pseudo_non_loc_dim_8,mo_num)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: tmp
integer, save :: ifirst = 0
if (ifirst == 0) then
pseudo_mo_term = 0.d0
ifirst = 1
endif
PROVIDE pseudo_ylm present_mos mo_pseudo_grid_scaled pseudo_non_loc_dim_count
dr_inv = dble(pseudo_grid_size)/pseudo_grid_rmax
do j=1,elec_num
tmp = 0.d0
kk=0
do k=1,nucl_num
r = nucl_elec_dist(k,j)
n = 1 + int(0.5+r*dr_inv)
if (n<=pseudo_grid_size) then
do ii=1,num_present_mos
i = present_mos(ii)
!DIR$ LOOP COUNT(4)
do l=kk+1,kk+pseudo_non_loc_dim_count(k)
tmp(l,i) = mo_pseudo_grid_scaled (l,i,n) * pseudo_ylm(l,j)
enddo
enddo
endif
kk = kk+pseudo_non_loc_dim_count(k)
enddo
do ii=1,num_present_mos
i = present_mos(ii)
pseudo_mo_term(i,j) = 0.d0
do kk=1,pseudo_non_loc_dim
pseudo_mo_term(i,j) = pseudo_mo_term(i,j) + tmp(kk,i) * v_pseudo_non_local(kk,j)
enddo
enddo
enddo
END_PROVIDER
real function ylm(l,m,x,y,z,r_inv)
implicit none
include 'constants.F'
BEGIN_DOC
! Y(l,m) evaluated at x,y,z
END_DOC
integer, intent(in) :: l,m
real, intent(in) :: x,y,z,r_inv
ylm = 0.5 * one_over_sqpi
select case(l)
case(0)
continue
case(1)
select case (m)
case (-1)
ylm = ylm * sq3 * y * r_inv
case (0)
ylm = ylm * sq3 * z * r_inv
case (1)
ylm = ylm * sq3 * x * r_inv
end select
! case(2)
! select case (m)
! case(-2)
! ylm = ylm * sqrt(15.) * x * y * r_inv * r_inv
! case(-1)
! ylm = ylm * sqrt(15.) * y * z * r_inv * r_inv
! case(0)
! ylm = 0.5 * ylm * sqrt(15.) * (2.*z*z - x*x - y*y) * r_inv * r_inv
! case(1)
! ylm = ylm * sqrt(15.) * z * x * r_inv * r_inv
! case(2)
! ylm = 0.5 * ylm * sqrt(15.) * (x*x - y*y) * r_inv * r_inv
! end select
case default
stop 'problem in Ylm of pseudo'
end select
end
BEGIN_PROVIDER [ double precision, pseudo_ylm, (pseudo_non_loc_dim_8,elec_num) ]
implicit none
BEGIN_DOC
! Y(l,m) evaluated for every electron position centered on every nuclei
END_DOC
integer :: i,j,l,m,kk
real, external :: ylm
integer, save :: ifirst = 0
if (ifirst == 0) then
pseudo_ylm = 0.d0
ifirst = 1
endif
do j=1,elec_num
kk = 0
do i=1,nucl_num
do l=0,pseudo_lmax
do m=-l,l
kk = kk+1
pseudo_ylm(kk,j) = ylm(l,m, &
nucl_elec_dist_vec(1,i,j), &
nucl_elec_dist_vec(2,i,j), &
nucl_elec_dist_vec(3,i,j), &
nucl_elec_dist_inv(i,j))
enddo
enddo
enddo
enddo
END_PROVIDER