TABLE OF CONTENTS


ABINIT/macroave [ Programs ]

[ Top ] [ Programs ]

NAME

 macroave

FUNCTION

 **********************************************************************
 The MACROAVE program implements the macroscopic average technique,
 introduced by A. Baldereschi and coworkers
 (A. Baldereschi, S. Baroni, and R. Resta, Phys. Rev. Lett. 61, 734 (1988) [[cite:Baldereschi1988]]).
 This is an extremely powerful method that relates
 microscopic quantities, typical outputs of first-principles codes,
 with macroscopic magnitudes, needed to perform electrostatic analysis.
 Within this methodology, we will be able to wash out all the
 wiggles of the rapidly-varying functions of position (resembling
 the underlying atomic structure) of the microscopic quantities,
 blowing up only the macroscopic features.
 It can be used to compute band offsets, work functions, effective
 charges, and high frequency dielectric constants, among others
 interesting physical properties.
 Ref: L. Colombo, R. Resta and S. Baroni  Phys Rev B  44, 5572 (1991) [[cite:Colombo1991]].
 Coded by P. Ordejon and J. Junquera, April 1999.
 Modified by J. Junquera, November 2001.
 **********************************************************************

COPYRIGHT

 Copyright (C) 1999-2008 (P. Ordejon, J. Junquera, J. Soler, A. Garcia)
 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 routine)

OUTPUT

  (main routine)

PARENTS

CHILDREN

      abi_io_redirect,abimem_init,abinit_doctor,four1,hdr_fort_read,hdr_free
      hdr_ncread,iorho,macroav_spline,macroav_splint,thetaft,xmpi_end
      xmpi_init

SOURCE

 48 #if defined HAVE_CONFIG_H
 49 #include "config.h"
 50 #endif
 51 
 52 #include "abi_common.h"
 53 
 54 
 55 program macroave
 56 
 57  use defs_basis
 58  use defs_abitypes
 59  use m_xmpi
 60  use m_abicore
 61  use m_errors
 62  use m_nctk
 63 #ifdef HAVE_NETCDF
 64  use netcdf
 65 #endif
 66  use m_hdr
 67  use m_macroave
 68 
 69  use m_fstrings,        only : sjoin, strcat, endswith
 70  use m_io_tools,        only : file_exists, open_file
 71 
 72 !This section has been created automatically by the script Abilint (TD).
 73 !Do not modify the following lines by hand.
 74 #undef ABI_FUNC
 75 #define ABI_FUNC 'macroave'
 76 !End of the abilint section
 77 
 78  implicit none
 79 
 80 !Arguments -----------------------------------
 81 
 82 !Local variables-------------------------------
 83 !no_abirules
 84 ! --------- PARAMETERS -------------------------------------------------
 85 ! INTEGER   NP    : Parameter needed to define maximum number of points
 86 !              for FFT grid.
 87 !              Number of points = (2**NP)
 88 ! INTEGER   N     : Maximum number of complex data point for FFT.
 89 !              MUST be a power of 2
 90 ! REAL*8   HARTREE: Conversion factor from Hartrees to Rydbergs
 91 !              1 hartree = 2 Ry
 92 ! REAL*8   RYDBERG: Conversion factor from Rydbergs to eV
 93 !              1 Ry = 13.6058 eV
 94 ! ----------------------------------------------------------------------
 95  integer, parameter ::  np=12
 96  integer, parameter ::  n=2**np
 97  real(dp), parameter :: hartree=two
 98  real(dp), parameter :: rydberg=13.6058d0
 99 ! --------- VARIABLES --------------------------------------------------
