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-2024 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_multipoles
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27  use m_distribfft
28  use m_xmpi
29  use m_atomdata
30 
31  use defs_abitypes,    only : mpi_type
32  use m_io_tools, only : open_file
33  use m_geometry, only : xred2xcart
34  use m_mpinfo,   only : ptabs_fourdp
35 
36  implicit none
37 
38  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

SOURCE

 72 subroutine multipoles_fftr(arraysp,dipole,nfft,ngfft,nspden,rprimd,origin,&
 73                            distribfft,mpi_comm_grid)
 74 
 75 !Arguments ------------------------------------
 76 !scalars
 77  integer,intent(in) :: nfft,nspden
 78  integer,intent(in),optional :: mpi_comm_grid
 79  type(distribfft_type),intent(in),optional,target :: distribfft
 80 !arrays
 81  integer,intent(in) :: ngfft(18)
 82  real(dp),intent(in) :: arraysp(nfft,nspden),origin(3),rprimd(3,3)
 83  real(dp),intent(out) :: dipole(3,nspden)
 84 !Local variables-------------------------------
 85 !scalars
 86  integer,parameter :: ishift=5
 87  integer :: ierr,ifft1,ifft2,ifft3,ifft0,ifft,ispden,i1,i2,i3,n1,n2,n3
 88  integer :: me_fft,my_mpi_comm,nfftot
 89  logical :: fftgrid_found
 90  real(dp) :: invn1,invn2,invn3,invnfftot,wrapfft1,wrapfft2,wrapfft3
 91  character(len=500) :: msg
 92  type(distribfft_type),pointer :: my_distribfft
 93 !arrays
 94  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
 95  real(dp) :: dipole_tmp(3,nspden)
 96 
 97 ! *************************************************************************
 98 
 99 
100 !Several initializations
101  my_mpi_comm=xmpi_comm_self;if (present(mpi_comm_grid)) my_mpi_comm=mpi_comm_grid
102  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
103  nfftot=n1*n2*n3;invnfftot=one/(dble(nfftot))
104  invn1=one/real(n1,kind=dp);invn2=one/real(n2,kind=dp);invn3=one/real(n3,kind=dp)
105  me_fft=xmpi_comm_rank(my_mpi_comm)
106 
107  dipole(:,:)=zero
108 
109 !Get the distrib associated with the FFT grid
110  if (present(distribfft)) then
111    my_distribfft => distribfft
112  else
113    ABI_MALLOC(my_distribfft,)
114    call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp')
115  end if
116  fftgrid_found=.false.
117  if (n2 == my_distribfft%n2_coarse ) then
118    if (n3 == size(my_distribfft%tab_fftdp3_distrib)) then
119      fftn3_distrib => my_distribfft%tab_fftdp3_distrib
120      ffti3_local => my_distribfft%tab_fftdp3_local
121      fftgrid_found=.true.
122    end if
123  end if
124  if (n2 == my_distribfft%n2_fine ) then
125    if (n3 == size(my_distribfft%tab_fftdp3dg_distrib)) then
126      fftn3_distrib => my_distribfft%tab_fftdp3dg_distrib
127      ffti3_local => my_distribfft%tab_fftdp3dg_local
128      fftgrid_found = .true.
129    end if
130  end if
131  if (.not.(fftgrid_found)) then
132    msg='Unable to find an allocated distrib for the FFT grid!'
133    ABI_BUG(msg)
134  end if
135 
136 !Loop over FFT grid points
137 !$OMP PARALLEL PRIVATE(ifft,i1,i2,ifft1,ifft2,ispden,wrapfft1,wrapfft2)
138  do ifft3=1,n3
139    i3=mod(ifft3-1+ishift*n3,n3)
140    if(fftn3_distrib(1+i3)==me_fft) then
141      wrapfft3=mod(real(ifft3-1,kind=dp)*invn3-origin(3)+1.5_dp,one)-half
142      ifft0=1+n1*n2*(ffti3_local(1+i3)-1)
143 !$OMP SINGLE
144      dipole_tmp=zero
145 !$OMP END SINGLE
146 !$OMP DO COLLAPSE(2) REDUCTION(+:dipole_tmp)
147      do ifft2=1,n2
148        do ifft1=1,n1
149          i2=mod(ifft2-1+ishift*n2,n2)
150          i1=mod(ifft1-1+ishift*n1,n1)
151          wrapfft2=mod(real(ifft2-1,kind=dp)*invn2-origin(2)+1.5_dp,one)-half
152          wrapfft1=mod(real(ifft1-1,kind=dp)*invn1-origin(1)+1.5_dp,one)-half
153          ifft=ifft0+n1*i2+i1
154 
155 !        Computation of integral(s)
156          do ispden=1,nspden
157            dipole_tmp(1,ispden)=dipole_tmp(1,ispden)+wrapfft1*arraysp(ifft,ispden)
158            dipole_tmp(2,ispden)=dipole_tmp(2,ispden)+wrapfft2*arraysp(ifft,ispden)
159            dipole_tmp(3,ispden)=dipole_tmp(3,ispden)+wrapfft3*arraysp(ifft,ispden)
160          end do
161 
162        end do
163      end do
164 !$OMP END DO
165 !$OMP SINGLE
166      dipole=dipole+dipole_tmp
167 !$OMP END SINGLE
168    end if
169  end do
170 !$OMP END PARALLEL
171 
172 !MPI parallelization
173  if (xmpi_comm_size(my_mpi_comm)>1) then
174    call xmpi_sum(dipole,my_mpi_comm,ierr)
175  end if
176 
177 !From reduced to cartesian coordinates
178  do ispden=1,nspden
179    dipole(:,ispden)=matmul(rprimd,dipole(:,ispden))*invnfftot
180  end do
181 
182  if (.not.present(distribfft)) then
183    call destroy_distribfft(my_distribfft)
184    ABI_FREE(my_distribfft)
185  end if
186 
187 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

