TABLE OF CONTENTS


ABINIT/m_electronpositron [ Modules ]

[ Top ] [ Modules ]

NAME

  m_electronpositron

FUNCTION

  This module provides the definition of the electronpositron_type used
  used to store data for the electron-positron two-component DFT
  as methods to operate on it.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MT, GJ)
 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

23 #if defined HAVE_CONFIG_H
24 #include "config.h"
25 #endif
26 #include "abi_common.h"
27 
28 MODULE m_electronpositron
29 
30  use defs_basis
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_energies
35  use m_xmpi
36  use m_cgtools
37 
38  use m_pawtab,   only : pawtab_type
39  use m_paw_an,   only : paw_an_type
40  use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy
41  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
42  use m_mpinfo,   only : proc_distrb_cycle
43  use m_xcpositron, only : xcpositron
44  use m_drivexc,    only : mkdenpos
45  use m_xctk,       only : xcden
46  use m_fft,        only : fourdp
47 
48  implicit none
49 
50  private
51 
52 ! public constants
53  integer,public,parameter :: EP_NOTHING  =-1
54  integer,public,parameter :: EP_ELECTRON = 0
55  integer,public,parameter :: EP_POSITRON = 1

ABINIT/rhohxcpositron [ Functions ]

[ Top ] [ Functions ]

NAME

 rhohxcpositron

FUNCTION

 Calculate the electrons/positron correlation term for the positron

INPUTS

  gprimd(3,3)=dimensional reciprocal space primitive translations
  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
  nhat(nfft,nspden*usepaw)= -PAW only- compensation density
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  nspden=number of spin density components
  n3xccc=dimension of the xccc3d array (0 or nfft).
  paral_kgb=flag for (k,band,FFT) parallelism
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
  ucvol = unit cell volume (Bohr**3)
  usexcnhat= -PAW only- flag controling use of compensation density in Vxc
  usepaw=flag for PAW
  xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
  xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

OUTPUT

  electronpositron%e_xc=electron-positron XC energy
  electronpositron%e_xcdc=Double-counting electron-positron XC energy
  strsxc(6)= contribution of xc to stress tensor (hartree/bohr^3),
  vhartr(nfft)=Hartree potential (returned if option/=0 and option/=10)
  vxcapn=XC electron-positron XC potential for the positron
  vxcavg=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r].
  kxcapn(nfft,nkxc)=electron-positron XC kernel (returned only if nkxc/=0)

SIDE EFFECTS

  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation

PARENTS

      energy,rhotov,setvtr

CHILDREN

      mean_fftr,mkdenpos,xcden,xcpositron,xmpi_sum

