TABLE OF CONTENTS


ABINIT/cut3d [ Programs ]

[ Top ] [ Programs ]

NAME

 cut3d

FUNCTION

 Main routine for the analysis of the density and potential files,
 as well as other files with the ABINIT header.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (GMR, RC, LSI, XG, NCJ, JFB, MCote, LPizzagalli)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  (main program)

OUTPUT

  (main program)

NOTES

 natom = number of atoms in the unit cell
 nr1,nr2,nr3 = grid size (nr1 x nr2 x nr3 = filrho dimension)
 ntypat = number of atom types
 ucvol = unit cell volume (> 0)
 filrho = name of the density file (binary or netcdf)

PARENTS

CHILDREN

      abi_io_redirect,abimem_init,abinit_doctor,crystal_free,crystal_from_hdr
      cut3d_hirsh,cut3d_lineint,cut3d_planeint,cut3d_pointint,cut3d_rrho
      cut3d_volumeint,cut3d_wffile,destroy_mpi_enreg,fftdatar_write
      flush_unit,hdr_echo,hdr_free,hdr_read_from_fname,herald
      init_distribfft_seq,initmpi_seq,metric,ngfft_seq,timein,wrtout,xmpi_end
      xmpi_init,xred2xcart

SOURCE

 42 #if defined HAVE_CONFIG_H
 43 #include "config.h"
 44 #endif
 45 
 46 #include "abi_common.h"
 47 
 48 program cut3d
 49 
 50  use defs_basis
 51  use defs_abitypes
 52  use m_errors
 53  use m_build_info
 54  use m_xmpi
 55  use m_nctk
 56  use m_abicore
 57 #ifdef HAVE_NETCDF
 58  use netcdf
 59 #endif
 60 #if defined FC_NAG
 61  use f90_unix_proc
 62 #endif
 63  use m_hdr
 64  use m_cut3d
 65  use m_crystal
 66  use m_crystal_io
 67 
 68  use m_specialmsg,      only : specialmsg_getcount, herald
 69  use m_fstrings,        only : endswith, sjoin, itoa
 70  use m_time,            only : timein
 71  use m_geometry,        only : xred2xcart, metric
 72  use m_mpinfo,          only : destroy_mpi_enreg, initmpi_seq
 73  use m_fftcore,         only : ngfft_seq
 74  use m_distribfft,      only : init_distribfft_seq
 75  use m_ioarr,           only : fftdatar_write
 76  use m_io_tools,        only : flush_unit, file_exists, open_file, is_open, get_unit, read_string
 77 
 78 !This section has been created automatically by the script Abilint (TD).
 79 !Do not modify the following lines by hand.
 80 #undef ABI_FUNC
 81 #define ABI_FUNC 'cut3d'
 82 !End of the abilint section
 83 
 84  implicit none
 85 
 86 !Local variables-------------------------------
 87  character(len=1) :: outputchar,blank=' '
 88 !scalars
 89  integer,parameter :: paral_kgb0=0,mfiles=10,exchn2n3d0=0
 90  integer :: fform0,gridshift1,gridshift2,gridshift3,i1,i2,i3
 91  integer :: iatom,ifiles,ii,ii1,ii2,ii3,index,iprompt,ir1,ir2,ir3,ispden,cplex
 92  integer :: itask,jfiles,natom,nfiles,nr1,nr2,unt,comm,iomode,nprocs,my_rank
 93  integer :: nr3,nr1_stored,nr2_stored,nr3_stored,nrws,nspden,nspden_stored,ntypat,timrev,nfft
 94  real(dp) :: dotdenpot,maxmz,normz,sumdenpot,ucvol,xm,xnow,xp,ym,ynow,yp,zm,znow,zp,tcpui,twalli
 95  character(len=24) :: codename
 96  character(len=fnlen) :: filnam,filrho,filrho_tmp
 97  character(len=nctk_slen) :: varname
 98  type(hdr_type) :: hdr
 99  type(abifile_t) :: abifile
