TABLE OF CONTENTS

ABINIT/multipoles_fftr [ Functions ]

[ Top ] [ 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

COPYRIGHT

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

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

 38 #if defined HAVE_CONFIG_H
 39 #include "config.h"
 40 #endif
 41 
 42 #include "abi_common.h"
 43 
 44 subroutine multipoles_fftr(arraysp,dipole,nfft,ngfft,nspden,rprimd,origin,&
 45 &                          distribfft,mpi_comm_grid)
 46 
 47  use m_profiling_abi
 48  use defs_basis
 49  use m_errors
 50  use defs_abitypes
 51  use m_distribfft
 52  use m_xmpi
 53 
 54 !This section has been created automatically by the script Abilint (TD).
 55 !Do not modify the following lines by hand.
 56 #undef ABI_FUNC
 57 #define ABI_FUNC 'multipoles_fftr'
 58 !End of the abilint section
 59 
 60  implicit none
 61 
 62 !Arguments ------------------------------------
 63 !scalars
 64  integer,intent(in) :: nfft,nspden
 65  integer,intent(in),optional :: mpi_comm_grid
 66  type(distribfft_type),intent(in),optional,target :: distribfft
 67 !arrays
 68  integer,intent(in) :: ngfft(18)
 69  real(dp),intent(in) :: arraysp(nfft,nspden),origin(3),rprimd(3,3)
 70  real(dp),intent(out) :: dipole(3,nspden)
 71 !Local variables-------------------------------
 72 !scalars
 73  integer,parameter :: ishift=5
 74  integer :: ierr,ifft1,ifft2,ifft3,ifft0,ifft,ispden,i1,i2,i3,n1,n2,n3
 75  integer :: me_fft,my_mpi_comm,nfftot
 76  logical :: fftgrid_found
 77  real(dp) :: invn1,invn2,invn3,invnfftot,wrapfft1,wrapfft2,wrapfft3
 78  character(len=500) :: msg
 79  type(distribfft_type),pointer :: my_distribfft
 80 !arrays
 81  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
 82  real(dp) :: dipole_tmp(3,nspden)
 83 
 84 ! *************************************************************************
 85 
 86 
 87 !Several initializations
 88  my_mpi_comm=xmpi_comm_self;if (present(mpi_comm_grid)) my_mpi_comm=mpi_comm_grid
 89  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
 90  nfftot=n1*n2*n3;invnfftot=one/(dble(nfftot))
 91  invn1=one/real(n1,kind=dp);invn2=one/real(n2,kind=dp);invn3=one/real(n3,kind=dp)
 92  me_fft=xmpi_comm_rank(my_mpi_comm)
 93 
 94  dipole(:,:)=zero
 95 
 96 !Get the distrib associated with the FFT grid
 97  if (present(distribfft)) then
 98    my_distribfft => distribfft
 99  else
100    ABI_DATATYPE_ALLOCATE(my_distribfft,)
101    call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp')
102  end if
103  fftgrid_found=.false.
104  if (n2 == my_distribfft%n2_coarse ) then
105    if (n3 == size(my_distribfft%tab_fftdp3_distrib)) then
106      fftn3_distrib => my_distribfft%tab_fftdp3_distrib
107      ffti3_local => my_distribfft%tab_fftdp3_local
108      fftgrid_found=.true.
109    end if
110  end if
111  if (n2 == my_distribfft%n2_fine ) then
112    if (n3 == size(my_distribfft%tab_fftdp3dg_distrib)) then
113      fftn3_distrib => my_distribfft%tab_fftdp3dg_distrib
114      ffti3_local => my_distribfft%tab_fftdp3dg_local
115      fftgrid_found = .true.
116    end if
117  end if
118  if (.not.(fftgrid_found)) then
119    msg='Unable to find an allocated distrib for the FFT grid!'
120    MSG_BUG(msg)
121  end if
122 
123 !Loop over FFT grid points
124 !$OMP PARALLEL PRIVATE(ifft,i1,i2,ifft1,ifft2,ispden,wrapfft1,wrapfft2)
125  do ifft3=1,n3
126    i3=mod(ifft3-1+ishift*n3,n3)
127    if(fftn3_distrib(1+i3)==me_fft) then
128      wrapfft3=mod(real(ifft3-1,kind=dp)*invn3-origin(3)+1.5_dp,one)-half
129      ifft0=1+n1*n2*(ffti3_local(1+i3)-1)
130 !$OMP SINGLE
131      dipole_tmp=zero
132 !$OMP END SINGLE
133 !$OMP DO COLLAPSE(2) REDUCTION(+:dipole_tmp)
134      do ifft2=1,n2
135        do ifft1=1,n1
136          i2=mod(ifft2-1+ishift*n2,n2)
137          i1=mod(ifft1-1+ishift*n1,n1)
138          wrapfft2=mod(real(ifft2-1,kind=dp)*invn2-origin(2)+1.5_dp,one)-half
139          wrapfft1=mod(real(ifft1-1,kind=dp)*invn1-origin(1)+1.5_dp,one)-half
140          ifft=ifft0+n1*i2+i1
141 
142 !        Computation of integral(s)
143          do ispden=1,nspden
144            dipole_tmp(1,ispden)=dipole_tmp(1,ispden)+wrapfft1*arraysp(ifft,ispden)
145            dipole_tmp(2,ispden)=dipole_tmp(2,ispden)+wrapfft2*arraysp(ifft,ispden)
146            dipole_tmp(3,ispden)=dipole_tmp(3,ispden)+wrapfft3*arraysp(ifft,ispden)
147          end do
148 
149        end do
150      end do
151 !$OMP END DO
152 !$OMP SINGLE
153      dipole=dipole+dipole_tmp
154 !$OMP END SINGLE
155    end if
156  end do
157 !$OMP END PARALLEL
158 
159 !MPI parallelization
160  if (xmpi_comm_size(my_mpi_comm)>1) then
161    call xmpi_sum(dipole,my_mpi_comm,ierr)
162  end if
163 
164 !From reduced to cartesian coordinates
165  do ispden=1,nspden
166    dipole(:,ispden)=matmul(rprimd,dipole(:,ispden))*invnfftot
167  end do
168 
169  if (.not.present(distribfft)) then
170    call destroy_distribfft(my_distribfft)
171    ABI_DATATYPE_DEALLOCATE(my_distribfft)
172  end if
173 
174 end subroutine multipoles_fftr