SOURCE

 857 subroutine rhohxcpositron(electronpositron,gprimd,kxcapn,mpi_enreg,nfft,ngfft,nhat,nkxc,nspden,n3xccc,&
 858 &                         paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxcapn,vxcavg,xccc3d,xc_denpos)
 859 
 860 
 861 !This section has been created automatically by the script Abilint (TD).
 862 !Do not modify the following lines by hand.
 863 #undef ABI_FUNC
 864 #define ABI_FUNC 'rhohxcpositron'
 865 !End of the abilint section
 866 
 867  implicit none
 868 
 869 !Arguments ------------------------------------
 870 !scalars
 871  integer,intent(in) :: nfft,nkxc,nspden,n3xccc,paral_kgb,usexcnhat,usepaw
 872  real(dp),intent(in) :: ucvol,xc_denpos
 873  real(dp),intent(out) :: vxcavg
 874  type(electronpositron_type),pointer :: electronpositron
 875 !arrays
 876  integer,intent(in) :: ngfft(18)
 877  real(dp),intent(in) :: gprimd(3,3)
 878  real(dp),intent(in) :: nhat(nfft,nspden*usepaw),rhor(nfft,nspden),xccc3d(n3xccc)
 879  real(dp),intent(out) :: kxcapn(nfft,nkxc),strsxc(6),vhartr(nfft),vxcapn(nfft,nspden)
 880  type(MPI_type),intent(in) :: mpi_enreg
 881 
 882 !Local variables-------------------------------
 883 !scalars
 884  integer :: cplex,ierr,ifft,ishift,iwarn,iwarnp,nfftot,ngr,ngrad,nspden_ep
 885  real(dp) :: exc,excdc,strdiag
 886  character(len=500) :: message
 887 !arrays
 888  real(dp),parameter :: qphon(3)=(/0._dp,0._dp,0._dp/)
 889  real(dp) :: vxcavg_tmp(1)
 890  real(dp),allocatable :: fxcapn(:),grho2apn(:),rhoe(:,:,:),rhop(:,:),rhotote(:),vxc_ep(:),vxcgr_ep(:)
 891 
 892 ! *************************************************************************
 893 
 894  if (electronpositron_calctype(electronpositron)/=1) then
 895    message = 'Only electronpositron%calctype=1 allowed !'
 896    MSG_BUG(message)
 897  end if
 898 
 899  if (nkxc>3) then
 900    message = 'nkxc>3 (Kxc for GGA) not yet implemented !'
 901    MSG_ERROR(message)
 902  end if
 903 
 904 !Hartree potential of the positron is zero
 905  vhartr=zero
 906 
 907 !Some allocations/inits
 908  ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
 909  ngr=0;if (ngrad==2) ngr=nfft
 910  ABI_ALLOCATE(fxcapn,(nfft))
 911  ABI_ALLOCATE(grho2apn,(ngr))
 912  nspden_ep=1;cplex=1;ishift=0
 913  iwarn=0;iwarnp=1
 914 
 915 !Compute total electronic density
 916  ABI_ALLOCATE(rhotote,(nfft))
 917  rhotote(:)=electronpositron%rhor_ep(:,1)
 918  if (n3xccc>0) rhotote(:)=rhotote(:)+xccc3d(:)
 919  if (usepaw==1.and.usexcnhat==0) rhotote(:)=rhotote(:)-electronpositron%nhat_ep(:,1)
 920 
 921 !Extra total electron/positron densities; compute gradients for GGA
 922  ABI_ALLOCATE(rhoe,(nfft,nspden_ep,ngrad**2))
 923  ABI_ALLOCATE(rhop,(nfft,nspden_ep))
 924  call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,paral_kgb,qphon,rhotote,rhoe)
 925  if (ngrad==2) grho2apn(:)=rhoe(:,1,2)**2+rhoe(:,1,3)**2+rhoe(:,1,4)**2
 926  rhop(:,1)=rhor(:,1);if (usepaw==1.and.usexcnhat==0) rhop(:,1)=rhop(:,1)-nhat(:,1)
 927  ABI_DEALLOCATE(rhotote)
 928 
 929 !Make the densities positive
 930  call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe(:,1,1),xc_denpos)
 931  if (.not.electronpositron%posdensity0_limit) then
 932    call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,xc_denpos)
 933  end if
 934 
 935 !Compute electron-positron Vxc_pos, Vxc_el, Fxc, Kxc, ...
 936  ABI_ALLOCATE(vxc_ep,(nfft))
 937  ABI_ALLOCATE(vxcgr_ep,(ngr))
 938  if (nkxc==0) then
 939    call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,&
 940 &   rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn)
 941  else
 942    call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,&
 943 &   rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn,dvxce=kxcapn)
 944  end if
 945  ABI_DEALLOCATE(rhoe)
 946  ABI_DEALLOCATE(vxc_ep)
 947  ABI_DEALLOCATE(vxcgr_ep)
 948  ABI_DEALLOCATE(grho2apn)
 949 
 950 !Store Vxc and Kxc according to spin components
 951  if (nspden>=2) vxcapn(:,2)=vxcapn(:,1)
 952  if (nspden==4) vxcapn(:,3:4)=zero
 953  if (nkxc==3) then
 954    kxcapn(:,1)=two*kxcapn(:,1)
 955    kxcapn(:,2)=kxcapn(:,1)
 956    kxcapn(:,3)=kxcapn(:,1)
 957  end if
 958 
 959 !Compute XC energies and contribution to stress tensor
 960  electronpositron%e_xc  =zero
 961  electronpositron%e_xcdc=zero
 962  strdiag=zero
 963  nfftot=PRODUCT(ngfft(1:3))
 964  do ifft=1,nfft
 965    electronpositron%e_xc  =electronpositron%e_xc  +fxcapn(ifft)
 966    electronpositron%e_xcdc=electronpositron%e_xcdc+vxcapn(ifft,1)*rhor(ifft,1)
 967 !  strdiag=strdiag+fxcapn(ifft)   ! Already stored in rhotoxc !
 968    strdiag=strdiag-vxcapn(ifft,1)*rhop(ifft,1)
 969  end do
 970  if (usepaw==1.and.usexcnhat==0) then
 971    do ifft=1,nfft
 972      electronpositron%e_xcdc=electronpositron%e_xcdc-vxcapn(ifft,1)*nhat(ifft,1)
 973    end do
 974  end if
 975  electronpositron%e_xc  =electronpositron%e_xc  *ucvol/dble(nfftot)
 976  electronpositron%e_xcdc=electronpositron%e_xcdc*ucvol/dble(nfftot)
 977  strdiag=strdiag/dble(nfftot)
 978  ABI_DEALLOCATE(fxcapn)
 979  ABI_DEALLOCATE(rhop)
 980 
 981 !Reduction in case of parallelism
 982  if(mpi_enreg%paral_kgb==1)then
 983    if(paral_kgb/=0)then
 984      exc=electronpositron%e_xc;excdc=electronpositron%e_xcdc
 985      call xmpi_sum(exc  ,mpi_enreg%comm_fft,ierr)
 986      call xmpi_sum(excdc,mpi_enreg%comm_fft,ierr)
 987      electronpositron%e_xc=exc;electronpositron%e_xcdc=excdc
 988      call xmpi_sum(strsxc,mpi_enreg%comm_fft,ierr)
 989    end if
 990  end if
 991 
 992 !Store stress tensor
 993  strsxc(1:3)=strdiag
 994  strsxc(4:6)=zero
 995 
 996 !Compute vxcavg
 997  call mean_fftr(vxcapn(:,1),vxcavg_tmp,nfft,nfftot,1)
 998  vxcavg=vxcavg_tmp(1)
 999 
1000 end subroutine rhohxcpositron

m_electronpositron/destroy_electronpositron [ Functions ]

[ Top ] [ m_electronpositron ] [ Functions ]

NAME

  destroy_electronpositron

FUNCTION

  Clean and destroy electronpositron datastructure

SIDE EFFECTS

  electronpositron=<type(electronpositron_type)>=electronpositron datastructure

