TABLE OF CONTENTS

ABINIT/multipoles_out [ Functions ]

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

COPYRIGHT

 Copyright (C) 2010-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

  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

 42 #if defined HAVE_CONFIG_H
 43 #include "config.h"
 44 #endif
 45 
 46 subroutine multipoles_out(rhor,mpi_enreg,natom,nfft,ngfft,nspden,&
 47 &                         ntypat,rprimd,typat,ucvol,unit_out,xred,ziontypat)
 48 
 49  use m_profiling_abi
 50  use defs_basis
 51  use defs_abitypes
 52 
 53 !This section has been created automatically by the script Abilint (TD).
 54 !Do not modify the following lines by hand.
 55 #undef ABI_FUNC
 56 #define ABI_FUNC 'multipoles_out'
 57  use interfaces_14_hidewrite
 58  use interfaces_53_spacepar
 59 !End of the abilint section
 60 
 61  implicit none
 62 
 63 !Arguments ------------------------------------
 64 !scalars
 65  integer,intent(in) :: natom,nfft,nspden,ntypat,unit_out
 66  real(dp), intent(in) :: ucvol
 67  type(MPI_type),intent(in) :: mpi_enreg
 68 !arrays
 69  integer, intent(in) :: ngfft(18),typat(natom)
 70  real(dp),intent(in) :: rhor(nfft,nspden),rprimd(3,3),xred(3,natom),ziontypat(ntypat)
 71 !Local variables ------------------------------
 72 !scalars
 73  integer :: iatom,nspden_updn
 74  real(dp) :: ziontotal
 75  character(len=500) :: message
 76 !arrays
 77  real(dp) :: center_of_charge(3),dipole_el(3,2),dipole_ions_cart(3)
 78  real(dp) :: dipole_ions_red(3),dipole_tot(3),tmp(3)
 79 
 80 ! *************************************************************************
 81 
 82 !Separate spins only for nspden=2
 83  nspden_updn=merge(1,2,nspden/=2)
 84 
 85 !Title
 86 
 87 !Get nuclear part of dipole
 88  dipole_ions_red(:) = zero ; ziontotal = zero
 89  do iatom = 1, natom
 90    dipole_ions_red(:) = dipole_ions_red(:) + xred(:,iatom)*ziontypat(typat(iatom))
 91    ziontotal = ziontotal + ziontypat(typat(iatom))
 92  end do
 93  dipole_ions_cart(:) = matmul(rprimd,dipole_ions_red(:))
 94 
 95 !Find coordinates of center of charge on FFT grid
 96  center_of_charge(1:3) = dipole_ions_red(1:3)/ziontotal
 97 
 98 !Get electronic part of dipole with respect to center of charge of ions (in cart. coord.)
 99  dipole_el = zero
100  call multipoles_fftr(rhor(:,1:nspden_updn),dipole_el(:,1:nspden_updn),nfft,ngfft,nspden_updn,&
101 & rprimd,center_of_charge,distribfft=mpi_enreg%distribfft,mpi_comm_grid=mpi_enreg%comm_fft)
102  dipole_el(1:3,1:nspden_updn)=dipole_el(1:3,1:nspden_updn)*ucvol
103 
104 !Take into account storage of rhor (up+dn,up)
105  if (nspden==2) then
106    tmp(1:3)=dipole_el(1:3,1)
107    dipole_el(1:3,1)=dipole_el(1:3,2)
108    dipole_el(1:3,2)=tmp(1:3)-dipole_el(1:3,2)
109  end if
110 
111 !Compute total dipole
112 !NOTE: wrt center of charge, dipole_ions is 0
113  dipole_tot(1) = - sum(dipole_el(1,1:nspden_updn))
114  dipole_tot(2) = - sum(dipole_el(2,1:nspden_updn))
115  dipole_tot(3) = - sum(dipole_el(3,1:nspden_updn))
116 
117 !Output
118  write (message, '(2a)') ch10,' ----- Electric nuclear dipole wrt the center of ionic charge ----- '
119  call wrtout(unit_out, message, 'COLL')
120  write (message, '(a,3(1x,ES12.5))') &
121 & ' Center of charge for ionic distribution (red. coord.): ',center_of_charge(1:3)
122  call wrtout(unit_out, message, 'COLL')
123  write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
124 & ' Ionic dipole (cart. coord.)     = ',dipole_ions_cart, ' (a.u.)',ch10, &
125 & '                                 = ',dipole_ions_cart/dipole_moment_debye,' (D)'
126  call wrtout(unit_out, message, 'COLL')
127  if (nspden/=2) then
128    !This is compatible with nspden=4
129    write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
130 &   ' Electronic dipole (cart. coord.)= ',dipole_el(:,1),' (a.u.)',ch10,&
131 &   '                                 = ',dipole_el(:,1)/dipole_moment_debye,' (D)'
132  else
133    write (message, '(3a,3(1x,E16.6),a,3(2a,3(1x,E16.6),a))') ' -----',ch10,&
134 &   ' Electronic dipole (cart. coord.)= ',dipole_el(:,1),' up (a.u.)',ch10,&
135 &   '                                 = ',dipole_el(:,2),' dn (a.u.)',ch10,&
136 &   '                                 = ',dipole_el(:,1)/dipole_moment_debye,' up (D)',ch10,&
137 &   '                                 = ',dipole_el(:,2)/dipole_moment_debye,' dn (D)'
138  end if
139  call wrtout(unit_out, message, 'COLL')
140  write (message, '(3a,3(1x,E16.6),3a,3(1x,E16.6),a)') ' -----',ch10,&
141 & ' Total dipole (cart. coord.)     = ',dipole_tot,' (a.u.)',ch10,&
142 & '                                 = ',dipole_tot/dipole_moment_debye,' (D)'
143  call wrtout(unit_out, message, 'COLL')
144  call wrtout(unit_out, ' ', 'COLL')
145 
146 end subroutine multipoles_out