Line data Source code
1 : MODULE m_checkdop
2 : CONTAINS
3 24 : SUBROUTINE checkdop(p,np,n,na,ivac,iflag,jsp,atoms,sphhar,stars,sym,&
4 : vacuum,cell,potden,potdenIm)
5 : ! ************************************************************
6 : ! subroutines checks the continuity of coulomb *
7 : ! potential or valence charge density *
8 : ! across the m.t. and vacuum boundaries *
9 : ! c.l.fu *
10 : ! (unifies vcheck and cdnchk g.b. *
11 : ! YM: this routine doesn't really work in the vacuum in 1D case yet
12 : ! ************************************************************
13 :
14 : USE m_types
15 : USE m_constants
16 : USE m_juDFT
17 : USE m_starf, ONLY : starf2,starf3
18 : USE m_angle
19 : USE m_ylm
20 : USE m_fitchk
21 :
22 : IMPLICIT NONE
23 :
24 : ! ..
25 : ! .. Scalar Arguments ..
26 :
27 : type(t_sphhar),intent(in) :: sphhar
28 : TYPE(t_stars),INTENT(IN) :: stars
29 : TYPE(t_atoms),INTENT(IN) :: atoms
30 : TYPE(t_sym),INTENT(IN) :: sym
31 : TYPE(t_vacuum),INTENT(IN) :: vacuum
32 : TYPE(t_cell),INTENT(IN) :: cell
33 : TYPE(t_potden),INTENT(IN) :: potden
34 : TYPE(t_potden),INTENT(IN), OPTIONAL :: potdenIm
35 : INTEGER, INTENT (IN) :: iflag,ivac,n,na,np,jsp
36 : !-odim
37 : !+odim
38 : ! .. Array Arguments ..
39 : REAL, INTENT (IN) :: p(:,:)!(3,DIMENSION%nspd)
40 : ! ..
41 : ! .. Local Scalars ..
42 : REAL av,dms,rms,s,ir2,help,phi,helpIm
43 : INTEGER i,j,k,lh,mem,nd,lm,ll1,nopa ,gz,m
44 : COMPLEX ic
45 : LOGICAL l_cdn, l_dfpt
46 : ! ..
47 : ! .. Local Arrays ..
48 24 : COMPLEX sf2(stars%ng2),sf3(stars%ng3),ylm( (atoms%lmaxd+1)**2 )
49 24 : REAL rcc(3),v1(SIZE(p,2)),v2(SIZE(p,2)),x(3),ri(3), v1Im(SIZE(p,2)), v2Im(SIZE(p,2))
50 : l_dfpt = .FALSE.
51 24 : l_dfpt = PRESENT(potdenIm)
52 24 : l_cdn = .FALSE. ! By default we assume that the input is a potential.
53 24 : IF (potden%potdenType.LE.0) CALL juDFT_error('unknown potden type', calledby='checkdop')
54 24 : IF (potden%potdenType.GT.1000) l_cdn = .TRUE. ! potdenTypes > 1000 are reserved for densities
55 :
56 :
57 : ! ..
58 : ! ..
59 : #ifdef __TOS_BGQ__
60 : RETURN
61 : #endif
62 24 : ic = CMPLX(0.,1.)
63 24 : IF (.NOT.iflag.LT.0) THEN
64 : ! ---> Interstitial part
65 0 : DO j = 1,np
66 0 : rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
67 0 : IF (l_dfpt) THEN
68 0 : CALL starf3(sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,p(:,j),sym%invtab,sf3,stars%center) ! Stars shifted by q
69 : ELSE
70 0 : CALL starf3(sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,p(:,j),sym%invtab,sf3)
71 : END IF
72 : !
73 0 : v1(j) = 0.0
74 0 : v1Im(j)=0.0
75 : If (l_dfpt) v1Im(j) = 0.0
76 0 : DO k = 1,stars%ng3
77 0 : v1(j) = v1(j) + REAL(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
78 0 : IF (l_dfpt) v1Im(j) = v1Im(j) + AIMAG(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
79 : ENDDO
80 : ENDDO
81 : ! ---> vacuum part
82 0 : IF (l_cdn) THEN
83 0 : IF (l_dfpt) THEN
84 0 : WRITE (oUnit,FMT=9005) ivac
85 : ELSE
86 0 : WRITE (oUnit,FMT=9000) ivac
87 : END IF
88 : ELSE
89 0 : IF (l_dfpt) THEN
90 0 : WRITE (oUnit,FMT=8005) ivac
91 : ELSE
92 0 : WRITE (oUnit,FMT=8000) ivac
93 : END IF
94 : ENDIF
95 0 : DO j = 1,np
96 0 : IF (l_dfpt) THEN
97 0 : CALL starf2(sym%nop2,stars%ng2,stars%kv2,sym%mrot,sym%symor,sym%tau,p(1:3,j),sym%invtab,sf2,stars%center)
98 : ELSE
99 0 : CALL starf2(sym%nop2,stars%ng2,stars%kv2,sym%mrot,sym%symor,sym%tau,p(1:3,j),sym%invtab,sf2)
100 : END IF
101 0 : v2(j)=0.0
102 0 : v2Im(j)=0.0
103 0 : IF ((norm2(stars%center)<1e-8)) THEN
104 0 : v2(j) = REAL(potden%vac(1,1,ivac,jsp))
105 0 : IF (l_dfpt) v2Im(j)=AIMAG(potden%vac(1,1,ivac,jsp)) !(l_dfpt Gamma Point)
106 0 : DO k = 2,stars%ng2
107 0 : v2(j) = v2(j) + REAL(potden%vac(1,k,ivac,jsp)*sf2(k))*stars%nstr2(k)
108 0 : IF (l_dfpt) v2Im(j) = v2Im(j) + AIMAG(potden%vac(1,k,ivac,jsp)*sf2(k))*stars%nstr2(k)
109 : END DO
110 : ELSE ! (l_dfpt)
111 0 : DO k = 1,stars%ng2
112 0 : v2(j) = v2(j) + REAL(potden%vac(1,k,ivac,jsp)*sf2(k))*stars%nstr2(k)
113 0 : v2Im(j) = v2Im(j) + AIMAG(potden%vac(1,k,ivac,jsp)*sf2(k))*stars%nstr2(k)
114 : ENDDO
115 : END IF
116 0 : rcc=MATMUL(cell%amat,p(:,j))
117 0 : IF (l_dfpt) THEN
118 0 : WRITE (oUnit,FMT=8025) (p(i,j),i=1,3),rcc, v1(j),v1Im(j),v2(j),v2Im(j)
119 : ELSE
120 0 : WRITE (oUnit,FMT=8020) (p(i,j),i=1,3),rcc, v1(j),v2(j)
121 : END IF
122 : END DO
123 0 : CALL fitchk(v1(:np),v2(:np),av,rms,dms)
124 0 : WRITE (oUnit,FMT=8030) av,rms,dms
125 0 : IF (l_dfpt) THEN
126 0 : CALL fitchk(v1Im(:np),v2Im(:np),av,rms,dms)
127 0 : WRITE (oUnit,*) "Imaginary Part:"
128 0 : WRITE (oUnit,FMT=8030) av,rms,dms
129 : END IF
130 0 : RETURN
131 : ENDIF
132 : ! ----> interstitial part
133 8424 : DO j = 1,np
134 134400 : rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
135 8400 : IF (l_dfpt) THEN
136 0 : CALL starf3(sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,sf3,stars%center)
137 : ELSE
138 8400 : CALL starf3(sym%nop,stars%ng3,sym%symor,stars%kv3,sym%mrot,sym%tau,rcc,sym%invtab,sf3)
139 : END IF
140 : !
141 8400 : v1(j) = 0.0
142 8400 : IF (l_dfpt) v1Im(j) = 0.0
143 5594424 : DO k = 1,stars%ng3
144 5586000 : v1(j) = v1(j) + REAL(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
145 5594400 : IF (l_dfpt) v1Im(j) = v1Im(j) + AIMAG(potden%pw(k,jsp)*sf3(k))*stars%nstr(k)
146 : ENDDO
147 : ENDDO
148 : ! ----> m.t. part
149 24 : IF (l_cdn) THEN
150 8 : IF (l_dfpt) THEN
151 0 : WRITE (oUnit,FMT=9015) n
152 : ELSE
153 8 : WRITE (oUnit,FMT=9010) n
154 : END IF
155 : ELSE
156 16 : IF (l_dfpt) THEN
157 0 : WRITE (oUnit,FMT=8015) n
158 : ELSE
159 16 : WRITE (oUnit,FMT=8010) n
160 : END IF
161 : ENDIF
162 24 : ir2 = 1.0
163 24 : IF (l_cdn) ir2 = 1.0 / ( atoms%rmt(n)*atoms%rmt(n) )
164 24 : nd = sym%ntypsy(na)
165 24 : nopa = sym%ngopr(na)
166 8424 : DO j = 1,np
167 33600 : DO i = 1,3
168 33600 : x(i) = p(i,j) - atoms%pos(i,na)
169 : ENDDO
170 : ! new
171 8400 : IF (nopa.NE.1) THEN
172 0 : rcc=MATMUL(cell%bmat,x)/tpi_const
173 :
174 0 : DO i = 1,3 ! rotate into representative
175 0 : ri(i) = 0.
176 0 : DO k = 1,3
177 0 : ri(i) = ri(i) + sym%mrot(i,k,nopa)*rcc(k)
178 : ENDDO
179 : ENDDO
180 0 : x=MATMUL(cell%amat,ri)
181 :
182 : END IF
183 : ! new
184 8400 : CALL ylm4(atoms%lmax(n),x,ylm)
185 8400 : help = 0.0
186 8400 : helpIm= 0.0
187 142800 : DO lh = 0,sphhar%nlh(nd)
188 134400 : s = 0.0
189 134400 : ll1 = sphhar%llh(lh,nd) * ( sphhar%llh(lh,nd) + 1 ) + 1
190 344400 : DO mem = 1,sphhar%nmem(lh,nd)
191 210000 : lm = ll1 + sphhar%mlh(mem,lh,nd)
192 344400 : s = s + REAL( sphhar%clnu(mem,lh,nd)* ylm(lm) )
193 : ENDDO
194 134400 : help = help + potden%mt(atoms%jri(n),lh,n,jsp) * s
195 142800 : IF (l_dfpt) helpIm=helpIm + potdenIm%mt(atoms%jri(n),lh,n,jsp) * s
196 : ENDDO
197 8400 : v2(j) = help * ir2
198 8400 : IF (l_dfpt) v2Im(j) = helpIm * ir2
199 8424 : IF (j.LE.8) THEN
200 3072 : rcc=MATMUL(cell%bmat,p(:,j))/tpi_const
201 :
202 192 : IF (l_dfpt) THEN
203 0 : WRITE (oUnit,FMT=8025) rcc, (p(i,j),i=1,3),v1(j),v1Im(j),v2(j),v2Im(j)
204 : ELSE
205 192 : WRITE (oUnit,FMT=8020) rcc, (p(i,j),i=1,3),v1(j),v2(j)
206 : END IF
207 : END IF
208 : END DO
209 24 : CALL fitchk(v1(:np),v2(:np),av,rms,dms)
210 24 : WRITE (oUnit,FMT=8030) av,rms,dms
211 24 : IF (l_dfpt) THEN
212 0 : CALL fitchk(v1Im(:np),v2Im(:np),av,rms,dms)
213 0 : WRITE (oUnit,*) "Imaginary Part:"
214 0 : WRITE (oUnit,FMT=8030) av,rms,dms
215 : END IF
216 : 8000 FORMAT (/,' int.-vac. boundary (potential): ivac=',i2,/,t10,&
217 : & 'int-coord',t36,'cart-coord',t57,' inter. ',t69,' vacuum ')
218 : 8005 FORMAT (/,' int.-vac. boundary (potential): ivac=',i2,/,t10,&
219 : & 'int-coord',t36,'cart-coord',t64,' inter. ',t87,' vacuum ')
220 : 8010 FORMAT (/,' int.-m.t. boundary (potential): atom type=',i2,/,&
221 : & t10,'int-coord',t36,'cart-coord',t57,' inter. ',t69,&
222 : & ' m. t. ')
223 : 8015 FORMAT (/,' int.-m.t. boundary (potential): atom type=',i2,/,&
224 : & t10,'int-coord',t36,'cart-coord',t64,' inter. ',t87,&
225 : & ' m. t. ')
226 : 8020 FORMAT (1x,2 (3f8.3,2x),2f12.6)
227 : 8025 FORMAT (1x,2 (3f8.3,2x),4f12.6) !;_dfpt
228 : 8030 FORMAT (/,10x,'average value = ',f12.6,/,t11,'rms,dmx=',2f10.3,&
229 : & ' per cent')
230 : 9000 FORMAT (/,' int.-vac. boundary (density): ivac=',i2,/,t10,&
231 : & 'int-coord',t36,'cart-coord',t57,' inter. ',t69,' vacuum ')
232 : 9005 FORMAT (/,' int.-vac. boundary (density): ivac=',i2,/,t10,&
233 : & 'int-coord',t36,'cart-coord',t64,' inter. ',t87,' vacuum ')
234 : 9010 FORMAT (/,' int.-m.t. boundary (density): atom type=',i2,/,t10,&
235 : & 'int-coord',t36,'cart-coord',t57,' inter. ',t69,' m. t. ')
236 : 9015 FORMAT (/,' int.-m.t. boundary (density): atom type=',i2,/,t10,&
237 : & 'int-coord',t36,'cart-coord',t64,' inter. ',t87,' m. t. ')
238 : RETURN
239 : END SUBROUTINE checkdop
240 : END MODULE m_checkdop
|