PARENTS

      gstate

CHILDREN

      energies_copy,fourdp,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawrhoij_alloc,pawrhoij_copy,pawrhoij_free

SOURCE

381 subroutine destroy_electronpositron(electronpositron)
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 'destroy_electronpositron'
388 !End of the abilint section
389 
390  implicit none
391 
392 !Arguments ------------------------------------
393 !scalars
394  type(electronpositron_type),pointer :: electronpositron
395 
396 !************************************************************************
397 
398  !@electronpositron_type
399 
400  if (associated(electronpositron)) then
401 
402   if (allocated(electronpositron%cg_ep))        then
403     ABI_DEALLOCATE(electronpositron%cg_ep)
404   end if
405   if (allocated(electronpositron%eigen_ep))     then
406     ABI_DEALLOCATE(electronpositron%eigen_ep)
407   end if
408   if (allocated(electronpositron%occ_ep))       then
409     ABI_DEALLOCATE(electronpositron%occ_ep)
410   end if
411   if (allocated(electronpositron%rhor_ep))      then
412     ABI_DEALLOCATE(electronpositron%rhor_ep)
413   end if
414   if (allocated(electronpositron%nhat_ep))      then
415     ABI_DEALLOCATE(electronpositron%nhat_ep)
416   end if
417   if (allocated(electronpositron%vha_ep))       then
418     ABI_DEALLOCATE(electronpositron%vha_ep)
419   end if
420   if (allocated(electronpositron%lmselect_ep))  then
421     ABI_DEALLOCATE(electronpositron%lmselect_ep)
422   end if
423   if (allocated(electronpositron%fred_ep))      then
424     ABI_DEALLOCATE(electronpositron%fred_ep)
425   end if
426   if (allocated(electronpositron%stress_ep))    then
427     ABI_DEALLOCATE(electronpositron%stress_ep)
428   end if
429 
430   if (electronpositron%has_pawrhoij_ep/=0) then
431    call pawrhoij_free(electronpositron%pawrhoij_ep)
432   end if
433   if (allocated(electronpositron%pawrhoij_ep))  then
434     ABI_DATATYPE_DEALLOCATE(electronpositron%pawrhoij_ep)
435   end if
436 
437   if (electronpositron%dimcprj/=0) then
438    call pawcprj_free(electronpositron%cprj_ep)
439   end if
440   if (allocated(electronpositron%cprj_ep))  then
441     ABI_DATATYPE_DEALLOCATE(electronpositron%cprj_ep)
442   end if
443 
444   electronpositron%calctype       =0
445   electronpositron%particle       =-1
446   electronpositron%dimcg          =0
447   electronpositron%dimcprj        =0
448   electronpositron%dimeigen       =0
449   electronpositron%dimocc         =0
450   electronpositron%has_pawrhoij_ep=0
451   electronpositron%has_pos_ham    =0
452   electronpositron%istep          =0
453   electronpositron%istep_scf      =0
454 
455   electronpositron%posdensity0_limit=.false.
456   electronpositron%scf_converged=.false.
457 
458   ABI_DATATYPE_DEALLOCATE(electronpositron)
459 
460  end if
461 
462 end subroutine destroy_electronpositron

m_electronpositron/electronpositron_calctype [ Functions ]

[ Top ] [ m_electronpositron ] [ Functions ]

NAME

  electronpositron_calctype

FUNCTION

  Returns the value of the calculation type from an electronpositron
  structure (can be eventually unassociated)

INPUTS

  electronpositron=<type(electronpositron_type)>=electronpositron datastructure

PARENTS

CHILDREN

SOURCE

784 integer function electronpositron_calctype(electronpositron)
785 
786 
787 !This section has been created automatically by the script Abilint (TD).
788 !Do not modify the following lines by hand.
789 #undef ABI_FUNC
790 #define ABI_FUNC 'electronpositron_calctype'
791 !End of the abilint section
792 
793  implicit none
794 
795 !Arguments ------------------------------------
796 !scalars
797  type(electronpositron_type),pointer :: electronpositron
798 
799 !************************************************************************
800 
801  if (associated(electronpositron)) then
802   electronpositron_calctype=electronpositron%calctype
803  else
804   electronpositron_calctype=0
805  end if
806 
807 
808 end function electronpositron_calctype

m_electronpositron/electronpositron_type [ Types ]

[ Top ] [ m_electronpositron ] [ Types ]

NAME

FUNCTION

NOTES

