Branch data 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
|
