-
Notifications
You must be signed in to change notification settings - Fork 2
/
pot_ao_pseudo_ints.irp.f
131 lines (119 loc) · 3.95 KB
/
pot_ao_pseudo_ints.irp.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
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
! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
!
! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
END_DOC
! l,m : Y(l,m) parameters
! c(3) : pseudopotential center
! a(3) : Atomic Orbital center
! n_a(3) : Powers of x,y,z in the Atomic Orbital
! g_a : Atomic Orbital exponent
! r : Distance between the Atomic Orbital center and the considered point
double precision, external :: ylm_orb
integer :: n_a(3)
double precision :: a(3), c(3), g_a
integer :: i,j,k,l,m,n,p
double precision :: dr, Ulc
double precision :: y
double precision, allocatable :: r(:)
allocate (r(pseudo_grid_size))
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
do j=2,pseudo_grid_size
r(j) = r(j-1) + dr
enddo
ao_pseudo_grid = 0.d0
do j=1,pseudo_grid_size
do k=1,nucl_num
c(1:3) = nucl_coord(k,1:3)
do l=0,pseudo_lmax
do i=1,ao_num
a(1:3) = nucl_coord(ao_nucl(i),1:3)
n_a(1:3) = ao_power(i,1:3)
do p=1,ao_prim_num(i)
g_a = ao_expo_ordered_transp(p,i)
do m=-l,l
y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
ao_pseudo_grid(i,m,l,k,j) = ao_pseudo_grid(i,m,l,k,j) + &
ao_coef_normalized_ordered_transp(p,i)*y
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(r)
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
! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \phi_i^{A} (r-r_A) d\Omega_C
!
! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
END_DOC
! l,m : Y(l,m) parameters
! c(3) : pseudopotential center
! a(3) : Atomic Orbital center
! n_a(3) : Powers of x,y,z in the Atomic Orbital
! g_a : Atomic Orbital exponent
! r : Distance between the Atomic Orbital center and the considered point
double precision, external :: ylm_orb
integer :: n_a(3)
double precision :: a(3), c(3), g_a
integer :: i,j,k,l,m,n,p
double precision :: dr, Ulc
double precision :: y
double precision, allocatable :: r(:)
allocate (r(pseudo_grid_size))
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
r(1) = 0.d0
do j=2,pseudo_grid_size
r(j) = r(j-1) + dr
enddo
mo_pseudo_grid = 0.d0
do n=1,pseudo_grid_size
do k=1,nucl_num
do l=0,pseudo_lmax
do m=-l,l
do i=1,ao_num
do j=1,mo_num
if (dabs(ao_pseudo_grid(i,m,l,k,n)) < 1.e-12) then
cycle
endif
if (dabs(mo_coef(i,j)) < 1.e-8) then
cycle
endif
mo_pseudo_grid(j,m,l,k,n) = mo_pseudo_grid(j,m,l,k,n) + &
ao_pseudo_grid(i,m,l,k,n) * mo_coef(i,j)
enddo
enddo
enddo
enddo
enddo
enddo
deallocate(r)
END_PROVIDER
double precision function test_pseudo_grid_ao(i,j)
implicit none
integer, intent(in) :: i,j
integer :: k,l,m,n
double precision :: r, dr,u
dr = pseudo_grid_rmax/dble(pseudo_grid_size)
test_pseudo_grid_ao = 0.d0
r = 0.d0
do k=1,pseudo_grid_size
do n=1,nucl_num
do l = 0,pseudo_lmax
u = pseudo_v_kl(n,l,1) * exp(-pseudo_dz_kl(n,l,1)*r*r)* r*r*dr
do m=-l,l
test_pseudo_grid_ao += ao_pseudo_grid(i,m,l,n,k) * ao_pseudo_grid(j,m,l,n,k) * u
enddo
enddo
enddo
r = r+dr
enddo
end