Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_outcdn
7 :
8 : CONTAINS
9 :
10 0 : SUBROUTINE outcdn(p, n, na, iv, iflag, jsp, l_potential, stars, vacuum, &
11 : sphhar, atoms, sym, cell, potDen, xdnout)
12 : USE m_types
13 : USE m_constants
14 : USE m_angle
15 : USE m_starf, ONLY : starf2,starf3
16 : USE m_ylm
17 :
18 : !--------------------------------------------------------------------------
19 : ! Calculates the charge density at a given point p(i=1,3).
20 : !--------------------------------------------------------------------------
21 :
22 : IMPLICIT NONE
23 :
24 : TYPE(t_stars),INTENT(IN) :: stars
25 : TYPE(t_vacuum),INTENT(IN) :: vacuum
26 : TYPE(t_sphhar),INTENT(IN) :: sphhar
27 : TYPE(t_atoms),INTENT(IN) :: atoms
28 : TYPE(t_sym),INTENT(IN) :: sym
29 : TYPE(t_cell),INTENT(IN) :: cell
30 :
31 : TYPE(t_potden),INTENT(IN) :: potDen
32 :
33 :
34 : ! Scalar Arguments
35 : INTEGER, INTENT (IN) :: iflag,jsp,n,na,iv
36 : REAL, INTENT (OUT) :: xdnout
37 : ! -odim
38 : ! +odim
39 :
40 : ! Array Arguments
41 : REAL, INTENT (INOUT) :: p(3)
42 :
43 : ! Logical argument
44 : LOGICAL, INTENT (IN) :: l_potential
45 :
46 : ! Local scalars
47 : REAL delta,s,sx,xd1,xd2,xx1,xx2,rrr,phi
48 : INTEGER i,j,jp3,jr,k,lh,mem,nd,nopa,ivac,ll1,lm ,gzi,m
49 :
50 : ! Local arrays
51 0 : COMPLEX, ALLOCATABLE :: sf2(:),sf3(:),ylm(:)
52 : REAL rcc(3),x(3)
53 :
54 0 : ALLOCATE( sf2(stars%ng2),sf3(stars%ng3),ylm((atoms%lmaxd+1)**2))
55 0 : ivac=iv
56 :
57 0 : IF (iflag.NE.1) THEN
58 0 : IF (iflag.NE.0) THEN
59 : ! Interstitial part:
60 0 : rcc=matmul(cell%bmat,p)/tpi_const
61 : CALL starf3(sym%nop, stars%ng3, sym%symor, stars%kv3, sym%mrot, &
62 0 : sym%tau, rcc, sym%invtab, sf3)
63 :
64 0 : xdnout=dot_product(real(potDen%pw(:,jsp)*sf3(:)),stars%nstr)
65 0 : RETURN
66 :
67 : ENDIF
68 :
69 : ! Vacuum part:
70 0 : xdnout = 0.
71 :
72 :
73 :
74 0 : IF (p(3).LT.0.0) THEN
75 0 : ivac = vacuum%nvac
76 0 : IF (sym%invs) THEN
77 0 : p(1:2) = -p(1:2)
78 : END IF
79 0 : p(3) = abs(p(3))
80 : END IF
81 0 : rcc=matmul(cell%bmat,p)/tpi_const
82 : CALL starf2(sym%nop2, stars%ng2, stars%kv2, sym%mrot, sym%symor, &
83 0 : sym%tau,rcc,sym%invtab,sf2)
84 :
85 0 : jp3 = (p(3)-cell%z1)/vacuum%delz
86 0 : delta = (p(3)-cell%z1)/vacuum%delz - jp3
87 : ! We count 0 as point 1.
88 0 : jp3 = jp3 + 1
89 0 : IF (jp3.LT.vacuum%nmz) THEN
90 : xdnout = REAL(potDen%vac(jp3,1,ivac,jsp)) + &
91 : delta*(REAL(potDen%vac(jp3+1,1,ivac,jsp)) - &
92 0 : REAL(potDen%vac(jp3,1,ivac,jsp)))
93 0 : IF (jp3.LT.vacuum%nmzxy) THEN
94 : xx1 = 0.
95 : xx2 = 0.
96 0 : DO k = 2,stars%ng2
97 : xx1 = xx1 + REAL(potDen%vac(jp3,k,ivac,jsp)*sf2(k))* &
98 0 : stars%nstr2(k)
99 : xx2 = xx2 + REAL(potDen%vac(jp3+1,k,ivac,jsp)*sf2(k))* &
100 0 : stars%nstr2(k)
101 : ENDDO
102 0 : xdnout = xdnout + xx1 + delta* (xx2-xx1)
103 : END IF
104 : ELSE
105 0 : xdnout = 0.0
106 : END IF
107 : ! Vacuum part finished.
108 :
109 0 : RETURN
110 : ENDIF
111 : ! MT part:
112 :
113 0 : nd = sym%ntypsy(na)
114 0 : nopa = sym%ngopr(na)
115 0 : sx = 0.0
116 0 : DO i = 1,3
117 0 : x(i) = p(i) - atoms%pos(i,na)
118 0 : sx = sx + x(i)*x(i)
119 : END DO
120 0 : sx = sqrt(sx)
121 0 : IF (nopa.NE.1) THEN
122 : ! Switch to internal units.
123 0 : rcc=matmul(cell%bmat,x)/tpi_const
124 : ! Rotate into representative.
125 0 : DO i = 1,3
126 0 : p(i) = 0.
127 0 : DO j = 1,3
128 :
129 0 : p(i) = p(i) + sym%mrot(i,j,nopa)*rcc(j)
130 :
131 : END DO
132 : END DO
133 : ! Switch back to cartesian units.
134 0 : x=matmul(cell%amat,p)
135 : END IF
136 0 : DO j = atoms%jri(n),2,-1
137 0 : IF (sx.GE.atoms%rmsh(j,n)) EXIT
138 : ENDDO
139 0 : jr = j
140 0 : CALL ylm4(atoms%lmax(n),x,ylm)
141 0 : xd1 = 0.0
142 0 : xd2 = 0.0
143 0 : DO lh = 0, sphhar%nlh(nd)
144 0 : ll1 = sphhar%llh(lh,nd) * ( sphhar%llh(lh,nd) + 1 ) + 1
145 0 : s = 0.0
146 0 : DO mem = 1,sphhar%nmem(lh,nd)
147 0 : lm = ll1 + sphhar%mlh(mem,lh,nd)
148 0 : s = s + real( sphhar%clnu(mem,lh,nd)*ylm(lm) )
149 : ENDDO
150 0 : IF (l_potential) THEN
151 0 : xd1 = xd1 + potDen%mt(jr,lh,n,jsp)*s
152 : ELSE
153 0 : xd1 = xd1 + potDen%mt(jr,lh,n,jsp)*s/(atoms%rmsh(jr,n)**2)
154 : END IF
155 0 : IF (jr.EQ.atoms%jri(n)) CYCLE
156 0 : IF (l_potential) THEN
157 0 : xd2 = xd2 + potDen%mt(jr+1,lh,n,jsp)*s
158 : ELSE
159 0 : xd2 = xd2 + potDen%mt(jr+1,lh,n,jsp)*s/(atoms%rmsh(jr+1,n)**2)
160 : END IF
161 : ENDDO
162 0 : IF (jr.EQ.atoms%jri(n)) THEN
163 0 : xdnout = xd1
164 : ELSE
165 : xdnout = xd1 + (xd2-xd1) * (sx-atoms%rmsh(jr,n))/ &
166 0 : (atoms%rmsh(jr+1,n)-atoms%rmsh(jr,n))
167 : END IF
168 : 8000 FORMAT (2f10.6)
169 :
170 : RETURN
171 0 : END SUBROUTINE outcdn
172 : END MODULE m_outcdn
|
