TABLE OF CONTENTS


ABINIT/m_multipoles [ Modules ]

[ Top ] [ Modules ]

NAME

  m_multipoles

FUNCTION

  Compute spatial multipole moments of input array on FFT grid

COPYRIGHT

  Copyright (C) 2003-2018 ABINIT group (MJV, MT, XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_multipoles
28 
29  use defs_basis
30  use defs_abitypes
31  use m_errors
32  use m_abicore
33  use m_distribfft
34  use m_xmpi
35  use m_atomdata
36 
37  use m_io_tools, only : open_file
38  use m_geometry, only : xred2xcart
39  use m_mpinfo,   only : ptabs_fourdp
40 
41  implicit none
42 
43  private

m_multipoles/multipoles_fftr [ Functions ]

[ Top ] [ m_multipoles ] [ Functions ]

NAME

 multipoles_fftr

FUNCTION

  Compute spatial multipole moments of input array on FFT grid
  Namely, the electrical dipole, quadrupole, etc... of the electron density

INPUTS

  arraysp(nfft,nspden)=the array whose average has to be computed
  [distribfft<type(distribfft_type)>]= -optional- infos related to FFT parallelism
  [mpi_comm_grid]= -optional- MPI communicator over the grid
  origin(3)=vector defining the origin of the dipole (point of observation, in reduced coordinates)
  nfft=number of FFT points stored by one proc
  nfftot=total number of FFT points
  ngfft =number of subdivisions along each lattice vector
  nspden=number of spin-density components
  rprimd = dimensionful lattice vectors

OUTPUT

  dipole(nspden)=mean value of the dipole of input array, for each nspden component

PARENTS

      multipoles_out

CHILDREN

      destroy_distribfft,init_distribfft_seq,xmpi_sum

SOURCE

 83 subroutine multipoles_fftr(arraysp,dipole,nfft,ngfft,nspden,rprimd,origin,&
 84 &                          distribfft,mpi_comm_grid)
 85 
 86 
 87 !This section has been created automatically by the script Abilint (TD).
 88 !Do not modify the following lines by hand.
 89 #undef ABI_FUNC
 90 #define ABI_FUNC 'multipoles_fftr'
 91 !End of the abilint section
 92 
 93  implicit none
 94 
 95 !Arguments ------------------------------------
 96 !scalars
 97  integer,intent(in) :: nfft,nspden
 98  integer,intent(in),optional :: mpi_comm_grid
 99  type(distribfft_type),intent(in),optional,target :: distribfft
100 !arrays
101  integer,intent(in) :: ngfft(18)
102  real(dp),intent(in) :: arraysp(nfft,nspden),origin(3),rprimd(3,3)
103  real(dp),intent(out) :: dipole(3,nspden)
104 !Local variables-------------------------------
105 !scalars
106  integer,parameter :: ishift=5
107  integer :: ierr,ifft1,ifft2,ifft3,ifft0,ifft,ispden,i1,i2,i3,n1,n2,n3
108  integer :: me_fft,my_mpi_comm,nfftot
109  logical :: fftgrid_found
110  real(dp) :: invn1,invn2,invn3,invnfftot,wrapfft1,wrapfft2,wrapfft3
111  character(len=500) :: msg
112  type(distribfft_type),pointer :: my_distribfft
113 !arrays
114  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
115  real(dp) :: dipole_tmp(3,nspden)
116 
117 ! *************************************************************************
118 
119 
120 !Several initializations
121  my_mpi_comm=xmpi_comm_self;if (present(mpi_comm_grid)) my_mpi_comm=mpi_comm_grid
122  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
123  nfftot=n1*n2*n3;invnfftot=one/(dble(nfftot))
124  invn1=one/real(n1,kind=dp);invn2=one/real(n2,kind=dp);invn3=one/real(n3,kind=dp)
125  me_fft=xmpi_comm_rank(my_mpi_comm)
126 
127  dipole(:,:)=zero
128 
129 !Get the distrib associated with the FFT grid
130  if (present(distribfft)) then
131    my_distribfft => distribfft
132  else
133    ABI_DATATYPE_ALLOCATE(my_distribfft,)
134    call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp')
135  end if
136  fftgrid_found=.false.
137  if (n2 == my_distribfft%n2_coarse ) then
138    if (n3 == size(my_distribfft%tab_fftdp3_distrib)) then
139      fftn3_distrib => my_distribfft%tab_fftdp3_distrib
140      ffti3_local => my_distribfft%tab_fftdp3_local
141      fftgrid_found=.true.
142    end if
143  end if
144  if (n2 == my_distribfft%n2_fine ) then
145    if (n3 == size(my_distribfft%tab_fftdp3dg_distrib)) then
146      fftn3_distrib => my_distribfft%tab_fftdp3dg_distrib
147      ffti3_local => my_distribfft%tab_fftdp3dg_local
148      fftgrid_found = .true.
149    end if
150  end if
151  if (.not.(fftgrid_found)) then
152    msg='Unable to find an allocated distrib for the FFT grid!'
153    MSG_BUG(msg)
154  end if
155 
156 !Loop over FFT grid points
157 !$OMP PARALLEL PRIVATE(ifft,i1,i2,ifft1,ifft2,ispden,wrapfft1,wrapfft2)
158  do ifft3=1,n3
159    i3=mod(ifft3-1+ishift*n3,n3)
160    if(fftn3_distrib(1+i3)==me_fft) then
161      wrapfft3=mod(real(ifft3-1,kind=dp)*invn3-origin(3)+1.5_dp,one)-half
162      ifft0=1+n1*n2*(ffti3_local(1+i3)-1)
163 !$OMP SINGLE
164      dipole_tmp=zero
165 !$OMP END SINGLE
166 !$OMP DO COLLAPSE(2) REDUCTION(+:dipole_tmp)
167      do ifft2=1,n2
168        do ifft1=1,n1
169          i2=mod(ifft2-1+ishift*n2,n2)
170          i1=mod(ifft1-1+ishift*n1,n1)
171          wrapfft2=mod(real(ifft2-1,kind=dp)*invn2-origin(2)+1.5_dp,one)-half
172          wrapfft1=mod(real(ifft1-1,kind=dp)*invn1-origin(1)+1.5_dp,one)-half
173          ifft=ifft0+n1*i2+i1
174 
175 !        Computation of integral(s)
176          do ispden=1,nspden
177            dipole_tmp(1,ispden)=dipole_tmp(1,ispden)+wrapfft1*arraysp(ifft,ispden)
178            dipole_tmp(2,ispden)=dipole_tmp(2,ispden)+wrapfft2*arraysp(ifft,ispden)
179            dipole_tmp(3,ispden)=dipole_tmp(3,ispden)+wrapfft3*arraysp(ifft,ispden)
180          end do
181 
182        end do
183      end do
184 !$OMP END DO
185 !$OMP SINGLE
186      dipole=dipole+dipole_tmp
187 !$OMP END SINGLE
188    end if
189  end do
190 !$OMP END PARALLEL
191 
192 !MPI parallelization
193  if (xmpi_comm_size(my_mpi_comm)>1) then
194    call xmpi_sum(dipole,my_mpi_comm,ierr)
195  end if
196 
197 !From reduced to cartesian coordinates
198  do ispden=1,nspden
199    dipole(:,ispden)=matmul(rprimd,dipole(:,ispden))*invnfftot
200  end do
201 
202  if (.not.present(distribfft)) then
203    call destroy_distribfft(my_distribfft)
204    ABI_DATATYPE_DEALLOCATE(my_distribfft)
205  end if
206 
207 end subroutine multipoles_fftr

m_multipoles/multipoles_out [ Functions ]

[ Top ] [ m_multipoles ] [ Functions ]

NAME

 multipoles_out

FUNCTION

  Output multipole moments of input array on FFT grid, calculated with multipoles_fftr
  Namely, the electrical dipole, quadrupole, etc... of the electron density

INPUTS

  rhor(nfft,nspden)=electronic density
  mpi_enreg=information about MPI parallelization
  natom=number of atoms
  nfft=number of FFT points stored by one proc
  ngfft =number of subdivisions along each lattice vector
  nspden=number of spin-density components
  ntypat=number of atom types
  rprimd(3,3)=dimensionful lattice vectors
  typat(ntypat)=types of atoms
  ucvol=unit cell volume
  unit_out=file unit to print out
  ziontypat(ntypat)=ionic charge of each type of atom

OUTPUT

NOTES

PARENTS

      outscfcv

CHILDREN

      multipoles_fftr,wrtout

SOURCE

244 subroutine multipoles_out(rhor,mpi_enreg,natom,nfft,ngfft,nspden,&
245 &                         ntypat,rprimd,typat,ucvol,unit_out,xred,ziontypat)
246 
247 
248 !This section has been created automatically by the script Abilint (TD).
249 !Do not modify the following lines by hand.
250 #undef ABI_FUNC
251 #define ABI_FUNC 'multipoles_out'
252 !End of the abilint section
253 
254  implicit none
255 
256 !Arguments ------------------------------------
257 !scalars
258  integer,intent(in) :: natom,nfft,nspden,ntypat,unit_out
259  real(dp), intent(in) :: ucvol
260  type(MPI_type),intent(in) :: mpi_enreg
261 !arrays
262  integer, intent(in) :: ngfft(18),typat(natom)
263  real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),xred(3,natom),ziontypat(ntypat)
264 !Local variables ------------------------------
265 !scalars
266  integer :: iatom,nspden_updn
267  real(dp) :: ziontotal
268  character(len=500) :: message
269 !arrays
270  real(dp) :: center_of_charge(3),dipole_el(3,2),dipole_ions_cart(3)
271  real(dp) :: dipole_ions_red(3),dipole_tot(3),tmp(3)
272 
273 ! *************************************************************************
274 
275 !Separate spins only for nspden=2
276  nspden_updn=merge(1,2,nspden/=2)
277 
278 !Title
279 
280 !Get nuclear part of dipole
281  dipole_ions_red(:) = zero ; ziontotal = zero
282  do iatom = 1, natom
283    dipole_ions_red(:) = dipole_ions_red(:) + xred(:,iatom)*ziontypat(typat(iatom))
284    ziontotal = ziontotal + ziontypat(typat(iatom))
285  end do
286  dipole_ions_cart(:) = matmul(rprimd,dipole_ions_red(:))
287 
288 !Find coordinates of center of charge on FFT grid
289  center_of_charge(1:3) = dipole_ions_red(1:3)/ziontotal
290 
291 !Get electronic part of dipole with respect to center of charge of ions (in cart. coord.)
292  dipole_el = zero
293  call multipoles_fftr(rhor(:,1:nspden_updn),dipole_el(:,1:nspden_updn),nfft,ngfft,nspden_updn,&
294 & rprimd,center_of_charge,distribfft=mpi_enreg%distribfft,mpi_comm_grid=mpi_enreg%comm_fft)
295  dipole_el(1:3,1:nspden_updn)=dipole_el(1:3,1:nspden_updn)*ucvol
296 
297 !Take into account storage of rhor (up+dn,up)
298  if (nspden==2) then
299    tmp(1:3)=dipole_el(1:3,1)
300    dipole_el(1:3,1)=dipole_el(1:3,2)
301    dipole_el(1:3,2)=tmp(1:3)-dipole_el(1:3,2)
302  end if
303 
304 !Compute total dipole
305 !NOTE: wrt center of charge, dipole_ions is 0
306  dipole_tot(1) = - sum(dipole_el(1,1:nspden_updn))
307  dipole_tot(2) = - sum(dipole_el(2,1:nspden_updn))
308  dipole_tot(3) = - sum(dipole_el(3,1:nspden_updn))
309 
310 !Output
311  write (message, '(2a)') ch10,' ----- Electric nuclear dipole wrt the center of ionic charge ----- '
312  call wrtout(unit_out, message, 'COLL')
313  write (message, '(a,3(1x,ES12.5))') &
314 & ' Center of charge for ionic distribution (red. coord.): ',center_of_charge(1:3)
315  call wrtout(unit_out, message, 'COLL')
316  write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
317 & ' Ionic dipole (cart. coord.)     = ',dipole_ions_cart, ' (a.u.)',ch10, &
318 & '                                 = ',dipole_ions_cart/dipole_moment_debye,' (D)'
319  call wrtout(unit_out, message, 'COLL')
320  if (nspden/=2) then
321    !This is compatible with nspden=4
322    write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
323 &   ' Electronic dipole (cart. coord.)= ',dipole_el(:,1),' (a.u.)',ch10,&
324 &   '                                 = ',dipole_el(:,1)/dipole_moment_debye,' (D)'
325  else
326    write (message, '(3a,3(1x,E16.6),a,3(2a,3(1x,E16.6),a))') ' -----',ch10,&
327 &   ' Electronic dipole (cart. coord.)= ',dipole_el(:,1),' up (a.u.)',ch10,&
328 &   '                                 = ',dipole_el(:,2),' dn (a.u.)',ch10,&
329 &   '                                 = ',dipole_el(:,1)/dipole_moment_debye,' up (D)',ch10,&
330 &   '                                 = ',dipole_el(:,2)/dipole_moment_debye,' dn (D)'
331  end if
332  call wrtout(unit_out, message, 'COLL')
333  write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
334 & ' Total dipole (cart. coord.)     = ',dipole_tot,' (a.u.)',ch10,&
335 & '                                 = ',dipole_tot/dipole_moment_debye,' (D)'
336  call wrtout(unit_out, message, 'COLL')
337  call wrtout(unit_out, ' ', 'COLL')
338 
339 end subroutine multipoles_out

