outcdn.f90 6.46 KB
 Markus Betzinger committed Apr 26, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 `````` MODULE m_outcdn use m_constants USE m_types ! ******************************************************** ! calculates the charge density at given point p(i=1,3) ! ******************************************************** CONTAINS SUBROUTINE outcdn(& & p,n,na,iv,iflag,jsp,sliceplot,stars,& & vacuum,sphhar,atoms,sym,cell,oneD,& & qpw,rhtxy,rho,rht,& & xdnout) ! `````` 14 `````` use m_constants `````` Markus Betzinger committed Apr 26, 2016 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 `````` USE m_angle USE m_starf, ONLY : starf2,starf3 USE m_ylm IMPLICIT NONE ! TYPE(t_sliceplot),INTENT(IN) :: sliceplot 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 ! .. Scalar Arguments .. INTEGER, INTENT (IN) :: iflag,jsp,n,na,iv REAL, INTENT (OUT) :: xdnout !-odim !+odim ! .. ! .. Array Arguments .. `````` Daniel Wortmann committed Nov 06, 2018 37 38 39 40 `````` COMPLEX, INTENT (IN) :: qpw(:,:) !(stars%ng3,input%jspins) COMPLEX, INTENT (IN) :: rhtxy(:,:,:,:) !(vacuum%nmzxyd,oneD%odi%n2d-1,2,input%jspins) REAL, INTENT (IN) :: rho(:,0:,:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins) REAL, INTENT (IN) :: rht(:,:,:) !(vacuum%nmzd,2,input%jspins) `````` Markus Betzinger committed Apr 26, 2016 41 42 43 44 45 46 47 `````` 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 ! .. ! .. Local Arrays .. `````` Daniel Wortmann committed Feb 13, 2017 48 `````` COMPLEX sf2(stars%ng2),sf3(stars%ng3),ylm((atoms%lmaxd+1)**2) `````` Markus Betzinger committed Apr 26, 2016 49 `````` REAL rcc(3),x(3) `````` Miriam Hinzen committed Oct 02, 2018 50 `````` `````` 51 `````` `````` Markus Betzinger committed Apr 26, 2016 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 `````` 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) ! xdnout=dot_product(real(qpw(:,jsp)*sf3(:)),stars%nstr) 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 xdnout = rht(jp3,ivac,jsp) + delta*& & (rht(jp3+1,ivac,jsp)-rht(jp3,ivac,jsp)) 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) xx1 = xx1 + real(rhtxy(jp3,k-1,ivac,jsp)*& `````` 86 `````` & exp(ImagUnit*m*phi)*exp(ImagUnit*gzi*cell%bmat(3,3)*p(3)))*& `````` Markus Betzinger committed Apr 26, 2016 87 88 `````` & oneD%odi%nst2(k) xx2 = xx2 + real(rhtxy(jp3+1,k-1,ivac,jsp)*& `````` 89 `````` & exp(ImagUnit*m*phi)*exp(ImagUnit*gzi*cell%bmat(3,3)*p(3)))*& `````` Markus Betzinger committed Apr 26, 2016 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 132 133 134 135 136 137 138 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 197 198 199 200 201 202 203 204 205 206 207 `````` & 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 xdnout = rht(jp3,ivac,jsp) + delta*& & (rht(jp3+1,ivac,jsp)-rht(jp3,ivac,jsp)) IF (jp3.LT.vacuum%nmzxy) THEN xx1 = 0. xx2 = 0. DO k = 2,stars%ng2 xx1 = xx1 + real(rhtxy(jp3,k-1,ivac,jsp)*sf2(k))*stars%nstr2(k) xx2 = xx2 + real(rhtxy(jp3+1,k-1,ivac,jsp)*sf2(k))*& & 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 IF (sliceplot%plpot) THEN xd1 = xd1 + rho(jr,lh,n,jsp)*s ELSE xd1 = xd1 + rho(jr,lh,n,jsp)*s/ (atoms%rmsh(jr,n)*atoms%rmsh(jr,n)) END IF IF (jr.EQ.atoms%jri(n)) CYCLE IF (sliceplot%plpot) THEN xd2 = xd2 + rho(jr+1,lh,n,jsp)*s ELSE xd2 = xd2 + rho(jr+1,lh,n,jsp)*s/& & (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``````