SOURCE

218 subroutine multipoles_out(rhor,mpi_enreg,natom,nfft,ngfft,nspden,&
219                           ntypat,rprimd,typat,ucvol,unit_out,xred,ziontypat)
220 
221 !Arguments ------------------------------------
222 !scalars
223  integer,intent(in) :: natom,nfft,nspden,ntypat,unit_out
224  real(dp), intent(in) :: ucvol
225  type(MPI_type),intent(in) :: mpi_enreg
226 !arrays
227  integer, intent(in) :: ngfft(18),typat(natom)
228  real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),xred(3,natom),ziontypat(ntypat)
229 !Local variables ------------------------------
230 !scalars
231  integer :: iatom,nspden_updn
232  real(dp) :: ziontotal
233  character(len=500) :: message
234 !arrays
235  real(dp) :: center_of_charge(3),dipole_el(3,2),dipole_ions_cart(3)
236  real(dp) :: dipole_ions_red(3),dipole_tot(3),tmp(3)
237 
238 ! *************************************************************************
239 
240 !Separate spins only for nspden=2
241  nspden_updn=merge(1,2,nspden/=2)
242 
243 !Title
244 
245 !Get nuclear part of dipole
246  dipole_ions_red(:) = zero ; ziontotal = zero
247  do iatom = 1, natom
248    dipole_ions_red(:) = dipole_ions_red(:) + xred(:,iatom)*ziontypat(typat(iatom))
249    ziontotal = ziontotal + ziontypat(typat(iatom))
250  end do
251  dipole_ions_cart(:) = matmul(rprimd,dipole_ions_red(:))
252 
253 !Find coordinates of center of charge on FFT grid
254  center_of_charge(1:3) = dipole_ions_red(1:3)/ziontotal
255 
256 !Get electronic part of dipole with respect to center of charge of ions (in cart. coord.)
257  dipole_el = zero
258  call multipoles_fftr(rhor(:,1:nspden_updn),dipole_el(:,1:nspden_updn),nfft,ngfft,nspden_updn,&
259 & rprimd,center_of_charge,distribfft=mpi_enreg%distribfft,mpi_comm_grid=mpi_enreg%comm_fft)
260  dipole_el(1:3,1:nspden_updn)=dipole_el(1:3,1:nspden_updn)*ucvol
261 
262 !Take into account storage of rhor (up+dn,up)
263  if (nspden==2) then
264    tmp(1:3)=dipole_el(1:3,1)
265    dipole_el(1:3,1)=dipole_el(1:3,2)
266    dipole_el(1:3,2)=tmp(1:3)-dipole_el(1:3,2)
267  end if
268 
269 !Compute total dipole
270 !NOTE: wrt center of charge, dipole_ions is 0
271  dipole_tot(1) = - sum(dipole_el(1,1:nspden_updn))
272  dipole_tot(2) = - sum(dipole_el(2,1:nspden_updn))
273  dipole_tot(3) = - sum(dipole_el(3,1:nspden_updn))
274 
275 !Output
276  write (message, '(2a)') ch10,' ----- Electric nuclear dipole wrt the center of ionic charge ----- '
277  call wrtout(unit_out, message, 'COLL')
278  write (message, '(a,3(1x,ES12.5))') &
279 & ' Center of charge for ionic distribution (red. coord.): ',center_of_charge(1:3)
280  call wrtout(unit_out, message, 'COLL')
281  write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
282 & ' Ionic dipole (cart. coord.)     = ',dipole_ions_cart, ' (a.u.)',ch10, &
283 & '                                 = ',dipole_ions_cart/dipole_moment_debye,' (D)'
284  call wrtout(unit_out, message, 'COLL')
285  if (nspden/=2) then
286    !This is compatible with nspden=4
287    write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
288 &   ' Electronic dipole (cart. coord.)= ',dipole_el(:,1),' (a.u.)',ch10,&
289 &   '                                 = ',dipole_el(:,1)/dipole_moment_debye,' (D)'
290  else
291    write (message, '(3a,3(1x,E16.6),a,3(2a,3(1x,E16.6),a))') ' -----',ch10,&
292 &   ' Electronic dipole (cart. coord.)= ',dipole_el(:,1),' up (a.u.)',ch10,&
293 &   '                                 = ',dipole_el(:,2),' dn (a.u.)',ch10,&
294 &   '                                 = ',dipole_el(:,1)/dipole_moment_debye,' up (D)',ch10,&
295 &   '                                 = ',dipole_el(:,2)/dipole_moment_debye,' dn (D)'
296  end if
297  call wrtout(unit_out, message, 'COLL')
298  write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
299 & ' Total dipole (cart. coord.)     = ',dipole_tot,' (a.u.)',ch10,&
300 & '                                 = ',dipole_tot/dipole_moment_debye,' (D)'
301  call wrtout(unit_out, message, 'COLL')
302  call wrtout(unit_out, ' ', 'COLL')
303 
304 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