SOURCE

 67  type, public :: electronpositron_type
 68 
 69 ! Integer scalars
 70   integer :: calctype        ! type of electron-positron calculation:
 71                              !   0: no calculation
 72                              !   1: positron in the electrons potential
 73                              !   2: electrons in the positron potential
 74   integer :: particle        ! current particle stored in electronpositron%xxx_ep arrays
 75                              !                 -1: no particle, 0: electron, 1: positron
 76   integer :: dimcg           ! Dimension of cg array dimcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
 77   integer :: dimcprj         ! Dimension of cprj array dimcprj=dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol*usecprj
 78   integer :: dimeigen        ! Dimension of eigen array dimeigen=dtset%mband*dtset%nkpt*dtset%nsppol
 79   integer :: dimocc          ! Dimension of occ array dimocc=dtset%mband*dtset%nkpt*dtset%nsppol
 80   integer :: has_pawrhoij_ep ! flag for pawrhoij_ep (0: not allocated, 1: allocated, 2: computed)
 81   integer :: has_pos_ham     ! flag: 1 if current Hamiltonian in memory (vtrial, vpsp, vhartr, vxc, paw_ij%dij)
 82 !                                    is the positronic hamiltonian, 0 is it is the electronic one
 83   integer :: ixcpositron     ! XC type for electron-positron correlation
 84   integer :: istep           ! Current index of TC-DFT SCF step
 85   integer :: istep_scf       ! Current index of DFT SCF step  in current electron/positron minimization
 86   integer :: lmmax           ! Max. number of (l,m) moments over all types of atom
 87   integer :: natom           ! Number of atoms
 88   integer :: nfft            ! Number of points in FFT grid
 89   integer :: nspden          ! Number of spin density components
 90   integer :: nstep           ! Max. number of steps for the TC-DFT SCF cycle
 91 
 92 ! Logical scalars
 93   logical :: posdensity0_limit ! True if we are in the zero positron density limit
 94   logical :: scf_converged     ! True if the SCF cycle is converged for a positronic/electronic GS calculation
 95 
 96 ! Real(dp) scalars
 97   real(dp) :: e_hartree      !  Hartree electron-positron interaction energy
 98   real(dp) :: e_xc           !  XC electron-positron interaction energy
 99   real(dp) :: e_xcdc         !  Double-counting XC electron-positron interaction energy
100   real(dp) :: e_paw          !  PAW electron-positron interaction energy
101   real(dp) :: e_pawdc        !  Double-counting PAW electron-positron interaction energy
102   real(dp) :: e0             !  Energy only due to particle(s) currently evolving
103                                   !   calctype=1, energy due to positron  only
104                                   !   calctype=2, energy due to electrons only
105   real(dp) :: etotal_prev    !  Total energy of the previous GS calculation
106   real(dp) :: lambda         ! Electron-positron annihilation rate
107   real(dp) :: lifetime       ! Positron lifetime
108   real(dp) :: maxfor_prev    ! Max. force of the previous GS calculation
109   real(dp) :: posocc         ! Occupation number for the positron
110   real(dp) :: postoldfe      ! Tolerance on total energy for the TC-DFT SCF cycle
111   real(dp) :: postoldff      ! Tolerance on max. force for the TC-DFT SCF cycle
112 
113 ! Other scalars
114   type(energies_type) :: energies_ep  !  Energies of the previous electronic/positronic SCF step
115 
116 ! Logical pointers
117   logical, allocatable :: lmselect_ep(:,:)
118 !  lmselect_ep(lmmax,my_natom)
119 !  flags selecting the non-zero LM-moments of on-site densities
120 
121 ! Real(dp) pointers
122   real(dp), allocatable :: cg_ep(:,:)
123 !  cg_ep(2,dimcg)
124 !  if typecalc=1: electronic wavefunctions
125 !  if typecalc=2: positronic wavefunctions
126 
127   real(dp), allocatable :: eigen_ep(:)
128 !  eigen(dimeigen)
129 !  if typecalc=1: electronic eigen energies
130 !  if typecalc=2: positronic eigen energies
131 
132   real(dp), allocatable :: fred_ep(:,:)
133 !  fred_ep(3,natom)
134 !  if typecalc=1: forces only due to electrons
135 !  if typecalc=2: forces only due to positron
136 
137   real(dp), allocatable :: nhat_ep(:,:)
138 !  nhat_ep(nfft,nspden)
139 !  if typecalc=1: electronic compensation charge density in real space
140 !  if typecalc=2: positronic compensation charge density in real space
141 
142   real(dp), allocatable :: occ_ep(:)
143 !  occ(dimocc)
144 !  if typecalc=1: electronic occupations
145 !  if typecalc=2: positronic occupations
146 
147   real(dp), allocatable :: rhor_ep(:,:)
148 !  rhor_ep(nfft,nspden)
149 !  if typecalc=1: electronic density in real space
150 !  if typecalc=2: positronic density in real space
151 
152   real(dp), allocatable :: stress_ep(:)
153 !  stress_ep(6)
154 !  if typecalc=1: stresses only due to electrons
155 !  if typecalc=2: stresses only due to positron
156 
157   real(dp), allocatable :: vha_ep(:)
158 !  vha_ep(nfft)
159 !  if typecalc=1: electronic Hartree potential
160 !  if typecalc=2: positronic Hartree potential
161 
162 ! Other pointers
163   type(pawrhoij_type), allocatable :: pawrhoij_ep(:)
164 !  pawrhoij_ep(natom)
165 !  Relevant only if PAW
166 !  if typecalc=1: electronic PAW occupation matrix associated with rhor_ep
167 !  if typecalc=2: positronic PAW occupation matrix associated with rhor_ep
168 
169   type(pawcprj_type), allocatable :: cprj_ep(:,:)
170 !  cprj_ep(natom,dimcprj)
171 !  Relevant only if PAW
172 !  if typecalc=1: electronic WF projected on nl projectors <p_i|Cnk>
173 !  if typecalc=2: positronic WF projected on nl projectors <p_i|Cnk>
174 
175  end type electronpositron_type
176 
177 ! public procedures
178  public :: init_electronpositron
179  public :: destroy_electronpositron
180  public :: exchange_electronpositron
181  public :: electronpositron_calctype
182  public :: rhohxcpositron
183 
184 CONTAINS
185 
186 !===========================================================

