Line data Source code
1 : MODULE m_cdntot
2 : #ifdef CPP_MPI
3 : use mpi
4 : #endif
5 : ! ********************************************************
6 : ! calculate the total charge density in the interstial.,
7 : ! vacuum, and mt regions c.l.fu
8 : ! ********************************************************
9 : CONTAINS
10 1531 : SUBROUTINE integrate_cdn(stars,nococonv,atoms,sym,vacuum,input,cell , integrand, &
11 1531 : q, qis, qmt, qvac, qtot, qistot, fmpi)
12 : ! if called with fmpi variable, distribute the calculation of the pwint
13 : ! over fmpi processes in fmpi%mpi_comm
14 : USE m_intgr, ONLY : intgr3
15 : USE m_constants
16 : USE m_qsf
17 : USE m_pwint
18 : USE m_types
19 : USE m_juDFT
20 : IMPLICIT NONE
21 : TYPE(t_stars),INTENT(IN) :: stars
22 : TYPE(t_nococonv),INTENT(IN):: nococonv
23 : TYPE(t_atoms),INTENT(IN) :: atoms
24 : TYPE(t_sym),INTENT(IN) :: sym
25 : TYPE(t_vacuum),INTENT(IN) :: vacuum
26 : TYPE(t_input),INTENT(IN) :: input
27 : TYPE(t_cell),INTENT(IN) :: cell
28 :
29 : TYPE(t_potden),INTENT(IN) :: integrand
30 : REAL, INTENT(OUT) :: q(input%jspins), qis(input%jspins), qmt(atoms%ntype,input%jspins),&
31 : qvac(2,input%jspins), qtot, qistot
32 : TYPE(t_mpi),INTENT(IN),OPTIONAL :: fmpi
33 : INTEGER :: jsp, j, ivac, nz, n, irank, nsize, intstart, intstop, chunk_size, leftover
34 1531 : REAL :: q2(vacuum%nmz), w(4), rht1(vacuum%nmzd,2,input%jspins)
35 : REAL :: sum_over_ng3
36 1531 : COMPLEX,ALLOCATABLE :: x(:) !(1:stars%ng3), may be distributed over fmpi ranks
37 : COMPLEX :: w_off
38 : #ifdef CPP_MPI
39 : INTEGER ierr
40 : #endif
41 1531 : IF (PRESENT(fmpi)) THEN
42 1352 : irank = fmpi%irank
43 1352 : nsize = fmpi%isize
44 : ELSE
45 : irank = 0
46 : nsize = 1
47 : ENDIF
48 :
49 1531 : qtot = 0.0
50 1531 : qistot = 0.0
51 8869 : qvac=0.0
52 3977 : q=0.0
53 3977 : qis=0.0
54 8209 : qmt=0.0
55 1531 : IF (irank.EQ.0) THEN
56 : ! -----mt charge
57 2365 : DO n = 1,atoms%ntype
58 4045 : DO jsp=1,size(integrand%mt,4)
59 4045 : CALL intgr3(integrand%mt(:,0,n,jsp),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),w(jsp))
60 : enddo
61 1510 : if (size(integrand%pw,2)>2) THEN
62 : !this is a noco-calculation
63 307 : if (size(integrand%mt,4)==4) THEN
64 78 : w_off=cmplx(w(3),w(4))
65 : else
66 229 : w_off=0.0
67 : endif
68 : !rotate into global frame
69 307 : call nococonv%rotdenmat(n,w(1),w(2),w_off,toGlobal=.true.)
70 : endif
71 4744 : DO jsp=1,input%jspins
72 2379 : qmt(n, jsp) = w(jsp)*sfp_const
73 3889 : q(jsp) = q(jsp) + atoms%neq(n)*qmt(n,jsp)
74 : ENDDO
75 : ENDDO
76 855 : IF (input%film) THEN
77 : ! -----vacuum region
78 258 : DO jsp = 1,input%jspins
79 455 : DO ivac = 1,vacuum%nvac
80 49447 : DO nz = 1,vacuum%nmz
81 49447 : rht1(nz,ivac,jsp) = REAL(integrand%vac(nz,1,ivac,jsp))
82 : END DO
83 197 : CALL qsf(vacuum%delz,rht1(1,ivac,jsp),q2,vacuum%nmz,0)
84 197 : qvac(ivac,jsp) = q2(1)*cell%area
85 :
86 346 : q(jsp) = q(jsp) + qvac(ivac,jsp)*2./real(vacuum%nvac)
87 : ENDDO
88 : enddo
89 : END IF
90 : END IF ! irank = 0
91 :
92 3977 : DO jsp = 1,input%jspins
93 : ! -----is region
94 2446 : chunk_size = stars%ng3/nsize
95 2446 : leftover = stars%ng3 - chunk_size*nsize
96 2446 : IF ( leftover > irank ) THEN
97 904 : chunk_size = chunk_size + 1
98 904 : intstart = irank * chunk_size + 1
99 : ELSE
100 1542 : intstart = leftover * (chunk_size+1) + (irank - leftover) * chunk_size + 1
101 : ENDIF
102 2446 : intstop = intstart + chunk_size -1
103 7338 : ALLOCATE(x(intstart:intstop))
104 2446 : CALL pwint_all(stars,atoms,sym ,cell,intstart,intstop,x)
105 2446 : sum_over_ng3 = 0.0
106 3458706 : DO j = intstart,intstop
107 3458706 : sum_over_ng3 = sum_over_ng3 + integrand%pw(j,jsp)*x(j)*stars%nstr(j)
108 : ENDDO
109 2446 : DEALLOCATE(x)
110 : #ifdef CPP_MPI
111 2446 : IF (PRESENT(fmpi)) THEN
112 2140 : CALL MPI_reduce(sum_over_ng3,qis(jsp),1,MPI_DOUBLE_PRECISION,MPI_SUM,0,fmpi%mpi_comm,ierr)
113 : ELSE
114 306 : qis(jsp) = sum_over_ng3
115 : ENDIF
116 : #else
117 : qis(jsp) = sum_over_ng3
118 : #endif
119 2446 : qistot = qistot + qis(jsp)
120 2446 : q(jsp) = q(jsp) + qis(jsp)
121 3977 : qtot = qtot + q(jsp)
122 : END DO ! loop over spins
123 1531 : END SUBROUTINE integrate_cdn
124 :
125 0 : SUBROUTINE integrate_realspace(xcpot, atoms, sym, sphhar, input, &
126 0 : stars, cell, vacuum, noco, mt, is, hint)
127 : use m_types
128 : use m_mt_tofrom_grid
129 : use m_pw_tofrom_grid
130 : use m_constants
131 : implicit none
132 : CLASS(t_xcpot), INTENT(inout) :: xcpot
133 : TYPE(t_atoms),INTENT(IN) :: atoms
134 : TYPE(t_sym), INTENT(in) :: sym
135 : TYPE(t_sphhar), INTENT(IN) :: sphhar
136 : TYPE(t_input), INTENT(IN) :: input
137 : TYPE(t_stars), INTENT(IN) :: stars
138 : TYPE(t_cell), INTENT(IN) :: cell
139 :
140 : TYPE(t_vacuum), INTENT(in) :: vacuum
141 : TYPE(t_noco), INTENT(in) :: noco
142 : real, intent(inout) :: mt(:,:,:), is(:,:)
143 : character(len=*), intent(in), optional :: hint
144 : integer :: n_atm, i
145 :
146 0 : TYPE(t_potden) :: tmp_potden
147 0 : REAL :: q(input%jspins), qis(input%jspins), &
148 0 : qmt(atoms%ntype,input%jspins), qvac(2,input%jspins),&
149 : qtot, qistot
150 :
151 0 : call tmp_potden%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
152 0 : call init_mt_grid(input%jspins, atoms, sphhar, xcpot%needs_grad(), sym)
153 0 : do n_atm =1,atoms%ntype
154 : call mt_from_grid(atoms, sym, sphhar, n_atm, input%jspins, mt(:,:,n_atm), &
155 0 : tmp_potden%mt(:,0:,n_atm,:))
156 :
157 0 : do i=1,atoms%jri(n_atm)
158 0 : tmp_potden%mt(i,:,n_atm,:) = tmp_potden%mt(i,:,n_atm,:) * atoms%rmsh(i,n_atm)**2
159 : enddo
160 : enddo
161 0 : call finish_mt_grid()
162 :
163 0 : call init_pw_grid(stars, sym, cell,xcpot)
164 0 : call pw_from_grid( stars, is, tmp_potden%pw) !THIS CODE SEEMS TO BE BROKEN!!
165 0 : call finish_pw_grid()
166 :
167 0 : call judft_error("Bug, integrate_realspace in cdntot")
168 : !call integrate_cdn(stars,atoms,sym,vacuum,input,cell , tmp_potden, &
169 : ! q, qis, qmt, qvac, qtot, qistot)
170 :
171 0 : call print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
172 0 : END SUBROUTINE integrate_realspace
173 :
174 1531 : SUBROUTINE cdntot(stars,nococonv,atoms,sym,vacuum,input,cell ,&
175 : den,l_printData,qtot,qistot,fmpi,l_par)
176 :
177 : USE m_types
178 : USE m_juDFT
179 : IMPLICIT NONE
180 :
181 : ! .. Scalar Arguments ..
182 : TYPE(t_stars),INTENT(IN) :: stars
183 : TYPE(t_nococonv),INTENT(IN):: nococonv
184 : TYPE(t_atoms),INTENT(IN) :: atoms
185 : TYPE(t_sym),INTENT(IN) :: sym
186 : TYPE(t_vacuum),INTENT(IN) :: vacuum
187 : TYPE(t_input),INTENT(IN) :: input
188 :
189 : TYPE(t_cell),INTENT(IN) :: cell
190 : TYPE(t_potden),INTENT(IN) :: den
191 : LOGICAL,INTENT(IN) :: l_printData,l_par
192 : REAL,INTENT(OUT) :: qtot,qistot
193 : TYPE(t_mpi),INTENT(IN) :: fmpi
194 :
195 : ! .. Local Scalars ..
196 1531 : REAL q(input%jspins),qis(input%jspins),w,mtCharge
197 : ! ..
198 : ! .. Local Arrays ..
199 1531 : REAL qmt(atoms%ntype,input%jspins),qvac(2,input%jspins)
200 :
201 1531 : CALL timestart("cdntot")
202 1531 : IF (l_par) THEN
203 : CALL integrate_cdn(stars,nococonv,atoms,sym,vacuum,input,cell , den, &
204 1352 : q, qis, qmt, qvac, qtot, qistot, fmpi)
205 : ELSE
206 : CALL integrate_cdn(stars,nococonv,atoms,sym,vacuum,input,cell , den, &
207 179 : q, qis, qmt, qvac, qtot, qistot)
208 : ENDIF
209 :
210 1531 : IF (fmpi%irank.EQ.0) CALL cdntot_writings(atoms,vacuum,input,l_printData,q,qis,qmt,qvac,qtot)
211 1531 : CALL timestop("cdntot")
212 1531 : END SUBROUTINE cdntot
213 :
214 776 : SUBROUTINE cdntot_writings(atoms,vacuum,input,l_printData,q,qis,qmt,qvac,qtot)
215 :
216 : USE m_constants
217 : USE m_types
218 : USE m_juDFT
219 : USE m_xmlOutput
220 : IMPLICIT NONE
221 :
222 : ! .. Scalar Arguments ..
223 : TYPE(t_atoms),INTENT(IN) :: atoms
224 : TYPE(t_vacuum),INTENT(IN) :: vacuum
225 : TYPE(t_input),INTENT(IN) :: input
226 : LOGICAL,INTENT(IN) :: l_printData
227 : REAL,INTENT(IN) :: q(input%jspins),qis(input%jspins)
228 : REAL,INTENT(IN) :: qmt(atoms%ntype,input%jspins),qvac(2,input%jspins)
229 : REAL,INTENT(IN) :: qtot
230 :
231 : ! .. Local Scalars ..
232 : REAL mtCharge
233 : INTEGER i,jsp,n
234 : ! ..
235 : ! .. Local Arrays ..
236 776 : INTEGER, ALLOCATABLE :: lengths(:,:)
237 : CHARACTER(LEN=20) :: attributes(6), names(6)
238 :
239 :
240 776 : IF (input%film) THEN
241 392 : ALLOCATE(lengths(4+vacuum%nvac,2))
242 : ELSE
243 678 : ALLOCATE(lengths(4,2))
244 : END IF
245 :
246 2020 : DO jsp = 1,input%jspins
247 1244 : WRITE (oUnit,FMT=8000) jsp,q(jsp),qis(jsp), (qmt(n,jsp),n=1,atoms%ntype)
248 1420 : IF (input%film) WRITE (oUnit,FMT=8010) (i,qvac(i,jsp),i=1,vacuum%nvac)
249 3405 : mtCharge = SUM(qmt(1:atoms%ntype,jsp) * atoms%neq(1:atoms%ntype))
250 1244 : names(1) = 'spin' ; WRITE(attributes(1),'(i0)') jsp ; lengths(1,1)=4 ; lengths(1,2)=1
251 1244 : names(2) = 'total' ; WRITE(attributes(2),'(f14.7)') q(jsp) ; lengths(2,1)=5 ; lengths(2,2)=14
252 1244 : names(3) = 'interstitial' ; WRITE(attributes(3),'(f14.7)') qis(jsp) ; lengths(3,1)=12 ; lengths(3,2)=14
253 1244 : names(4) = 'mtSpheres' ; WRITE(attributes(4),'(f14.7)') mtCharge ; lengths(4,1)=9 ; lengths(4,2)=14
254 2020 : IF(l_printData) THEN
255 1100 : IF(input%film) THEN
256 270 : DO i = 1, vacuum%nvac
257 155 : WRITE(names(4+i),'(a6,i0)') 'vacuum', i
258 155 : WRITE(attributes(4+i),'(f14.7)') qvac(i,jsp)
259 155 : lengths(4+i,1)=7
260 270 : lengths(4+i,2)=14
261 : END DO
262 : CALL writeXMLElementFormPoly('spinDependentCharge',names(1:4+vacuum%nvac),&
263 115 : attributes(1:4+vacuum%nvac),lengths)
264 : ELSE
265 985 : CALL writeXMLElementFormPoly('spinDependentCharge',names(1:4),attributes(1:4),lengths)
266 : END IF
267 : END IF
268 : END DO ! loop over spins
269 776 : WRITE (oUnit,FMT=8020) qtot
270 776 : IF(l_printData) THEN
271 2073 : CALL writeXMLElementFormPoly('totalCharge',(/'value'/),(/qtot/),reshape((/5,20/),(/1,2/)))
272 : END IF
273 : 8000 FORMAT (/,10x,'total charge for spin',i3,'=',f12.6,/,10x,&
274 : 'interst. charge = ',f12.6,/,&
275 : (10x,'mt charge= ',4f12.6,/))
276 : 8010 FORMAT (10x,'vacuum ',i2,' charge= ',f12.6)
277 : 8020 FORMAT (/,10x,'total charge =',f12.6)
278 :
279 776 : END SUBROUTINE cdntot_writings
280 :
281 0 : SUBROUTINE print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
282 : use ieee_arithmetic
283 : implicit none
284 : REAL, INTENT(in) :: q(:), qis(:), qmt(:,:), qvac(:,:), qtot, qistot
285 : character(len=*), intent(in), optional :: hint
286 : integer :: n_mt
287 :
288 :
289 0 : if(present(hint)) write (*,*) "DEN of ", hint
290 0 : write (*,*) "q = ", q
291 0 : write (*,*) "qis = ", qis
292 :
293 0 : write (*,*) "qmt"
294 0 : do n_mt = 1,size(qmt, dim=1)
295 0 : write (*,*) "mt = ", n_mt, qmt(n_mt,:)
296 : enddo
297 :
298 : if(.not. any(ieee_is_nan(qvac))) then
299 0 : write (*, *) "qvac", qvac
300 : endif
301 0 : write (*, *) "qtot", qtot
302 0 : write (*, *) "qis_tot", qistot
303 0 : write (*, *) "-------------------------"
304 0 : END SUBROUTINE print_cdn_inte
305 : END MODULE m_cdntot
|