m_multipoles/out1dm [ Functions ]

[ Top ] [ m_multipoles ] [ Functions ]

NAME

 out1dm

FUNCTION

 Output the 1 dimensional mean of potential and density on the three reduced axis.

INPUTS

  fnameabo_app_1dm=name of the file in which the data is written, appended with _1DM
  natom=number of atoms in unit cell
  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nspden=number of spin-density components
  ntypat=number of types of atoms in unit cell.
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
  rprimd(3,3)=real space dimensional primitive translations (bohr)
  typat(natom)=type integer for each atom in cell
  ucvol=unit cell volume in bohr**3.
  vtrial(nfft,nspden)=INPUT Vtrial(r).
  xred(3,natom)=reduced coordinates of atoms
  znucl(ntypat)=real(dp), nuclear number of atom type

OUTPUT

SIDE EFFECTS

  data written in file fnameabo_app_1dm

NOTES

PARENTS

      outscfcv

CHILDREN

      atomdata_from_znucl,ptabs_fourdp,wrtout,xmpi_sum,xred2xcart

SOURCE

380 subroutine out1dm(fnameabo_app_1dm,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,&
381 &  rhor,rprimd,typat,ucvol,vtrial,xred,znucl)
382 
383 
384 !This section has been created automatically by the script Abilint (TD).
385 !Do not modify the following lines by hand.
386 #undef ABI_FUNC
387 #define ABI_FUNC 'out1dm'
388 !End of the abilint section
389 
390  implicit none
391 
392 !Arguments ------------------------------------
393 !scalars
394  integer,intent(in) :: natom,nfft,nspden,ntypat
395  real(dp),intent(in) :: ucvol
396  character(len=fnlen),intent(in) :: fnameabo_app_1dm
397  type(MPI_type),intent(in) :: mpi_enreg
398 !arrays
399  integer,intent(in) :: ngfft(18),typat(natom)
400  real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),vtrial(nfft,nspden)
401  real(dp),intent(in) :: znucl(ntypat),xred(3,natom)
402 
403 !Local variables-------------------------------
404 !scalars
405  integer :: ia,ib,idim,ierr,ifft,islice,ispden,na,nb,ndig,nslice,nu,temp_unit
406  integer :: comm_fft, nproc_fft, me_fft
407  real(dp) :: global_den,global_pot
408  character(len=2) :: symbol
409  character(len=500) :: message
410  type(atomdata_t) :: atom
411 !arrays
412  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
413  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
414  real(dp),allocatable :: lin_den(:),mean_pot(:),reduced_coord(:),xcart(:,:)
415  character(len=8),allocatable :: iden(:)
416 
417 ! *************************************************************************
418 
419 !Initialize the file
420  write(message, '(a,a)' ) ' io1dm : about to open file ',fnameabo_app_1dm
421  call wrtout(std_out,message,'COLL')
422  call wrtout(ab_out,message,'COLL')
423 
424  comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft
425 
426  if (me_fft == 0) then
427    if (open_file(fnameabo_app_1dm,message,newunit=temp_unit,status='unknown',form='formatted') /= 0) then
428      MSG_ERROR(message)
429    end if
430    rewind(temp_unit)
431  end if
432 
433  write(message, '(a,a)' ) ch10,'# ABINIT package : 1DM file '
434  if (me_fft == 0) write(temp_unit,'(a)') message
435 
436  write(message, '(a,a)' )ch10,'# Primitive vectors of the periodic cell (bohr)'
437  if (me_fft == 0) write(temp_unit,'(a)') message
438  do nu=1,3
439    write(message, '(1x,a,i1,a,3f10.5)' ) '#  R(',nu,')=',rprimd(:,nu)
440    if (me_fft == 0) write(temp_unit,'(a)') message
441  end do
442 
443  write(message, '(a,a)' ) ch10,&
444 & '# Atom list        Reduced coordinates          Cartesian coordinates (bohr)'
445  if (me_fft == 0) write(temp_unit,'(a)') message
446 
447 !Set up a list of character identifiers for all atoms : iden(ia)
448  ABI_ALLOCATE(iden,(natom))
449  do ia=1,natom
450    call atomdata_from_znucl(atom,znucl(typat(ia)))
451    symbol = atom%symbol
452 
453    ndig=int(log10(dble(ia)+0.5_dp))+1
454    if(ndig==1) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')   '
455    if(ndig==2) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')  '
456    if(ndig==3) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,') '
457    if(ndig==4) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')'
458    if(ndig>4)then
459      close(temp_unit)
460      write(message, '(a,i0)' )&
461 &     ' out1dm : cannot handle more than 9999 atoms, while natom=',natom
462      MSG_ERROR(message)
463    end if
464  end do
465 
466 !Compute cartesian coordinates, and print reduced and cartesian coordinates
467  ABI_ALLOCATE(xcart,(3,natom))
468  call xred2xcart(natom,rprimd,xcart,xred)
469  do ia=1,natom
470    write(message, '(a,a,3f10.5,a,3f10.5)' ) &
471 &   '#   ',iden(ia),xred(1:3,ia),'    ',xcart(1:3,ia)
472    if (me_fft == 0) write(temp_unit,'(a)') message
473  end do
474  ABI_DEALLOCATE(iden)
475  ABI_DEALLOCATE(xcart)
476 
477 !Get the distrib associated with this fft_grid
478  call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
479 
480  do idim=1,3
481 
482    nslice=ngfft(idim)
483 
484 !  Dummy initialisation of na, nb and ifft.
485    na=0 ; nb=0; ifft=0
486    select case(idim)
487    case(1)
488      na=ngfft(2) ; nb=ngfft(3)
489    case(2)
490      na=ngfft(1) ; nb=ngfft(3)
491    case(3)
492      na=ngfft(1) ; nb=ngfft(2)
493    end select
494 
495    ABI_ALLOCATE( reduced_coord,(nslice))
496    ABI_ALLOCATE(mean_pot,(nslice))
497    ABI_ALLOCATE(lin_den,(nslice))
498 
499    do ispden=1,nspden
500 
501      if(ispden==1)then
502        write(message, '(a,a,a)' ) ch10,'#===========',&
503 &       '====================================================================='
504        if (me_fft == 0) write(temp_unit,'(a)') message
505      end if
506 
507      select case(idim)
508      case(1)
509        write(message, '(a)' )'# Projection along the first dimension '
510      case(2)
511        write(message, '(a)' )'# Projection along the second dimension '
512      case(3)
513        write(message, '(a)' )'# Projection along the third dimension '
514      end select
515      if (me_fft == 0) write(temp_unit,'(a)') message
516 
517      if(nspden==2)then
518        select case(ispden)
519        case(1)
520          write(message, '(a)' )'# Spin up '
521        case(2)
522          write(message, '(a)' )'# Spin down '
523        end select
524        if (me_fft == 0) write(temp_unit,'(a)') message
525      else if (nspden == 4) then
526        select case(ispden)
527        case(1)
528          write(message, '(a)' )'# Spinor up up '
529        case(2)
530          write(message, '(a)' )'# Spinor down down'
531        case(3)
532          write(message, '(a)' )'# Spinor up down'
533        case(4)
534          write(message, '(a)' )'# Spinor down up'
535        end select
536        if (me_fft == 0) write(temp_unit,'(a)') message
537      end if
538 
539      write(message, '(2a)' ) ch10,&
540 &     '#     Red. coord. Mean KS potential  Linear density  '
541      if (me_fft == 0) write(temp_unit,'(a)') message
542 
543      write(message, '(a)' ) &
544 &     '#                  (Hartree unit)   (electron/red. unit)'
545      if (me_fft == 0) write(temp_unit,'(a)') message
546 
547      global_pot=zero
548      global_den=zero
549      mean_pot(:)=zero
550      lin_den(:)=zero
551      do islice=1,nslice
552        reduced_coord(islice)=(islice-1)/dble(nslice)
553      end do
554      if (idim == 1) then
555        do islice=1,nslice
556          do ib=1,nb
557 ! skip z values not on this processor
558            if (fftn3_distrib(ib)/=mpi_enreg%me_fft) cycle
559            do ia=1,na
560              ifft = islice + nslice*( ia    -1 + na    *(ffti3_local(ib)    -1) )
561              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
562              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
563            end do
564          end do
565        end do
566      else if (idim == 2) then
567        do islice=1,nslice
568          do ib=1,nb
569 ! skip z values not on this processor
570            if (fftn3_distrib(ib)/=mpi_enreg%me_fft) cycle
571            do ia=1,na
572              ifft = ia     + na    *( islice-1 + nslice*(ffti3_local(ib)    -1) )
573              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
574              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
575            end do
576          end do
577        end do
578      else if (idim == 3) then
579 ! this full z slice is mine in parallel case
580        do islice=1,nslice
581          if (fftn3_distrib(islice)/=mpi_enreg%me_fft) cycle
582          do ib=1,nb
583            do ia=1,na
584              ifft = ia     + na    *( ib    -1 + nb    *(ffti3_local(islice)-1) )
585              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
586              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
587            end do
588          end do
589        end do
590      end if
591      mean_pot(:)=mean_pot(:)/dble(na*nb)
592      lin_den(:)=lin_den(:)/dble(na*nb)*ucvol
593      global_pot=sum(mean_pot)/dble(nslice)
594      global_den=sum(lin_den)/dble(nslice)
595 
596 !    FFT parallelization
597      if(mpi_enreg%nproc_fft>1)then
598        call xmpi_sum(mean_pot  ,mpi_enreg%comm_fft,ierr)
599        call xmpi_sum(global_pot,mpi_enreg%comm_fft,ierr)
600        call xmpi_sum(lin_den   ,mpi_enreg%comm_fft,ierr)
601        call xmpi_sum(global_den,mpi_enreg%comm_fft,ierr)
602      end if
603 
604      do islice=1,ngfft(idim)
605        write(message, '(f10.4,es20.6,es16.6)' )&
606 &       reduced_coord(islice),mean_pot(islice),lin_den(islice)
607        if (me_fft == 0) write(temp_unit,'(a)') message
608      end do
609 
610      write(message, '(a,a,es15.6,es16.6,a)' ) ch10,&
611 &     '# Cell mean       :',global_pot,global_den, ch10
612      if (me_fft == 0) write(temp_unit,'(a)') message
613 
614 
615 !    End of the loop on spins
616    end do
617 
618    ABI_DEALLOCATE(reduced_coord)
619    ABI_DEALLOCATE(mean_pot)
620    ABI_DEALLOCATE(lin_den)
621 
622 !  End of the loops on the three dimensions
623  end do
624 
625  if (me_fft == 0) then
626    write(message, '(a,a,a)' ) ch10,'#===========',&
627 &   '====================================================================='
628    call wrtout(temp_unit,message,'COLL')
629    close(temp_unit)
630  end if
631 
632 end subroutine out1dm