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-2022 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)

SOURCE

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