Line data Source code
1 : !!****m* ABINIT/m_pptools
2 : !! NAME
3 : !! m_pptools
4 : !!
5 : !! FUNCTION
6 : !! Helper functions used for post-processing.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2002-2026 ABINIT group (MG, ZL, MJV, BXu)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_pptools
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_krank
29 :
30 : use m_io_tools, only : open_file
31 : use m_fstrings, only : sjoin, itoa
32 : use m_numeric_tools, only : wrap2_pmhalf
33 :
34 : implicit none
35 :
36 : private
37 :
38 : public :: prmat ! print real(dp) matrices in an attractive format.
39 : public :: write_xsf ! Open file and Write a generic array in the XSF format (XCrysden format)
40 : public :: printxsf ! Write a generic array in the XSF format (XCrysden format)
41 : public :: print_fofr_ri ! Print the [real, imaginary] part of an array
42 : public :: print_fofr_xyzri ! Print the Cartesian coordinates and the [real,imaginary] part of an array
43 : public :: print_fofr_cube ! Print ||fofr|| in CUBE format.
44 : public :: printbxsf ! Print band structure energies in XCrysDen format.
45 : public :: printvtk ! Print band structure energies and velocities in VTK format.
46 :
47 : CONTAINS !===========================================================
48 : !!***
49 :
50 : !!****f* m_pptools/prmat
51 : !! NAME
52 : !! prmat
53 : !!
54 : !! FUNCTION
55 : !! This subroutine prints real*8 matrices in an attractive format.
56 : !!
57 : !! INPUTS
58 : !! mat(mi,nj)= matrix to be printed
59 : !! mi = no rows of mat
60 : !! ni = no rows to print
61 : !! nj = no colums of mat
62 : !! unitm = unit to print to, if not provided std_out is chosen
63 : !!
64 : !! OUTPUT
65 : !! (only writing)
66 : !!
67 : !! SOURCE
68 :
69 624 : subroutine prmat(mat, ni, nj, mi, unitm)
70 :
71 : !Arguments ------------------------------------
72 : !scalars
73 : integer,intent(in) :: mi,ni,nj
74 : integer,intent(in), optional :: unitm
75 : !arrays
76 : real(dp),intent(in) :: mat(mi,nj)
77 :
78 : !Local variables-------------------------------
79 : !scalars
80 : character(len=1000) :: msg
81 : integer,parameter :: nline=10
82 : integer :: ii,jj,jstart,jstop,unitn
83 : ! *************************************************************************
84 :
85 624 : unitn = std_out; if (present(unitm)) unitn = unitm
86 :
87 1248 : do jstart = 1, nj, nline
88 624 : jstop = min(nj, jstart+nline-1)
89 2452 : write(msg, '(3x,10(i4,8x))' ) (jj,jj=jstart,jstop)
90 1248 : call wrtout(unitn,msg)
91 : end do
92 :
93 1872 : do ii = 1,ni
94 3120 : do jstart= 1, nj, nline
95 1248 : jstop = min(nj, jstart+nline-1)
96 2496 : if (jstart==1) then
97 1248 : write(msg, '(i3,1p,10e12.4)' ) ii, (mat(ii,jj),jj=jstart,jstop)
98 1248 : call wrtout(unitn,msg)
99 : else
100 0 : write(msg, '(3x,1p,10e12.4)' ) (mat(ii,jj),jj=jstart,jstop)
101 0 : call wrtout(unitn,msg)
102 : end if
103 : end do
104 : end do
105 :
106 624 : end subroutine prmat
107 : !!***
108 :
109 : !----------------------------------------------------------------------
110 :
111 : !!****f* m_pptools/write_xsf
112 : !! NAME
113 : !! write_xsf
114 : !!
115 : !! FUNCTION
116 : !! Open file and write a generic array in the XSF format (XCrysden format)
117 : !! See printxsf for the meaning of the arguments.
118 : !!
119 : !! OUTPUT
120 : !! Only write
121 : !!
122 : !! SOURCE
123 :
124 8 : subroutine write_xsf(filepath, n1, n2, n3, datagrid, basis, origin, natom, ntypat, typat, xcart, znucl, realrecip)
125 :
126 : !Arguments ------------------------------------
127 : !scalars
128 : character(len=*),intent(in) :: filepath
129 : integer,intent(in) :: n1, n2, n3, realrecip, natom, ntypat
130 : !arrays
131 : integer,intent(in) :: typat(natom)
132 : real(dp),intent(in) :: basis(3,3), datagrid(n1*n2*n3), origin(3), xcart(3,natom), znucl(ntypat)
133 :
134 : !Local variables-------------------------------
135 : integer :: ount
136 : character(len=500) :: msg
137 : ! *************************************************************************
138 :
139 8 : if (open_file(filepath, msg, newunit=ount, form='formatted', status='unknown', action="write") /= 0) then
140 0 : ABI_ERROR(msg)
141 : end if
142 :
143 8 : call printxsf(n1, n2, n3, datagrid, basis, origin, natom, ntypat, typat, xcart, znucl, ount, realrecip)
144 8 : close(ount)
145 :
146 8 : end subroutine write_xsf
147 : !!***
148 :
149 : !----------------------------------------------------------------------
150 :
151 : !!****f* m_pptools/printxsf
152 : !! NAME
153 : !! printxsf
154 : !!
155 : !! FUNCTION
156 : !! Write a generic array in the XSF format (XCrysden format)
157 : !!
158 : !! INPUTS
159 : !! n1, n2, n3=grid size along x, y, z
160 : !! datagrid(n1*n2*n3) = datagrid values stored using the fortran convention
161 : !! basis(3,3) = basis vectors of the direct real lattice or of the reciprocal lattice (fortran convention)
162 : !! (Bohr units if realrecip=0, Bohr^-1 if realrecip=1, see below)
163 : !! origin(3) = origin of the grid
164 : !! ount = unit number of the output file (already open by the caller, not closed here!)
165 : !! realrecip = 0 for a plot in real space
166 : !! 1 for a plot in reciprocal space
167 : !!
168 : !! OUTPUT
169 : !! Only write
170 : !!
171 : !! SOURCE
172 :
173 20 : subroutine printxsf(n1, n2, n3, datagrid, basis, origin, natom, ntypat, typat, xcart, znucl, ount, realrecip)
174 :
175 : !Arguments ------------------------------------
176 : !scalars
177 : integer,intent(in) :: n1, n2, n3, ount, realrecip, natom, ntypat
178 : !arrays
179 : integer,intent(in) :: typat(natom)
180 : real(dp),intent(in) :: basis(3,3), datagrid(n1*n2*n3), origin(3), xcart(3,natom), znucl(ntypat)
181 :
182 : !Local variables-------------------------------
183 : !scalars
184 : integer :: ix,iy,iz,nslice,nsym,iatom
185 : real(dp) :: fact
186 : !arrays
187 40 : real(dp) :: tau(3,natom)
188 : ! *************************************************************************
189 :
190 20 : if (all(realrecip /= [0, 1])) then
191 0 : ABI_BUG(sjoin('The argument realrecip should be 0 or 1, received:', itoa(realrecip)))
192 : end if
193 :
194 : ! conversion between ABINIT default units and XCrysden units
195 20 : fact = Bohr_Ang; if (realrecip == 1) fact=one/fact ! since we are in reciprocal space
196 :
197 : ! TODO insert crystalline structure and dummy atoms in case of reciprocal space need to convert basis too
198 20 : write(ount,'(1X,A)') 'DIM-GROUP'
199 20 : write(ount,*) '3 1'
200 20 : write(ount,'(1X,A)') 'PRIMVEC'
201 80 : do iy = 1,3
202 260 : write(ount,'(3(ES17.10,2X))') (Bohr_Ang*basis(ix,iy), ix=1,3)
203 : end do
204 :
205 : ! generate translated coordinates to fit origin shift
206 9144 : do iatom = 1,natom
207 36516 : tau(:,iatom) = xcart(:,iatom) - origin(:)
208 : end do
209 :
210 20 : write(ount,'(1X,A)') 'PRIMCOORD'
211 20 : write(ount,*) natom, ' 1'
212 9144 : do iatom = 1,natom
213 9124 : write(ount,'(i9,3(3X,ES17.10))') NINT(znucl(typat(iatom))), & ! WARNING alchemy not supported by XCrysden
214 18268 : Bohr_Ang*tau(1,iatom), Bohr_Ang*tau(2,iatom), Bohr_Ang*tau(3,iatom)
215 : end do
216 20 : write(ount,'(1X,A)') 'ATOMS'
217 9144 : do iatom = 1,natom
218 9124 : write(ount,'(i9,3(3X,ES17.10))') NINT(znucl(typat(iatom))), & ! WARNING alchemy not supported by XCrysden
219 18268 : Bohr_Ang*tau(1,iatom), Bohr_Ang*tau(2,iatom), Bohr_Ang*tau(3,iatom)
220 : end do
221 :
222 20 : write(ount,'(a)')' BEGIN_BLOCK_DATAGRID3D'
223 20 : write(ount,'(a)')' datagrid'
224 20 : write(ount,'(a)')' DATAGRID_3D_DENSITY'
225 : ! NOTE: XCrysden uses aperiodical data grid
226 20 : write(ount,*)n1+1,n2+1,n3+1
227 20 : write(ount,*)origin
228 80 : write(ount,*)basis(:,1)*fact
229 80 : write(ount,*)basis(:,2)*fact
230 80 : write(ount,*)basis(:,3)*fact
231 :
232 20 : nslice=1
233 626 : do iz=1,n3
234 27624 : do iy=1,n2
235 27018 : write(ount,'(8es16.8)') datagrid(1+n1*(nslice-1):n1+n1*(nslice-1)),datagrid(1+n1*(nslice-1))
236 27624 : nslice = nslice+1
237 : end do
238 606 : nsym=nslice-n2
239 626 : write(ount,'(8es16.8)') datagrid(1+n1*(nsym-1):n1+n1*(nsym-1)),datagrid(1+n1*(nsym-1))
240 : end do
241 :
242 : ! Now write upper plane
243 : nslice = 1
244 626 : do iy=1,n2
245 606 : write(ount,'(8es16.8)') datagrid(1+n1*(nslice-1):n1+n1*(nslice-1)),datagrid(1+n1*(nslice-1))
246 626 : nslice=nslice+1
247 : end do
248 :
249 20 : nsym = nslice-n2
250 20 : write(ount,'(8es16.8)') datagrid(1+n1*(nsym-1):n1+n1*(nsym-1)),datagrid(1+n1*(nsym-1))
251 :
252 20 : write(ount,'(a)')' END_DATAGRID_3D'
253 20 : write(ount,'(a)')' END_BLOCK_DATAGRID3D'
254 :
255 20 : end subroutine printxsf
256 : !!***
257 :
258 : !----------------------------------------------------------------------
259 :
260 : !!****f* m_pptools/print_fofr_ri
261 : !! NAME
262 : !! print_fofr_ri
263 : !!
264 : !! FUNCTION
265 : !! Print the [real,imaginary] part of fofr on unit unit
266 : !!
267 : !! INPUTS
268 : !! ri_mode =
269 : !! "RI" if both real and imag part are wanted
270 : !! "R" for real part
271 : !! "I" for imaginary part
272 : !! nx,ny,nz,ldx,ldy,ldz = Logical and physical dimensions of the array.
273 : !! fofr(2,ldx,ldy,ldz) = Input data
274 : !! [unit] = Fortran unit number. Default: std_out
275 : !!
276 : !! OUTPUT
277 : !! Only writing
278 : !!
279 : !! SOURCE
280 :
281 2 : subroutine print_fofr_ri(ri_mode,nx,ny,nz,ldx,ldy,ldz,fofr,unit)
282 :
283 : !Arguments -----------------------------------------------
284 : !scalars
285 : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz
286 : integer,optional,intent(in) :: unit
287 : character(len=*),intent(in) :: ri_mode
288 : !arrays
289 : real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
290 :
291 : !Local variables-------------------------------
292 : integer :: ount,ix,iy,iz
293 : ! *************************************************************************
294 :
295 2 : ount = std_out; if (PRESENT(unit)) ount = unit
296 :
297 1 : SELECT CASE (ri_mode)
298 : CASE ("RI","ri")
299 13 : do iz=1,nz
300 157 : do iy=1,ny
301 1884 : do ix=1,nx
302 1872 : write(ount,'(2f20.16)') fofr(:,ix,iy,iz)
303 : end do
304 : end do
305 : end do
306 :
307 : CASE ("R","r")
308 25 : do iz=1,nz
309 409 : do iy=1,ny
310 4248 : do ix=1,nx
311 4224 : write(ount,'(f20.16)') fofr(1,ix,iy,iz)
312 : end do
313 : end do
314 : end do
315 :
316 : CASE ("I","i")
317 0 : do iz=1,nz
318 0 : do iy=1,ny
319 0 : do ix=1,nx
320 0 : write(ount,'(f20.16)') fofr(2,ix,iy,iz)
321 : end do
322 : end do
323 : end do
324 :
325 : CASE DEFAULT
326 2 : ABI_ERROR("Wrong ri_mode")
327 : END SELECT
328 :
329 2 : end subroutine print_fofr_ri
330 : !!***
331 :
332 : !----------------------------------------------------------------------
333 :
334 : !!****f* m_pptools/print_fofr_xyzri
335 : !! NAME
336 : !! print_fofr_xyzri
337 : !!
338 : !! FUNCTION
339 : !! Print the Cartesian coordinates and the [real,imaginary] part of fofr on unit unit
340 : !!
341 : !! INPUTS
342 : !! ri_mode =
343 : !! "RI" if both real and imag part are wanted
344 : !! "R" for real part
345 : !! "I" for imaginary part
346 : !! nx,ny,nz,ldx,ldy,ldz = Logical and physical dimensions of the array.
347 : !! fofr(2,ldx,ldy,ldz) = Input data
348 : !! rprimd(3,3)=Lattive vectors in Bohr
349 : !! [conv_fact] = Conversion factor for rprimd (rprimd is multiplied by conv_fact). Default is one
350 : !! [unit] = Fortran unit number. Default: std_out
351 : !!
352 : !! OUTPUT
353 : !! Only writing
354 : !!
355 : !! SOURCE
356 :
357 1 : subroutine print_fofr_xyzri(ri_mode,nx,ny,nz,ldx,ldy,ldz,fofr,rprimd,conv_fact,unit)
358 :
359 : !Arguments -----------------------------------------------
360 : !scalars
361 : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz
362 : integer,optional,intent(in) :: unit
363 : real(dp),optional,intent(in) :: conv_fact
364 : character(len=*),intent(in) :: ri_mode
365 : !arrays
366 : real(dp),intent(in) :: rprimd(3,3)
367 : real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
368 :
369 : !Local variables-------------------------------
370 : !scalars
371 : integer :: ount,ix,iy,iz
372 : real(dp) :: xnow,ynow,znow,my_cfact
373 : ! *************************************************************************
374 :
375 1 : ount = std_out; if (PRESENT(unit)) ount = unit
376 1 : my_cfact = one; if (PRESENT(conv_fact)) my_cfact = conv_fact
377 :
378 1 : SELECT CASE (ri_mode)
379 : CASE ("RI","ri")
380 13 : do iz=1,nz
381 157 : do iy=1,ny
382 1884 : do ix=1,nx
383 1728 : xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
384 1728 : ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
385 1728 : znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
386 1872 : write(ount,'(3f16.10,2f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(:,ix,iy,iz)
387 : end do
388 : end do
389 : end do
390 :
391 : CASE ("R","r")
392 0 : do iz=1,nz
393 0 : do iy=1,ny
394 0 : do ix=1,nx
395 0 : xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
396 0 : ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
397 0 : znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
398 0 : write(ount,'(3f16.10,f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(1,ix,iy,iz)
399 : end do
400 : end do
401 : end do
402 :
403 : CASE ("I","i")
404 0 : do iz=1,nz
405 0 : do iy=1,ny
406 0 : do ix=1,nx
407 0 : xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
408 0 : ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
409 0 : znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
410 0 : write(ount,'(3f16.10,f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(2,ix,iy,iz)
411 : end do
412 : end do
413 : end do
414 :
415 : CASE DEFAULT
416 1 : ABI_ERROR("Wrong ri_mode")
417 : END SELECT
418 :
419 1 : end subroutine print_fofr_xyzri
420 : !!***
421 :
422 : !----------------------------------------------------------------------
423 :
424 : !!****f* m_pptools/print_fofr_cube
425 : !! NAME
426 : !! print_fofr_cube
427 : !!
428 : !! FUNCTION
429 : !! Print array fofr in the cube file format
430 : !!
431 : !! INPUTS
432 : !! nx,ny,nz,ldx,ldy,ldz = Logical and physical dimensions of the array.
433 : !! fofr(2,ldx,ldy,ldz) = Input data
434 : !! rprimd(3,3)=Lattive vectors in Bohr
435 : !! [unit] = Fortran unit number. Default: std_out
436 : !!
437 : !! OUTPUT
438 : !! Only writing
439 : !!
440 : !! SOURCE
441 :
442 :
443 1 : subroutine print_fofr_cube(nx,ny,nz,ldx,ldy,ldz,fofr,rprimd,natom,znucl_atom,xcart,unit)
444 :
445 : !Arguments -----------------------------------------------
446 : !scalars
447 : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,natom
448 : integer,optional,intent(in) :: unit
449 : !arrays
450 : integer,intent(in) :: znucl_atom(natom)
451 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
452 : real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
453 :
454 : !Local variables-------------------------------
455 : !scalars
456 : integer,parameter :: cplx=2
457 : integer :: ount,ix,iy,iz,iatom
458 : ! *************************************************************************
459 :
460 1 : ount = std_out; if (PRESENT(unit)) ount = unit
461 :
462 : ! EXAMPLE FROM THE WEB
463 : ! CPMD CUBE FILE.
464 : ! OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z
465 : ! 3 0.000000 0.000000 0.000000
466 : ! 40 0.283459 0.000000 0.000000
467 : ! 40 0.000000 0.283459 0.000000
468 : ! 40 0.000000 0.000000 0.283459
469 : ! 8 0.000000 5.570575 5.669178 5.593517
470 : ! 1 0.000000 5.562867 5.669178 7.428055
471 : ! 1 0.000000 7.340606 5.669178 5.111259
472 : ! -0.25568E-04 0.59213E-05 0.81068E-05 0.10868E-04 0.11313E-04 0.35999E-05
473 :
474 1 : write(ount,'(a)') 'ABINIT generated cube file'
475 1 : write(ount,'(a)') 'from cut3d tool'
476 :
477 1 : write(ount,'(i9,3(1x,f12.6))') natom,0.,0.,0.
478 4 : write(ount,'(i9,3(1x,f12.6))') nx,(rprimd(iy,1)/nx, iy=1,3)
479 4 : write(ount,'(i9,3(1x,f12.6))') ny,(rprimd(iy,2)/ny, iy=1,3)
480 4 : write(ount,'(i9,3(1x,f12.6))') nz,(rprimd(iy,3)/nz, iy=1,3)
481 :
482 2 : do iatom=1,natom
483 2 : write(ount,'(i9,4(3X,ES17.10))') znucl_atom(iatom),0.d0,xcart(1:3,iatom)
484 : end do
485 :
486 : ! Note C ordering of the indexes
487 : if (cplx==2) then
488 17 : do ix=1,nx
489 273 : do iy=1,ny
490 4368 : do iz=1,nz
491 4352 : write(ount,'(6(f12.6,2x))') sqrt(fofr(1,ix,iy,iz)**2 + fofr(2,ix,iy,iz)**2 )
492 : end do
493 : end do
494 : end do
495 : else
496 : do ix=1,nx
497 : do iy=1,ny
498 : do iz=1,nz
499 : write(ount,'(6(f12.6,2x))') fofr(1,ix,iy,iz)
500 : end do
501 : end do
502 : end do
503 : end if
504 :
505 1 : end subroutine print_fofr_cube
506 : !!***
507 :
508 : !----------------------------------------------------------------------
509 :
510 : !!****f* m_pptools/printbxsf
511 : !! NAME
512 : !! printbxsf
513 : !!
514 : !! FUNCTION
515 : !! Print band structure energies in XCrysDen format.
516 : !!
517 : !! INPUTS
518 : !! eigen(mband,nkpt,nsppol) = eigenvalues in hartree
519 : !! ewind = energy window around the fermi level.
520 : !! if ewind /= 0 ==> a band is considered in the plot of FSurf
521 : !! only if it is inside [ ef-ewind, ef+ewind ] for some k point
522 : !! if ewind == 0 ==> all bands will be keept in the _BXSF file
523 : !! fermie = Fermi energy (Hartree)
524 : !! gprimd(3,3) = dimensional primitive translations for reciprocal space (bohr^-1)
525 : !! kptrlatt(3,3) = reciprocal of lattice vectors for full kpoint grid
526 : !! mband = maximum number of bands
527 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
528 : !! shiftk(3,nshiftk) =shift vector for k point grid
529 : !! fname = filename for the fortran file
530 : !! symafm(nsym)=(Anti)ferromagnetic symmetries.
531 : !! use_afm=.TRUE. if (anti)ferromagnetic symmetries are used.
532 : !!
533 : !! OUTPUT
534 : !! ierr=Status error.
535 : !! BXSF file.
536 : !!
537 : !! SOURCE
538 :
539 3 : subroutine printbxsf(eigen,ewind,fermie,gprimd,kptrlatt,mband,&
540 3 : nkptirred,kptirred,nsym,use_afm,symrec,symafm,use_tr,nsppol,shiftk,nshiftk,fname,ierr)
541 :
542 : !Arguments ------------------------------------
543 : !scalars
544 : integer,intent(in) :: mband,nkptirred,nshiftk,nsppol,nsym
545 : integer,intent(out) :: ierr
546 : real(dp),intent(in) :: ewind,fermie
547 : logical,intent(in) :: use_afm,use_tr
548 : character(len=*),intent(in) :: fname
549 : !arrays
550 : integer,intent(in) :: kptrlatt(3,3),symafm(nsym),symrec(3,3,nsym)
551 : real(dp),intent(in) :: eigen(mband,nkptirred,nsppol),gprimd(3,3)
552 : real(dp),intent(in) :: kptirred(3,nkptirred),shiftk(3,nshiftk)
553 :
554 : !Local variables-------------------------------
555 : !scalars
556 : integer,parameter :: enough = 50
557 : integer :: iband,ik1,ik2,ik3,ikgrid,ikpt,indx
558 : integer :: isppol,isym,maxband,minband,nk1,nk2,nk3,nkptfull,ubxsf,timrev
559 : integer :: symkptrank, nsymfm, isymfm
560 : real(dp) :: ene
561 : character(len=500) :: msg
562 3 : type(krank_t) :: krank
563 : !arrays
564 3 : integer,allocatable :: fulltoirred(:),symrecfm(:,:,:)
565 : real(dp) :: kptgrid(3),gmet(3,3)
566 : ! *************************************************************************
567 :
568 3 : ierr = 0
569 :
570 : ! Error if klatt is no simple orthogonal lattice (in red space)
571 : ! for generalization to MP grids, need new version of XCrysDen
572 : if (kptrlatt(1,2)/=0 .or. kptrlatt(1,3)/=0 .or. kptrlatt(2,1)/=0 .or. &
573 3 : kptrlatt(2,3)/=0 .or. kptrlatt(3,1)/=0 .or. kptrlatt(3,2)/=0 ) then
574 : write(msg,'(3a)')&
575 0 : 'kptrlatt should be diagonal, for the FS calculation ',ch10,&
576 0 : 'Action: use an orthogonal k-grid for the GS calculation '
577 0 : ABI_COMMENT(msg)
578 0 : ierr = ierr + 1
579 : end if
580 :
581 : ! Error if there are not at least 2 kpts in each direction:
582 : ! kptrank will fail for the intermediate points below
583 3 : if (abs(kptrlatt(1,1)) < 2 .or. abs(kptrlatt(2,2)) < 2 .or. abs(kptrlatt(3,3)) < 2) then
584 : write(msg,'(3a)')&
585 0 : 'You need at least 2 points in each direction in k space to output BXSF files ',ch10,&
586 0 : 'Action: use an augmented k-grid for the GS calculation (at least 2x2x2) '
587 0 : ABI_COMMENT(msg)
588 0 : ierr = ierr + 1
589 : end if
590 :
591 15 : if (ANY(ABS(shiftk) > tol10)) then
592 : write(msg,'(3a)')&
593 0 : 'Origin of the k-grid should be (0,0,0) for the FS calculation ',ch10,&
594 0 : 'Action: use a non-shifted k-grid for the GS calculation. Returning '
595 0 : ABI_COMMENT(msg)
596 0 : ierr = ierr + 1
597 : end if
598 :
599 3 : if (ierr /= 0) return
600 :
601 : ! Compute reciprocal space metric.
602 3 : gmet = MATMUL(TRANSPOSE(gprimd), gprimd)
603 :
604 3 : if (use_afm) then
605 0 : nsymfm = 0
606 0 : do isym = 1, nsym
607 0 : if (symafm(isym) == 1) nsymfm = nsymfm+1
608 : end do
609 0 : ABI_MALLOC(symrecfm, (3,3,nsymfm))
610 0 : isymfm = 0
611 0 : do isym = 1, nsym
612 0 : if (symafm(isym) == 1) then
613 0 : isymfm = isymfm + 1
614 0 : symrecfm(:,:,isymfm) = symrec(:,:,isym)
615 : end if
616 : end do
617 : else
618 3 : nsymfm = nsym
619 9 : ABI_MALLOC(symrecfm, (3,3,nsymfm))
620 1566 : symrecfm = symrec
621 : end if
622 :
623 : ! Xcrysden uses aperiodic data-grid (images are included in the grid)
624 3 : nk1 = kptrlatt(1,1); nk2 = kptrlatt(2,2); nk3 = kptrlatt(3,3)
625 3 : nkptfull = (nk1+1) * (nk2+1) * (nk3+1)
626 :
627 9 : ABI_MALLOC(fulltoirred, (nkptfull))
628 3 : timrev = 0; if (use_tr) timrev=1
629 :
630 3 : call krank%init(nkptirred, kptirred, nsym=nsymfm, symrec=symrecfm, time_reversal=use_tr)
631 :
632 : ! Xcrysden employs the C-ordering for the Fermi Surface (x-y-z)
633 3 : ikgrid=0
634 38 : do ik1=0,nk1
635 489 : do ik2=0,nk2
636 6857 : do ik3=0,nk3
637 :
638 6371 : ikgrid = ikgrid+1
639 6371 : kptgrid(1) = DBLE(ik1)/kptrlatt(1,1)
640 6371 : kptgrid(2) = DBLE(ik2)/kptrlatt(2,2)
641 6371 : kptgrid(3) = DBLE(ik3)/kptrlatt(3,3)
642 :
643 : ! Find correspondence between the Xcrysden grid and the IBZ
644 6371 : symkptrank = krank%get_rank(kptgrid)
645 6371 : fulltoirred(ikgrid) = krank%invrank(symkptrank)
646 :
647 6822 : if (fulltoirred(ikgrid) < 1) then
648 0 : if (ierr <= enough) then
649 : write(msg,'(a,3es16.8,2a,i0,2a)')&
650 0 : 'kpt = ',kptgrid,ch10,' with rank ', symkptrank, ch10,&
651 0 : 'has no symmetric among the k-points used in the GS calculation '
652 0 : ABI_WARNING(msg)
653 : end if
654 0 : ierr = ierr + 1
655 : end if
656 :
657 : end do !ik1
658 : end do !ik2
659 : end do !ik3
660 :
661 3 : call krank%free()
662 :
663 3 : ABI_CHECK(ierr == 0, "See above warnings")
664 :
665 3 : if (abs(ewind) < tol12 ) then
666 : ! Keep all bands.
667 3 : minband=1
668 3 : maxband=mband
669 : else
670 : ! Select a subset of bands.
671 0 : minband = mband
672 0 : maxband = 0
673 0 : ene=abs(ewind)
674 0 : do isppol=1,nsppol
675 0 : do iband=1,mband
676 0 : if(minval(eigen(iband,:,isppol))-fermie < -ene) minband = iband
677 : end do
678 0 : do iband=mband,1,-1
679 0 : if (maxval(eigen(iband,:,isppol))-fermie > ene) maxband = iband
680 : end do
681 : end do ! isppol
682 :
683 : end if ! abs(energy_window)
684 :
685 : ! Dump results to file
686 3 : if (open_file(fname,msg, newunit=ubxsf, status='unknown', action="write", form='formatted') /= 0 ) then
687 0 : ABI_WARNING(msg)
688 0 : ierr=ierr +1; RETURN
689 : end if
690 :
691 : ! Write header
692 3 : write(ubxsf,*)' BEGIN_INFO'
693 3 : write(ubxsf,*)' #'
694 3 : write(ubxsf,*)' # this is a Band-XCRYSDEN-Structure-File for Visualization of Fermi Surface'
695 3 : write(ubxsf,*)' # generated by the ABINIT package'
696 3 : write(ubxsf,*)' #'
697 3 : write(ubxsf,*)' # bands between ',minband,' and ',maxband
698 3 : write(ubxsf,*)' #'
699 3 : if (nsppol == 2 ) then
700 0 : write(ubxsf,*)' # NOTE: the first band is relative to spin-up electrons,'
701 0 : write(ubxsf,*)' # the second band to spin-down and so on .. '
702 0 : write(ubxsf,*)' #'
703 : end if
704 3 : write(ubxsf,*)' # Launch as: xcrysden --bxsf '
705 3 : write(ubxsf,*)' #'
706 3 : write(ubxsf,'(a,es16.8)')' Fermi Energy: ',fermie
707 3 : write(ubxsf,*)' END_INFO'
708 3 : write(ubxsf,*)' '
709 3 : write(ubxsf,*)' BEGIN_BLOCK_BANDGRID_3D'
710 3 : write(ubxsf,*)' band_energies'
711 3 : write(ubxsf,*)' BEGIN_BANDGRID_3D'
712 :
713 3 : write(ubxsf,*)' ',(maxband-minband+1)*nsppol
714 3 : write(ubxsf,*)' ',nk1+1,nk2+1,nk3+1
715 3 : write(ubxsf,*)' ',shiftk(:,1)
716 : ! Angstrom units are used in the BXSF format
717 12 : write(ubxsf,*)' ',gprimd(:,1)/Bohr_Ang
718 12 : write(ubxsf,*)' ',gprimd(:,2)/Bohr_Ang
719 12 : write(ubxsf,*)' ',gprimd(:,3)/Bohr_Ang
720 :
721 : ! print out data for all relevant bands and full kpt grid (redundant, yes)
722 : ! for each kpt in full zone, find equivalent irred kpt and print eigenval
723 3 : indx = 0
724 22 : do iband=minband,maxband
725 41 : do isppol=1,nsppol
726 19 : write(ubxsf,*)' BAND: ',indx+minband
727 26422 : write(ubxsf,'(7(es16.8))')(eigen(iband,fulltoirred(ikpt),isppol),ikpt=1,nkptfull)
728 38 : indx=indx+1
729 : end do
730 : end do
731 :
732 3 : write(ubxsf,*)' END_BANDGRID_3D'
733 3 : write(ubxsf,*)' END_BLOCK_BANDGRID_3D'
734 3 : close(ubxsf)
735 :
736 3 : ABI_FREE(fulltoirred)
737 3 : ABI_FREE(symrecfm)
738 :
739 3 : end subroutine printbxsf
740 : !!***
741 :
742 : !!****f* m_pptools/printvtk
743 : !! NAME
744 : !! printvtk
745 : !!
746 : !! FUNCTION
747 : !! Print band structure energies and velocities in VTK format.
748 : !!
749 : !! INPUTS
750 : !! eigen(mband,nkpt,nsppol) = eigenvalues in hartree
751 : !! ewind = energy window around the fermi level.
752 : !! if ewind /= 0 ==> a band is considered in the plot of FSurf
753 : !! only if it is inside [ ef-ewind, ef+ewind ] for some k point
754 : !! if ewind == 0 ==> all bands will be keept in the _BXSF file
755 : !! fermie = Fermi energy (Hartree)
756 : !! gprimd(3,3) = dimensional primitive translations for reciprocal space (bohr^-1)
757 : !! kptrlatt(3,3) = reciprocal of lattice vectors for full kpoint grid
758 : !! mband = maximum number of bands
759 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
760 : !! shiftk(3,nshiftk) =shift vector for k point grid
761 : !! fname = filename for the fortran file
762 : !! symafm(nsym)=(Anti)ferromagnetic symmetries.
763 : !! use_afm=.TRUE. if (anti)ferromagnetic symmetries are used.
764 : !!
765 : !! OUTPUT
766 : !! ierr=Status error.
767 : !! BXSF file.
768 : !!
769 : !! SOURCE
770 :
771 0 : subroutine printvtk(eigen,v_surf,ewind,fermie,gprimd,kptrlatt,mband,&
772 0 : & nkptirred,kptirred,nsym,use_afm,symrec,symafm,use_tr,nsppol,shiftk,nshiftk,fname,ierr)
773 :
774 :
775 : !Arguments ------------------------------------
776 : !scalars
777 : integer,intent(in) :: mband,nkptirred,nshiftk,nsppol,nsym
778 : integer,intent(out) :: ierr
779 : real(dp),intent(in) :: ewind,fermie
780 : logical,intent(in) :: use_afm,use_tr
781 : character(len=*),intent(in) :: fname
782 : !arrays
783 : integer,intent(in) :: kptrlatt(3,3),symafm(nsym),symrec(3,3,nsym)
784 : real(dp),intent(in) :: eigen(mband,nkptirred,nsppol),gprimd(3,3)
785 : real(dp),intent(in) :: v_surf(mband,kptrlatt(1,1)+1,kptrlatt(2,2)+1,kptrlatt(3,3)+1,3,nsppol)
786 : real(dp),intent(in) :: kptirred(3,nkptirred),shiftk(3,nshiftk)
787 :
788 : !Local variables-------------------------------
789 : !scalars
790 : integer :: iband,ikgrid,ikpt1,indx
791 : integer :: ikpt,jkpt,kkpt,ikpt_fine, ik1, ik2, ik3
792 : integer :: isppol,isym,itim,maxband,minband,nk1,nk2,nk3,nkptfull,uvtk,timrev
793 : real(dp) :: ene,res,ss,timsign
794 : logical :: found
795 : character(len=500) :: msg, format_str
796 : !arrays
797 0 : integer,allocatable :: fulltoirred(:)
798 : real(dp) :: kconv(3),kpt(3),kptgrid(3),kptsym(3)
799 : ! *************************************************************************
800 :
801 0 : ierr=0
802 :
803 : !Error if klatt is no simple orthogonal lattice (in red space)
804 : !for generalization to MP grids, need new version of XCrysDen
805 :
806 : if (kptrlatt(1,2)/=0 .or. kptrlatt(1,3)/=0 .or. kptrlatt(2,1)/=0 .or. &
807 0 : kptrlatt(2,3)/=0 .or. kptrlatt(3,1)/=0 .or. kptrlatt(3,2)/=0 ) then
808 : write(msg,'(3a)')&
809 0 : 'kptrlatt should be diagonal, for the FS calculation ',ch10,&
810 0 : 'action: use an orthogonal k-grid for the GS calculation '
811 0 : ABI_COMMENT(msg)
812 0 : ierr=ierr+1
813 : end if
814 :
815 0 : if (ANY(ABS(shiftk(:,:))>tol10)) then
816 : write(msg,'(3a)')&
817 0 : 'Origin of the k-grid should be (0,0,0) for the FS calculation ',ch10,&
818 0 : 'Action: use a non-shifted k-grid for the GS calculation. Returning '
819 0 : ABI_COMMENT(msg)
820 0 : ierr=ierr+1
821 : end if
822 :
823 0 : if (ierr/=0) RETURN
824 :
825 : ! Xcrysden uses aperiodical data-grid
826 0 : nk1 = kptrlatt(1,1)
827 0 : nk2 = kptrlatt(2,2)
828 0 : nk3 = kptrlatt(3,3)
829 0 : nkptfull=(nk1+1)*(nk2+1)*(nk3+1)
830 :
831 0 : ABI_MALLOC(fulltoirred,(nkptfull))
832 0 : timrev=0; if (use_tr) timrev=1
833 :
834 : !Xcrysden employs C-ordering for the Fermi Surface.
835 0 : ierr = 0
836 0 : ikgrid=0
837 0 : do ik1=0,nk1
838 0 : do ik2=0,nk2
839 0 : do ik3=0,nk3
840 :
841 0 : ikgrid=ikgrid+1
842 0 : kptgrid(1)=DBLE(ik1)/kptrlatt(1,1)
843 0 : kptgrid(2)=DBLE(ik2)/kptrlatt(2,2)
844 0 : kptgrid(3)=DBLE(ik3)/kptrlatt(3,3)
845 0 : call wrap2_pmhalf(kptgrid(1),kpt(1),res)
846 0 : call wrap2_pmhalf(kptgrid(2),kpt(2),res)
847 0 : call wrap2_pmhalf(kptgrid(3),kpt(3),res)
848 :
849 : ! === Find correspondence between the Xcrysden grid and the IBZ ===
850 : ! If AFM case, use only Ferromagetic symmetries.
851 0 : found=.FALSE.
852 0 : irred: do ikpt1=1,nkptirred
853 0 : do itim=0,timrev
854 0 : do isym=1,nsym
855 0 : if (use_afm.and.symafm(isym)==-1) CYCLE
856 0 : timsign = one-two*itim
857 : kptsym(:) = timsign*(symrec(:,1,isym)*kptirred(1,ikpt1) + &
858 : symrec(:,2,isym)*kptirred(2,ikpt1) + &
859 0 : symrec(:,3,isym)*kptirred(3,ikpt1))
860 0 : call wrap2_pmhalf(kptsym(1),kconv(1),res)
861 0 : call wrap2_pmhalf(kptsym(2),kconv(2),res)
862 0 : call wrap2_pmhalf(kptsym(3),kconv(3),res)
863 : ! is kconv equivalent to kpt?
864 0 : ss= (kpt(1)-kconv(1))**2 + (kpt(2)-kconv(2))**2 + (kpt(3)-kconv(3))**2
865 0 : if (ss < tol6) then
866 0 : found=.TRUE.
867 0 : fulltoirred(ikgrid)=ikpt1
868 : exit irred
869 : end if
870 :
871 : end do !itim
872 : end do !isym
873 : end do irred
874 :
875 0 : if (.not.found) then
876 : write(msg,'(a,3es16.8,2a)')&
877 0 : ' kpt = ',kpt,ch10,' has no symmetric among the irred k-points used in the GS calculation '
878 0 : ierr=ierr+1
879 0 : ABI_ERROR(msg)
880 : end if
881 :
882 : end do !ik1
883 : end do !ik2
884 : end do !ik3
885 :
886 0 : if (ierr/=0) then
887 0 : ABI_FREE(fulltoirred)
888 0 : RETURN
889 : end if
890 :
891 0 : if (abs(ewind) < tol12 ) then
892 : ! Keep all bands.
893 : minband=1
894 : maxband=mband
895 : else
896 : ! Select a subset of bands.
897 : minband = mband
898 : maxband = 0
899 : ene=abs(ewind)
900 0 : do isppol=1,nsppol
901 0 : do iband=1,mband
902 0 : if(minval(eigen(iband,:,isppol))-fermie < -ene) then
903 0 : minband = iband
904 : end if
905 : end do
906 0 : do iband=mband,1,-1
907 0 : if (maxval(eigen(iband,:,isppol))-fermie > ene) then
908 0 : maxband = iband
909 : end if
910 : end do
911 : end do ! isppol
912 :
913 : end if ! abs(energy_window)
914 :
915 : ! Dump the results on file ===
916 0 : if (open_file(fname,msg,newunit=uvtk,status='unknown',form='formatted') /= 0) then
917 0 : ABI_FREE(fulltoirred)
918 0 : ABI_WARNING(msg)
919 0 : ierr=ierr +1; RETURN
920 : end if
921 :
922 : ! write header
923 0 : write(uvtk,"(a)") '# vtk DataFile Version 2.0'
924 0 : write(uvtk,"(a)") 'Eigen values for the Fermi surface'
925 0 : write(uvtk,"(a)") 'ASCII'
926 0 : write(uvtk,*) ''
927 0 : write(uvtk,"(a)") 'DATASET STRUCTURED_GRID'
928 0 : write(uvtk,"(a,3i6)") 'DIMENSIONS', nk1+1,nk2+1,nk3+1
929 0 : write(uvtk,"(a,i6,a)") 'POINTS',nkptfull,' float'
930 :
931 0 : do ik3 = 0, nk3
932 0 : do ik2 = 0, nk2
933 0 : do ik1 = 0, nk1
934 : write(uvtk,'(3es16.8)') dble(ik1)/nk1*gprimd(1,1)+ &
935 : dble(ik2)/nk2*gprimd(1,2)+ &
936 0 : dble(ik3)/nk3*gprimd(1,3), &
937 : dble(ik1)/nk1*gprimd(2,1)+ &
938 : dble(ik2)/nk2*gprimd(2,2)+ &
939 0 : dble(ik3)/nk3*gprimd(2,3), &
940 : dble(ik1)/nk1*gprimd(3,1)+ &
941 : dble(ik2)/nk2*gprimd(3,2)+ &
942 0 : dble(ik3)/nk3*gprimd(3,3)
943 : end do
944 : end do
945 : end do
946 :
947 : !print out data for all relevant bands and full kpt grid (redundant, yes)
948 : !for each kpt in full zone, find equivalent irred kpt and print eigenval
949 0 : write(uvtk,*) ''
950 0 : write(uvtk,"(a,i6)") 'POINT_DATA',nkptfull
951 0 : indx=0
952 0 : do iband=minband,maxband
953 0 : do isppol=1,nsppol
954 0 : if (minband+indx < 10) then
955 0 : format_str="(a14,i1,1X,a)"
956 : else
957 0 : format_str="(a14,i2,1X,a)"
958 : end if
959 0 : write(uvtk,format_str) 'SCALARS eigval', minband+indx, 'float 1'
960 0 : write(uvtk,"(a)") 'LOOKUP_TABLE default'
961 0 : write(uvtk,*) ' '
962 0 : do kkpt = nk3/2+1, nk3+nk3/2+1
963 0 : do jkpt = nk2/2+1, nk2+nk2/2+1
964 0 : do ikpt = nk1/2+1, nk1+nk1/2+1
965 0 : ik1 = ikpt
966 0 : ik2 = jkpt
967 0 : ik3 = kkpt
968 0 : if (ikpt > nk1+1) ik1 = ikpt - nk1
969 0 : if (jkpt > nk2+1) ik2 = jkpt - nk2
970 0 : if (kkpt > nk3+1) ik3 = kkpt - nk3
971 : ! get the index with zyx order
972 0 : ikpt_fine = (ik1-1)*(nk2+1)*(nk3+1) + (ik2-1)*(nk3+1) + ik3
973 0 : write(uvtk,'(es16.8)') eigen(iband,fulltoirred(ikpt_fine),isppol)
974 : end do
975 : end do
976 : end do
977 0 : indx=indx+1
978 : end do
979 : end do
980 :
981 0 : write(uvtk,*) ''
982 0 : indx=0
983 0 : do iband=minband,maxband
984 0 : do isppol=1,nsppol
985 0 : if (minband+indx < 10) then
986 0 : format_str="(a10,i1,1X,a)"
987 : else
988 0 : format_str="(a10,i2,1X,a)"
989 : end if
990 0 : write(uvtk,format_str) 'SCALARS ve', minband+indx, 'float'
991 0 : write(uvtk,"(a)") 'LOOKUP_TABLE default'
992 0 : write(uvtk,*) ' '
993 0 : do kkpt = nk3/2+1, nk3+nk3/2+1
994 0 : do jkpt = nk2/2+1, nk2+nk2/2+1
995 0 : do ikpt = nk1/2+1, nk1+nk1/2+1
996 0 : ik1 = ikpt
997 0 : ik2 = jkpt
998 0 : ik3 = kkpt
999 0 : if (ikpt > nk1+1) ik1 = ikpt - nk1
1000 0 : if (jkpt > nk2+1) ik2 = jkpt - nk2
1001 0 : if (kkpt > nk3+1) ik3 = kkpt - nk3
1002 : ! write(uvtk,'(3i6,3es16.8)') ik1,ik2,ik3,v_surf(iband,ik1,ik2,ik3,1,isppol), &
1003 : ! & v_surf(iband,ik1,ik2,ik3,2,isppol), &
1004 : ! & v_surf(iband,ik1,ik2,ik3,3,isppol)
1005 : write(uvtk,'(es16.8)') sqrt(v_surf(iband,ik1,ik2,ik3,1,isppol)*v_surf(iband,ik1,ik2,ik3,1,isppol)+ &
1006 : & v_surf(iband,ik1,ik2,ik3,2,isppol)*v_surf(iband,ik1,ik2,ik3,2,isppol)+ &
1007 0 : & v_surf(iband,ik1,ik2,ik3,3,isppol)*v_surf(iband,ik1,ik2,ik3,3,isppol))
1008 : end do
1009 : end do
1010 : end do
1011 0 : write(uvtk,format_str) 'SCALARS vx', minband+indx, 'float'
1012 0 : write(uvtk,"(a)") 'LOOKUP_TABLE default'
1013 0 : write(uvtk,*) ' '
1014 0 : do kkpt = nk3/2+1, nk3+nk3/2+1
1015 0 : do jkpt = nk2/2+1, nk2+nk2/2+1
1016 0 : do ikpt = nk1/2+1, nk1+nk1/2+1
1017 0 : ik1 = ikpt
1018 0 : ik2 = jkpt
1019 0 : ik3 = kkpt
1020 0 : if (ikpt > nk1+1) ik1 = ikpt - nk1
1021 0 : if (jkpt > nk2+1) ik2 = jkpt - nk2
1022 0 : if (kkpt > nk3+1) ik3 = kkpt - nk3
1023 0 : write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,1,isppol)
1024 : end do
1025 : end do
1026 : end do
1027 0 : write(uvtk,format_str) 'SCALARS vy', minband+indx, 'float'
1028 0 : write(uvtk,"(a)") 'LOOKUP_TABLE default'
1029 0 : write(uvtk,*) ' '
1030 0 : do kkpt = nk3/2+1, nk3+nk3/2+1
1031 0 : do jkpt = nk2/2+1, nk2+nk2/2+1
1032 0 : do ikpt = nk1/2+1, nk1+nk1/2+1
1033 0 : ik1 = ikpt
1034 0 : ik2 = jkpt
1035 0 : ik3 = kkpt
1036 0 : if (ikpt > nk1+1) ik1 = ikpt - nk1
1037 0 : if (jkpt > nk2+1) ik2 = jkpt - nk2
1038 0 : if (kkpt > nk3+1) ik3 = kkpt - nk3
1039 0 : write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,2,isppol)
1040 : end do
1041 : end do
1042 : end do
1043 0 : write(uvtk,format_str) 'SCALARS vz', minband+indx, 'float'
1044 0 : write(uvtk,"(a)") 'LOOKUP_TABLE default'
1045 0 : write(uvtk,*) ' '
1046 0 : do kkpt = nk3/2+1, nk3+nk3/2+1
1047 0 : do jkpt = nk2/2+1, nk2+nk2/2+1
1048 0 : do ikpt = nk1/2+1, nk1+nk1/2+1
1049 0 : ik1 = ikpt
1050 0 : ik2 = jkpt
1051 0 : ik3 = kkpt
1052 0 : if (ikpt > nk1+1) ik1 = ikpt - nk1
1053 0 : if (jkpt > nk2+1) ik2 = jkpt - nk2
1054 0 : if (kkpt > nk3+1) ik3 = kkpt - nk3
1055 0 : write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,3,isppol)
1056 : end do
1057 : end do
1058 : end do
1059 0 : indx=indx+1
1060 : end do
1061 : end do
1062 :
1063 0 : close (uvtk)
1064 0 : ABI_FREE(fulltoirred)
1065 :
1066 : end subroutine printvtk
1067 : !!***
1068 :
1069 : END MODULE m_pptools
|