m_electronpositron/exchange_electronpositron [ Functions ]

[ Top ] [ m_electronpositron ] [ Functions ]

NAME

  exchange_electronpositron

FUNCTION

  Invert electron and positron quantities between an electronpositron datastructure
  and current evoving variables
  Example: exchange electronpositron%rhor_ep and rhor

INPUTS

  dtset <type(dataset_type)>=all input variables for this dataset
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current proc
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT
  npwarr(nkpt)=number of planewaves in basis at this k point
  usecprj= 1 if cprj array is stored in memory

SIDE EFFECTS

  cg(2,mcg)=wavefunctions
  cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
                             cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
  electronpositron=<type(electronpositron_type)>=electronpositron datastructure
  energies <type(energies_type)>=all part of total energy.
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  fred(3,natom)=forces in reduced coordinates
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  occ(mband*nkpt*nsppol)=occupation number for each band at each k point
  paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh
  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  rhog(2,nfft)=Fourier transform of total electron/positron density
  rhor(nfft,nspden)=total electron/positron density (el/bohr**3)
  stress(6)=components of the stress tensor (hartree/bohr^3) for the
  vhartr(nfftf)=array for holding Hartree potential

PARENTS

      afterscfloop

CHILDREN

      energies_copy,fourdp,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawrhoij_alloc,pawrhoij_copy,pawrhoij_free

SOURCE

