outcdn.f90 6.81 KB
Newer Older
 Markus Betzinger committed Apr 26, 2016 1 2 3 4 5 6 7 `````` MODULE m_outcdn use m_constants USE m_types ! ******************************************************** ! calculates the charge density at given point p(i=1,3) ! ******************************************************** CONTAINS `````` Robin Hilgers committed Sep 26, 2019 8 9 10 11 12 13 14 `````` !SUBROUTINE outcdn(& !& p,n,na,iv,iflag,jsp,l_potential,stars,& !& vacuum,sphhar,atoms,sym,cell,oneD,& !& qpw,rhtxy,rho,rht,& !& xdnout) SUBROUTINE outcdn(& & p,n,na,iv,iflag,jsp,l_potential,stars,& `````` Markus Betzinger committed Apr 26, 2016 15 `````` & vacuum,sphhar,atoms,sym,cell,oneD,& `````` Robin Hilgers committed Sep 26, 2019 16 `````` & potDen,& `````` Markus Betzinger committed Apr 26, 2016 17 18 `````` & xdnout) ! `````` 19 `````` use m_constants `````` Markus Betzinger committed Apr 26, 2016 20 21 22 23 24 `````` USE m_angle USE m_starf, ONLY : starf2,starf3 USE m_ylm IMPLICIT NONE ! `````` Robin Hilgers committed Sep 26, 2019 25 ``````! TYPE(t_sliceplot),INTENT(IN) :: sliceplot TODO:Remove `````` Markus Betzinger committed Apr 26, 2016 26 27 28 29 30 31 32 `````` TYPE(t_stars),INTENT(IN) :: stars TYPE(t_vacuum),INTENT(IN) :: vacuum TYPE(t_sphhar),INTENT(IN) :: sphhar TYPE(t_atoms),INTENT(IN) :: atoms TYPE(t_sym),INTENT(IN) :: sym TYPE(t_cell),INTENT(IN) :: cell TYPE(t_oneD),INTENT(IN) :: oneD `````` Robin Hilgers committed Sep 26, 2019 33 `````` TYPE(t_potden),INTENT(IN) :: potDen `````` Markus Betzinger committed Apr 26, 2016 34 35 36 37 38 39 40 41 42 `````` ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: iflag,jsp,n,na,iv REAL, INTENT (OUT) :: xdnout !-odim !+odim ! .. ! .. Array Arguments .. `````` Robin Hilgers committed Sep 26, 2019 43 `````` `````` Markus Betzinger committed Apr 26, 2016 44 45 46 47 48 `````` REAL, INTENT (INOUT) :: p(3) ! .. ! .. Local Scalars .. REAL delta,s,sx,xd1,xd2,xx1,xx2,rrr,phi INTEGER i,j,jp3,jr,k,lh,mem,nd,nopa,ivac,ll1,lm ,gzi,m `````` Robin Hilgers committed Sep 26, 2019 49 50 51 52 `````` ! .. Logical Argumens .. LOGICAL, INTENT (IN) :: l_potential `````` Markus Betzinger committed Apr 26, 2016 53 54 ``````! .. ! .. Local Arrays .. `````` Daniel Wortmann committed Feb 13, 2017 55 `````` COMPLEX sf2(stars%ng2),sf3(stars%ng3),ylm((atoms%lmaxd+1)**2) `````` Markus Betzinger committed Apr 26, 2016 56 `````` REAL rcc(3),x(3) `````` Robin Hilgers committed Sep 26, 2019 57 58 `````` !Temp. Variable Assignment TODO: Can be removed as soon as the potDen datatype has been implemented in the whole file. `````` Robin Hilgers committed Sep 26, 2019 59 60 61 62 `````` !qpw=potDen%pw !rhtxz=potDen%vacxy !rho=potDen%mt !rht=potDen%vacz `````` Robin Hilgers committed Sep 26, 2019 63 `````` `````` 64 `````` `````` Markus Betzinger committed Apr 26, 2016 65 66 67 68 69 70 71 72 73 74 75 `````` ivac=iv if (iflag.ne.1) THEN if (iflag.ne.0) THEN ! ---> interstitial part !CALL cotra1(p(1),rcc,cell%bmat) rcc=matmul(cell%bmat,p)/tpi_const CALL starf3(& & sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,& & sf3) ! `````` Robin Hilgers committed Sep 26, 2019 76 `````` xdnout=dot_product(real(potDen%pw(:,jsp)*sf3(:)),stars%nstr) `````` Markus Betzinger committed Apr 26, 2016 77 78 79 80 81 82 83 84 85 86 87 88 89 `````` RETURN ! ---> vacuum part ENDIF xdnout = 0. !-odim IF (oneD%odi%d1) THEN rrr = sqrt( p(1)**2 + p(2)**2 ) phi = angle(p(1),p(2)) jp3 = (rrr-cell%z1)/vacuum%delz delta = (rrr-cell%z1)/vacuum%delz - jp3 !* we count 0 as point 1 jp3 = jp3 + 1 IF (jp3.LT.vacuum%nmz) THEN `````` Robin Hilgers committed Sep 26, 2019 90 91 `````` xdnout = potDen%vacz(jp3,ivac,jsp) + delta*& & (potDen%vacz(jp3+1,ivac,jsp)-potDen%vacz(jp3,ivac,jsp)) `````` Markus Betzinger committed Apr 26, 2016 92 93 94 95 96 97 `````` IF (jp3.LT.vacuum%nmzxy) THEN xx1 = 0. xx2 = 0. DO k = 2,oneD%odi%nq2 m = oneD%odi%kv(2,k) gzi = oneD%odi%kv(1,k) `````` Robin Hilgers committed Sep 26, 2019 98 `````` xx1 = xx1 + real(potDen%vacxy(jp3,k-1,ivac,jsp)*& `````` 99 `````` & exp(ImagUnit*m*phi)*exp(ImagUnit*gzi*cell%bmat(3,3)*p(3)))*& `````` Markus Betzinger committed Apr 26, 2016 100 `````` & oneD%odi%nst2(k) `````` Robin Hilgers committed Sep 26, 2019 101 `````` xx2 = xx2 + real(potDen%vacxy(jp3+1,k-1,ivac,jsp)*& `````` 102 `````` & exp(ImagUnit*m*phi)*exp(ImagUnit*gzi*cell%bmat(3,3)*p(3)))*& `````` Markus Betzinger committed Apr 26, 2016 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 `````` & oneD%odi%nst2(k) ENDDO xdnout = xdnout + xx1 + delta* (xx2-xx1) END IF ELSE xdnout = 0.0 END IF ELSE !+odim IF (p(3).LT.0.0) THEN ivac = vacuum%nvac IF (sym%invs) THEN p(1:2) = -p(1:2) END IF p(3) = abs(p(3)) END IF !CALL cotra1(p,rcc,cell%bmat) rcc=matmul(cell%bmat,p)/tpi_const CALL starf2(& & sym%nop2,stars%ng2,stars%kv3,sym%mrot,sym%symor,sym%tau,rcc,sym%invtab,& & sf2) ! jp3 = (p(3)-cell%z1)/vacuum%delz delta = (p(3)-cell%z1)/vacuum%delz - jp3 !* we count 0 as point 1 jp3 = jp3 + 1 IF (jp3.LT.vacuum%nmz) THEN `````` Robin Hilgers committed Sep 26, 2019 131 132 `````` xdnout = potDen%vacz(jp3,ivac,jsp) + delta*& & (potDen%vacz(jp3+1,ivac,jsp)-potDen%vacz(jp3,ivac,jsp)) `````` Markus Betzinger committed Apr 26, 2016 133 134 135 136 `````` IF (jp3.LT.vacuum%nmzxy) THEN xx1 = 0. xx2 = 0. DO k = 2,stars%ng2 `````` Robin Hilgers committed Sep 26, 2019 137 138 `````` xx1 = xx1 + real(potDen%vacxy(jp3,k-1,ivac,jsp)*sf2(k))*stars%nstr2(k) xx2 = xx2 + real(potDen%vacxy(jp3+1,k-1,ivac,jsp)*sf2(k))*& `````` Markus Betzinger committed Apr 26, 2016 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 `````` & stars%nstr2(k) enddo xdnout = xdnout + xx1 + delta* (xx2-xx1) END IF ELSE xdnout = 0.0 END IF !----> vacuum finishes ENDIF RETURN ENDIF ! ----> m.t. part nd = atoms%ntypsy(na) nopa = atoms%ngopr(na) IF (oneD%odi%d1) nopa = oneD%ods%ngopr(na) sx = 0.0 DO i = 1,3 x(i) = p(i) - atoms%pos(i,na) sx = sx + x(i)*x(i) enddo sx = sqrt(sx) IF (nopa.NE.1) THEN !... switch to internal units !CALL cotra1(x,rcc,cell%bmat) rcc=matmul(cell%bmat,x)/tpi_const !... rotate into representative DO i = 1,3 p(i) = 0. DO j = 1,3 IF (.NOT.oneD%odi%d1) THEN p(i) = p(i) + sym%mrot(i,j,nopa)*rcc(j) ELSE p(i) = p(i) + oneD%ods%mrot(i,j,nopa)*rcc(j) END IF enddo enddo !... switch back to cartesian units !CALL cotra0(p,x,cell%amat) x=matmul(cell%amat,p) END IF DO j = atoms%jri(n),2,-1 IF (sx.GE.atoms%rmsh(j,n)) EXIT ENDDO jr = j CALL ylm4(& & atoms%lmax(n),x,& & ylm) xd1 = 0.0 xd2 = 0.0 DO lh = 0, sphhar%nlh(nd) ll1 = sphhar%llh(lh,nd) * ( sphhar%llh(lh,nd) + 1 ) + 1 s = 0.0 DO mem = 1,sphhar%nmem(lh,nd) lm = ll1 + sphhar%mlh(mem,lh,nd) s = s + real( sphhar%clnu(mem,lh,nd)*ylm(lm) ) ENDDO `````` Robin Hilgers committed Sep 26, 2019 197 `````` IF (l_potential) THEN `````` Robin Hilgers committed Sep 26, 2019 198 `````` xd1 = xd1 + potDen%mt(jr,lh,n,jsp)*s `````` Markus Betzinger committed Apr 26, 2016 199 `````` ELSE `````` Robin Hilgers committed Sep 26, 2019 200 `````` xd1 = xd1 + potDen%mt(jr,lh,n,jsp)*s/ (atoms%rmsh(jr,n)*atoms%rmsh(jr,n)) `````` Markus Betzinger committed Apr 26, 2016 201 202 `````` END IF IF (jr.EQ.atoms%jri(n)) CYCLE `````` Robin Hilgers committed Sep 26, 2019 203 `````` IF (l_potential) THEN `````` Robin Hilgers committed Sep 26, 2019 204 `````` xd2 = xd2 + potDen%mt(jr+1,lh,n,jsp)*s `````` Markus Betzinger committed Apr 26, 2016 205 `````` ELSE `````` Robin Hilgers committed Sep 26, 2019 206 `````` xd2 = xd2 + potDen%mt(jr+1,lh,n,jsp)*s/& `````` Markus Betzinger committed Apr 26, 2016 207 208 209 210 211 212 213 214 215 216 217 218 219 220 `````` & (atoms%rmsh(jr+1,n)*atoms%rmsh(jr+1,n)) END IF ENDDO IF (jr.EQ.atoms%jri(n)) THEN xdnout = xd1 ELSE xdnout = xd1 + (xd2-xd1) *& & (sx-atoms%rmsh(jr,n)) / (atoms%rmsh(jr+1,n)-atoms%rmsh(jr,n)) END IF 8000 FORMAT (2f10.6) ! RETURN END SUBROUTINE outcdn END MODULE m_outcdn``````