Line data Source code
1 : MODULE m_nabla
2 : use m_juDFT
3 :
4 : CONTAINS
5 :
6 0 : SUBROUTINE nabla(
7 : > ispecies,number_of_j1,grid_size,delta_x,
8 0 : > nstd,ntypd,j1,l1,lmax,ms,ri,psi,phi,dphi,
9 0 : < psi_phi)
10 : !----------------------------------------------------------------
11 : !
12 : ! ispecies ... number of the atom (itype)
13 : ! number_of_j1 ... number of core wavefunction
14 : ! grid_size ... rumber of radial gridpoints (jri)
15 : ! delta_x ... logarithmic increment of grid (dx)
16 : ! ri(grid_size)... radial mesh
17 : ! nstd ... number of corelevels (dimension)
18 : ! ntypd ... number of atoms (dimension)
19 : ! j1,l1 ... quantum numbers of the core wavefunction
20 : ! ms ... -1/2 or +1/2 for spin 1 or 2
21 : ! lmax ... parameter, 3 (s,p,d,f)
22 : ! psi(r) ... core wavefunction
23 : ! phi(r,l) ... valence wavefunction
24 : ! dphi(r,l) ... radial derivative of valence wavefunction
25 : !
26 : !----------------------------------------------------------------
27 : USE m_constants
28 : USE m_clebsch
29 : USE m_intgr, ONLY : intgr3
30 :
31 : IMPLICIT NONE
32 :
33 : INTEGER, INTENT(IN) :: ispecies, number_of_j1, grid_size
34 : INTEGER, INTENT(IN) :: l1, lmax, nstd, ntypd
35 : REAL, INTENT(IN) :: delta_x, j1, ms
36 : REAL, INTENT(IN) :: psi(grid_size), ri(grid_size)
37 : REAL, INTENT(IN) :: phi(grid_size,0:lmax)
38 : REAL, INTENT(IN) ::dphi(grid_size,0:lmax)
39 : COMPLEX, INTENT(INOUT):: psi_phi(:,:,:)!nstd,(lmax+1)**2,3*ntypd)
40 :
41 : INTEGER :: m1, l2, m2, index, alloc_error, lmn1, lmn2
42 : REAL :: result, result1, total_result, spin
43 : REAL :: mu, cnst_one_over_sqrt_two, cnst_zero
44 : COMPLEX :: cnst_i
45 : REAL, DIMENSION(:), POINTER :: f
46 0 : spin = 0.50
47 0 : cnst_one_over_sqrt_two = 1.0/sqrt(2.0)
48 0 : cnst_i = cmplx(0.0,1.0)
49 0 : cnst_zero = 0.0
50 :
51 0 : NULLIFY(f)
52 :
53 : IF ( ASSOCIATED(f) ) THEN
54 : WRITE(oUnit,*)'nabla: f association status:',ASSOCIATED(f)
55 : STOP
56 : ENDIF
57 :
58 0 : ALLOCATE ( f(grid_size),STAT=alloc_error )
59 0 : IF (alloc_error /= 0) CALL juDFT_error("Couldn't allocate f",
60 0 : + calledby ="nabla")
61 :
62 0 : lmn1 = 2 * (number_of_j1 - 1) * l1
63 0 : mu = -j1
64 :
65 0 : DO WHILE (mu <= j1)
66 0 : lmn1 = lmn1 + 1
67 0 : m1 = INT(mu - ms)
68 0 : lmn2 = 0
69 0 : DO l2 = 0, lmax
70 0 : DO m2 = -l2, l2
71 0 : lmn2 = lmn2 + 1
72 0 : IF(l1 == l2 + 1)THEN
73 0 : total_result = 0.00
74 : result = 0.00
75 : result1 = 0.00
76 : !
77 : ! (l+1)/srqt[(2l+1)(2l+3)] < phi_core | (d phi_valence / dr ) >
78 : !
79 0 : f(:) = psi(:) * dphi(:,l2) ! assumed to be already multiplied with * ri(:) * ri(:)
80 :
81 0 : CALL intgr3(f,ri,delta_x,grid_size,result)
82 :
83 : result = result * (l2 + 1.00) /
84 0 : + sqrt((2.00*l2 +1.00)*(2.*l2+3.00))
85 :
86 : !
87 : ! - l(l+1)/srqt[(2l+1)(2l+3)] < phi_core | (1/r) phi_valence >
88 : !
89 0 : f(:) = psi(:) * phi(:,l2) ! assumed to be already multiplied with 1G* ri(:)
90 0 : CALL intgr3(f,ri,delta_x,grid_size,result1)
91 : result1 = - result1 * l2 * (l2 + 1.0) /
92 0 : + sqrt((2.0 * l2 + 1.00) * (2. * l2 + 3.00))
93 : !
94 : ! Sum up and decorate with Clebsch-Gordon coefficients
95 : !
96 0 : result = result + result1
97 0 : total_result = result*clebsch(real(l1),spin,mu-ms,ms,j1,mu)
98 :
99 0 : index = (ispecies - 1) * 3 + 1
100 : psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,1,m1) * ! left polarization
101 0 : + total_result / cgc(l2,1,l1,0,0,0)
102 :
103 0 : index = index + 1
104 : psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,-1,m1) * ! right polarization
105 0 : + total_result / cgc(l2,1,l1,0,0,0)
106 :
107 0 : index = index + 1
108 : psi_phi(lmn1,lmn2,index)= cgc(l2,1,l1,m2,0,m1) * ! z-polarization
109 0 : + total_result / cgc(l2,1,l1,0,0,0)
110 :
111 0 : ELSEIF(l1== l2-1)THEN
112 :
113 : result = cnst_zero
114 : result1 = cnst_zero
115 : !
116 : ! l/srqt[(2l-1)(2l+1)] < phi_core | (d phi_valence / dr ) >
117 : !
118 0 : f(:) = psi(:)* dphi(:,l2) * ri(:) * ri(:)
119 0 : CALL intgr3(f,ri,delta_x,grid_size,result)
120 0 : result = result * l2 / sqrt((2.0*l2 - 1.0)*(2.0*l2 + 1.0))
121 : !
122 : ! l(l+1)/srqt[(2l-1)(2l+1)] < phi_core | (1/r) phi_valence >
123 : !
124 0 : f(:) = psi(:)* phi(:,l2) * ri(:)
125 0 : CALL intgr3(f,ri,delta_x,grid_size,result1)
126 : result1 = result1 * l2 * (l2 + 1.0) /
127 0 : + sqrt((2.00 * l2 - 1.00) * (2.00 * l2 + 1.00))
128 : !
129 : ! Sum up and decorate with Clebsch-Gordon coefficients
130 : !
131 0 : result = result + result1
132 0 : total_result= result*clebsch(real(l1),spin,mu-ms,ms,j1,mu)
133 :
134 :
135 0 : index = (ispecies - 1) * 3 + 1
136 : psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,1,m1) * ! left polarization
137 0 : + total_result / cgc(l2,1,l1,0,0,0)
138 0 : index = index + 1
139 : psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,-1,m1) * ! right polarization
140 0 : + total_result / cgc(l2,1,l1,0,0,0)
141 0 : index = index + 1
142 : psi_phi(lmn1,lmn2,index) = cgc(l2,1,l1,m2,0,m1) * ! z-polarization
143 0 : + total_result / cgc(l2,1,l1,0,0,0)
144 : ENDIF
145 : ENDDO
146 : ENDDO
147 0 : mu = mu + 1.00
148 : ENDDO
149 0 : DEALLOCATE(f,STAT=alloc_error)
150 : IF (alloc_error /= 0) CALL juDFT_error("Couldn't deallocate f"
151 0 : + ,calledby ="nabla")
152 0 : END SUBROUTINE nabla
153 :
154 0 : FUNCTION cgc(l1,l2,l3,m1,m2,m3)
155 :
156 : IMPLICIT NONE
157 : INTEGER :: l1, l2, l3, m1, m2, m3
158 : REAL :: two_l1p1, two_l1p2, l1pm3, l1pm3p1, l1mm3p1, l1mm3, cgc
159 :
160 0 : IF (m3 /= m1 + m2) THEN
161 0 : cgc = 0.0
162 : RETURN
163 : END IF
164 : ! gb m3 = m1 + m2
165 0 : two_l1p1 = 2 * l1 + 1
166 0 : two_l1p2 = 2 * l1 + 2
167 0 : l1pm3 = l1 + m3
168 0 : l1pm3p1 = l1 + m3 + 1
169 0 : l1mm3p1 = l1 - m3 + 1
170 0 : l1mm3 = l1 - m3
171 0 : cgc = 0.0
172 0 : IF (l3 == l1 + 1) THEN
173 0 : IF (m2 == 1) then
174 0 : cgc = sqrt( (l1pm3 * l1pm3p1) / (two_l1p1 * two_l1p2))
175 0 : ELSEIF (m2 == 0) THEN
176 0 : cgc = sqrt( (l1mm3p1 * l1pm3p1) / (two_l1p1 * (l1 + 1)))
177 0 : ELSEIF (m2 == -1) THEN
178 0 : cgc = sqrt( (l1mm3 * l1mm3p1) / (two_l1p1 * two_l1p2))
179 : END IF
180 0 : ELSE IF(l3 == l1 -1) THEN
181 0 : IF (m2 == 1) then
182 0 : cgc = sqrt( (l1mm3 * l1mm3p1) / (2.0 * l1 * two_l1p1))
183 0 : ELSEIF (m2 == 0) THEN
184 0 : cgc = -sqrt( (l1mm3 * l1pm3) / (l1 * (two_l1p1)))
185 0 : ELSEIF (m2 == -1) THEN
186 0 : cgc = sqrt( (l1pm3p1 * l1pm3) / (2.0 * l1 * two_l1p1))
187 : END IF
188 : END IF
189 : END FUNCTION cgc
190 :
191 :
192 : END MODULE m_nabla
|