512 subroutine exchange_electronpositron(cg,cprj,dtset,eigen,electronpositron,energies,fred,mcg,mcprj,&
513 &                                    mpi_enreg,my_natom,nfft,ngfft,nhat,npwarr,occ,paw_an,pawrhoij,&
514 &                                    rhog,rhor,stress,usecprj,vhartr)
515 
516 
517 !This section has been created automatically by the script Abilint (TD).
518 !Do not modify the following lines by hand.
519 #undef ABI_FUNC
520 #define ABI_FUNC 'exchange_electronpositron'
521 !End of the abilint section
522 
523  implicit none
524 
525 !Arguments ------------------------------------
526 !scalars
527  integer,intent(in) :: mcg,mcprj,my_natom,nfft,usecprj
528  type(dataset_type),intent(in) :: dtset
529  type(electronpositron_type),pointer :: electronpositron
530  type(energies_type),intent(inout) :: energies
531  type(MPI_type),intent(in) :: mpi_enreg
532 !arrays
533  integer,intent(in) :: ngfft(18),npwarr(dtset%nkpt)
534  real(dp),intent(inout) :: cg(2,mcg)
535  real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
536  real(dp),intent(inout) :: fred(3,dtset%natom),nhat(nfft,dtset%nspden)
537  real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
538  real(dp), intent(inout) :: rhog(2,nfft),rhor(nfft,dtset%nspden)
539  real(dp),intent(inout) :: stress(6),vhartr(nfft)
540  type(pawcprj_type) :: cprj(dtset%natom,mcprj*usecprj)
541  type(paw_an_type),intent(inout) :: paw_an(my_natom*dtset%usepaw)
542  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*dtset%usepaw)
543 
544 !Local variables-------------------------------
545 !scalars
546  integer :: comm,iatom,ib,ibsp,icg,icgb,ifft,ii,ilm,ikpt
547  integer :: ispden,isppol,ispinor,me,my_nspinor,nband_k,npw_k,sz1,sz2,sz3
548  logical :: ltmp
549  real(dp) :: rtmp
550  type(energies_type) :: energies_tmp
551 !arrays
552  integer,allocatable :: nlmn(:),typ(:)
553  real(dp) :: ctmp(2)
554  type(pawcprj_type),allocatable :: cprj_tmp(:,:)
555  type(pawrhoij_type),allocatable :: pawrhoij_tmp(:)
556 
557 !*********************************************************************
558 
559  if (associated(electronpositron)) then
560   if (electronpositron%particle/=EP_NOTHING) then
561 
562 !  Type of particle stored
563    if (electronpositron%particle==EP_ELECTRON) then
564      electronpositron%particle=EP_POSITRON
565    else if (electronpositron%particle==EP_POSITRON) then
566      electronpositron%particle=EP_ELECTRON
567    end if
568 
569 !  Energies
570    ctmp(1)=energies%e_electronpositron
571 !  ctmp(2)=energies%edc_electronpositron
572    call energies_copy(electronpositron%energies_ep,energies_tmp)
573    call energies_copy(energies,electronpositron%energies_ep)
574    call energies_copy(energies_tmp,energies)
575    energies%e_electronpositron=ctmp(1)
576 !  energies%edc_electronpositron=ctmp(2)
577    energies%e0_electronpositron=electronpositron%e0
578    electronpositron%e0=electronpositron%energies_ep%e0_electronpositron
579 
580 !  Density and PAW occupation matrix
581    do ispden=1,dtset%nspden
582      do ifft=1,nfft
583        rtmp=rhor(ifft,ispden)
584        rhor(ifft,ispden)=electronpositron%rhor_ep(ifft,ispden)
585        electronpositron%rhor_ep(ifft,ispden)=rtmp
586      end do
587      if (allocated(electronpositron%nhat_ep).and.size(nhat,2)>0) then
588        do ifft=1,nfft
589          rtmp=nhat(ifft,ispden)
590          nhat(ifft,ispden)=electronpositron%nhat_ep(ifft,ispden)
591          electronpositron%nhat_ep(ifft,ispden)=rtmp
592        end do
593      end if
594    end do
595    call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0)
596    if (dtset%usepaw==1.and.my_natom>0) then
597     if (electronpositron%has_pawrhoij_ep==1) then
598       ABI_DATATYPE_ALLOCATE(pawrhoij_tmp,(my_natom))
599       ABI_ALLOCATE(typ,(my_natom))
600       ABI_ALLOCATE(nlmn,(my_natom))
601       do iatom=1,my_natom
602         typ(iatom)=iatom
603         nlmn(iatom)=pawrhoij(iatom)%lmn_size
604       end do
605 !     Be careful: parallelism over atoms is ignored...
606       call pawrhoij_alloc(pawrhoij_tmp,pawrhoij(1)%cplex,pawrhoij(1)%nspden,&
607 &                      pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,typ, &
608 &                      lmnsize=nlmn,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
609 &                      use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres)
610       ABI_DEALLOCATE(typ)
611       ABI_DEALLOCATE(nlmn)
612       call pawrhoij_copy(pawrhoij,pawrhoij_tmp)
613       call pawrhoij_copy(electronpositron%pawrhoij_ep,pawrhoij)
614       call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep)
615       if (pawrhoij_tmp(1)%ngrhoij>0.and.pawrhoij(1)%ngrhoij==0) then
616         do iatom=1,my_natom
617           sz1=pawrhoij_tmp(iatom)%ngrhoij
618           sz2=pawrhoij_tmp(iatom)%cplex*pawrhoij_tmp(iatom)%lmn2_size
619           sz3=pawrhoij_tmp(iatom)%nspden
620           ABI_ALLOCATE(pawrhoij(iatom)%grhoij,(sz1,sz2,sz3))
621           pawrhoij(iatom)%grhoij(:,:,:)=pawrhoij_tmp(iatom)%grhoij(:,:,:)
622         end do
623       end if
624       if (pawrhoij_tmp(1)%use_rhoijres>0.and.pawrhoij(1)%use_rhoijres==0) then
625         do iatom=1,my_natom
626           sz1=pawrhoij_tmp(iatom)%cplex*pawrhoij_tmp(iatom)%lmn2_size
627           sz2=pawrhoij_tmp(iatom)%nspden
628           ABI_ALLOCATE(pawrhoij(iatom)%rhoijres,(sz1,sz2))
629           pawrhoij(iatom)%rhoijres(:,:)=pawrhoij_tmp(iatom)%rhoijres(:,:)
630         end do
631       end if
632       if (pawrhoij_tmp(1)%use_rhoij_>0.and.pawrhoij(1)%use_rhoij_==0) then
633         do iatom=1,my_natom
634           sz1=pawrhoij_tmp(iatom)%cplex*pawrhoij_tmp(iatom)%lmn2_size
635           sz2=pawrhoij_tmp(iatom)%nspden
636           ABI_ALLOCATE(pawrhoij(iatom)%rhoij_,(sz1,sz2))
637           pawrhoij(iatom)%rhoij_(:,:)=pawrhoij_tmp(iatom)%rhoij_(:,:)
638         end do
639       end if
640       if (pawrhoij_tmp(1)%lmnmix_sz>0.and.pawrhoij(1)%lmnmix_sz==0) then
641         do iatom=1,my_natom
642           ABI_ALLOCATE(pawrhoij(iatom)%kpawmix,(pawrhoij_tmp(iatom)%lmnmix_sz))
643           pawrhoij(iatom)%kpawmix(:)=pawrhoij_tmp(iatom)%kpawmix(:)
644         end do
645       end if
646       call pawrhoij_free(pawrhoij_tmp)
647       ABI_DATATYPE_DEALLOCATE(pawrhoij_tmp)
648     else
649       do iatom=1,my_natom
650         pawrhoij(iatom)%rhoijp=zero
651       end do
652     end if
653    end if
654 
655 !  Hartree potential
656    do ifft=1,nfft
657     rtmp=vhartr(ifft)
658     vhartr(ifft)=electronpositron%vha_ep(ifft)
659     electronpositron%vha_ep(ifft)=rtmp
660    end do
661 
662 !  PAW LM-moment selection flags
663    if (dtset%usepaw==1.and.my_natom>0) then
664     if (electronpositron%lmmax>0) then
665      do iatom=1,my_natom
666       do ilm=1,paw_an(iatom)%lm_size
667        ltmp=electronpositron%lmselect_ep(ilm,iatom)
668        electronpositron%lmselect_ep(ilm,iatom)=paw_an(iatom)%lmselect(ilm)
669        paw_an(iatom)%lmselect(ilm)=ltmp
670       end do
671      end do
672     else
673      do iatom=1,my_natom
674       paw_an(iatom)%lmselect(:)=.true.
675      end do
676     end if
677    end if
678 
679 !  Wave-functions
680    if (electronpositron%dimcg>0) then
681     do ii=1,electronpositron%dimcg
682      ctmp(1:2)=electronpositron%cg_ep(1:2,ii)
683      electronpositron%cg_ep(1:2,ii)=cg(1:2,ii)
684      cg(1:2,ii)=ctmp(1:2)
685     end do
686    else
687     icg=0
688     my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
689     comm=mpi_enreg%comm_cell
690     me=xmpi_comm_rank(comm)
691     do isppol=1,dtset%nsppol
692      do ikpt=1,dtset%nkpt
693       npw_k=npwarr(ikpt);nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
694       if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
695       icgb=icg;ibsp=0
696       do ib=1,nband_k
697        cg(:,icgb+1:icgb+my_nspinor*npw_k)=zero
698        do ispinor=1,my_nspinor
699         ibsp=ibsp+1;if (ibsp<my_nspinor*npw_k) cg(1,icgb+ibsp)=one
700        end do
701        icgb=icgb+my_nspinor*npw_k
702       end do
703       if (dtset%mkmem/=0) icg=icg+my_nspinor*npw_k*nband_k
704      end do
705     end do
706    end if
707    if (dtset%usepaw==1) then
708     if(electronpositron%dimcprj>0) then
709      ABI_ALLOCATE(nlmn,(dtset%natom))
710      ABI_DATATYPE_ALLOCATE(cprj_tmp,(dtset%natom,electronpositron%dimcprj))
711      do iatom=1,dtset%natom;nlmn(iatom)=cprj(iatom,1)%nlmn;end do
712      call pawcprj_alloc(cprj_tmp,cprj(1,1)%ncpgr,nlmn)
713      ABI_DEALLOCATE(nlmn)
714      call pawcprj_copy(electronpositron%cprj_ep,cprj_tmp)
715      call pawcprj_copy(cprj,electronpositron%cprj_ep)
716      call pawcprj_copy(cprj_tmp,cprj)
717      call pawcprj_free(cprj_tmp)
718      ABI_DATATYPE_DEALLOCATE(cprj_tmp)
719     else
720 !TO BE ACTIVATED WHEN cprj IS PRESENT
721 !    call pawcprj_set_zero(cprj)
722     end if
723    end if
724 
725 !  Eigenvalues
726    if (electronpositron%dimeigen>0) then
727     do ii=1,electronpositron%dimeigen
728      rtmp=eigen(ii)
729      eigen(ii)=electronpositron%eigen_ep(ii)
730      electronpositron%eigen_ep(ii)=rtmp
731     end do
732    else
733     eigen(:)=9.99999_dp
734    end if
735 
736 !  Occupations
737    if (electronpositron%dimocc>0) then
738     do ii=1,electronpositron%dimocc
739      rtmp=occ(ii)
740      occ(ii)=electronpositron%occ_ep(ii)
741      electronpositron%occ_ep(ii)=rtmp
742     end do
743    else
744     occ(:)=9.99999_dp
745    end if
746 
747 !  Forces
748    if (allocated(electronpositron%fred_ep)) then
749     do iatom=1,dtset%natom
750      electronpositron%fred_ep(1:3,iatom)=fred(1:3,iatom)-electronpositron%fred_ep(1:3,iatom)
751     end do
752    end if
753 
754 !  Stresses
755    if (allocated(electronpositron%stress_ep)) then
756     electronpositron%stress_ep(1:6)=stress(1:6)-electronpositron%stress_ep(1:6)
757    end if
758 
759   end if
760  end if
761 
762 end subroutine exchange_electronpositron