100  type(MPI_type) :: mpi_enreg
101  type(crystal_t) :: cryst
102 !arrays
103  integer, allocatable :: isdenpot(:)
104  integer :: ngfft(18)
105  real(dp) :: rprimd(3,3),shift_tau(3),tsec(2)
106  real(dp) :: xcart2(3),gmet(3,3),gprimd(3,3),rmet(3,3)
107  real(dp),allocatable :: grid(:,:,:),grid_full(:,:,:,:),grid_full_stored(:,:,:,:,:),gridtt(:,:,:)
108  real(dp),allocatable :: tau2(:,:),xcart(:,:),xred(:,:),rhomacu(:,:),gridmz(:,:,:),gridmy(:,:,:),gridmx(:,:,:)
109  character(len=fnlen),allocatable :: filrho_stored(:)
110  character(len=500) :: message
111 
112 !******************************************************************
113 
114 !Change communicator for I/O (mandatory!)
115  call abi_io_redirect(new_io_comm=xmpi_world)
116 
117 !Initialize MPI
118  call xmpi_init()
119  comm = xmpi_world
120  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
121  ABI_CHECK(nprocs == 1, "cut3d not programmed for parallel execution")
122 
123 !Initialize memory profiling if it is activated
124 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
125 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
126 #ifdef HAVE_MEM_PROFILING
127  call abimem_init(0)
128 #endif
129 
130  call timein(tcpui,twalli)
131 
132 !Default for sequential use
133 !Other values of mpi_enreg are dataset dependent, and should NOT be initialized inside cut3d.F90.
134  call initmpi_seq(mpi_enreg)
135 
136  codename='CUT3D '//repeat(' ',18)
137  call herald(codename,abinit_version,std_out)
138 
139 !BIG LOOP on files
140  ABI_MALLOC(isdenpot,(mfiles))
141  isdenpot=0
142  ABI_MALLOC(filrho_stored,(mfiles))
143  iomode = IO_MODE_FORTRAN
144 
145  do ifiles=1,mfiles
146 
147 !  Get name of density file
148    write(std_out,*)
149    write(std_out,*) ' What is the name of the 3D function (density, potential or wavef) file ?'
150    if (read_string(filrho, unit=std_in) /= 0) then
151      MSG_ERROR("Fatal error!")
152    end if
153    filrho_tmp=adjustl(filrho)
154    do ii=1,len_trim(filrho_tmp)
155      if(filrho_tmp(ii:ii)==blank)then
156        filrho=trim(filrho_tmp(1:ii-1))
157        exit
158      end if
159    end do
160    write(std_out,*) ' => Your 3D function file is: ',trim(filrho)
161    write(std_out,*)
162    ! Checking the existence of data file
163    if (nctk_try_fort_or_ncfile(filrho, message) /= 0) then
164      MSG_ERROR(message)
165    end if
166 
167 !  Treat the different cases: formatted or unformatted
168    iomode = IO_MODE_FORTRAN; if (endswith(filrho, ".nc")) iomode = IO_MODE_ETSF
169    if (iomode == IO_MODE_FORTRAN) then
170      write(std_out,"(a)") '- Your file contains unformatted binary header + 3D data'
171    else
172      write(std_out,"(a)") '- Your file contains ETSF data'
173    end if
174 
175    ! Read the header and extract dimensions.
176    write(std_out,*)
177    call hdr_read_from_fname(hdr, filrho, fform0, comm)
178    ABI_CHECK(fform0 /= 0, "hdr_read returned fform = 0")
179    abifile = abifile_from_fform(fform0)
180    ABI_CHECK(abifile%fform /= 0, "Cannot detect abifile from fform")
181 
182 !  Echo part of the header
183    call hdr_echo(hdr, fform0, 4)
184 
185    nr1=hdr%ngfft(1); nr2=hdr%ngfft(2); nr3=hdr%ngfft(3)
186    natom=hdr%natom
187    nspden=hdr%nspden
188    ntypat=hdr%ntypat
189    rprimd(:,:)=hdr%rprimd(:,:)
190 
191 !  Need to know natom in order to allocate xcart
192    ABI_MALLOC(xcart,(3,natom))
193    ABI_MALLOC(xred,(3,natom))
194    xred(:,:)=hdr%xred(:,:)
195    call xred2xcart(natom,rprimd,xcart,xred)
196 
197    ispden=0
198    if (abifile%class == "density" .or. abifile%class == "potential") then
199      if(nspden/=1)then
200        write(std_out,'(a)' )' '
201        write(std_out,'(a)' )' * This file contains more than one spin component,'
202        write(std_out,'(a,i3,a)' )'  (indeed, nspden=',nspden,' )'
203        write(std_out,'(a)' )'  Some of the tasks that you will define later will concern all spin components.'
204        write(std_out,'(a)' )'  Others tasks might require you to have chosen among the following:'
205      end if
206      if(nspden==2)then
207        write(std_out,'(a)' )'   ispden= 0 ==> Total density'
208        write(std_out,'(a)' )'   ispden= 1 ==> spin-up density'
209        write(std_out,'(a)' )'   ispden= 2 ==> spin-down density'
210        write(std_out,'(a)' )'   ispden= 3 ==> spin-polarization (or magnetization) density'
211        write(std_out,'(a)' )'                 spin up - spin down difference.'
212      end if
213      if(nspden==4)then
214        write(std_out,'(a)' )'   ispden= 0 ==> Total density'
215        write(std_out,'(a)' )'   ispden= 1 ==> magnetization in the x direction'
216        write(std_out,'(a)' )'   ispden= 2 ==> magnetization in the y direction'
217        write(std_out,'(a)' )'   ispden= 3 ==> magnetization in the z direction'
218        write(std_out,'(a)' )'   ispden= 4 might be used to plot the magnetization (3D) in the XCrysDen format,'
219      end if
220      if(nspden/=1)then
221        write(std_out,*)'  Please define ispden:'
222        read(std_in,*)ispden
223        write(std_out,'(a,i3)' )' You entered ispden=',ispden
224      end if
225    end if
226 
227    write(std_out,*)
228    write(std_out,*) '==========================================================='
229    write(std_out,*)
230 
231 !  Echo the value of different input parameters
232    write(std_out,*)'ECHO important input variables ...'
233    write(std_out,*)
234    write(std_out,*) ' Dimensional primitive vectors (ABINIT equivalent: rprimd):'
235    write(std_out,'(3es16.6)' ) rprimd(1:3,1)
236    write(std_out,'(3es16.6)' ) rprimd(1:3,2)
237    write(std_out,'(3es16.6)' ) rprimd(1:3,3)
238 
239 !  Compute ucvol and test the non-collinearity of rprimd vectors.
240    call metric(gmet,gprimd,dev_null,rmet,rprimd,ucvol)
241 
242    write(std_out,'(a,3i5)' ) '  Grid density (ABINIT equivalent: ngfft): ',nr1,nr2,nr3
243    write(std_out,*) ' Number of atoms       :',natom
244    write(std_out,*) ' Number of atomic types:',ntypat
245 
246    write(std_out,*)
247    write(std_out,*) '  #    Atomic positions (cartesian coordinates - Bohr)'
248    do iatom=1,natom
249      write(std_out,'(i4,3es16.6)' )iatom,xcart(1:3,iatom)
250    end do
251    write(std_out,*)
252 
253 !  ------------------------------------------------------------------------
254 !  Branching: either WF file, or DEN/POT file.
255 
256    if (abifile%class == "wf_planewave") then
257      write(std_out,*)' This file is a WF file. '
258      isdenpot(ifiles)=0
259      iprompt = 0 ! this needs to be initialized, as it is used after the loop on files...
260 
261      call cut3d_wffile(filrho,hdr%ecut_eff,exchn2n3d0,hdr%istwfk,hdr%kptns,natom,hdr%nband,hdr%nkpt,hdr%npwarr,&
262 &     nr1,nr2,nr3,hdr%nspinor,hdr%nsppol,ntypat,paral_kgb0,rprimd,xcart,hdr%typat,hdr%znucltypat)
263      call hdr_free(hdr)
264 
265 !    -------------------------------------------------------------------------
266 ! This is a DEN/POT file
267    else if (abifile%class == "density" .or. abifile%class == "potential") then
268 
269 !    This should become a subroutine
270      write(std_out,*)' This file is a Density or Potential file '
271      isdenpot(ifiles)=1
272 
273 !    Read the function on the 3D grid
274      ABI_MALLOC(grid,(nr1,nr2,nr3))
275      ABI_MALLOC(grid_full,(nr1,nr2,nr3,nspden))
276      ABI_MALLOC(gridtt,(nr1,nr2,nr3))
277      ABI_MALLOC(gridmx,(nr1,nr2,nr3))
278      ABI_MALLOC(gridmy,(nr1,nr2,nr3))
279      ABI_MALLOC(gridmz,(nr1,nr2,nr3))
280 
281      varname = varname_from_fname(filrho)
282      if (iomode == IO_MODE_ETSF) then
283        call wrtout(std_out, sjoin("- Reading netcdf variable: ", varname))
284      end if
285 
286      call cut3d_rrho(filrho,varname,iomode,grid_full,nr1,nr2,nr3,nspden)
287 
288 !    Do not forget that the first sub-array of a density file is the total density,
289 !    while the first sub-array of a potential file is the spin-up potential
290      if (abifile%class == "density") then
291 
292 !      gridtt= grid --> Total density or potential.
293 !      gridmx= grid --> spin-Up density, or magnetization density in X direction.
294 !      gridmy= grid --> spin-Down density, or magnetization density in Y direction.
295 !      gridmz= grid --> spin-polarization density (Magnetization),
296 !      or magnetization density in Z direction.
297        gridtt(:,:,:)=grid_full(:,:,:,1)
298        if(nspden==2)then
299          gridmx = grid_full(:,:,:,2)
300          gridmy = grid_full(:,:,:,1)-grid_full(:,:,:,2)
301          gridmz = -grid_full(:,:,:,1)+two*grid_full(:,:,:,2)
302        else if(nspden==4)then
303          gridmx = grid_full(:,:,:,2)
304          gridmy = grid_full(:,:,:,3)
305          gridmz = grid_full(:,:,:,4)
306        end if
307 
308        if(nspden==1)then
309          grid = grid_full(:,:,:,1)
310        else
311          if(ispden==0)then
312            grid = gridtt
313          else if(ispden==1)then
314            grid = gridmx
315          else if(ispden==2)then
316            grid = gridmy
317          else if(ispden==3)then
318            grid = gridmz
319 !          if(ispden==0)then
320 !          grid(:,:,:)=grid_full(:,:,:,1)
321 !          else if(ispden==1)then
322 !          grid(:,:,:)=grid_full(:,:,:,2)
323 !          else if(ispden==2)then
324 !          grid(:,:,:)=grid_full(:,:,:,1)-grid_full(:,:,:,2)
325 !          else if(ispden==-1)then
326 !          grid(:,:,:)=-grid_full(:,:,:,1)+two*grid_full(:,:,:,2)
327          else if(ispden==4)then
328            write(std_out,*) ' '
329          else
330            MSG_ERROR(sjoin('bad ispden value = ',itoa(ispden)))
331          end if
332        end if
333 
334      else if (abifile%class == "potential") then   ! Potential case
335        if(ispden==0)then
336          grid(:,:,:)=grid_full(:,:,:,1)
337        else if(ispden==1 .or. ispden==2)then
338          grid(:,:,:)=grid_full(:,:,:,ispden)
339        else
340          MSG_ERROR(sjoin('bad ispden value = ',itoa(ispden)))
341        end if
342        gridtt = grid
343      end if
344 
345      write(std_out,*)
346      write(std_out,*) ' 3D function was read. Ready for further treatment.'
347      write(std_out,*)
348      write(std_out,*) '==========================================================='
349      write(std_out,*)
350 
351 !    ------------------------------------------------------------------------
352 
353 !    At this moment all the input is done
354 !    The code knows the geometry of the system,
355 !    and the data file (electron density, potential, etc).
356 !    It will further calculate the electron density by interpolation in
357 !    a point, along a line or in a plane.
358 
359      do
360        do
361          write(std_out,*) ' What is your choice ? Type:'
362          write(std_out,*) '  0 => exit'
363          write(std_out,*) '  1 => point  (interpolation of data for a single point)'
364          write(std_out,*) '  2 => line   (interpolation of data along a line)'
365          write(std_out,*) '  3 => plane  (interpolation of data in a plane)'
366          write(std_out,*) '  4 => volume (interpolation of data in a volume)'
367          write(std_out,*) '  5 => 3D formatted data (output the bare 3D data - one column)'
368          write(std_out,*) '  6 => 3D indexed data (bare 3D data, preceeded by 3D index)'
369          write(std_out,*) '  7 => 3D Molekel formatted data '
370          write(std_out,*) '  8 => 3D data with coordinates (tecplot ASCII format)'
371          write(std_out,*) '  9 => output .xsf file for XCrysDen'
372          write(std_out,*) ' 11 => compute atomic charge using the Hirshfeld method'
373          write(std_out,*) ' 14 => Gaussian/cube wavefunction module'
374          write(std_out,*) ' 15 => Write data to netcdf file'
375          read(std_in,*) itask
376          write(std_out,'(a,a,i2,a)' ) ch10,' Your choice is ',itask,ch10
377 
378          if ((5 <= itask .and. itask <= 9) .or. any(itask == [14, 15]) )then
379            write(std_out,*) ch10,'  Enter the name of an output file:'
380            if (read_string(filnam, unit=std_in) /= 0) then
381              MSG_ERROR("Fatal error!")
382            end if
383            write(std_out,*) '  The name of your file is: ',trim(filnam)
384          end if
385 
386          select case(itask)
387 
388          case(1) ! point calculation
389            call cut3d_pointint(gridtt,gridmx,gridmy,gridmz,nr1,nr2,nr3,nspden,rprimd)
390            exit
391 
392          case(2) ! line calculation
393            call cut3d_lineint(gridtt,gridmx,gridmy,gridmz,nr1,nr2,nr3,nspden,rprimd)
394            exit
395 
396          case(3) ! plane calculation
397            call cut3d_planeint(gridtt,gridmx,gridmy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,xcart)
398            exit
399 
400          case(4) ! volume calculation
401            write(std_out,*) ' Enter volume calculation'
402            call cut3d_volumeint(gridtt,gridmx,gridmy,gridmz,natom,nr1,nr2,nr3,nspden,rprimd,xcart)
403            exit
404 
405          case(5)
406            ! Rewrite the data on a formatted file, just in one (or four) column(s)
407            if (open_file(filnam,message, newunit=unt, status='unknown', action="write") /= 0) then
408              MSG_ERROR(message)
409            end if
410 
411            if(nspden==1)then
412              do i3=1,nr3
413                do i2=1,nr2
414                  do i1=1,nr1
415                    write(unt,'(4(es22.12))') grid(i1,i2,i3)
416                  end do
417                end do
418              end do
419            else
420              do i3=1,nr3
421                do i2=1,nr2
422                  do i1=1,nr1
423                    write(unt,'(4(es22.12))') gridtt(i1,i2,i3), gridmx(i1,i2,i3), gridmy(i1,i2,i3), gridmz(i1,i2,i3)
424                  end do
425                end do
426              end do
427            end if
428            close(unt)
429            exit
430 
431          case(6)
432 !            Rewrite the data on a formatted file, 3D index + density
433            if (open_file(filnam,message, newunit=unt, status='unknown', action="write") /= 0) then
434              MSG_ERROR(message)
435            end if
436 
437            if(nspden==1)then
438              write(unt,*)'   i1    i2    i3      data '
439              do i3=1,nr3
440                do i2=1,nr2
441                  do i1=1,nr1
442                    write(unt,'(3i6,4(es24.14))') i1,i2,i3,grid(i1,i2,i3)
443                  end do
444                end do
445              end do
446            else
447              if(nspden==2)then
448                write(unt,*)'   i1    i2    i3     non-spin-polarized spin up  spin down  difference  '
449              else if(nspden==4)then
450                write(unt,*)'   i1    i2    i3     non-spin-polarized   x       y      z   '
451              end if
452              do i3=1,nr3
453                do i2=1,nr2
454                  do i1=1,nr1
455                    write(unt,'(3i6,4(es24.14))') i1,i2,i3,gridtt(i1,i2,i3),gridmx(i1,i2,i3),gridmy(i1,i2,i3),gridmz(i1,i2,i3)
456                  end do
457                end do
458              end do
459            end if ! nspden
460            close(unt)
461            exit
462 
463          case(7)
464            if (open_file(filnam,message, newunit=unt, form='unformatted', action="write") /= 0) then
465              MSG_ERROR(message)
466            end if
467 
468            xm=0 ; xp=rprimd(1,1)*Bohr_Ang
469            ym=0 ; yp=rprimd(2,2)*Bohr_Ang
470            zm=0 ; zp=rprimd(3,3)*Bohr_Ang
471            write(std_out,'(/,a,/)' )' Extremas (x,y,z) of the cube in which the molecule is placed, in Angstroms'
472            write(std_out,'(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp
473            write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1,nr2,nr3
474            write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', nr1*nr2*nr3
475            write(unt) xm,xp,ym,yp,zm,zp,nr1,nr2,nr3
476            ABI_MALLOC(rhomacu,(nr1,nr2))
477            do i3=1,nr3
478              do i2=1,nr2
479                do i1=1,nr1
480                  rhomacu(i1,i2)=grid(i1,i2,i3)
481                end do
482              end do
483              write(unt) rhomacu(:,:)
484            end do
485            close(unt)
486            exit
487 
488          case (8)
489            if (open_file(filnam, message, newunit=unt, form='formatted', action="write") /= 0) then
490              MSG_ERROR(message)
491            end if
492 
493            write(std_out,'(/,a,/)' )' Extremas (x,y,z) of the cube in which the molecule is placed, in Angstroms'
494            write(std_out,'(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp
495            write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1,nr2,nr3
496            write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', nr1*nr2*nr3
497            write(unt,'(a)') 'TITLE = "  " '
498            write(unt,'(a)') 'VARIABLES = "X"  "Y"  "Z" (all three in Angstrom)  "DENSITY or POTENTIAL" (atomic units) '
499            write(unt,'(3(a,i6),a)') 'ZONE I=',nr1, ' J=', nr2, ' K=', nr3, ' F=POINT'
500            do i3=1,nr3
501              do i2=1,nr2
502                do i1=1,nr1
503                  xnow = rprimd(1,1)*(i1-1)/nr1 + rprimd(1,2)*(i2-1)/nr2 + rprimd(1,3)*(i3-1)/nr3
504                  ynow = rprimd(2,1)*(i1-1)/nr1 + rprimd(2,2)*(i2-1)/nr2 + rprimd(2,3)*(i3-1)/nr3
505                  znow = rprimd(3,1)*(i1-1)/nr1 + rprimd(3,2)*(i2-1)/nr2 + rprimd(3,3)*(i3-1)/nr3
506                  write(unt,'(4es22.15)') Bohr_Ang*xnow, Bohr_Ang*ynow, Bohr_Ang*znow, grid (i1,i2,i3)
507                end do
508              end do
509            end do
510            close(unt)
511            exit
512 
513          case (9)
514            if (open_file(filnam, message, newunit=unt, form='formatted', action="write") /= 0) then
515              MSG_ERROR(message)
516            end if
517            xm=0 ; xp=rprimd(1,1)*Bohr_Ang
518            ym=0 ; yp=rprimd(2,2)*Bohr_Ang
519            zm=0 ; zp=rprimd(3,3)*Bohr_Ang
520            write(std_out,'(/,a,/)' )' Extremas (x,y,z) of the cube in which the molecule is placed, in Angstroms'
521            write(std_out,'(5x,6f10.5)' ) xm,xp,ym,yp,zm,zp
522            write(std_out,'(/,a,2x,3i5)' )' Number of points per side: ',nr1+1,nr2+1,nr3+1
523            write(std_out,'(/,a,2x,i10,//)' )' Total number of points:', (nr1+1)*(nr2+1)*(nr3+1)
524            write(std_out,*) '  znucl = ', hdr%znucltypat, ' type = ', hdr%typat, ' ntypat = ', ntypat
525 
526            gridshift1 = 0
527            gridshift2 = 0
528            gridshift3 = 0
529            write(std_out,*) 'Do you want to shift the grid along the x,y or z axis (y/n)?'
530            write(std_out,*)
531            shift_tau(:) = zero
532            read(std_in,"(a)") outputchar
533            if (outputchar == 'y' .or. outputchar == 'Y') then
534              write(std_out,*) 'Give the three shifts (x,y,z < ',nr1,nr2,nr3,'):'
535              write(std_out,*)
536              read(std_in,*) gridshift1, gridshift2, gridshift3
537              shift_tau(:) = gridshift1*rprimd(:,1)/(nr1+1) + gridshift2*rprimd(:,2)/(nr2+1) + gridshift3*rprimd(:,3)/(nr3+1)
538            end if
539 !
540 !            Generate translated coordinates to match density shift
541 !
542            ABI_MALLOC(tau2,(3,natom))
543            do iatom = 1,natom
544              tau2(:,iatom) = xcart(:,iatom) - shift_tau(:)
545            end do
546 !            ################################################################### (LD)
547 !            Option only available for "xcrysden" format as documented at the beginning
548            if (ispden==4) then
549 !              It is necessary to know previously how many atoms will be used.
550 !              in order to plot the necessary magnetization arrows only.
551              write(std_out,*)'Is it possible to decrease the number of arrows in order to improve the'
552              write(std_out,*)'visualization in the screen, and decrease the size of the xcrysden output file.'
553              write(std_out,*)'How many arrows would you like to skip? 0 = take all. 1 = skip every other point...'
554              read (std_in,*) nrws
555              nrws=nrws+1
556              index=natom
557              maxmz=0.0
558              do i1=1,nr1,nrws
559                do i2=1,nr2,nrws
560                  do i3=1,nr3,nrws
561                    normz=gridmx(i1,i2,i3)**2+gridmy(i1,i2,i3)**2+gridmz(i1,i2,i3)**2
562                    if(normz > maxmz) maxmz=normz
563                  end do
564                end do
565              end do
566              if(abs(maxmz)<tol10)then
567                MSG_ERROR('At least, one of the components must differ from zero.')
568              end if
569              do i1=1,nr1,nrws
570                do i2=1,nr2,nrws
571                  do i3=1,nr3,nrws
572                    normz=gridmx(i1,i2,i3)**2+gridmy(i1,i2,i3)**2+gridmz(i1,i2,i3)**2
573                    if(0.1*maxmz <= normz) index=index+1
574                  end do
575                end do
576              end do
577 
578              write(unt,'(1X,A)') 'CRYSTAL'
579              write(unt,'(1X,A)') 'PRIMVEC'
580              do i1 = 1,3
581                write(unt,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3)
582              end do
583              write(unt,'(1X,A)') 'PRIMCOORD'
584              write(unt,*) index, '1'
585 
586              ! write out atom types and positions
587              do iatom = 1,natom
588                write(unt,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),Bohr_Ang*tau2(1:3,iatom)
589              end do
590 
591              ! write out magnetization vectors.
592              ! xcrysden consider these as X (dummy) atoms.
593              do i1=1,nr1,nrws
594                do i2=1,nr2,nrws
595                  do i3=1,nr3,nrws
596                    normz=gridmx(i1,i2,i3)**2+gridmy(i1,i2,i3)**2+gridmz(i1,i2,i3)**2
597                    if(0.1*maxmz <= normz) then
598                      xcart2 = matmul (rprimd, (/(i1-one)/nr1, (i2-one)/nr2, (i3-one)/nr3/))
599                      write(unt,'(A,1X,6(ES17.10,2X))')'X',&
600                      Bohr_Ang*(xcart2(1)-shift_tau(1)),&
601                      Bohr_Ang*(xcart2(2)-shift_tau(2)),&
602                      Bohr_Ang*(xcart2(3)-shift_tau(3)),&
603                      gridmx(i1,i2,i3),&
604                      gridmy(i1,i2,i3),&
605                      gridmz(i1,i2,i3)
606                    end if
607                  end do
608                end do
609              end do
610            else
611 !              ################################################################### (LD)
612 !
613 !              normal case: output density or potential (scalar field)
614              write(unt,'(1X,A)')  'DIM-GROUP'
615              write(unt,*) '3  1'
616              write(unt,'(1X,A)') 'PRIMVEC'
617              do i1 = 1,3
618                write(unt,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3)
619              end do
620              write(unt,'(1X,A)') 'PRIMCOORD'
621              write(unt,*) natom, ' 1'
622              do iatom = 1,natom
623                write(unt,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),Bohr_Ang*tau2(1:3,iatom)
624              end do
625              write(unt,'(1X,A)') 'ATOMS'
626              do iatom = 1,natom
627                write(unt,'(i9,3(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),Bohr_Ang*tau2(1:3,iatom)
628              end do
629 !              write(31,'(1X,A)') 'FRAMES'
630              write(unt,'(1X,A)') 'BEGIN_BLOCK_DATAGRID3D'
631              write(unt,*) 'datagrids'
632              write(unt,'(1X,A)') 'DATAGRID_3D_DENSITY'
633              write(unt,*) nr1+1,nr2+1,nr3+1
634              write(unt,*) '0.0 0.0 0.0 '
635              do i1 = 1,3
636                write(unt,'(3(ES17.10,2X))') (Bohr_Ang*rprimd(i2,i1), i2=1,3)
637              end do
638 
639              index = 0
640              do ir3=gridshift3+1,nr3+1
641                ii3=mod(ir3-1,nr3) + 1
642                do ir2=gridshift2+1,nr2+1
643                  ii2=mod(ir2-1,nr2) + 1
644                  do ir1=gridshift1+1,nr1+1
645                    ii1=mod(ir1-1,nr1) + 1
646                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
647                    index = index+1
648                    if (mod (index,6) == 0) write (unt,*)
649                  end do
650                  do ir1=1,gridshift1
651                    ii1=mod(ir1-1,nr1) + 1
652                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
653                    index = index+1
654                    if (mod (index,6) == 0) write (unt,*)
655                  end do
656                end do
657                do ir2=1,gridshift2
658                  ii2=mod(ir2-1,nr2) + 1
659                  do ir1=gridshift1+1,nr1+1
660                    ii1=mod(ir1-1,nr1) + 1
661                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
662                    index = index+1
663                    if (mod (index,6) == 0) write (unt,*)
664                  end do
665                  do ir1=1,gridshift1
666                    ii1=mod(ir1-1,nr1) + 1
667                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
668                    index = index+1
669                    if (mod (index,6) == 0) write (unt,*)
670                  end do
671                end do
672              end do
673              do ir3=1,gridshift3
674                ii3=mod(ir3-1,nr3) + 1
675                do ir2=gridshift2+1,nr2+1
676                  ii2=mod(ir2-1,nr2) + 1
677                  do ir1=gridshift1+1,nr1+1
678                    ii1=mod(ir1-1,nr1) + 1
679                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
680                    index = index+1
681                    if (mod (index,6) == 0) write (unt,*)
682                  end do
683                  do ir1=1,gridshift1
684                    ii1=mod(ir1-1,nr1) + 1
685                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
686                    index = index+1
687                    if (mod (index,6) == 0) write (unt,*)
688                  end do
689                end do
690                do ir2=1,gridshift2
691                  ii2=mod(ir2-1,nr2) + 1
692                  do ir1=gridshift1+1,nr1+1
693                    ii1=mod(ir1-1,nr1) + 1
694                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
695                    index = index+1
696                    if (mod (index,6) == 0) write (unt,*)
697                  end do
698                  do ir1=1,gridshift1
699                    ii1=mod(ir1-1,nr1) + 1
700                    write(unt,'(e20.5,2x)',ADVANCE='NO') grid(ii1,ii2,ii3)
701                    index = index+1
702                    if (mod (index,6) == 0) write (unt,*)
703                  end do
704                end do
705              end do
706              write (unt,*)
707              write(unt,'(1X,A)') 'END_DATAGRID_3D'
708              write(unt,'(1X,A)') 'END_BLOCK_DATAGRID3D'
709 
710            end if
711 
712            close(unt)
713            exit
714 
715          case(11)
716            call cut3d_hirsh(grid,natom,nr1,nr2,nr3,ntypat,rprimd,xcart,hdr%typat,hdr%zionpsp,hdr%znucltypat)
717            exit
718 
719          case(14) ! CUBE file format from GAUSSIAN
720            write(std_out,*)
721            write(std_out,*) 'Output a cube file of 3D volumetric data'
722            write(std_out,*)
723 
724 !            EXAMPLE FROM THE WEB
725 !            CPMD CUBE FILE.
726 !            OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z
727 !            3    0.000000    0.000000    0.000000
728 !            40    0.283459    0.000000    0.000000
729 !            40    0.000000    0.283459    0.000000
730 !            40    0.000000    0.000000    0.283459
731 !            8    0.000000    5.570575    5.669178    5.593517
732 !            1    0.000000    5.562867    5.669178    7.428055
733 !            1    0.000000    7.340606    5.669178    5.111259
734 !            -0.25568E-04  0.59213E-05  0.81068E-05  0.10868E-04  0.11313E-04  0.35999E-05
735 
736            if (open_file(filnam,message,newunit=unt, status='unknown', form='formatted', action="write") /= 0) then
737              MSG_ERROR(message)
738            end if
739 
740            !%% call print_fofr_cube(nr1,nr2,n3,nr1,nr2,nr3,fofr,rprimd,natom,znucl_atom,xcart,unit=unt)
741            write(unt,'(a)') 'ABINIT generated cube file'
742            write(unt,'(a)') 'from cut3d tool'
743 
744            write(unt,'(i9,3(1x,f12.6))') natom,0.,0.,0.
745            write(unt,'(i9,3(1x,f12.6))') nr1,(rprimd(ir2,1)/nr1, ir2=1,3)
746            write(unt,'(i9,3(1x,f12.6))') nr2,(rprimd(ir2,2)/nr2, ir2=1,3)
747            write(unt,'(i9,3(1x,f12.6))') nr3,(rprimd(ir2,3)/nr3, ir2=1,3)
748 
749            do iatom=1,natom
750              write(unt,'(i9,4(3X,ES17.10))') nint(hdr%znucltypat(hdr%typat(iatom))),0.d0, &
751 &             xcart(1,iatom),xcart(2,iatom),xcart(3,iatom)
752            end do
753 
754 !            C ordering of the indexes
755            do i1=1,nr1
756              do i2=1,nr2
757                do i3=1,nr3
758                  write(unt,'(6(f12.6,2x))') grid(i1,i2,i3)
759                end do
760              end do
761            end do
762 
763            close(unt)
764            exit
765 
766          case (15)
767            ! Write netcdf file.
768            timrev = 2; if (any(hdr%kptopt == [3, 4])) timrev = 1
769            call crystal_from_hdr(cryst, hdr, timrev)
770            call ngfft_seq(ngfft, [nr1, nr2, nr3])
771            ngfft(4:6) = ngfft(1:3)
772            nfft = product(ngfft(1:3))
773            cplex = 1
774            call init_distribfft_seq(mpi_enreg%distribfft, 'c', ngfft(2), ngfft(3), 'all')
775            call init_distribfft_seq(mpi_enreg%distribfft, 'f', ngfft(2), ngfft(3), 'all')
776 
777            call fftdatar_write(varname,filnam,IO_MODE_ETSF,hdr,cryst,ngfft,cplex,nfft,nspden,grid_full,mpi_enreg)
778            call crystal_free(cryst)
779 
780          case(0)
781            write(std_out,*)' Exit requested by user'
782            exit
783 
784          case default
785            MSG_ERROR(sjoin("Wrong task:", itoa(itask)))
786          end select
787        end do
788 
789        write(std_out,*) ' Task ',itask,' has been done !'
790        write(std_out,*)
791        write(std_out,'(a)') ' More analysis of the 3D file ? ( 0=no ; 1=default=yes ; 2= treat another file - restricted usage)'
792        read(std_in,*) iprompt
793        if(iprompt/=1) then
794          call hdr_free(hdr)
795          exit
796        else
797          cycle
798        end if
799      end do
800 
801    else
802      MSG_ERROR(sjoin("Don't know how to handle file class ", abifile%class))
803    end if ! WF file or DEN/POT file
804 
805 !  A maximum number of files had been previously specified, but set the actual number of files
806 !  to 1 if one does not read at least one other.
807    if(ifiles==1)then
808      nfiles=1
809      if(iprompt==2)nfiles=mfiles
810 
811 !    A data structure for storing the important information should be created ...
812 !    Here, one supposes that the files are compatible ...
813      if(isdenpot(ifiles)==1)then
814        ABI_MALLOC(grid_full_stored,(nr1,nr2,nr3,nspden,nfiles))
815        nr1_stored=nr1
816        nr2_stored=nr2
817        nr3_stored=nr3
818        nspden_stored=nspden
819      else if(isdenpot(ifiles)/=1 .and. iprompt==2)then
820        MSG_ERROR("in case of storage mode, the first file must be a density/potential file.")
821      end if
822    end if
823 
824    if(isdenpot(ifiles)==1) grid_full_stored(:,:,:,:,ifiles)=grid_full(:,:,:,:)
825    if(isdenpot(ifiles)==1) filrho_stored(ifiles)=filrho
826 
827    if(allocated(xcart)) then
828      ABI_FREE(xcart)
829    end if
830    if(allocated(xred)) then
831      ABI_FREE(xred)
832    end if
833    if(allocated(grid)) then
834      ABI_FREE(grid)
835    end if
836    if(allocated(grid_full)) then
837      ABI_FREE(grid_full)
838    end if
839    if(allocated(gridtt)) then
840      ABI_FREE(gridtt)
841    end if
842    if(allocated(gridmx)) then
843      ABI_FREE(gridmx)
844    end if
845    if(allocated(gridmy)) then
846      ABI_FREE(gridmy)
847    end if
848    if(allocated(gridmz)) then
849      ABI_FREE(gridmz)
850    end if
851    if(allocated(rhomacu)) then
852      ABI_FREE(rhomacu)
853    end if
854    if(allocated(tau2)) then
855      ABI_FREE(tau2)
856    end if
857 
858    if(iprompt/=2) then
859      exit
860    end if
861 
862  end do ! End big loop on files
863 
864 !Will provide different information on the density and potential files
865  do ifiles=1,nfiles
866    if(isdenpot(ifiles)==1)then
867      write(std_out,*)
868      write(std_out,*) ' Provide some global information about the density and/or potential file(s)'
869      exit
870    end if
871  end do
872  do ifiles=1,nfiles
873    if(isdenpot(ifiles)==1)then
874      write(std_out,*)
875      write(std_out, '(a,i5,3a)' ) '-  File number ',ifiles,', with name "',trim(filrho_stored(ifiles)),'"'
876      write(std_out, '(a,i12,a,es14.6)' ) '  Number of grid points =',nr1*nr2*nr3,' ; Volume of real space cell (Bohr^3)=',ucvol
877      do ispden=1,nspden
878        sumdenpot=sum(grid_full_stored(:,:,:,ispden,ifiles))
879        write(std_out, '(a,i5,3a)' ) '   Spin-component number ',ispden
880        write(std_out, '(a,3es16.6)' ) '      Sum of values, mean, mean times cell volume=',&
881 &       sumdenpot,sumdenpot/real(nr1*nr2*nr3),sumdenpot*ucvol/real(nr1*nr2*nr3)
882      end do
883    end if
884  end do
885 
886  if(nspden==1)then
887 !  At present, only nspden=1 is correctly implemented, due to specificities of the treatment of the spin-density
888    do ifiles=1,nfiles
889      if(isdenpot(ifiles)==1)then
890        write(std_out,*)
891        write(std_out,'(a)') ' Provide some global joint information about the stored density and potential file(s)'
892        exit
893      end if
894    end do
895    do ifiles=1,nfiles
896      if(isdenpot(ifiles)==1)then
897        do jfiles=ifiles,nfiles
898          if(isdenpot(jfiles)==1)then
899            write(std_out,*)
900            write(std_out, '(a,2i5)' )'  File numbers: ',ifiles,jfiles
901            do ispden=1,nspden
902              dotdenpot=zero
903              do ir1=1,nr1
904                do ir2=1,nr2
905                  do ir3=1,nr3
906                    dotdenpot=dotdenpot+grid_full_stored(ir1,ir2,ir3,ispden,ifiles)*grid_full_stored(ir1,ir2,ir3,ispden,jfiles)
907                  end do
908                end do
909              end do
910              write(std_out, '(a,i5,3a)' ) '   Spin-component number ',ispden
911              write(std_out, '(a,3es16.6)' ) '      Dot product of values, mean, mean times cell volume=',&
912 !            write(std_out, '(a,3es20.10)' ) '      Dot product of values, mean, mean times cell volume=',&
913 &             dotdenpot,dotdenpot/real(nr1*nr2*nr3),dotdenpot*ucvol/real(nr1*nr2*nr3)
914            end do
915          end if
916        end do
917      end if
918    end do
919  end if
920 
921  ABI_FREE(filrho_stored)
922 
923  if(allocated(grid_full_stored)) then
924    ABI_FREE(grid_full_stored)
925  end if
926 
927  if(allocated(isdenpot)) then
928    ABI_FREE(isdenpot)
929  end if
930 
931  call timein(tsec(1),tsec(2))
932  tsec(1)=tsec(1)-tcpui
933  tsec(2)=tsec(2)-twalli
934 
935  write(std_out, '(3a,f13.1,a,f13.1)' )'-',ch10,'- Proc.   0 individual time (sec): cpu=',tsec(1),'  wall=',tsec(2)
936 
937  write(std_out,*)
938  write(std_out,*) ' Thank you for using me'
939  write(std_out,*)
940 
941  call flush_unit(std_out)
942  call destroy_mpi_enreg(mpi_enreg)
943  call abinit_doctor("__cut3d")
944  call xmpi_end()
945 
946  end program cut3d