SOURCE

339 subroutine out1dm(fnameabo_app_1dm,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,&
340                   rhor,rprimd,typat,ucvol,vtrial,xred,znucl)
341 
342 !Arguments ------------------------------------
343 !scalars
344  integer,intent(in) :: natom,nfft,nspden,ntypat
345  real(dp),intent(in) :: ucvol
346  character(len=fnlen),intent(in) :: fnameabo_app_1dm
347  type(MPI_type),intent(in) :: mpi_enreg
348 !arrays
349  integer,intent(in) :: ngfft(18),typat(natom)
350  real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),vtrial(nfft,nspden)
351  real(dp),intent(in) :: znucl(ntypat),xred(3,natom)
352 
353 !Local variables-------------------------------
354 !scalars
355  integer :: ia,ib,idim,ierr,ifft,islice,ispden,na,nb,ndig,nslice,nu,temp_unit
356  integer :: comm_fft, nproc_fft, me_fft
357  real(dp) :: global_den,global_pot
358  character(len=2) :: symbol
359  character(len=500) :: message
360  type(atomdata_t) :: atom
361 !arrays
362  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
363  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
364  real(dp),allocatable :: lin_den(:),mean_pot(:),reduced_coord(:),xcart(:,:)
365  character(len=8),allocatable :: iden(:)
366 
367 ! *************************************************************************
368 
369 !Initialize the file
370  write(message, '(a,a)' ) ' io1dm : about to open file ',trim(fnameabo_app_1dm)
371  call wrtout(std_out,message,'COLL')
372  call wrtout(ab_out,message,'COLL')
373 
374  comm_fft = mpi_enreg%comm_fft; nproc_fft = xmpi_comm_size(comm_fft); me_fft = mpi_enreg%me_fft
375 
376  if (me_fft == 0) then
377    if (open_file(fnameabo_app_1dm,message,newunit=temp_unit,status='unknown',form='formatted') /= 0) then
378      ABI_ERROR(message)
379    end if
380    rewind(temp_unit)
381  end if
382 
383  write(message, '(a,a)' ) ch10,'# ABINIT package : 1DM file '
384  if (me_fft == 0) write(temp_unit,'(a)') message
385 
386  write(message, '(a,a)' )ch10,'# Primitive vectors of the periodic cell (bohr)'
387  if (me_fft == 0) write(temp_unit,'(a)') message
388  do nu=1,3
389    write(message, '(1x,a,i1,a,3f10.5)' ) '#  R(',nu,')=',rprimd(:,nu)
390    if (me_fft == 0) write(temp_unit,'(a)') message
391  end do
392 
393  write(message, '(a,a)' ) ch10,&
394 & '# Atom list        Reduced coordinates          Cartesian coordinates (bohr)'
395  if (me_fft == 0) write(temp_unit,'(a)') message
396 
397 !Set up a list of character identifiers for all atoms : iden(ia)
398  ABI_MALLOC(iden,(natom))
399  do ia=1,natom
400    call atomdata_from_znucl(atom,znucl(typat(ia)))
401    symbol = atom%symbol
402 
403    ndig=int(log10(dble(ia)+0.5_dp))+1
404    if(ndig==1) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')   '
405    if(ndig==2) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')  '
406    if(ndig==3) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,') '
407    if(ndig==4) write(iden(ia), '(a,a,i1,a)' )symbol,'(',ia,')'
408    if(ndig>4)then
409      close(temp_unit)
410      write(message, '(a,i0)' )&
411 &     ' out1dm : cannot handle more than 9999 atoms, while natom=',natom
412      ABI_ERROR(message)
413    end if
414  end do
415 
416 !Compute cartesian coordinates, and print reduced and cartesian coordinates
417  ABI_MALLOC(xcart,(3,natom))
418  call xred2xcart(natom,rprimd,xcart,xred)
419  do ia=1,natom
420    write(message, '(a,a,3f10.5,a,3f10.5)' ) &
421 &   '#   ',iden(ia),xred(1:3,ia),'    ',xcart(1:3,ia)
422    if (me_fft == 0) write(temp_unit,'(a)') message
423  end do
424  ABI_FREE(iden)
425  ABI_FREE(xcart)
426 
427 !Get the distrib associated with this fft_grid
428  call ptabs_fourdp(mpi_enreg,ngfft(2),ngfft(3),fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
429 
430  do idim=1,3
431 
432    nslice=ngfft(idim)
433 
434 !  Dummy initialisation of na, nb and ifft.
435    na=0 ; nb=0; ifft=0
436    select case(idim)
437    case(1)
438      na=ngfft(2) ; nb=ngfft(3)
439    case(2)
440      na=ngfft(1) ; nb=ngfft(3)
441    case(3)
442      na=ngfft(1) ; nb=ngfft(2)
443    end select
444 
445    ABI_MALLOC( reduced_coord,(nslice))
446    ABI_MALLOC(mean_pot,(nslice))
447    ABI_MALLOC(lin_den,(nslice))
448 
449    do ispden=1,nspden
450 
451      if(ispden==1)then
452        write(message, '(a,a,a)' ) ch10,'#===========',&
453 &       '====================================================================='
454        if (me_fft == 0) write(temp_unit,'(a)') message
455      end if
456 
457      select case(idim)
458      case(1)
459        write(message, '(a)' )'# Projection along the first dimension '
460      case(2)
461        write(message, '(a)' )'# Projection along the second dimension '
462      case(3)
463        write(message, '(a)' )'# Projection along the third dimension '
464      end select
465      if (me_fft == 0) write(temp_unit,'(a)') message
466 
467      if(nspden==2)then
468        select case(ispden)
469        case(1)
470          write(message, '(a)' )'# Spin up '
471        case(2)
472          write(message, '(a)' )'# Spin down '
473        end select
474        if (me_fft == 0) write(temp_unit,'(a)') message
475      else if (nspden == 4) then
476        select case(ispden)
477        case(1)
478          write(message, '(a)' )'# Spinor up up '
479        case(2)
480          write(message, '(a)' )'# Spinor down down'
481        case(3)
482          write(message, '(a)' )'# Spinor up down'
483        case(4)
484          write(message, '(a)' )'# Spinor down up'
485        end select
486        if (me_fft == 0) write(temp_unit,'(a)') message
487      end if
488 
489      write(message, '(2a)' ) ch10,&
490 &     '#     Red. coord. Mean KS potential  Linear density  '
491      if (me_fft == 0) write(temp_unit,'(a)') message
492 
493      write(message, '(a)' ) &
494 &     '#                  (Hartree unit)   (electron/red. unit)'
495      if (me_fft == 0) write(temp_unit,'(a)') message
496 
497      global_pot=zero
498      global_den=zero
499      mean_pot(:)=zero
500      lin_den(:)=zero
501      do islice=1,nslice
502        reduced_coord(islice)=(islice-1)/dble(nslice)
503      end do
504      if (idim == 1) then
505        do islice=1,nslice
506          do ib=1,nb
507 ! skip z values not on this processor
508            if (fftn3_distrib(ib)/=mpi_enreg%me_fft) cycle
509            do ia=1,na
510              ifft = islice + nslice*( ia    -1 + na    *(ffti3_local(ib)    -1) )
511              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
512              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
513            end do
514          end do
515        end do
516      else if (idim == 2) then
517        do islice=1,nslice
518          do ib=1,nb
519 ! skip z values not on this processor
520            if (fftn3_distrib(ib)/=mpi_enreg%me_fft) cycle
521            do ia=1,na
522              ifft = ia     + na    *( islice-1 + nslice*(ffti3_local(ib)    -1) )
523              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
524              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
525            end do
526          end do
527        end do
528      else if (idim == 3) then
529 ! this full z slice is mine in parallel case
530        do islice=1,nslice
531          if (fftn3_distrib(islice)/=mpi_enreg%me_fft) cycle
532          do ib=1,nb
533            do ia=1,na
534              ifft = ia     + na    *( ib    -1 + nb    *(ffti3_local(islice)-1) )
535              mean_pot(islice)=mean_pot(islice)+vtrial(ifft,ispden)
536              lin_den(islice)=lin_den(islice)+rhor(ifft,ispden)
537            end do
538          end do
539        end do
540      end if
541      mean_pot(:)=mean_pot(:)/dble(na*nb)
542      lin_den(:)=lin_den(:)/dble(na*nb)*ucvol
543      global_pot=sum(mean_pot)/dble(nslice)
544      global_den=sum(lin_den)/dble(nslice)
545 
546 !    FFT parallelization
547      if(mpi_enreg%nproc_fft>1)then
548        call xmpi_sum(mean_pot  ,mpi_enreg%comm_fft,ierr)
549        call xmpi_sum(global_pot,mpi_enreg%comm_fft,ierr)
550        call xmpi_sum(lin_den   ,mpi_enreg%comm_fft,ierr)
551        call xmpi_sum(global_den,mpi_enreg%comm_fft,ierr)
552      end if
553 
554      do islice=1,ngfft(idim)
555        write(message, '(f10.4,es20.6,es16.6)' )&
556 &       reduced_coord(islice),mean_pot(islice),lin_den(islice)
557        if (me_fft == 0) write(temp_unit,'(a)') message
558      end do
559 
560      write(message, '(a,a,es15.6,es16.6,a)' ) ch10,&
561 &     '# Cell mean       :',global_pot,global_den, ch10
562      if (me_fft == 0) write(temp_unit,'(a)') message
563 
564 
565 !    End of the loop on spins
566    end do
567 
568    ABI_FREE(reduced_coord)
569    ABI_FREE(mean_pot)
570    ABI_FREE(lin_den)
571 
572 !  End of the loops on the three dimensions
573  end do
574 
575  if (me_fft == 0) then
576    write(message, '(a,a,a)' ) ch10,'#===========',&
577 &   '====================================================================='
578    call wrtout(temp_unit,message,'COLL')
579    close(temp_unit)
580  end if
581 
582 end subroutine out1dm