m_electronpositron/init_electronpositron [ Functions ]

[ Top ] [ m_electronpositron ] [ Functions ]

NAME

  init_electronpositron

FUNCTION

  Init all scalars and pointers in the structure.

INPUTS

  ireadwf=if 1, read the wavefunction
  dtset <type(dataset_type)>=all input variables for this dataset
  mpi_enreg=information about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  pawrhoij(natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data

SIDE EFFECTS

  electronpositron=<type(electronpositron_type)>=electronpositron datastructure

PARENTS

      gstate

CHILDREN

      energies_copy,fourdp,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawrhoij_alloc,pawrhoij_copy,pawrhoij_free

SOURCE

216 subroutine init_electronpositron(ireadwf,dtset,electronpositron,mpi_enreg,nfft,pawrhoij,pawtab)
217 
218 
219 !This section has been created automatically by the script Abilint (TD).
220 !Do not modify the following lines by hand.
221 #undef ABI_FUNC
222 #define ABI_FUNC 'init_electronpositron'
223 !End of the abilint section
224 
225  implicit none
226 
227 !Arguments ------------------------------------
228 !scalars
229  integer,intent(in) :: ireadwf,nfft
230  type(dataset_type),intent(in) :: dtset
231  type(electronpositron_type),pointer :: electronpositron
232  type(MPI_type),intent(in) :: mpi_enreg
233 !arrays
234  type(pawrhoij_type), intent(in) :: pawrhoij(mpi_enreg%my_natom*dtset%usepaw)
235  type(pawtab_type),intent(in)  :: pawtab(dtset%ntypat*dtset%usepaw)
236 
237 !Local variables-------------------------------
238 !scalars
239  integer :: ii,my_nspinor,ncpgr,optfor,optstr
240  logical,parameter :: include_nhat_in_gamma=.false.
241 !arrays
242  integer,allocatable :: nlmn(:)
243 
244 !************************************************************************
245 
246  !@electronpositron_type
247 
248  if (dtset%positron/=0) then
249 
250   ABI_DATATYPE_ALLOCATE(electronpositron,)
251 
252   electronpositron%calctype=0
253   electronpositron%particle=-1
254 
255   electronpositron%ixcpositron=dtset%ixcpositron
256   electronpositron%natom=dtset%natom
257   electronpositron%nfft=nfft
258   electronpositron%nspden=dtset%nspden
259   electronpositron%istep=0
260   electronpositron%istep_scf=0
261 
262   electronpositron%posocc=dtset%posocc
263   electronpositron%nstep=dtset%posnstep
264   electronpositron%postoldfe=dtset%postoldfe
265   electronpositron%postoldff=dtset%postoldff
266   electronpositron%posdensity0_limit=(dtset%ixcpositron/=2)
267   electronpositron%scf_converged=.false.
268   electronpositron%has_pos_ham=0
269 
270   call energies_init(electronpositron%energies_ep)
271 
272   electronpositron%e_hartree  =zero
273   electronpositron%e_xc       =zero
274   electronpositron%e_xcdc     =zero
275   electronpositron%e_paw      =zero
276   electronpositron%e_pawdc    =zero
277   electronpositron%e0         =zero
278   electronpositron%etotal_prev=zero
279   electronpositron%maxfor_prev=zero
280 
281   electronpositron%lambda=zero
282   electronpositron%lifetime=zero
283 
284   ABI_ALLOCATE(electronpositron%rhor_ep,(nfft,dtset%nspden))
285   ABI_ALLOCATE(electronpositron%vha_ep,(nfft))
286   ABI_DATATYPE_ALLOCATE(electronpositron%pawrhoij_ep,(mpi_enreg%my_natom*dtset%usepaw))
287 
288   if (dtset%usepaw==1) then
289    electronpositron%has_pawrhoij_ep=1
290    if (mpi_enreg%my_natom>0) then
291     call pawrhoij_alloc(electronpositron%pawrhoij_ep,pawrhoij(1)%cplex,pawrhoij(1)%nspden,&
292 &                    pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,dtset%typat,&
293 &                    mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,&
294 &                    pawtab=pawtab,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,&
295 &                    use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres)
296    end if
297    electronpositron%lmmax=0
298    do ii=1,dtset%ntypat
299     electronpositron%lmmax=max(electronpositron%lmmax,pawtab(ii)%lcut_size**2)
300    end do
301    ABI_ALLOCATE(electronpositron%lmselect_ep,(electronpositron%lmmax,mpi_enreg%my_natom))
302    if (maxval(pawtab(1:dtset%ntypat)%usexcnhat)==0.or.(.not.include_nhat_in_gamma)) then
303      ABI_ALLOCATE(electronpositron%nhat_ep,(nfft,dtset%nspden))
304    end if
305   else
306    electronpositron%has_pawrhoij_ep=0
307    electronpositron%lmmax=0
308   end if
309 
310   if (dtset%positron<=-10.or.dtset%posdoppler>0) then
311    my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
312    electronpositron%dimcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
313    electronpositron%dimocc=dtset%mband*dtset%nkpt*dtset%nsppol
314    electronpositron%dimeigen=dtset%mband*dtset%nkpt*dtset%nsppol
315    ABI_ALLOCATE(electronpositron%cg_ep,(2,electronpositron%dimcg))
316    ABI_ALLOCATE(electronpositron%eigen_ep,(electronpositron%dimeigen))
317    ABI_ALLOCATE(electronpositron%occ_ep,(electronpositron%dimocc))
318    electronpositron%dimcprj=0
319 !  if (.false.) then !TEMPORARY: will be activated later
320    if (dtset%usepaw==1.and.dtset%pawusecp>0.and.dtset%posdoppler>0) then
321     electronpositron%dimcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
322     if (mpi_enreg%paral_kgb/=0) electronpositron%dimcprj=electronpositron%dimcprj/mpi_enreg%nproc_band
323     ABI_DATATYPE_ALLOCATE(electronpositron%cprj_ep,(dtset%natom,electronpositron%dimcprj))
324     ABI_ALLOCATE(nlmn,(dtset%natom))
325     ncpgr=0
326     do ii=1,dtset%natom;nlmn(ii)=pawtab(dtset%typat(ii))%lmn_size;end do
327     call pawcprj_alloc(electronpositron%cprj_ep,ncpgr,nlmn)
328     ABI_DEALLOCATE(nlmn)
329    else
330     ABI_DATATYPE_ALLOCATE(electronpositron%cprj_ep,(dtset%natom,electronpositron%dimcprj))
331    end if
332   else
333    electronpositron%dimcg   =0
334    electronpositron%dimcprj =0
335    electronpositron%dimeigen=0
336    electronpositron%dimocc  =0
337   end if
338 
339   optfor=0;optstr=0
340   if ((dtset%optforces>0.or.dtset%ionmov/=0.or.abs(dtset%toldff)>tiny(0._dp))) optfor=1
341   if (dtset%optstress>0.and.dtset%iscf>0.and.(dtset%nstep>0.or.ireadwf==1)) optstr=1
342 
343   if (optfor>0) then
344    ABI_ALLOCATE(electronpositron%fred_ep,(3,dtset%natom))
345    electronpositron%fred_ep(:,:)=zero
346   end if
347 
348   if (optstr>0) then
349    ABI_ALLOCATE(electronpositron%stress_ep,(6))
350    electronpositron%stress_ep(:)=zero
351   end if
352 
353  else !dtset%positron==0
354   nullify(electronpositron)
355  end if
356 
357 end subroutine init_electronpositron