Line data Source code
1 : MODULE m_potdis
2 : CONTAINS
3 0 : SUBROUTINE potdis(stars,vacuum,atoms,sphhar, input,cell,sym)
4 : !
5 : ! *****************************************************
6 : ! calculates the root-mean-square distance between input
7 : ! and output potentials (averaged over unit cell volume)
8 : ! based on code by c.l.fu
9 : ! *****************************************************
10 : !
11 : USE m_intgr, ONLY : intgr3, intgz0
12 : USE m_constants
13 : USE m_loddop
14 : USE m_cfft
15 : USE m_types
16 : IMPLICIT NONE
17 : TYPE(t_input),INTENT(IN) :: input
18 : TYPE(t_vacuum),INTENT(IN) :: vacuum
19 : TYPE(t_sym),INTENT(IN) :: sym
20 : TYPE(t_stars),INTENT(IN) :: stars
21 : TYPE(t_cell),INTENT(IN) :: cell
22 : TYPE(t_sphhar),INTENT(IN) :: sphhar
23 : TYPE(t_atoms),INTENT(IN) :: atoms
24 : ! .. Local Scalars ..
25 : REAL fact,rhs,sumis,sumz
26 : COMPLEX phase
27 : INTEGER i,i1,i2,i3,id2,id3,io,ip,iter,ivac,j,k1,k2,k3,lh,n,&
28 : nk12,npz,nt,num,na
29 : LOGICAL tail
30 : ! ..
31 : ! .. Local Arrays ..
32 0 : COMPLEX rhpw(stars%ng3,2,2),rhv1(vacuum%nmzxyd,stars%ng2-1,2,2,2),rhv2(vacuum%nmzd,stars%ng2,2,2,2)
33 0 : REAL rhsp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,2,2),rhv0(vacuum%nmzd,2,2,2)
34 0 : REAL dis(2),disz(vacuum%nmzd),rh(atoms%jmtd)
35 0 : REAL af3((2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1),input%jspins),&
36 0 : bf3((2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1),input%jspins), &
37 0 : bf2((2*stars%mx1+1)*(2*stars%mx2+1),input%jspins),&
38 0 : af2((2*stars%mx1+1)*(2*stars%mx2+1),input%jspins)
39 0 : REAL pdis(0:4,0:atoms%ntype,2)
40 : ! ..
41 : !
42 :
43 0 : dis = 0.
44 0 : pdis=0.0
45 :
46 0 : tail = .TRUE.
47 0 : npz = vacuum%nmz + 1
48 0 : nt = (2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1)
49 0 : nk12 = (2*stars%mx1+1)*(2*stars%mx2+1)
50 0 : fact = cell%omtil/REAL(nt)
51 : ! ---> reload potentials
52 0 : OPEN (9,file='nrp',form='unformatted',status='unknown')
53 0 : REWIND 9
54 0 : DO io = 1,2
55 : CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
56 0 : 9, iter,rhsp(:,0:,:,:,io),rhpw(:,:,io), rhv2(:,:,:,:,io))
57 : ENDDO
58 0 : CLOSE (9)
59 0 : IF (input%jspins.EQ.1) THEN
60 : ! ---> total potential difference
61 0 : CALL priv_cdndif(1,1,1,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
62 : ELSE
63 : ! ---> total input potential
64 0 : CALL priv_cdndif(1,1,2,1,1,1,1.0,1.0,input%film, rhsp,rhpw,rhv0,rhv1)
65 : ! ---> total output potential
66 0 : CALL priv_cdndif(1,1,2,2,2,2,1.0,1.0,input%film, rhsp,rhpw,rhv0,rhv1)
67 : ! ---> input 'spin' potential
68 0 : CALL priv_cdndif(2,1,2,1,1,1,1.0,-2.0,input%film, rhsp,rhpw,rhv0,rhv1)
69 : ! ---> output 'spin' potential
70 0 : CALL priv_cdndif(2,1,2,2,2,2,1.0,-2.0,input%film, rhsp,rhpw,rhv0,rhv1)
71 : ! ---> potential difference
72 0 : CALL priv_cdndif(1,1,1,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
73 : ! ---> spin potential difference
74 0 : CALL priv_cdndif(2,2,2,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
75 : END IF
76 0 : DO num = 1,input%jspins
77 : ! ----> m.t. part
78 0 : DO n = 1,atoms%ntype
79 0 : na = atoms%firstAtom(n)
80 0 : rh(:atoms%jri(n))=rhsp(:atoms%jri(n),0,n,num,1)*rhsp(:atoms%jri(n),0,n,num,1) *fpi_const
81 0 : CALL intgr3(rh,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),pdis(0,n,num))
82 0 : pdis(4,n,num)=pdis(0,n,num)
83 0 : DO lh = 1,sphhar%nlh(sym%ntypsy(na))
84 : rh(:atoms%jri(n))=rh(:atoms%jri(n))+rhsp(:atoms%jri(n),lh,n,num,1)*&
85 0 : rhsp(:atoms%jri(n),lh,n,num,1)* atoms%rmsh(:atoms%jri(n),n)*atoms%rmsh(:atoms%jri(n),n)
86 0 : CALL intgr3(rh,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),rhs)
87 0 : IF (lh<4) pdis(lh,n,num)=rhs
88 0 : pdis(4,n,num)=pdis(4,n,num)+rhs
89 : ENDDO
90 0 : pdis(:,n,num) = pdis(:,n,num)*atoms%neq(n)
91 0 : dis(num)=dis(num)+SUM(pdis(:,n,num))
92 0 : pdis(:,n,num) = SQRT(pdis(:,n,num)/atoms%volmts(n))*1000.
93 : ENDDO
94 : ! ----> interstitial part
95 : ! ---> create density in the real space
96 : i = 0
97 0 : DO i3 = 0,2*stars%mx3
98 0 : k3 = i3
99 0 : IF (k3.GT.stars%mx3) k3 = k3 - 2*stars%mx3-1
100 0 : DO i2 = 0,2*stars%mx2
101 0 : k2 = i2
102 0 : IF (k2.GT.stars%mx2) k2 = k2 - 2*stars%mx2-1
103 0 : DO i1 = 0,2*stars%mx1
104 0 : k1 = i1
105 0 : IF (k1.GT.stars%mx1) k1 = k1 - 2*stars%mx3-1
106 0 : i = i + 1
107 0 : id3 = stars%ig(k1,k2,k3)
108 0 : phase = stars%rgphs(k1,k2,k3)
109 0 : IF (id3.EQ.0) THEN
110 0 : af3(i,1) = 0.
111 0 : bf3(i,1) = 0.
112 : ELSE
113 : af3(i,1) = REAL(rhpw(id3,num,1))*REAL(phase) &
114 0 : - AIMAG(rhpw(id3,num,1))*AIMAG(phase)
115 : bf3(i,1) = AIMAG(rhpw(id3,num,1))*REAL(phase) &
116 0 : + REAL(rhpw(id3,num,1))*AIMAG(phase)
117 : END IF
118 : ENDDO
119 : ENDDO
120 : ENDDO
121 :
122 0 : CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx1*2+1,stars%mx1*2+1,1)
123 0 : CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx2*2+1,nk12,1)
124 0 : CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx3*2+1,nt,1)
125 : ! ---> form the dot product
126 0 : i = 0
127 0 : DO i3 = 0,stars%mx3*2
128 0 : DO i2 = 0,stars%mx2*2
129 0 : DO i1 = 0,stars%mx3*2
130 0 : i = i + 1
131 0 : af3(i,1) = af3(i,1)*af3(i,1)
132 0 : bf3(i,1) = 0.
133 : ENDDO
134 : ENDDO
135 : ENDDO
136 : ! ---> back fft
137 0 : CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx1*2+1,stars%mx1*2+1,-1)
138 0 : CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx2*2+1,nk12,-1)
139 0 : CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx3*2+1,nt,-1)
140 0 : sumis = 0.
141 0 : i = 0
142 0 : DO i3 = 0,stars%mx3*2
143 0 : k3 = i3
144 0 : IF (k3.GT.stars%mx3) k3 = k3 - stars%mx3*2-1
145 0 : DO i2 = 0,stars%mx2*2
146 0 : k2 = i2
147 0 : IF (k2.GT.stars%mx2) k2 = k2 - stars%mx2*2-1
148 0 : DO i1 = 0,stars%mx1*2
149 0 : k1 = i1
150 0 : IF (k1.GT.stars%mx1) k1 = k1 - stars%mx1*2-1
151 0 : i = i + 1
152 0 : phase = stars%rgphs(k1,k2,k3)
153 0 : id3 = stars%ig(k1,k2,k3)
154 0 : IF (id3.NE.0) THEN
155 0 : sumis = sumis + REAL(CMPLX(af3(i,1),bf3(i,1))* CONJG(stars%ustep(id3)))*phase*fact
156 : END IF
157 : ENDDO
158 : ENDDO
159 : ENDDO
160 0 : pdis(0,0,num)=SQRT(sumis/cell%volint)*1000.
161 0 : dis(num) = dis(num) + sumis
162 0 : IF (input%film) THEN
163 : ! ----> vacuum part
164 0 : DO ivac = 1,vacuum%nvac
165 0 : DO ip = 1,vacuum%nmzxy
166 : ! ---> create density in the real space
167 : i = 0
168 0 : DO i2 = 0,stars%mx2*2
169 0 : k2 = i2
170 0 : IF (k2.GT.stars%mx2) k2 = k2 - stars%mx2*2-1
171 0 : DO i1 = 0,stars%mx1*2
172 0 : k1 = i1
173 0 : IF (k1.GT.stars%mx1) k1 = k1 - stars%mx1*2-1
174 0 : i = i + 1
175 0 : id3 = stars%ig(k1,k2,0)
176 0 : IF (id3.EQ.0) THEN
177 0 : af2(i,1) = 0.
178 0 : bf2(i,1) = 0.
179 : ELSE
180 0 : id2 = stars%ig2(id3)
181 0 : phase = stars%rgphs(k1,k2,0)
182 0 : IF (id2.EQ.1) THEN
183 0 : af2(i,1) = rhv0(ip,ivac,num,1)
184 0 : bf2(i,1) = 0.
185 : ELSE
186 0 : af2(i,1) = REAL(rhv1(ip,id2-1,ivac,num,1))
187 0 : bf2(i,1) =AIMAG(rhv1(ip,id2-1,ivac,num,1))
188 : END IF
189 : END IF
190 : ENDDO
191 : ENDDO
192 0 : CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx1*2+1,stars%mx1*2+1,1)
193 0 : CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx2*2+1,nk12,1)
194 : ! ---> form dot product
195 0 : i = 0
196 0 : DO i2 = 0,stars%mx2*2
197 0 : DO i1 = 0,stars%mx1*2
198 0 : i = i + 1
199 0 : af2(i,1) = af2(i,1)*af2(i,1)
200 0 : bf2(i,1) = 0.
201 : ENDDO
202 : ENDDO
203 : ! ---> back fft
204 0 : CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx1*2+1,stars%mx1*2+1,-1)
205 0 : CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx2*2+1,nk12,-1)
206 0 : disz(npz-ip) = cell%area*af2(1,1)/REAL(nk12)
207 : ENDDO
208 : ! ---> beyond warping region
209 0 : DO ip = vacuum%nmzxy + 1,vacuum%nmz
210 0 : disz(npz-ip) = cell%area*rhv0(ip,ivac,num,1)
211 : ENDDO
212 0 : CALL intgz0(disz,vacuum%delz,vacuum%nmz,sumz,tail)
213 0 : dis(num) = dis(num) + 2.0/vacuum%nvac*sumz
214 : ENDDO
215 : END IF
216 : ENDDO
217 :
218 0 : dis(:) = SQRT(dis(:)/cell%vol)*1000.
219 :
220 0 : WRITE (oUnit,FMT=8000) iter,dis(1)
221 : 8000 FORMAT (/,'----> distance of the potential for it=',i3,':',f11.6, ' mhtr/bohr**3')
222 0 : WRITE(oUnit,*) "Details of potential differences for each atom type"
223 0 : WRITE(oUnit,*) "Atom: total difference: difference of first sphhar"
224 0 : DO n=1,atoms%ntype
225 0 : WRITE(oUnit,"(i5,' : ',f10.6,' : ',4f10.6)") n,pdis(4,n,1),pdis(0:3,n,1)
226 : ENDDO
227 0 : WRITE(oUnit,*) "Difference of interstitial:",pdis(0,0,1)
228 0 : IF (input%jspins.EQ.2) THEN
229 0 : WRITE (oUnit,FMT=8010) iter,dis(2)
230 : 8010 FORMAT (/,'----> distance of spin potential for it=',i3,':', f11.6,' mhtr/bohr**3')
231 0 : WRITE(oUnit,*) "Details of potential differences for each atom type"
232 0 : WRITE(oUnit,*) "Atom: total difference: difference of first sphhar"
233 0 : DO n=1,atoms%ntype
234 0 : WRITE(oUnit,"(i5,' : ',f10.6,' : ',4f10.6)") n,pdis(4,n,2),pdis(0:3,n,2)
235 : ENDDO
236 0 : WRITE(oUnit,*) "Difference of interstitial:",pdis(0,0,2)
237 : END IF
238 : CONTAINS
239 0 : SUBROUTINE priv_cdndif(m1,m2,m3,n1,n2,n3,s1,s2,film,rhsp,rhpw,rhv0,rhv1)
240 : IMPLICIT NONE
241 : INTEGER, INTENT (IN) :: m1,m2,m3,n1,n2,n3
242 : REAL, INTENT (IN) :: s1,s2
243 : COMPLEX, INTENT (INOUT) :: rhpw(:,:,:)
244 : COMPLEX, INTENT (INOUT) :: rhv1(:,:,:,:,:)
245 : REAL, INTENT (INOUT) :: rhsp(:,:,:,:,:)
246 : REAL, INTENT (INOUT) :: rhv0(:,:,:,:)
247 : LOGICAL,INTENT (IN) :: film
248 :
249 0 : rhsp(:,:,:,m1,n1)=s1*rhsp(:,:,:,m2,n2) + s2*rhsp(:,:,:,m3,n3)
250 0 : rhpw(:,m1,n1) =s1*rhpw(:,m2,n2) + s2*rhpw(:,m3,n3)
251 0 : IF (film) THEN
252 0 : rhv0(:,:,m1,n1) =s1*rhv0(:,:,m2,n2) + s2*rhv0(:,:,m3,n3)
253 0 : rhv1(:,:,:,m1,n1)=s1*rhv1(:,:,:,m2,n2) + s2*rhv1(:,:,:,m3,n3)
254 : ENDIF
255 0 : END SUBROUTINE priv_cdndif
256 :
257 : END SUBROUTINE potdis
258 : END MODULE m_potdis
|