100  integer :: i,ii,ij,ip,is,j,iomode
101  integer :: nconv,npoints,npt,nsm,nspin,nz,nspden
102  integer :: unit1,unit2,unit3,unit4,varid
103  integer :: mesh(3)
104  character(len=10) :: code,interp
105  character(len=15) :: SNAME,inpdata
106  character(len=fnlen) :: fnamerho,fnamedelv,fnameplave
107  character(len=500) :: msg
108  character(len=nctk_slen) :: varname
109  logical :: siesta,abinit,potential,charge,totalcharge
110  logical :: found,linear,splin
111  real,allocatable :: rhos(:,:)
112  real(dp),allocatable :: rho(:,:)
113  real(dp) :: cell(3,3),dcell(3,3)
114  real(dp) :: l,sur,ds,length,convfac,qren,qtot,lav1,lav2,vol
115  real(dp),allocatable :: z(:),rhoz(:),d2rhoz(:),drhodz(:)
116  real(dp) :: data(2*n),th(2*n),re(n),im(n)!,v(2*n)
117  real(dp) :: x,delta,yp1,ypn,phi
118  complex(dp) :: a,b,c
119 ! ABINIT variables
120  integer :: fform
121  integer :: comm,nproc,my_rank
122  type(hdr_type) :: hdr
123 ! end ABINIT variables
124 
125 !************************************************************************
126 !CHARACTER CODE       : First principles-code used to generate the
127 !electrostatic potential at the points of a grid
128 !in real space. It is read from file 'macroave.in'
129 !CHARACTER SNAME      : System Label
130 !(If code = ABINIT, then SNAME = FNAMERHO)
131 !CHARACTER INPDATA    : Calculate the band offset from the charge
132 !density or from the electrostatic potential?
133 !CHARACTER FNAMERHO   : Name of the file where the electrostatic
134 !potential at the mesh is stored
135 !CHARACTER FNAMEPLAVE : Name of the file where the planar average of the
136 !electronic potential will be stored
137 !CHARACTER FNAMEDELV  : Name of the file where the profile
138 !of the electrostatic potential will be stored
139 !LOGICAL   SIESTA     : Have you used SIESTA to get the electrostatic potential
140 !or the charge density?
141 !LOGICAL   ABINIT     : Have you used ABINIT to get the electrostatic potential
142 !or the charge density?
143 !LOGICAL   LINEAR     : Linear interpolation to get the charge
144 !density/potential in the FFT grid.
145 !LOGICAL   SPLIN      : Cubic spline interpolation to get the charge
146 !density/potential in the FFT grid.
147 !LOGICAL   POTENTIAL  : We are going to compute the band offset from
148 !the electrostatic potential
149 !LOGICAL   CHARGE     : We are going to compute the band offset from
150 !the charge density
151 !LOGICAL   FOUND      : Were data found? (only when task in iorho ='read')
152 !INTEGER   NATOMS     : Number of atoms in unit cell
153 !INTEGER   NSPIN      : Number of spin polarizations (1 or 2)
154 !INTEGER   MESH(3)    : Number of mesh divisions of each lattice vectors,
155 !INCLUDING subgrid
156 !INTEGER   NSM        : Number of sub-mesh points per mesh point
157 !(not used in this version)
158 !INTEGER   NPOINTS    : Number of mesh subdivisions in the normal direction
159 !to the interface
160 !INTEGER   NCONV      : Number of convolutions required to calculate the
161 !macroscopic average
162 !INTEGER   NPT        : Total number of mesh points (included subpoints)
163 !REAL*8    CELL(3,3)  : Unit cell lattice vectors (a.u.) CELL(IXYZ,IVECT)
164 !REAL*8    DS         : Differential area per point of the mesh
165 !REAL*8    SUR        : Area of a plane parallel to the interface
166 !REAL*8    LENGTH     : Distance between two planes parallel to the interface
167 !REAL*8    L          : Length of the cell in the direction nomal
168 !to the interface (a.u.)
169 !REAL*8    LAV1       : Linear period of the electrostatic potential
170 !in the bulklike region for the first material
171 !REAL*8    LAV2       : Linear period of the electrostatic potential
172 !in the bulklike region for the second material
173 !REAL*8    CONVFAC    : Conversion factor for the output units
174 !REAL*8    QTOT       : Total electronic charge in the unit cell
175 !REAL*8    QREN       : Total electronic charge calculated from
176 !the input data file
177 !REAL*4    Z(MESH(3)) : Z-coordinate of the planes where the elec. density
178 !is averaged
179 !REAL*4    RHO        : Electron density
180 !Notice single precision in this version
181 !REAL*8    RHOZ       : Planar average of the electron density
182 !REAL*8    DATA(2*N)  : Fourier coefficients of the planar average density
183 !REAL*8    TH(2*N)    : Fourier coefficients of the step functions
184 !REAL*8    V(2*N)     : Fourier coefficients of the potential
185 !REAL*4    VREEC(N)   : Real part of the electronic potential in real space
186 !REAL*4    VIMEC(N)   : Imag. part of the electronic potential in real space
187 !REAL*4    VRENC(N)   : Real part of the nuclear potential in real space
188 !REAL*4    VIMNC(N)   : Imag. part of the nuclear potential in real space
189 !*********************************************************************
190 
191 !Change communicator for I/O (mandatory!)
192  call abi_io_redirect(new_io_comm=xmpi_world)
193 
194 !Initialize MPI
195  call xmpi_init()
196  comm = xmpi_world
197  nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
198 
199 !Initialize memory profiling if it is activated
200 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
201 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
202 #ifdef HAVE_MEM_PROFILING
203  call abimem_init(0)
204 #endif
205 
206 !Reading input data from a file ---------------------------------------
207  if (open_file('macroave.in', msg, newunit=UNIT1, STATUS='OLD') /= 0) then
208    MSG_ERROR(msg)
209  end if
210 
211  READ(UNIT1,'(A)')CODE
212  READ(UNIT1,'(A)')INPDATA
213  READ(UNIT1,'(A)')SNAME
214  READ(UNIT1,*)NCONV
215  READ(UNIT1,*)LAV1
216  READ(UNIT1,*)LAV2
217  READ(UNIT1,*)QTOT
218  READ(UNIT1,*)INTERP
219  close(UNIT1)
220 
221 !Which code has been used to get the electrostatic potential? ---------
222  if ( CODE == 'siesta' .OR. CODE == 'SIESTA' .OR.&
223 & CODE == 'Siesta' ) then
224    SIESTA = .TRUE.
225    ABINIT = .FALSE.
226  else if ( CODE == 'abinit' .OR. CODE == 'ab-init' .OR.&
227 &   CODE == 'ABINIT' .OR. CODE == 'AB-INIT' .OR.&
228 &   CODE == 'Abinit' .OR. CODE == 'Ab-init') then
229    SIESTA = .FALSE.
230    ABINIT = .TRUE.
231  else
232    MSG_ERROR(sjoin('macroave: Unknown code: ', CODE))
233  end if
234 
235 !Are we going to compute the band offset from the charge density or
236 !from the electrostatic potential? ------------------------------------
237  if ( INPDATA == 'potential' .OR. INPDATA == 'POTENTIAL' .OR.&
238 & INPDATA == 'Potential' ) then
239    POTENTIAL   = .TRUE.
240    CHARGE      = .FALSE.
241    TOTALCHARGE = .FALSE.
242  else if ( INPDATA == 'charge' .OR. INPDATA == 'CHARGE' .OR.&
243 &   INPDATA == 'Charge' ) then
244    POTENTIAL   = .FALSE.
245    CHARGE      = .TRUE.
246    TOTALCHARGE = .FALSE.
247  else if ( INPDATA == 'totalcharge' .OR. &
248 &   INPDATA == 'Totalcharge' .OR.&
249 &   INPDATA == 'TotalCharge' .OR.&
250 &   INPDATA == 'TOTALCHARGE' ) then
251    POTENTIAL   = .FALSE.
252    CHARGE      = .FALSE.
253    TOTALCHARGE = .TRUE.
254  else
255    MSG_ERROR(sjoin('macroave: Unknown input data  ', INPDATA))
256  end if
257 
258 !What kind of interpolation will we use to get the charge density/
259 !potential in a FFT grid? ---------------------------------------------
260  if ( INTERP == 'linear' .OR. INTERP == 'Linear' .OR.&
261 & INTERP == 'LINEAR' ) then
262    LINEAR = .TRUE.
263    SPLIN  = .FALSE.
264  else if ( INTERP == 'spline' .OR. INTERP == 'Spline' .OR.&
265 &   INTERP == 'SPLINE' ) then
266    LINEAR = .FALSE.
267    SPLIN  = .TRUE.
268  end if
269 
270 !Reading charge density from a file -----------------------------------
271  if ( SIESTA ) then
272    if (POTENTIAL) then
273      FNAMERHO = strcat(trim(SNAME),'.VH')
274    elseif (CHARGE) then
275      FNAMERHO = strcat(trim(SNAME),'.RHO')
276    elseif (TOTALCHARGE) then
277      FNAMERHO = strcat(trim(SNAME),'.TOCH')
278    end if
279  else if ( ABINIT ) then
280    FNAMERHO = trim(SNAME)
281  end if
282 
283  if (SIESTA) then
284    NSM   = 1
285    NPT   = 0
286    NSPIN = 0
287    CALL IORHO( 'READ', trim(FNAMERHO), DCELL, MESH, NSM, NPT, NSPIN,&
288 &   RHOS, FOUND )
289    if (FOUND) then
290      ABI_ALLOCATE( RHOS,(NPT,NSPIN))
291      ABI_ALLOCATE( RHO,(NPT,NSPIN))
292      CALL IORHO( 'READ', trim(FNAMERHO), DCELL, MESH, NSM, NPT, NSPIN,&
293 &     RHOS, FOUND )
294      do I = 1, 3
295        do J = 1, 3
296          CELL(J,I) = DCELL(J,I)
297        end do
298      end do
299 !    Transform the density or the potential read from SIESTA
300 !    from a single precision variable to a double precision variable
301      do IS = 1, NSPIN
302        do IP = 1, NPT
303          RHO(IP,IS) = RHOS(IP,IS) * 1.0D0
304        end do
305      end do
306 
307    else
308      MSG_ERROR(sjoin('macroave: file not found: ', FNAMERHO))
309    end if
310 
311  else if (ABINIT) then
312 
313    if (nctk_try_fort_or_ncfile(FNAMERHO, msg) /= 0) then
314      MSG_ERROR(msg)
315    end if
316    iomode = IO_MODE_FORTRAN; if (endswith(FNAMERHO, ".nc")) iomode = IO_MODE_ETSF
317 
318    if (iomode == IO_MODE_FORTRAN) then
319      if (open_file(FNAMERHO, msg, newunit=unit2, form="unformatted", status="old") /= 0) then
320        MSG_ERROR(msg)
321      end if
322      call hdr_fort_read(hdr, unit2, fform)
323      ABI_CHECK(FFORM /= 0, "fform == 0")
324    else
325 #ifdef HAVE_NETCDF
326      NCF_CHECK(nctk_open_read(unit2, fnamerho, xmpi_comm_self))
327      call hdr_ncread(hdr, unit2, fform)
328 #else
329      MSG_ERROR("Netcdf support is missing!")
330 #endif
331    end if
332 
333 !  For debugging
334 !  call hdr_echo(hdr,fform,4,std_out)
335 
336    do I = 1, 3
337      MESH(I) = HDR%NGFFT(I)
338      do J = 1, 3
339        CELL(J,I) = HDR%RPRIMD(J,I)
340      end do
341    end do
342    NSPIN = HDR%NSPPOL
343    nspden = hdr%nspden
344    ABI_CHECK(hdr%nspinor == 1, "nspinor == 2 not coded")
345    call hdr_free(hdr)
346 
347    NPT = MESH(1) * MESH(2) * MESH(3)
348    ABI_ALLOCATE( RHO,(NPT,NSPIN))
349 
350    if (iomode == IO_MODE_FORTRAN) then
351      do IS = 1, NSPIN
352        READ(UNIT2) (RHO(IP,IS),IP=1,NPT)
353      end do
354      close(UNIT2)
355    else
356 #ifdef HAVE_NETCDF
357      varname = varname_from_fname(fnamerho)
358      NCF_CHECK(nf90_inq_varid(unit2, varname, varid))
359      ! [cplex, n1, n2, n3, nspden]
360      NCF_CHECK(nf90_get_var(unit2, varid, rho, start=[1,1,1,1,1], count=[1, mesh(1), mesh(2), mesh(3), nspden]))
361 #endif
362    end if
363 
364 !    Units for the potential in Ab-init are in Hartrees,
365 !    so we transform them into Ry. No transformation is
366 !    needed for the charge density
367 !    (it is directly read in electrons/bohr**3).
368 
369    if (POTENTIAL) then
370      do IS = 1, NSPIN
371        do IP = 1, NPT
372          RHO(IP,IS) = RHO(IP,IS) * HARTREE
373        end do
374      end do
375    end if
376  end if
377 
378 !Initialize some variables (we suppose cells with a c axis orthogonal to a and b) -------------
379 
380  L  = CELL(3,3)
381  SUR = SURPLA( CELL ) ! surface of unit cell in xy plane, perpendicular to z
382  VOL = VOLCEL( CELL )
383  DS = SUR/( MESH(1) * MESH(2) ) ! this seems adapted to arbitrary in-plane cells
384  LENGTH = L/DBLE(N)
385  NPOINTS = MESH(3)
386 
387  ABI_ALLOCATE(Z,(NPOINTS+1))
388  ABI_ALLOCATE(RHOZ,(NPOINTS+1))
389  ABI_ALLOCATE(D2RHOZ,(NPOINTS+1))
390  ABI_ALLOCATE(DRHODZ,(N))
391 
392  RHOZ(1:NPOINTS+1)   = 0.D0
393  D2RHOZ(1:NPOINTS+1) = 0.D0
394  DRHODZ(1:N)         = 0.D0
395 
396  if (POTENTIAL) then
397    CONVFAC = RYDBERG
398  else if (CHARGE) then
399    CONVFAC = 1.0D0
400  else if (TOTALCHARGE) then
401    CONVFAC = 1.0D0
402  end if
403 
404 
405 !Loop over all points and calculate the planar average ----------------
406 !Warning: The planar average is only done for the first component of
407 !RHO. Spin polarization is not implemented yet ---------------
408  do IP = 1, NPT
409    NZ = (IP-1) / (MESH(1)*MESH(2)) + 1
410    RHOZ(NZ) =  RHOZ(NZ) + RHO(IP,1)*DS
411  end do
412 
413  do IP = 1, NPOINTS
414    RHOZ(IP) = RHOZ(IP) / SUR
415  end do
416 
417  do IP  = 1, NPOINTS
418    Z(IP) = (IP-1)*CELL(3,3)/DBLE(NPOINTS)
419  end do
420 
421 !Calculate electrostatic potential or electronic charge density -------
422 !in fft grid, interpolating the planar average calculated before ------
423  if (SPLIN) then
424    Z(NPOINTS+1)    = L
425    RHOZ(NPOINTS+1) = RHOZ(1)
426    DELTA = L/DBLE(NPOINTS)
427    YP1 = ( RHOZ(2) - RHOZ(NPOINTS) )   / (2.0D0*DELTA)
428    YPN = YP1
429    CALL MACROAV_SPLINE(DELTA, RHOZ, NPOINTS+1, YP1, YPN, D2RHOZ)
430    I = 0
431    do II = 1, 2*N-1, 2
432      I = I + 1
433      X = (I-1)*L/DBLE(N)
434      CALL MACROAV_SPLINT( DELTA, RHOZ, D2RHOZ, NPOINTS+1, X, DATA(II), &
435 &     DRHODZ(I) )
436      DATA(II+1) = 0.D0
437    end do
438  else if (LINEAR) then
439    I = 0
440    do II = 1,2*N-1,2
441      I = I + 1
442      X = (I-1)*L/DBLE(N)
443      do IJ = 1, NPOINTS
444        if (X == Z(IJ)) then
445          DATA(II) = RHOZ(IJ)
446          DATA(II+1) = 0.D0
447          GOTO 20
448        end if
449        if (Z(IJ) > X) then
450          DATA(II) = RHOZ(IJ-1) +&
451 &         (X-Z(IJ-1))*(RHOZ(IJ)-RHOZ(IJ-1))/&
452 &         (Z(IJ)-Z(IJ-1))
453          DATA(II+1) = 0.D0
454          GOTO 20
455        end if
456      end do
457      DATA(II)=RHOZ(NPOINTS) +&
458 &     (X-Z(NPOINTS))*(RHOZ(1)-RHOZ(NPOINTS))/&
459 &     (Z(NPOINTS)-Z(NPOINTS-1))
460      DATA(II+1) = 0.D0
461      20       CONTINUE
462    end do
463  end if
464 
465 !Renormalize the charge density ---------------------------------------
466  if (CHARGE .OR. TOTALCHARGE) then
467    QREN = 0.D0
468    do IP = 1, 2*N-1, 2
469      QREN = QREN + DATA(IP)*LENGTH*SUR
470    end do
471    do IP = 1, 2*N-1, 2
472      if (CHARGE) then
473        DATA(IP) = DATA(IP) * QTOT/QREN
474      elseif(TOTALCHARGE) then
475        DATA(IP) = DATA(IP) - QREN/VOL
476      end if
477    end do
478    QREN = 0.D0
479    do IP = 1, 2*N-1, 2
480      QREN = QREN + DATA(IP)*LENGTH*SUR
481    end do
482 !  For debugging
483 !  write(std_out,*)' QREN = ', QREN
484  end if
485 !...
486 
487 !Print planar average of the electrostatic potential or ---------------
488 !the electronic charge density ----------------------------------------
489  FNAMEPLAVE = strcat(SNAME,'.PAV')
490  if (open_file(FNAMEPLAVE,msg,newunit=UNIT3,STATUS='UNKNOWN') /= 0) then
491    MSG_ERROR(msg)
492  end if
493  I = 0
494  do II = 1, 2*N-1, 2
495    I = I+1
496    X=(I-1)*L/DBLE(N)
497 !  WRITE(UNIT3,'(3F20.12)')X,
498 !  .           DATA(II)*CONVFAC,DATA(II+1)*CONVFAC
499    WRITE(UNIT3,'(2F20.12)')X,&
500 &   DATA(II)*CONVFAC
501  end do
502  close(UNIT3)
503 !...
504 
505 
506 !Calculate Fourier transform of the electrostatic potential or
507 !the electronic density -----------------------------------------------
508  CALL FOUR1(DATA,N,1)
509 !...
510 
511 !Calculate macroscopic average of the electrostatic potential or the
512 !electronic charge density taking the convolution with two step functions.
513 !In Fourier space, it is a product of the Fourier transform components -
514 !The decompositions in the sum over II is due to the special way in which
515 !the data are stored in subroutine four1( see Fig. 12.2.2, in
516 !'Numerical Recipes, The Art of Scientific Computing'
517 !by W.H. Press, S.A. Teukolsky, W.T. Veterling and B.P. Flannery,
518 !Cambridge U.P. 1987 and 1992.
519 
520 
521  CALL THETAFT(N,L,LAV1,TH)
522 
523  do II = 1, N+1, 2
524    A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
525    B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
526    C = A*B
527    DATA(II) = REAL(C)*L/DBLE(N)
528    DATA(II+1) = AIMAG(C)*L/DBLE(N)
529  end do
530 
531  do II = N+3, 2*N-1, 2
532    A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
533    B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
534    C = A*B
535    DATA(II) = REAL(C)*L/DBLE(N)
536    DATA(II+1) = AIMAG(C)*L/DBLE(N)
537  end do
538 
539 
540  if (NCONV == 2) then
541    CALL THETAFT(N,L,LAV2,TH)
542 
543    do II = 1, N+1, 2
544      A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
545      B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
546      C = A*B
547      DATA(II) = REAL(C)*L/DBLE(N)
548      DATA(II+1) = AIMAG(C)*L/DBLE(N)
549 !    if ( POISON ) then
550 !    IG = (II-1) / 2
551 !    GSQ= (2.D0*PI*IG/L)**2
552 !    if(GSQ > 0.D0) then
553 !    V(II) = DATA(II) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
554 !    V(II+1) = DATA(II+1) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
555 !    else
556 !    V(II) = 0.D0
557 !    V(II+1) = 0.D0
558 !    endif
559 !    endif
560    end do
561 
562    do II = N+3, 2*N-1, 2
563      A = DATA(II)*(1.D0,0.D0) + DATA(II+1)*(0.D0,1.D0)
564      B = TH(II)*(1.D0,0.D0) + TH(II+1)*(0.D0,1.D0)
565      C = A*B
566      DATA(II) = REAL(C)*L/DBLE(N)
567      DATA(II+1) = AIMAG(C)*L/DBLE(N)
568 !    if ( POISON ) then
569 !    IG = (-2*N+II-1) / 2
570 !    GSQ= (2.D0*PI*IG/L)**2
571 !    if(GSQ > 0.D0) then
572 !    V(II) = DATA(II) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
573 !    V(II+1) = DATA(II+1) * (4.D0*PI/GSQ) * HARTREE * RYDBERG
574 !    else
575 !    V(II) = 0.D0
576 !    V(II+1) = 0.D0
577 !    endif
578 !    endif
579    end do
580 
581  end if
582 !...
583 
584 !Transform average electronic density and potential to real space -----
585 !The decompositions in the sum over J is due to the special way in which
586 !the data are stored in subroutine four1( see Fig. 12.2.2, in
587 !'Numerical Recipes, The Art of Scientific Computing'
588 !by W.H. Press, S.A. Teukolsky, W.T. Veterling and B.P. Flannery,
589 !Cambridge U.P. 1987 and 1992.
590 
591  do II = 1, N
592    RE(II) = 0.D0
593    IM(II) = 0.D0
594 !  if ( POISON ) then
595 !  VREEC(II) = 0.D0
596 !  VIMEC(II) = 0.D0
597 !  endif
598    do J = 1, N+1, 2
599      PHI = -2.D0 * PI * (II-1) * ( (J-1)/2 ) / DBLE(N)
600      RE(II)=RE(II)+(1.D0/DBLE(N))*(DATA(J)*COS(PHI)&
601 &     -DATA(J+1)*SIN(PHI))
602      IM(II)=IM(II)+(1.D0/DBLE(N))*(DATA(J)*SIN(PHI)&
603 &     +DATA(J+1)*COS(PHI))
604 !    if ( POISON ) then
605 !    VREEC(II)=VREEC(II)+(1.D0/DBLE(N))*(V(J)*COS(PHI)
606 !    .                           -V(J+1)*SIN(PHI))
607 !    VIMEC(II)=VIMEC(II)+(1.D0/DBLE(N))*(V(J)*SIN(PHI)
608 !    .                           +V(J+1)*COS(PHI))
609 !    endif
610    end do
611 
612    do J = N+3, 2*N-1, 2
613      PHI = -2.0D0 * PI * (II-1) * ((-2*N+J-1)/2) / DBLE(N)
614      RE(II)=RE(II)+(1.D0/DBLE(N))*(DATA(J)*COS(PHI)&
615 &     -DATA(J+1)*SIN(PHI))
616      IM(II)=IM(II)+(1.D0/DBLE(N))*(DATA(J)*SIN(PHI)&
617 &     +DATA(J+1)*COS(PHI))
618 !    if ( POISON ) then
619 !    VREEC(II)=VREEC(II)+(1.D0/DBLE(N))*(V(J)*COS(PHI)
620 !    .                           -V(J+1)*SIN(PHI))
621 !    VIMEC(II)=VIMEC(II)+(1.D0/DBLE(N))*(V(J)*SIN(PHI)
622 !    .                           +V(J+1)*COS(PHI))
623 !    endif
624    end do
625  end do
626 !...
627 
628 !Print averaged electronic charge density and potential ---------------
629  FNAMEDELV = strcat( trim(SNAME),'.MAV')
630  if (open_file(FNAMEDELV, msg, newunit=UNIT4, STATUS='UNKNOWN') /= 0) then
631    MSG_ERROR(msg)
632  end if
633  do I = 1, N
634    X=(I-1)*L/DBLE(N)
635 !  WRITE(UNIT4,'(3F20.5)')X,
636 !  .           RE(I)*CONVFAC ,IM(I)*CONVFAC
637    WRITE(UNIT4,'(2F20.12)')X,&
638 &   RE(I)*CONVFAC
639  end do
640  close(UNIT4)
641 !...
642 
643 !Print electrostatic potential ----------------------------------------
644 !if (POISON) then
645 !FNAMEVEC = strcat( SNAME,'.VEC')
646 !if (open_file(FNAMEVEC, msg newunit=UNIT5, STATUS='UNKNOWN') /= 0) then
647 !  MSG_ERROR(msg)
648 !end if
649 !do I = 1, N
650 !X=(I-1)*L/DBLE(N)
651 !WRITE(UNIT5,'(3F20.12)')X, VREEC(I), VIMEC(I)
652 !enddo
653 !close(unit5)
654 !endif
655 !...
656 
657  ABI_DEALLOCATE(Z)
658  ABI_DEALLOCATE(RHOZ)
659  ABI_DEALLOCATE(DRHODZ)
660  ABI_DEALLOCATE(D2RHOZ)
661  if (allocated(rho)) then
662    ABI_FREE(rho)
663  end if
664 
665 !Write information on file about the memory before ending mpi module, if memory profiling is enabled
666  call abinit_doctor("__macroave")
667 
668  call xmpi_end()
669 
670  end program macroave