Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2019 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_plot
7 : USE m_types
8 : USE m_juDFT
9 : USE m_constants
10 :
11 : IMPLICIT NONE
12 :
13 : !-----------------------------------------------------------------------------
14 : ! A general purpose plotting routine for FLEUR.
15 : !
16 : ! Based on the older plotting routines pldngen.f90 and plotdop.f90 that were
17 : ! originally called by optional.F90 and are now used in the scf-loop instead.
18 : ! At the cost of reduced postprocess functionality, this allowed us to re-
19 : ! move I/O (using plot.hdf/text files) from the plotting routine completely.
20 : !
21 : !
22 : ! A. Neukirchen & R. Hilgers, September 2019
23 : !
24 : ! Added OMP+MPI Parallelization, R. Hilgers March 2020
25 : !
26 : ! Added 3D vectorplot option and plotting only MT option as well as
27 : ! single MT only plots, A. Neukirchen and R. Hilgers April 2020
28 : !
29 : !-----------------------------------------------------------------------------
30 :
31 : PUBLIC :: checkplotinp, vectorsplit, matrixsplit, savxsf, vectorplot, &
32 : matrixplot, makeplots, procplot, getMTSphere
33 :
34 : CONTAINS
35 :
36 0 : SUBROUTINE checkplotinp(fmpi)
37 :
38 : !--------------------------------------------------------------------------
39 : ! Checks for existing plot input. If an ancient plotin file is used, an
40 : ! error is called. If no usable plot_inp exists, a new one is generated.
41 : !--------------------------------------------------------------------------
42 : TYPE(t_mpi), INTENT(IN) :: fmpi
43 : LOGICAL :: oldform,newform
44 :
45 0 : oldform=.FALSE.
46 0 : IF(fmpi%irank.EQ.0) INQUIRE(file="plotin",exist = oldform)
47 :
48 0 : IF (oldform) THEN
49 0 : CALL juDFT_error("Use of plotin file no longer supported", calledby="plot")
50 : END IF
51 0 : newform=.FALSE.
52 0 : IF(fmpi%irank.EQ.0) INQUIRE(file="plot_inp",exist = newform)
53 0 : IF (newform) THEN
54 0 : CALL juDFT_error("Use of plot_inp file no longer supported", calledby="plot")
55 : END IF
56 :
57 :
58 0 : END SUBROUTINE checkplotinp
59 :
60 0 : SUBROUTINE vectorsplit(stars, atoms, sphhar, vacuum, input, factor, denmat, &
61 : cden,mden)
62 :
63 : !--------------------------------------------------------------------------
64 : ! Takes a spin-polarized density vector and rearranges it into two plottable
65 : ! seperate ones, i.e. (rho_up, rho_down) ---> n, m.
66 : !
67 : ! Can also be applied to a potential with an additional factor, i.e. while
68 : ! the densities n/m are defined as (rho_up+-rho_down), V_eff/B_eff are de-
69 : ! fined as (V_up+-V_down)/2.
70 : !--------------------------------------------------------------------------
71 :
72 : TYPE(t_stars), INTENT(IN) :: stars
73 : TYPE(t_atoms), INTENT(IN) :: atoms
74 : TYPE(t_sphhar), INTENT(IN) :: sphhar
75 : TYPE(t_vacuum), INTENT(IN) :: vacuum
76 : TYPE(t_input), INTENT(IN) :: input
77 : REAL, INTENT(IN) :: factor
78 : TYPE(t_potden), INTENT(IN) :: denmat
79 : TYPE(t_potden), INTENT(OUT) :: cden, mden
80 :
81 0 : IF (factor==1.0) THEN
82 : CALL cden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
83 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
84 : POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,&
85 0 : stars%ng2)
86 : CALL mden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
87 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
88 : POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,&
89 0 : stars%ng2)
90 : ELSE
91 : CALL cden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
92 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
93 : POTDEN_TYPE_POTTOT,vacuum%nmzd,&
94 0 : vacuum%nmzxyd,stars%ng2)
95 : CALL mden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
96 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
97 : POTDEN_TYPE_POTTOT,vacuum%nmzd,&
98 0 : vacuum%nmzxyd,stars%ng2)
99 : END IF
100 :
101 0 : cden%mt(:,0:,1:,1) = (denmat%mt(:,0:,1:,1)+denmat%mt(:,0:,1:,2))/factor
102 0 : cden%pw(1:,1) = (denmat%pw(1:,1)+denmat%pw(1:,2))/factor
103 0 : IF (input%film) THEN
104 0 : cden%vac(1:,1:,1:,1) = (denmat%vac(1:,1:,1:,1)+denmat%vac(1:,1:,1:,2))/factor
105 : END IF
106 :
107 0 : mden%mt(:,0:,1:,1) = (denmat%mt(:,0:,1:,1)-denmat%mt(:,0:,1:,2))/factor
108 0 : mden%pw(1:,1) = (denmat%pw(1:,1)-denmat%pw(1:,2))/factor
109 0 : IF (input%film) THEN
110 0 : mden%vac(1:,1:,1:,1) = (denmat%vac(1:,1:,1:,1)-denmat%vac(1:,1:,1:,2))/factor
111 : END IF
112 :
113 0 : END SUBROUTINE vectorsplit
114 :
115 64 : SUBROUTINE matrixsplit(sym,stars, atoms, sphhar, vacuum, input, noco,nococonv, factor, &
116 : denmat, cden, mxden, myden, mzden)
117 : USE m_fft2d
118 : USE m_fft3d
119 : USE m_rotdenmat
120 :
121 : !--------------------------------------------------------------------------
122 : ! Takes a 2x2 density matrix and rearranges it into four plottable seperate
123 : ! ones, i.e. ((rho_11, rho_12),(rho_21, rho_22)) ---> n, mx, my, mz.
124 : !
125 : ! This is an adaptation of the old pldngen.f90 by P. Kurz.
126 : !
127 : ! Can also be applied to a potential with additional factors, i.e. while
128 : ! the densities n/m are defined as (rho_up+-rho_down), V_eff/B_eff are de-
129 : ! fined as (V_up+-V_down)/2 and while rho_21 is of the form a-i*b, for the
130 : ! potential we have a+i*b.
131 : !
132 : ! TODO: Find out, whether the latter form of rho is still valid. Having the
133 : ! additive rho_21=m_x+i*m_y would be nicer and more convenient/consistent.
134 : !--------------------------------------------------------------------------
135 :
136 : IMPLICIT NONE
137 : TYPE(t_sym), INTENT(IN) :: sym
138 : TYPE(t_stars), INTENT(IN) :: stars
139 : TYPE(t_atoms), INTENT(IN) :: atoms
140 : TYPE(t_sphhar), INTENT(IN) :: sphhar
141 : TYPE(t_vacuum), INTENT(IN) :: vacuum
142 : TYPE(t_input), INTENT(IN) :: input
143 : TYPE(t_noco), INTENT(IN) :: noco
144 : TYPE(t_nococonv),INTENT(IN) :: nococonv
145 : REAL, INTENT(IN) :: factor
146 : TYPE(t_potden), INTENT(IN) :: denmat
147 : TYPE(t_potden), INTENT(OUT) :: cden, mxden, myden, mzden
148 :
149 : ! Local scalars: iteration indices, matrix elements etc.
150 : INTEGER iden, ivac, ifft2, ifft3, imz, ityp, iri, ilh, imesh
151 : REAL cdnup, cdndown, chden, mgden, theta, phi, rhotot, mx, my, mz
152 : REAL zero, rziw, vz_r, vz_i, rho_11, rho_22, rho_21r, rho_21i, cdn11, cdn22
153 : COMPLEX czero, cdn21
154 :
155 : ! Local arrays: densities in real space and off-diagonal elements.
156 64 : REAL, ALLOCATABLE :: rvacxy(:,:,:,:), ris(:,:), ris2(:,:), fftwork(:)
157 64 : REAL, ALLOCATABLE :: rho(:,:,:,:), rht(:,:,:)
158 64 : COMPLEX, ALLOCATABLE :: qpw(:,:), qpww(:,:), rhtxy(:,:,:,:)
159 64 : COMPLEX, ALLOCATABLE :: cdom(:), cdomw(:), cdomvz(:,:), cdomvxy(:,:,:)
160 : complex :: mat(2,2)
161 64 : COMPLEX :: rhfull(stars%ng2)
162 :
163 64 : zero = 0.0; czero = CMPLX(0.0,0.0)
164 64 : ifft3 = 27*stars%mx1*stars%mx2*stars%mx3
165 64 : ifft2 = 9*stars%mx1*stars%mx2
166 :
167 : ! Allocation of arrays and initialization of those that make up the real
168 : ! space density matrix.
169 0 : ALLOCATE (rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,4), qpw(stars%ng3,4), &
170 0 : qpww(stars%ng3,4), ris2(0:27*stars%mx1*stars%mx2*stars%mx3-1,4), &
171 0 : rht(vacuum%nmzd,2,4), rhtxy(vacuum%nmzxyd,stars%ng2-1,2,4), &
172 0 : cdomvz(vacuum%nmzd,2), cdomvxy(vacuum%nmzxyd,stars%ng2-1,2), &
173 : cdom(stars%ng3), fftwork(0:27*stars%mx1*stars%mx2*stars%mx3-1), &
174 0 : ris(0:27*stars%mx1*stars%mx2*stars%mx3-1,4), cdomw(stars%ng3), &
175 2304 : rvacxy(0:9*stars%mx1*stars%mx2-1,vacuum%nmzxyd,2,4))
176 :
177 18504526 : rho(:,:,:,:) = zero; qpw(:,:) = czero; cdom(:) = czero
178 :
179 64 : IF (input%film) THEN
180 0 : cdomvz(:,:) = czero; rhtxy(:,:,:,:) = czero
181 0 : cdomvxy(:,:,:) = czero; rht(:,:,:) = zero
182 : END IF
183 :
184 8471968 : rho(:,0:,1:,:input%jspins) = denmat%mt(:,0:,1:,:input%jspins)
185 624300 : qpw(1:,:input%jspins) = denmat%pw(1:,:input%jspins)
186 64 : IF (ALLOCATED(denmat%pw_w)) THEN
187 12290 : qpww(1:,:input%jspins) = denmat%pw_w(1:,:input%jspins)
188 : END IF
189 64 : IF (input%film) THEN
190 0 : rht(1:,1:,:input%jspins) = REAL(denmat%vac(1:,1,1:,:input%jspins))
191 0 : rhtxy(1:,1:,1:,:input%jspins) = denmat%vac(1:vacuum%nmzxyd,2:,1:,:input%jspins)
192 : END IF
193 64 : IF(noco%l_noco) THEN
194 281305 : cdom = denmat%pw(:,3)
195 45 : IF (ALLOCATED(denmat%pw_w)) THEN
196 6146 : cdomw = denmat%pw_w(:,3)
197 : END IF
198 45 : IF (input%film) THEN
199 0 : cdomvz(:,:) = denmat%vac(:,1,:,3)
200 0 : cdomvxy = denmat%vac(:vacuum%nmzxyd,2:,:,3)
201 : END IF
202 : END IF
203 :
204 : ! Calculate the charge and magnetization densities in the muffin tins.
205 :
206 138 : DO ityp = 1,atoms%ntype
207 74 : theta = nococonv%beta(ityp)
208 74 : phi = nococonv%alph(ityp)
209 6132 : DO ilh = 0,sphhar%nlh(sym%ntypsy(atoms%firstAtom(ityp)))
210 : !$OMP PARALLEL DO DEFAULT(shared)& !Note: The default(shared) is because explicitly defining nococonv as shared is problematic for some compilers
211 : !$OMP SHARED(rho,ityp,ilh,denmat)&
212 6068 : !$OMP PRIVATE(iri,cdnup,cdndown,chden,mgden,cdn11,cdn21,cdn22,mat)
213 : DO iri = 1,atoms%jri(ityp)
214 : #if 1==1
215 : mat(1,1)=rho(iri,ilh,ityp,1)
216 : mat(2,2)=rho(iri,ilh,ityp,2)
217 : if (size(denmat%mt,4)>2) Then
218 : mat(1,2)=CMPLX(denmat%mt(iri,ilh,ityp,3), denmat%mt(iri,ilh,ityp,4))
219 : mat(2,1)=conjg(mat(1,2))
220 : else
221 : mat(1,2)=0;mat(2,1)=0
222 : endif
223 :
224 : call nococonv%rotdenmat(ityp,mat,toGlobal=.true.)
225 :
226 : rho(iri,ilh,ityp,:) =nococonv%denmat_to_mag(mat)
227 : #endif
228 : #if 1==2
229 : IF (SIZE(denmat%mt,4).LE.2) THEN
230 : cdnup = rho(iri,ilh,ityp,1)
231 : cdndown = rho(iri,ilh,ityp,2)
232 :
233 : chden = cdnup + cdndown
234 : mgden = cdnup - cdndown
235 : rho(iri,ilh,ityp,1) = chden
236 : rho(iri,ilh,ityp,2) = mgden*COS(phi)*SIN(theta)
237 : rho(iri,ilh,ityp,3) = mgden*SIN(phi)*SIN(theta)
238 : rho(iri,ilh,ityp,4) = mgden*COS(theta)
239 : ELSE
240 : cdn11 = rho(iri,ilh,ityp,1)
241 : cdn22 = rho(iri,ilh,ityp,2)
242 : cdn21 = CMPLX(denmat%mt(iri,ilh,ityp,3), denmat%mt(iri,ilh,ityp,4))
243 : IF(factor.NE.1.0) cdn21 = CMPLX(denmat%mt(iri,ilh,ityp,3), -denmat%mt(iri,ilh,ityp,4))
244 :
245 : CALL rot_den_mat(nococonv%alph(ityp), nococonv%beta(ityp), &
246 : cdn11,cdn22,cdn21)
247 :
248 : rho(iri,ilh,ityp,1) = cdn11 + cdn22
249 : rho(iri,ilh,ityp,2) = 2.0*REAL(cdn21)
250 : ! Note: The missing minus sign in the following line is a discrepancy
251 : ! from the other regions (IR, vac). But it was like that in version v0.26.
252 : rho(iri,ilh,ityp,3) = 2.0*AIMAG(cdn21)
253 : IF(factor.NE.1.0) rho(iri,ilh,ityp,3) = -2.0*AIMAG(cdn21)
254 : rho(iri,ilh,ityp,4) = cdn11 - cdn22
255 : END IF
256 : #endif
257 : END DO
258 : !$OMP END PARALLEL DO
259 : END DO
260 : END DO
261 :
262 : ! Fourier transform the diagonal part of the density matrix in the
263 : ! interstitial (qpw) to real space (ris).
264 192 : DO iden = 1,2
265 128 : CALL fft3d(ris(0:,iden),fftwork,qpw(1,iden),stars,1)
266 192 : IF (ALLOCATED(denmat%pw_w)) THEN
267 4 : CALL fft3d(ris2(0:,iden),fftwork,qpww(1,iden),stars,1)
268 : END IF
269 : END DO
270 :
271 : ! Also do that for the off-diagonal part. Real part goes into index
272 : ! 3 and imaginary part into index 4.
273 64 : CALL fft3d(ris(0:,3),ris(0:,4),cdom(1),stars,1)
274 64 : IF (ALLOCATED(denmat%pw_w)) THEN
275 2 : CALL fft3d(ris2(0:,3),ris2(0:,4),cdomw(1),stars,1)
276 : END IF
277 :
278 : ! Calculate the charge and magnetization densities in the interstitial.
279 : !$OMP PARALLEL DO DEFAULT(none)&
280 : !$OMP SHARED(ris,ris2,denmat,ifft3)&
281 64 : !$OMP PRIVATE(imesh,rho_11,rho_22,rho_21r,rhotot,rho_21i,mx,my,mz)
282 : DO imesh = 0,ifft3-1
283 : rho_11 = ris(imesh,1)
284 : rho_22 = ris(imesh,2)
285 : rho_21r = ris(imesh,3)
286 : rho_21i = ris(imesh,4)
287 : rhotot = rho_11 + rho_22
288 : mx = 2*rho_21r
289 : my = -2*rho_21i ! TODO: This is a magic minus. It should be 2*rho_21i.
290 : mz = (rho_11-rho_22)
291 : ris(imesh,1) = rhotot
292 : ris(imesh,2) = mx
293 : ris(imesh,3) = my
294 : ris(imesh,4) = mz
295 :
296 : IF (ALLOCATED(denmat%pw_w)) THEN
297 : rho_11 = ris2(imesh,1)
298 : rho_22 = ris2(imesh,2)
299 : rho_21r = ris2(imesh,3)
300 : rho_21i = ris2(imesh,4)
301 : rhotot = rho_11 + rho_22
302 : mx = 2*rho_21r
303 : my = -2*rho_21i ! TODO: This is a magic minus. It should be 2*rho_21i.
304 : mz = (rho_11-rho_22)
305 : ris2(imesh,1) = rhotot
306 : ris2(imesh,2) = mx
307 : ris2(imesh,3) = my
308 : ris2(imesh,4) = mz
309 : END IF
310 : END DO
311 : !$OMP END PARALLEL DO
312 :
313 : ! Invert the transformation to put the four densities back into
314 : ! reciprocal space.
315 320 : DO iden = 1,4
316 10191784 : fftwork=zero
317 256 : CALL fft3d(ris(0:,iden),fftwork,qpw(1,iden),stars,-1)
318 10191784 : fftwork=zero
319 320 : IF (ALLOCATED(denmat%pw_w)) THEN
320 8 : CALL fft3d(ris2(0:,iden),fftwork,qpww(1,iden),stars,-1)
321 : END IF
322 : END DO
323 :
324 : ! As above, but for the vacuum components in xy- and z-direction.
325 64 : IF (input%film) THEN
326 0 : DO iden = 1,2
327 0 : DO ivac = 1,vacuum%nvac
328 0 : DO imz = 1,vacuum%nmzxyd
329 0 : rziw = 0.0
330 0 : rhfull(1) = rht(imz,ivac,iden)
331 0 : rhfull(2:)= rhtxy(imz,:,ivac,iden)
332 : CALL fft2d(stars, rvacxy(0,imz,ivac,iden), fftwork, &
333 : rhfull, &
334 0 : 1)
335 : END DO
336 : END DO
337 : END DO
338 :
339 0 : DO ivac = 1,vacuum%nvac
340 0 : DO imz = 1,vacuum%nmzxyd
341 0 : rziw = 0.0
342 0 : vz_r = REAL(cdomvz(imz,ivac))
343 0 : vz_i = AIMAG(cdomvz(imz,ivac))
344 0 : rhfull(1) = cdomvz(imz,ivac)
345 0 : rhfull(2:)= cdomvxy(imz,:,ivac)
346 : CALL fft2d(stars, rvacxy(0,imz,ivac,3), rvacxy(0,imz,ivac,4), &
347 0 : rhfull, 1)
348 : END DO
349 : END DO
350 :
351 0 : DO ivac = 1,vacuum%nvac
352 0 : DO imz = 1,vacuum%nmzxyd
353 : !$OMP PARALLEL DO DEFAULT(none)&
354 : !$OMP SHARED(rvacxy,ivac,imz,ifft2)&
355 0 : !$OMP PRIVATE(imesh,rho_11,rho_22,rhotot,rho_21r,rho_21i,mx,my,mz)
356 : DO imesh = 0,ifft2-1
357 : rho_11 = rvacxy(imesh,imz,ivac,1)
358 : rho_22 = rvacxy(imesh,imz,ivac,2)
359 : rho_21r = rvacxy(imesh,imz,ivac,3)
360 : rho_21i = rvacxy(imesh,imz,ivac,4)
361 : rhotot = rho_11 + rho_22
362 : mx = 2*rho_21r
363 : my = -2*rho_21i
364 : mz = (rho_11-rho_22)
365 :
366 : rvacxy(imesh,imz,ivac,1) = rhotot
367 : rvacxy(imesh,imz,ivac,2) = mx
368 : rvacxy(imesh,imz,ivac,3) = my
369 : rvacxy(imesh,imz,ivac,4) = mz
370 : END DO
371 : !$OMP END PARALLEL DO
372 :
373 : END DO
374 : !$OMP PARALLEL DO DEFAULT(none)&
375 : !$OMP SHARED(rht,ivac,cdomvz,vacuum)&
376 0 : !$OMP PRIVATE(imz,rho_11,rho_22,rho_21r,rho_21i,mx,my,mz,rhotot)
377 : DO imz = vacuum%nmzxyd+1,vacuum%nmzd
378 : rho_11 = rht(imz,ivac,1)
379 : rho_22 = rht(imz,ivac,2)
380 : rho_21r = REAL(cdomvz(imz,ivac))
381 : rho_21i = AIMAG(cdomvz(imz,ivac))
382 : rhotot = rho_11 + rho_22
383 : mx = 2*rho_21r
384 : my = -2*rho_21i
385 : mz = (rho_11-rho_22)
386 :
387 : rht(imz,ivac,1) = rhotot
388 : rht(imz,ivac,2) = mx
389 : rht(imz,ivac,3) = my
390 : rht(imz,ivac,4) = mz
391 : END DO
392 : !$OMP END PARALLEL DO
393 : END DO
394 :
395 0 : DO iden = 1,4
396 0 : DO ivac = 1,vacuum%nvac
397 0 : DO imz = 1,vacuum%nmzxyd
398 0 : fftwork=zero
399 0 : rhfull(1) = rht(imz,ivac,iden)
400 0 : rhfull(2:)= rhtxy(imz,:,ivac,iden)
401 : CALL fft2d(stars, rvacxy(0,imz,ivac,iden), fftwork, &
402 : rhfull, &
403 0 : -1)
404 : END DO
405 : END DO
406 : END DO
407 : END IF
408 :
409 :
410 :
411 : ! Initialize and save the four output densities.
412 64 : IF (factor==1.0) THEN
413 : CALL cden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
414 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
415 : POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,&
416 62 : stars%ng2)
417 : CALL mxden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
418 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
419 : POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,&
420 62 : stars%ng2)
421 : CALL myden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
422 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
423 : POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,&
424 62 : stars%ng2)
425 : CALL mzden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
426 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
427 : POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,&
428 62 : stars%ng2)
429 :
430 3990354 : cden%mt(:,0:,1:,1) = rho(:,0:,1:,1)
431 305974 : cden%pw(1:,1) = qpw(1:,1)
432 62 : IF (input%film) THEN
433 0 : cden%vac(1:,1,1:,1) = rht(1:,1:,1)
434 0 : cden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,1)
435 : END IF
436 :
437 3990354 : mxden%mt(:,0:,1:,1) = rho(:,0:,1:,2)
438 305974 : mxden%pw(1:,1) = qpw(1:,2)
439 62 : IF (input%film) THEN
440 0 : mxden%vac(1:,1,1:,1) = rht(1:,1:,2)
441 0 : mxden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,2)
442 : END IF
443 :
444 3990354 : myden%mt(:,0:,1:,1) = rho(:,0:,1:,3)
445 305974 : myden%pw(1:,1) = qpw(1:,3)
446 62 : IF (input%film) THEN
447 0 : myden%vac(1:,1,1:,1) = rht(1:,1:,3)
448 0 : myden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,3)
449 : END IF
450 :
451 3990354 : mzden%mt(:,0:,1:,1) = rho(:,0:,1:,4)
452 305974 : mzden%pw(1:,1) = qpw(1:,4)
453 62 : IF (input%film) THEN
454 0 : mzden%vac(1:,1,1:,1) = rht(1:,1:,4)
455 0 : mzden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,4)
456 : END IF
457 :
458 62 : IF (ALLOCATED(denmat%pw_w)) THEN
459 0 : ALLOCATE (cden%pw_w, mold=cden%pw)
460 0 : ALLOCATE (mxden%pw_w, mold=mxden%pw)
461 0 : ALLOCATE (myden%pw_w, mold=myden%pw)
462 0 : ALLOCATE (mzden%pw_w, mold=mzden%pw)
463 0 : cden%pw_w(1:,1) = qpww(1:,1)
464 0 : mxden%pw_w(1:,1) = qpww(1:,2)
465 0 : myden%pw_w(1:,1) = qpww(1:,3)
466 0 : mzden%pw_w(1:,1) = qpww(1:,4)
467 : END IF
468 :
469 : ELSE
470 : CALL cden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
471 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
472 : POTDEN_TYPE_POTTOT,vacuum%nmzd,&
473 2 : vacuum%nmzxyd,stars%ng2)
474 : CALL mxden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
475 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
476 : POTDEN_TYPE_POTTOT,vacuum%nmzd,&
477 2 : vacuum%nmzxyd,stars%ng2)
478 : CALL myden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
479 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
480 : POTDEN_TYPE_POTTOT,vacuum%nmzd,&
481 2 : vacuum%nmzxyd,stars%ng2)
482 : CALL mzden%init_potden_simple(stars%ng3,atoms%jmtd,atoms%msh,sphhar%nlhd,&
483 : atoms%ntype,atoms%n_u,atoms%n_vPairs,1,.FALSE.,.FALSE.,&
484 : POTDEN_TYPE_POTTOT,vacuum%nmzd,&
485 2 : vacuum%nmzxyd,stars%ng2)
486 :
487 : ! Correction for the case of plotting the total potential.
488 : ! Needed due to the different definitons of density/potential matrices in
489 : ! FLEUR:
490 : !
491 : ! rhoMat = 0.5*((n +m_z,m_x+i*m_y),(m_x-i*m_y,n -m_z))
492 : ! vMat = ((V_eff+B_z,B_x-i*B_y),(B_x+i*m_y,V_eff-B_z))
493 :
494 982394 : rho(:,0:,1:,:) = rho(:,0:,1:,:)/2.0
495 24578 : qpw(1:,:) = qpw(1:,:)/2.0
496 4026 : rht(1:,1:,:) = rht(1:,1:,:)/2.0
497 26 : rhtxy(1:,1:,1:,:) = rhtxy(1:,1:,1:,:)/2.0
498 :
499 245598 : rho(:,0:,1:,3) = -rho(:,0:,1:,3)
500 6144 : qpw(1:,3) = -qpw(1:,3) ! TODO: This is a consequence of the magic minus.
501 1006 : rht(1:,1:,3) = -rht(1:,1:,3)
502 6 : rhtxy(1:,1:,1:,3) = -rhtxy(1:,1:,1:,3)
503 :
504 245598 : cden%mt(:,0:,1:,1) = rho(:,0:,1:,1)
505 6144 : cden%pw(1:,1) = qpw(1:,1)
506 2 : IF (input%film) THEN
507 0 : cden%vac(1:,1,1:,1) = rht(1:,1:,1)
508 0 : cden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,1)
509 : END IF
510 :
511 245598 : mxden%mt(:,0:,1:,1) = rho(:,0:,1:,2)
512 6144 : mxden%pw(1:,1) = qpw(1:,2)
513 2 : IF (input%film) THEN
514 0 : mxden%vac(1:,1,1:,1) = rht(1:,1:,2)
515 0 : mxden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,2)
516 : END IF
517 :
518 245598 : myden%mt(:,0:,1:,1) = rho(:,0:,1:,3)
519 6144 : myden%pw(1:,1) = qpw(1:,3)
520 2 : IF (input%film) THEN
521 0 : myden%vac(1:,1,1:,1) = rht(1:,1:,3)
522 0 : myden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,3)
523 : END IF
524 :
525 245598 : mzden%mt(:,0:,1:,1) = rho(:,0:,1:,4)
526 6144 : mzden%pw(1:,1) = qpw(1:,4)
527 2 : IF (input%film) THEN
528 0 : mzden%vac(1:,1,1:,1) = rht(1:,1:,4)
529 0 : mzden%vac(1:vacuum%nmzxyd,2:,1:,1) = rhtxy(1:,1:,1:,4)
530 : END IF
531 :
532 2 : IF (ALLOCATED(denmat%pw_w)) THEN
533 8 : ALLOCATE (cden%pw_w, mold=cden%pw)
534 8 : ALLOCATE (mxden%pw_w, mold=mxden%pw)
535 8 : ALLOCATE (myden%pw_w, mold=myden%pw)
536 8 : ALLOCATE (mzden%pw_w, mold=mzden%pw)
537 24578 : qpww(1:,:) = qpww(1:,:)/2.0
538 6144 : qpww(1:,3) = -qpww(1:,3) ! TODO: This is a consequence of the magic minus.
539 6144 : cden%pw_w(1:,1) = qpww(1:,1)
540 6144 : mxden%pw_w(1:,1) = qpww(1:,2)
541 6144 : myden%pw_w(1:,1) = qpww(1:,3)
542 6144 : mzden%pw_w(1:,1) = qpww(1:,4)
543 : END IF
544 :
545 : END IF
546 :
547 0 : DEALLOCATE (rho, qpw, rht, rhtxy, cdomvz, cdomvxy, &
548 64 : cdom, fftwork, ris, rvacxy)
549 :
550 64 : END SUBROUTINE matrixsplit
551 :
552 0 : SUBROUTINE savxsf(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, &
553 : noco, nococonv,score, potnorm, denName, denf, denA1, denA2, denA3)
554 : USE m_outcdn
555 : USE m_xsf_io
556 : #ifdef CPP_MPI
557 : USE mpi
558 : #endif
559 : USE m_polangle
560 : ! Takes one/several t_potden variable(s), i.e. scalar fields in MT-sphere/
561 : ! plane wave representation and makes it/them into plottable .xsf file(s)
562 : ! according to a scheme given in plot_inp.
563 : !
564 : ! This is an adaptation of the old plotdop.f90 by P. Kurz.
565 : !
566 : ! Naming convention:
567 : ! The output .xsf files will have names composed of the plotted density
568 : ! (c.f. identifier) and an optional tag. For a scalar density, only said
569 : ! density is plotted. For a spin-polarized density with an up and down
570 : ! component, f denotes the sum of both densities and g denotes their
571 : ! difference u-d possibly with additional modifying factors as explained
572 : ! above in vector-/matrixsplit. f and g are to be understood as scalar
573 : ! fields f(r_vec) and g(r_vec). For a density matrix, f is still the
574 : ! scalar part. Aditionally, three vector components A_[1-3] arise.
575 : !
576 : ! E.g.: scalar density rho (n(r_vec))
577 : ! ---> rho.xsf (n(r_vec))
578 : ! spin-polarized density rho (n_up(r_vec),n_down(r_vec))
579 : ! ---> rho_f.xsf, rho_g.xsf [n(r_vec),m(r_vec)]
580 : ! matrix density rho ((n_11(r_vec),n_12(r_vec)),
581 : ! (n_21(r_vec),n_22(r_vec)))
582 : ! ---> rho_f.xsf, rho_A1.xsf, rho_A2.xsf, rho_A3.xsf
583 : ! [n(r_vec),m_vec(r_vec)]
584 : type(t_sliceplot), INTENT(IN) :: sliceplot
585 : TYPE(t_stars), INTENT(IN) :: stars
586 : TYPE(t_atoms), INTENT(IN) :: atoms
587 : TYPE(t_sphhar), INTENT(IN) :: sphhar
588 : TYPE(t_vacuum), INTENT(IN) :: vacuum
589 : TYPE(t_mpi), INTENT(IN) :: fmpi
590 : TYPE(t_input), INTENT(IN) :: input
591 :
592 : TYPE(t_sym), INTENT(IN) :: sym
593 : TYPE(t_cell), INTENT(IN) :: cell
594 : TYPE(t_noco), INTENT(IN) :: noco
595 : TYPE(t_nococonv), INTENT(IN) :: nococonv
596 : LOGICAL, INTENT(IN) :: score, potnorm
597 : CHARACTER(len=20), INTENT(IN) :: denName
598 : TYPE(t_potden), INTENT(IN) :: denf
599 : TYPE(t_potden), OPTIONAL, INTENT(IN) :: denA1
600 : TYPE(t_potden), OPTIONAL, INTENT(IN) :: denA2
601 : TYPE(t_potden), OPTIONAL, INTENT(IN) :: denA3
602 :
603 : REAL :: tec, qint, phi0, angss
604 : INTEGER :: i, j, ix, iy, iz, na, nplo, iv, iflag, nfile
605 : INTEGER :: nplot, nt, jm, jspin, numInDen, numOutFiles
606 : LOGICAL :: twodim, cartesian, xsf, polar,unwind
607 :
608 0 : TYPE(t_potden), ALLOCATABLE :: den(:)
609 0 : REAL, ALLOCATABLE :: xdnout(:)
610 0 : REAL, ALLOCATABLE :: tempResults(:,:,:,:)
611 0 : REAL, ALLOCATABLE :: points(:,:,:,:)
612 0 : REAL, ALLOCATABLE :: tempVecs(:,:,:,:)
613 : REAL :: pt(3), vec1(3), vec2(3), vec3(3), &
614 : zero(3), help(3), qssc(3), point(3)
615 : INTEGER :: grid(3),k,strt,fin
616 0 : REAL :: rhocc(atoms%jmtd)
617 0 : CHARACTER (len=20), ALLOCATABLE :: outFilenames(:)
618 : CHARACTER (len=30) :: filename
619 : CHARACTER (len=7) :: textline
620 :
621 : REAL, PARAMETER :: eps = 1.0e-15
622 :
623 : !NAMELIST /plot/twodim,cartesian,unwind,vec1,vec2,vec3,grid,zero,phi0,filename
624 :
625 : integer:: ierr
626 :
627 0 : unwind = .FALSE.
628 0 : nfile = 120
629 :
630 0 : IF (PRESENT(denA2)) THEN
631 0 : ALLOCATE(den(4))
632 0 : den(1) = denf
633 0 : den(2) = denA1
634 0 : den(3) = denA2
635 0 : den(4) = denA3
636 0 : numInDen = 4
637 0 : numOutFiles = 4
638 0 : ELSE IF (PRESENT(denA1)) THEN
639 0 : ALLOCATE(den(2))
640 0 : den(1) = denf
641 0 : den(2) = denA1
642 0 : numInDen = 2
643 0 : numOutFiles = 2
644 : ELSE
645 0 : ALLOCATE(den(1))
646 0 : den(1) = denf
647 0 : numInDen = 1
648 0 : numOutFiles = 1
649 : END IF
650 :
651 0 : polar = sliceplot%polar
652 0 : xsf=sliceplot%format==PLOT_XSF_FORMAT
653 :
654 0 : IF((polar).AND.(.NOT.noco%l_noco)) THEN
655 0 : CALL juDFT_warn("l_noco=F and making polar plots is not compatible.",calledby="plot.f90")
656 : END IF
657 :
658 0 : IF (polar.AND.(numOutFiles==4)) THEN
659 0 : numOutFiles = 7
660 : END IF
661 :
662 0 : ALLOCATE(outFilenames(numOutFiles))
663 0 : ALLOCATE(xdnout(numOutFiles))
664 :
665 0 : DO i = 1, numInDen
666 :
667 : ! TODO: Does this work as intended?
668 0 : IF ((numInDen.NE.4).AND.(score).AND.(fmpi%irank.EQ.0)) THEN
669 0 : OPEN (17,file='cdnc',form='unformatted',status='old')
670 0 : REWIND 17
671 0 : DO jspin = 1, input%jspins
672 0 : DO nt = 1, atoms%ntype
673 0 : jm = atoms%jri(nt)
674 0 : READ (17) (rhocc(j),j=1,jm)
675 0 : DO j = 1, jm
676 0 : den(i)%mt(j,0,nt,jspin) = den(i)%mt(j,0,nt,jspin) - rhocc(j)/2.0/SQRT(pi_const)
677 : END DO
678 0 : READ (17) tec
679 : END DO
680 0 : READ (17) qint
681 0 : den(i)%pw(1,jspin) = den(i)%pw(1,jspin) - qint/cell%volint
682 : END DO
683 0 : CLOSE (17)
684 0 : ELSE IF (score) THEN
685 0 : CALL juDFT_error('Subtracting core charge in noco calculations not supported', calledby = 'plot')
686 : END IF
687 : END DO
688 :
689 0 : IF (noco%l_ss) THEN
690 0 : qssc = MATMUL(TRANSPOSE(cell%bmat),nococonv%qss)
691 : END IF
692 :
693 0 : IF (numOutFiles.EQ.1) THEN
694 0 : outFilenames(1) = TRIM(denName)
695 0 : ELSE IF (numOutFiles.EQ.2) THEN
696 0 : outFilenames(1) = TRIM(denName) // '_f'
697 0 : outFilenames(2) = TRIM(denName) // '_g'
698 : ELSE
699 0 : outFilenames(1) = TRIM(denName) // '_f'
700 0 : outFilenames(2) = TRIM(denName) // '_A1'
701 0 : outFilenames(3) = TRIM(denName) // '_A2'
702 0 : outFilenames(4) = TRIM(denName) // '_A3'
703 0 : IF (polar) THEN
704 0 : outFilenames(5) = TRIM(denName) // '_Aabs'
705 0 : outFilenames(6) = TRIM(denName) // '_Atheta'
706 0 : outFilenames(7) = TRIM(denName) // '_Aphi'
707 : END IF
708 : END IF
709 :
710 : ! If xsf is specified we create input files for xcrysden
711 0 : IF (xsf.AND.(fmpi%irank.EQ.0)) THEN
712 0 : DO i = 1, numOutFiles
713 0 : OPEN(nfile+i,file=TRIM(ADJUSTL(outFilenames(i)))//'.xsf',form='formatted')
714 0 : CALL xsf_WRITE_atoms(nfile+i,atoms,input%film ,cell%amat)
715 : END DO
716 : END IF
717 :
718 0 : IF (size(sliceplot%plot).EQ.0) THEN
719 0 : CALL juDFT_error('You set iplot/=0 without specifying a plot in the input.', calledby = 'plot')
720 : END IF
721 :
722 : ! Loop over all plots
723 0 : DO nplo = 1, size(sliceplot%plot)
724 0 : twodim=sliceplot%plot(nplo)%twodim
725 0 : cartesian=sliceplot%plot(nplo)%cartesian
726 0 : grid=sliceplot%plot(nplo)%grid
727 0 : IF(.NOT.(MODULO(grid(3),fmpi%isize)).EQ.0) CALL juDFT_error('Your grid z component doesnt fit the # of MPI processed. ',calledby='plot.F90')
728 0 : ALLOCATE(tempResults(0:grid(1)-1, 0:grid(2)-1,0:grid(3)-1,numOutFiles))
729 0 : ALLOCATE(points(0:grid(1)-1, 0:grid(2)-1,0:grid(3)-1,3))
730 0 : ALLOCATE(tempVecs(0:grid(1)-1, 0:grid(2)-1,0:grid(3)-1,6))
731 0 : vec1=sliceplot%plot(nplo)%vec1
732 0 : vec2=sliceplot%plot(nplo)%vec2
733 0 : vec3=sliceplot%plot(nplo)%vec3
734 0 : zero=sliceplot%plot(nplo)%zero
735 0 : filename=sliceplot%plot(nplo)%filename
736 :
737 0 : IF(.NOT.sliceplot%plot(nplo)%edgesAllSides) THEN
738 0 : vec1 = vec1 * float(grid(1)) / float(grid(1)+1)
739 0 : vec2 = vec2 * float(grid(2)) / float(grid(2)+1)
740 0 : IF(.NOT.input%film) vec3 = vec3 * float(grid(3)) / float(grid(3)+1)
741 : END IF
742 :
743 0 : IF (xsf.AND.sliceplot%plot(nplo)%vecField.AND.(fmpi%irank.EQ.0)) THEN
744 0 : OPEN(nfile+10,file=TRIM(denName)//'_A_vec_'//TRIM(filename)//'.xsf',form='formatted')
745 0 : CALL xsf_WRITE_atoms(nfile+10,atoms,input%film ,cell%amat)
746 : END IF
747 :
748 0 : IF (twodim.AND.ANY(grid(1:2)<1).AND.(fmpi%irank .EQ. 0)) &
749 0 : CALL juDFT_error("Illegal grid size in plot",calledby="plot")
750 0 : IF (.NOT.twodim.AND.ANY(grid<1).AND.(fmpi%irank .EQ. 0)) &
751 0 : CALL juDFT_error("Illegal grid size in plot",calledby="plot")
752 0 : IF (twodim) grid(3) = 1
753 :
754 : !calculate cartesian coordinates if needed
755 0 : IF (.NOT.cartesian) THEN
756 0 : vec1=matmul(cell%amat,vec1)
757 0 : vec2=matmul(cell%amat,vec2)
758 0 : vec3=matmul(cell%amat,vec3)
759 0 : zero=matmul(cell%amat,zero)
760 : END IF
761 :
762 : !Open the file
763 0 : IF (filename =="default".AND.(fmpi%irank .EQ. 0)) WRITE(filename,'(a,i2)') "plot",nplo
764 :
765 0 : IF (xsf) THEN
766 0 : DO i = 1, numOutFiles
767 0 : IF (fmpi%irank .EQ.0) CALL xsf_WRITE_header(nfile+i,twodim,filename,vec1,vec2,vec3,zero,grid)
768 : END DO
769 : ELSE
770 0 : IF (fmpi%irank .EQ. 0) OPEN (nfile,file = TRIM(ADJUSTL(denName))//'_'//filename,form='formatted')
771 : END IF
772 :
773 0 : IF ((.NOT.xsf).AND.fmpi%irank .EQ. 0) THEN
774 0 : IF (twodim) THEN
775 0 : IF (numOutFiles.EQ.1) THEN
776 0 : WRITE(nfile,'(3a15)') 'x','y','f'
777 0 : ELSE IF (numOutFiles.EQ.2) THEN
778 0 : WRITE(nfile,'(4a15)') 'x','y','f','g'
779 0 : ELSE IF (numOutFiles.EQ.4) THEN
780 0 : WRITE(nfile,'(6a15)') 'x','y','f','A1','A2','A3'
781 : ELSE
782 0 : WRITE(nfile,'(9a15)') 'x','y','f','A1','A2','A3','|A|','theta','phi'
783 : END IF
784 : ELSE
785 0 : IF (numOutFiles.EQ.1) THEN
786 0 : WRITE(nfile,'(4a15)') 'x','y','z','f'
787 0 : ELSE IF (numOutFiles.EQ.2) THEN
788 0 : WRITE(nfile,'(5a15)') 'x','y','z','f','g'
789 0 : ELSE IF (numOutFiles.EQ.4) THEN
790 0 : WRITE(nfile,'(7a15)') 'x','y','z','f','A1','A2','A3'
791 : ELSE
792 0 : WRITE(nfile,'(10a15)') 'x','y','z','f','A1','A2','A3','|A|','theta','phi'
793 : END IF
794 : END IF
795 : END IF
796 :
797 : #ifdef CPP_MPI
798 0 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr)
799 : #endif
800 0 : CALL timestart("loop over points")
801 : !print*, "loop over points", fmpi%irank
802 : !loop over all points
803 : !strt= fmpi%irank*(grid(3)-1)/fmpi%isize + 1
804 : !fin = ((fmpi%irank+1)*(grid(3)-1))/fmpi%isize
805 : !IF ( fmpi%irank == 0) strt = 0
806 0 : strt= fmpi%irank *grid(3)/fmpi%isize ! Proposed fix to issue #540.
807 0 : fin = (fmpi%irank+1)*grid(3)/fmpi%isize - 1 ! Continued
808 0 : DO iz = strt,fin
809 : !$OMP PARALLEL DO DEFAULT(none) &
810 : !$OMP& SHARED(iz,grid,zero,vec1,vec2,vec3,twodim,input,cell ,numInDen,potnorm) &
811 : !$OMP& SHARED(stars,vacuum,sphhar,atoms,sym,den,sliceplot,noco,points) &
812 : !$OMP& SHARED(unwind,polar,tempResults,tempVecs,nplo,qssc,xsf,NumOutFiles) &
813 0 : !$OMP& PRIVATE(ix,point,nt,na,pt,iv,iflag,xdnout,angss,help,phi0)
814 : DO iy = 0, grid(2)-1
815 : DO ix = 0, grid(1)-1
816 :
817 : point = zero + vec1*REAL(ix)/(grid(1)-1) +&
818 : vec2*REAL(iy)/(grid(2)-1)
819 : IF (.NOT.twodim) point = point + vec3*REAL(iz)/(grid(3)-1)
820 :
821 : ! Set region specific parameters for point
822 :
823 : ! Get MT sphere for point if point is in MT sphere
824 : CALL getMTSphere(input,cell,atoms ,point,nt,na,pt)
825 : IF (na.NE.0) THEN
826 : ! In MT sphere
827 : iv = 0
828 : iflag = 1
829 : ELSE IF (input%film.AND.ABS(point(3))>=cell%z1) THEN
830 : ! In vacuum in 2D system
831 : iv = 1
832 : iflag = 0
833 : pt(:) = point(:)
834 : ELSE
835 : ! In interstitial region
836 : iv = 0
837 : iflag = 2
838 : pt(:) = point(:)
839 : END IF
840 :
841 : DO i = 1, numInDen
842 : CALL outcdn(pt,nt,na,iv,iflag,1,potnorm,stars,&
843 : vacuum,sphhar,atoms,sym,cell ,&
844 : den(i),xdnout(i))
845 : END DO
846 :
847 : IF (sliceplot%plot(nplo)%onlyMT.AND.(iflag.NE.1)) THEN
848 : xdnout=0.0
849 : ELSE IF (sliceplot%plot(nplo)%typeMT.NE.0) THEN
850 : IF (sliceplot%plot(nplo)%typeMT.NE.nt) THEN
851 : xdnout=0.0
852 : END IF
853 : END IF
854 :
855 : IF (na.NE.0) THEN
856 : IF (noco%l_ss) THEN
857 : ! rotate magnetization "backward"
858 : angss = DOT_PRODUCT(qssc,pt-atoms%pos(:,na))
859 : help(1) = xdnout(2)
860 : help(2) = xdnout(3)
861 : xdnout(2) = +help(1)*COS(angss)+help(2)*SIN(angss)
862 : xdnout(3) = -help(1)*SIN(angss)+help(2)*COS(angss)
863 : ! xdnout(2)=0. ; xdnout(3)=0. ; xdnout(4)=0.
864 : END IF
865 : END IF
866 :
867 : IF (noco%l_ss .AND. (.NOT. unwind)) THEN
868 : ! rotate magnetization
869 : angss = DOT_PRODUCT(qssc,point)
870 : help(1) = xdnout(2)
871 : help(2) = xdnout(3)
872 : xdnout(2) = +help(1)*COS(angss) -help(2)*SIN(angss)
873 : xdnout(3) = +help(1)*SIN(angss) +help(2)*COS(angss)
874 : END IF
875 :
876 : IF (polar) THEN
877 : CALL pol_angle(xdnout(2),xdnout(3),xdnout(4),xdnout(6),xdnout(7))
878 : xdnout(5) = SQRT(ABS(xdnout(2)**2+xdnout(3)**2+xdnout(4)**2))
879 : xdnout(6)= xdnout(6)/pi_const
880 : xdnout(7)= xdnout(7)/pi_const
881 : END IF ! (polar)
882 : IF (xsf) THEN
883 : DO i = 1, numOutFiles
884 : tempResults(ix,iy,iz,i)=xdnout(i)
885 : END DO
886 : IF ((size(xdnout).GE.4).AND.sliceplot%plot(nplo)%vecField) THEN
887 : tempVecs(ix,iy,iz,1:3)=point(:)/1.8897269
888 : tempVecs(ix,iy,iz,4:6)=xdnout(2:4)
889 : ELSE IF (sliceplot%plot(nplo)%vecField.AND.(.NOT.noco%l_noco)) THEN
890 : CALL juDFT_warn("l_noco=F and making vector plots is not compatible. Do a regular plot for a spin-polarized system please, because vectors pointing only into z are not that interesting!",calledby="plot.f90")
891 : END IF
892 : ELSE
893 : tempResults(ix,iy,iz,:)=xdnout(:)
894 : points(ix,iy,iz,:)=point(:)/1.8897269
895 : END IF
896 : END DO !x-loop
897 : END DO !y-loop
898 : !$OMP END PARALLEL DO
899 : END DO !z-loop
900 0 : CALL timestop("loop over points")
901 :
902 0 : CALL timestart("output")
903 : !Print out results of the different MPI processes in correct order.
904 0 : IF(fmpi%irank.EQ.0) THEN
905 0 : IF (xsf) THEN
906 0 : DO i = 1, numOutFiles
907 0 : CLOSE(nfile+i)
908 : END DO
909 0 : IF (sliceplot%plot(nplo)%vecField) THEN
910 0 : CLOSE(nfile+10)
911 : END IF
912 : ELSE
913 0 : CLOSE(nfile)
914 : END IF
915 : END IF
916 0 : DO k=0, (fmpi%isize-1)
917 : #ifdef CPP_MPI
918 0 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr)
919 : #endif
920 0 : IF(fmpi%irank.EQ.k) THEN
921 :
922 0 : IF (xsf) THEN
923 0 : DO i = 1, numOutFiles
924 0 : OPEN(nfile+i,file=TRIM(ADJUSTL(outFilenames(i)))//'.xsf',form='formatted',position="append", action="write")
925 0 : DO iz = strt, fin
926 0 : DO iy = 0, grid(2)-1
927 0 : DO ix = 0, grid(1)-1
928 0 : WRITE(nfile+i,*) tempResults(ix,iy,iz,i)
929 : END DO
930 : END DO
931 : END DO
932 0 : CLOSE(nfile+i)
933 : END DO
934 :
935 0 : IF (sliceplot%plot(nplo)%vecField) THEN
936 0 : OPEN(nfile+10,file=TRIM(denName)//'_A_vec_'//TRIM(filename)//'.xsf',form='formatted',position="append", action="write")
937 0 : DO iz = strt, fin
938 0 : DO iy = 0, grid(2)-1
939 0 : DO ix = 0, grid(1)-1
940 0 : WRITE(nfile+10,'(a,6f20.11)') 'X', tempVecs(ix,iy,iz,:)
941 : END DO
942 : END DO
943 : END DO
944 0 : CLOSE(nfile+10)
945 : END IF
946 : ELSE
947 0 : OPEN (nfile,file = TRIM(ADJUSTL(denName))//'_'//filename,form='formatted',position="append", action="write")
948 0 : DO iz = strt, fin
949 0 : DO iy = 0, grid(2)-1
950 0 : DO ix = 0, grid(1)-1
951 0 : WRITE(nfile,'(10e15.7)') points(ix,iy,iz,:) ,tempResults(ix,iy,iz,:)
952 : END DO
953 : END DO
954 : END DO
955 0 : CLOSE(nfile)
956 : END IF
957 :
958 : END IF
959 : #ifdef CPP_MPI
960 0 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr)
961 : #endif
962 : END DO
963 0 : CALL timestop("output")
964 :
965 :
966 0 : IF (xsf.AND.(fmpi%irank.EQ.0)) THEN
967 0 : DO i = 1, numOutFiles
968 0 : OPEN(nfile+i,file=TRIM(ADJUSTL(outFilenames(i)))//'.xsf',form='formatted',position="append", action="write")
969 0 : CALL xsf_WRITE_endblock(nfile+i,twodim)
970 0 : CLOSE(nfile+i)
971 : END DO
972 : END IF
973 :
974 0 : DEALLOCATE(tempResults)
975 0 : DEALLOCATE(points)
976 : END DO !nplo
977 :
978 0 : DEALLOCATE(xdnout, outFilenames)
979 :
980 0 : END SUBROUTINE savxsf
981 :
982 0 : SUBROUTINE vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, &
983 : noco,nococonv, factor, score, potnorm, denmat, denName)
984 : #ifdef CPP_MPI
985 : USE mpi
986 : #endif
987 : ! Takes a spin-polarized t_potden variable, i.e. a 2D vector in MT-sphere/
988 : ! plane wave representation, splits it into two spinless ones, which are
989 : ! then passed on to the savxsf routine to get 2 .xsf files out.
990 : TYPE(t_sliceplot), INTENT(IN) :: sliceplot
991 : TYPE(t_stars), INTENT(IN) :: stars
992 : TYPE(t_atoms), INTENT(IN) :: atoms
993 : TYPE(t_sphhar), INTENT(IN) :: sphhar
994 : TYPE(t_vacuum), INTENT(IN) :: vacuum
995 : TYPE(t_input), INTENT(IN) :: input
996 :
997 : TYPE(t_mpi), INTENT(IN) :: fmpi
998 : TYPE(t_sym), INTENT(IN) :: sym
999 : TYPE(t_cell), INTENT(IN) :: cell
1000 : TYPE(t_noco), INTENT(IN) :: noco
1001 : TYPE(t_nococonv), INTENT(IN) :: nococonv
1002 : REAL, INTENT(IN) :: factor
1003 : LOGICAL, INTENT(IN) :: score, potnorm
1004 : TYPE(t_potden), INTENT(IN) :: denmat
1005 : CHARACTER(len=20), INTENT(IN) :: denName
1006 :
1007 0 : TYPE(t_potden) :: cden, mden
1008 :
1009 : CALL vectorsplit(stars, atoms, sphhar, vacuum, input, factor, denmat, &
1010 0 : cden, mden)
1011 :
1012 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, noco,nococonv, &
1013 0 : score, potnorm, denName, cden, mden)
1014 :
1015 0 : END SUBROUTINE vectorplot
1016 :
1017 0 : SUBROUTINE matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, &
1018 : noco, nococonv,factor, score, potnorm, denmat, denName)
1019 : #ifdef CPP_MPI
1020 : USE mpi
1021 : #endif
1022 : ! Takes a 2x2 t_potden variable, i.e. a sum of Pauli matrices in MT-
1023 : ! sphere/ plane wave representation and splits it into four spinless ones,
1024 : ! which are then passed on to the savxsf routine to get 4 .xsf files out.
1025 : TYPE(t_sliceplot), INTENT(IN) :: sliceplot
1026 : TYPE(t_stars), INTENT(IN) :: stars
1027 : TYPE(t_atoms), INTENT(IN) :: atoms
1028 : TYPE(t_sphhar), INTENT(IN) :: sphhar
1029 : TYPE(t_vacuum), INTENT(IN) :: vacuum
1030 : TYPE(t_input), INTENT(IN) :: input
1031 : TYPE(t_mpi), INTENT(IN) :: fmpi
1032 :
1033 : TYPE(t_sym), INTENT(IN) :: sym
1034 : TYPE(t_cell), INTENT(IN) :: cell
1035 : TYPE(t_noco), INTENT(IN) :: noco
1036 : TYPE(t_nococonv), INTENT(IN) :: nococonv
1037 : REAL, INTENT(IN) :: factor
1038 : LOGICAL, INTENT(IN) :: score, potnorm
1039 : TYPE(t_potden), INTENT(IN) :: denmat
1040 : CHARACTER(len=20), INTENT(IN) :: denName
1041 :
1042 0 : TYPE(t_potden) :: cden, mxden, myden, mzden
1043 :
1044 : CALL matrixsplit(sym,stars, atoms, sphhar, vacuum, input, noco,nococonv, factor, &
1045 0 : denmat, cden, mxden, myden, mzden)
1046 :
1047 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, noco, nococonv,&
1048 0 : score, potnorm, denName, cden, mxden, myden, mzden)
1049 :
1050 0 : END SUBROUTINE matrixplot
1051 :
1052 0 : SUBROUTINE procplot(stars, atoms, sphhar, sliceplot,vacuum, input, fmpi , sym, cell, &
1053 : noco, nococonv,denmat, plot_const)
1054 :
1055 : ! According to iplot, we process which exact plots we make after we assured
1056 : ! that we do any. n-th digit (from the back) of iplot ==1 --> plot with
1057 : ! identifier n is done.
1058 : #ifdef CPP_MPI
1059 : USE mpi
1060 : #endif
1061 :
1062 : TYPE(t_stars), INTENT(IN) :: stars
1063 : TYPE(t_atoms), INTENT(IN) :: atoms
1064 : TYPE(t_sphhar), INTENT(IN) :: sphhar
1065 : TYPE(t_sliceplot), INTENT(IN) :: sliceplot
1066 : TYPE(t_vacuum), INTENT(IN) :: vacuum
1067 : TYPE(t_input), INTENT(IN) :: input
1068 :
1069 : TYPE(t_mpi), INTENT(IN) :: fmpi
1070 : TYPE(t_sym), INTENT(IN) :: sym
1071 : TYPE(t_cell), INTENT(IN) :: cell
1072 : TYPE(t_noco), INTENT(IN) :: noco
1073 : TYPE(t_nococonv), INTENT(IN) :: nococonv
1074 : TYPE(t_potden), INTENT(IN) :: denmat
1075 : INTEGER, INTENT(IN) :: plot_const
1076 :
1077 : INTEGER :: i
1078 : REAL :: factor
1079 : CHARACTER (len=20) :: denName
1080 : LOGICAL :: score, potnorm
1081 :
1082 : ! Plotting the input density matrix as n / n, m / n, mx, my, mz.
1083 : ! Plot identifier: PLOT_INPDEN = 1
1084 : ! Additive term for iplot: 2
1085 0 : IF (plot_const.EQ.1) THEN
1086 0 : factor = 1.0
1087 0 : IF(sliceplot%slice) denName='slice'
1088 0 : IF(.NOT.sliceplot%slice) denName = 'denIn'
1089 0 : score = .FALSE.
1090 0 : potnorm = .FALSE.
1091 0 : IF (input%jspins.EQ.2) THEN
1092 0 : IF (noco%l_noco) THEN
1093 : CALL matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input,fmpi , sym, &
1094 : cell, noco, nococonv,factor, score, potnorm, denmat, &
1095 0 : denName)
1096 : ELSE
1097 : CALL vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1098 : cell, noco, nococonv,factor, score, potnorm, denmat, &
1099 0 : denName)
1100 : END IF
1101 : ELSE
1102 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, &
1103 0 : noco, nococonv,score, potnorm, denName, denmat)
1104 : END IF
1105 : END IF
1106 :
1107 : ! Plotting the output density matrix as n / n, m / n, mx, my, mz.
1108 : ! Plot identifier: PLOT_OUTDEN_Y_CORE = 5
1109 : ! No core subtraction done!
1110 : ! Additive term for iplot: 32
1111 0 : IF (plot_const.EQ.5) THEN
1112 0 : factor = 1.0
1113 0 : denName = 'denOutWithCore'
1114 0 : score = .FALSE.
1115 0 : potnorm = .FALSE.
1116 0 : IF (input%jspins.EQ.2) THEN
1117 0 : IF (noco%l_noco) THEN
1118 : CALL matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1119 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1120 0 : denName)
1121 : ELSE
1122 : CALL vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1123 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1124 0 : denName)
1125 : END IF
1126 : ELSE
1127 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, &
1128 0 : noco, nococonv, score, potnorm, denName, denmat)
1129 : END IF
1130 : END IF
1131 :
1132 : !! Plotting the output density matrix as n / n, m / n, mx, my, mz.
1133 : !! Plot identifier: PLOT_OUTDEN_N_CORE = 3
1134 : !! Core subtraction done!
1135 : !! Additive term for iplot: 8
1136 : !IF (plot_const.EQ.3) THEN
1137 : ! factor = 1.0
1138 : ! denName = 'denOutNOCore'
1139 : ! score = .TRUE.
1140 : ! potnorm = .FALSE.
1141 : ! IF (input%jspins.EQ.2) THEN
1142 : ! IF (noco%l_noco) THEN
1143 : ! CALL matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1144 : ! cell, noco, nococonv, factor, score, potnorm, denmat, &
1145 : ! denName)
1146 : ! ELSE
1147 : ! CALL vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1148 : ! cell, noco, nococonv, factor, score, potnorm, denmat, &
1149 : ! denName)
1150 : ! END IF
1151 : ! ELSE
1152 : ! CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, &
1153 : ! noco, nococonv, score, potnorm, denName, denmat)
1154 : ! END IF
1155 : !END IF
1156 :
1157 : ! Plotting the mixed density matrix as n / n, m / n, mx, my, mz.
1158 : ! Plot identifier: PLOT_MIXDEN_Y_CORE = 6
1159 : ! No core subtraction done!
1160 : ! Additive term for iplot: 128
1161 0 : IF (plot_const.EQ.6) THEN
1162 0 : factor = 1.0
1163 0 : denName = 'denOutMixWithCore'
1164 0 : score = .FALSE.
1165 0 : potnorm = .FALSE.
1166 0 : IF (input%jspins.EQ.2) THEN
1167 0 : IF (noco%l_noco) THEN
1168 : CALL matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1169 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1170 0 : denName)
1171 : ELSE
1172 : CALL vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1173 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1174 0 : denName)
1175 : END IF
1176 : ELSE
1177 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input,fmpi , sym, cell, &
1178 0 : noco, nococonv, score, potnorm, denName, denmat)
1179 : END IF
1180 : END IF
1181 :
1182 : !! Plotting the mixed density matrix as n / n, m / n, mx, my, mz.
1183 : !! Plot identifier: PLOT_MIXDEN_N_CORE = 5
1184 : !! Core subtraction done!
1185 : !! Additive term for iplot: 32
1186 : !IF (plot_const.EQ.5) THEN
1187 : ! factor = 1.0
1188 : ! denName = 'denOutMixNoCore'
1189 : ! score = .TRUE.
1190 : ! potnorm = .FALSE.
1191 : ! IF (input%jspins.EQ.2) THEN
1192 : ! IF (noco%l_noco) THEN
1193 : ! CALL matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input,fmpi, sym, &
1194 : ! cell, noco, nococonv, factor, score, potnorm, denmat, &
1195 : ! denName)
1196 : ! ELSE
1197 : ! CALL vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1198 : ! cell, noco, nococonv, factor, score, potnorm, denmat, &
1199 : ! denName)
1200 : ! END IF
1201 : ! ELSE
1202 : ! CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input,fmpi , sym, cell, &
1203 : ! noco, nococonv, score, potnorm, denName, denmat)
1204 : ! END IF
1205 : !END IF
1206 :
1207 :
1208 : ! Plotting the total potential as vTot / v_eff, B_eff / v_eff,
1209 : ! B_xc_1, B_xc_2, B_xc_3.
1210 : ! Plot identifier: PLOT_POT_TOT = 2
1211 : ! No core subtraction done!
1212 : ! Additive term for iplot: 4
1213 0 : IF (plot_const.EQ.2) THEN
1214 0 : IF(any(noco%l_alignMT)) CALL juDFT_warn("l_RelaxMT=T and plotting potentials can lead to wrong potentials visualized inside the MT",calledby="plot.f90")
1215 0 : factor = 2.0
1216 0 : denName = 'vTot'
1217 0 : score = .FALSE.
1218 0 : potnorm = .TRUE.
1219 0 : IF (input%jspins.EQ.2) THEN
1220 0 : IF (noco%l_noco) THEN
1221 : CALL matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1222 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1223 0 : denName)
1224 : ELSE
1225 : CALL vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1226 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1227 0 : denName)
1228 : END IF
1229 : ELSE
1230 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input,fmpi , sym, cell, &
1231 0 : noco, nococonv, score, potnorm, denName, denmat)
1232 : END IF
1233 : END IF
1234 :
1235 : ! Plotting the Coulomb potential as vCoul.
1236 : ! Plot identifier: PLOT_POT_COU = 3
1237 : ! No core subtraction done!
1238 : ! Additive term for iplot: 8
1239 0 : IF (plot_const.EQ.3) THEN
1240 0 : IF(any(noco%l_alignMT)) CALL juDFT_warn("l_alignMT=T and plotting potentials can lead to wrong potentials visualized inside the MT",calledby="plot.f90")
1241 0 : factor = 1.0
1242 0 : denName = 'vCoul'
1243 0 : score = .FALSE.
1244 0 : potnorm = .TRUE.
1245 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, cell, &
1246 0 : noco, nococonv, score, potnorm, denName, denmat)
1247 : END IF
1248 :
1249 : ! Plotting the xc potential as vXc / v_Xc, B_Xc / v_Xc,
1250 : ! B_xc_1, B_xc_2, B_xc_3.
1251 : ! Plot identifier: PLOT_POT_VXC = 4
1252 : ! No core subtraction done!
1253 : ! Additive term for iplot: 16
1254 0 : IF (plot_const.EQ.4) THEN
1255 0 : IF(any(noco%l_alignMT)) CALL juDFT_warn("l_alignMT=T and plotting potentials can lead to wrong potentials visualized inside the MT",calledby="plot.f90")
1256 0 : factor = 2.0
1257 0 : denName = 'vXc'
1258 0 : score = .FALSE.
1259 0 : potnorm = .TRUE.
1260 0 : IF (input%jspins.EQ.2) THEN
1261 0 : IF (noco%l_noco) THEN
1262 : CALL matrixplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1263 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1264 0 : denName)
1265 : ELSE
1266 : CALL vectorplot(sliceplot,stars, atoms, sphhar, vacuum, input, fmpi , sym, &
1267 : cell, noco, nococonv, factor, score, potnorm, denmat, &
1268 0 : denName)
1269 : END IF
1270 : ELSE
1271 : CALL savxsf(sliceplot,stars, atoms, sphhar, vacuum, input,fmpi , sym, cell, &
1272 0 : noco, nococonv, score, potnorm, denName, denmat)
1273 : END IF
1274 : END IF
1275 :
1276 0 : END SUBROUTINE procplot
1277 :
1278 0 : SUBROUTINE makeplots(stars, atoms, sphhar, vacuum, input, fmpi, sym, cell, &
1279 : noco, nococonv,denmat, plot_const, sliceplot)
1280 : USE m_Relaxspinaxismagn
1281 : ! Checks, based on the iplot switch that is given in the input, whether or
1282 : ! not plots should be made. Before the plot command is processed, we check
1283 : ! whether the plot_inp is there or an oldform is given. Both are outdated.
1284 : ! If that is not the case, we start plotting.
1285 : #ifdef CPP_MPI
1286 : USE mpi
1287 : #endif
1288 :
1289 : TYPE(t_stars), INTENT(IN) :: stars
1290 : TYPE(t_atoms), INTENT(IN) :: atoms
1291 : TYPE(t_sphhar), INTENT(IN) :: sphhar
1292 : TYPE(t_vacuum), INTENT(IN) :: vacuum
1293 : TYPE(t_input), INTENT(IN) :: input
1294 : TYPE(t_mpi), INTENT(IN) :: fmpi
1295 :
1296 : TYPE(t_sym), INTENT(IN) :: sym
1297 : TYPE(t_cell), INTENT(IN) :: cell
1298 : TYPE(t_noco), INTENT(IN) :: noco
1299 : TYPE(t_nococonv), INTENT(INOUT) :: nococonv
1300 : TYPE(t_potden), INTENT(INOUT) :: denmat
1301 : INTEGER, INTENT(IN) :: plot_const
1302 : TYPE(t_sliceplot), INTENT(IN) :: sliceplot
1303 : LOGICAL :: allowplot
1304 : INTEGER :: ierr
1305 : #ifdef CPP_MPI
1306 0 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr)
1307 : #endif
1308 : ! The check is done via bitwise operations. If the i-th position of iplot
1309 : ! in binary representation (2^n == n-th position) has a 1, the correspon-
1310 : ! ding plot with number 2^n is plotted.
1311 : !
1312 : ! E.g.: If the plots with identifying constants 1,2 and 4 are to be plotted
1313 : ! and none else, iplot would need to be 2^1 + 2^2 + 2^3 = 2 + 4 + 8 = 14.
1314 : ! iplot=1 or any odd number will *always* plot all possible options.
1315 :
1316 0 : CALL timestart("Plotting iplot plots")
1317 0 : allowplot=BTEST(sliceplot%iplot,plot_const).OR.(MODULO(sliceplot%iplot,2).EQ.1)
1318 : IF (allowplot) THEN
1319 0 : CALL toGlobalSpinFrame(noco, nococonv, vacuum, sphhar, stars, sym, cell, input, atoms, Denmat,fmpi,.true.)
1320 0 : CALL checkplotinp(fmpi)
1321 : CALL procplot(stars, atoms, sphhar,sliceplot, vacuum, input,fmpi, sym, cell, &
1322 0 : noco, nococonv, denmat, plot_const)
1323 0 : CALL toLocalSpinFrame(fmpi,vacuum, sphhar, stars, sym, cell, noco, nococonv, input, atoms,.false., denmat,.true.)
1324 : END IF
1325 0 : CALL timestop("Plotting iplot plots")
1326 :
1327 0 : END SUBROUTINE makeplots
1328 :
1329 0 : SUBROUTINE getMTSphere(input, cell, atoms, point, iType, iAtom, pt)
1330 :
1331 : ! Old subroutine originally from plotdop.f90, which is needed in savxsf.
1332 :
1333 : TYPE(t_input), INTENT(IN) :: input
1334 : TYPE(t_cell), INTENT(IN) :: cell
1335 : TYPE(t_atoms), INTENT(IN) :: atoms
1336 :
1337 : INTEGER, INTENT(OUT) :: iType, iAtom
1338 : REAL, INTENT(OUT) :: pt(3)
1339 : REAL, INTENT(IN) :: point(3)
1340 :
1341 : INTEGER :: ii1, ii2, ii3, i1, i2, i3, nq
1342 : REAL :: s
1343 :
1344 0 : ii1 = 3
1345 0 : ii2 = 3
1346 0 : ii3 = 3
1347 0 : IF (input%film ) ii3 = 0
1348 :
1349 :
1350 0 : DO i1 = -ii1, ii1
1351 0 : DO i2 = -ii2, ii2
1352 0 : DO i3 = -ii3, ii3
1353 0 : pt = point+MATMUL(cell%amat,(/i1,i2,i3/))
1354 0 : iAtom = 0
1355 0 : DO iType = 1, atoms%ntype
1356 0 : DO nq = 1, atoms%neq(iType)
1357 0 : iAtom = iAtom + 1
1358 0 : s = SQRT(DOT_PRODUCT(atoms%pos(:,iAtom)-pt,atoms%pos(:,iAtom)-pt))
1359 0 : IF (s<atoms%rmsh(atoms%jri(iType),iType)) THEN
1360 : ! Return with the current iType, iAtom, pt
1361 0 : RETURN
1362 : END IF
1363 : END DO
1364 : END DO
1365 : END DO
1366 : END DO
1367 : END DO !i1
1368 :
1369 : ! If no MT sphere encloses the point return 0 for iType, iAtom
1370 0 : iType = 0
1371 0 : iAtom = 0
1372 0 : pt(:) = point(:)
1373 :
1374 : END SUBROUTINE getMTSphere
1375 :
1376 : END MODULE m_plot
|
