TABLE OF CONTENTS
- ABINIT/m_electronpositron
- ABINIT/rhohxcpositron
- m_electronpositron/destroy_electronpositron
- m_electronpositron/electronpositron_calctype
- m_electronpositron/electronpositron_type
- m_electronpositron/exchange_electronpositron
- m_electronpositron/init_electronpositron
ABINIT/m_electronpositron [ 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-2024 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 .
SOURCE
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 #include "abi_common.h" 22 23 MODULE m_electronpositron 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_energies 29 use m_xmpi 30 use m_cgtools 31 use m_dtset 32 33 use defs_abitypes, only : MPI_type 34 use m_pawtab, only : pawtab_type 35 use m_paw_an, only : paw_an_type 36 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy 37 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 38 use m_mpinfo, only : proc_distrb_cycle 39 use m_xcpositron, only : xcpositron 40 use m_drivexc, only : mkdenpos 41 use m_xctk, only : xcden 42 use m_fft, only : fourdp 43 44 implicit none 45 46 private 47 48 ! public constants 49 integer,public,parameter :: EP_NOTHING =-1 50 integer,public,parameter :: EP_ELECTRON = 0 51 integer,public,parameter :: EP_POSITRON = 1
ABINIT/rhohxcpositron [ 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
SOURCE
788 subroutine rhohxcpositron(electronpositron,gprimd,kxcapn,mpi_enreg,nfft,ngfft,nhat,nkxc,nspden,n3xccc,& 789 & paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxcapn,vxcavg,xccc3d,xc_denpos) 790 791 !Arguments ------------------------------------ 792 !scalars 793 integer,intent(in) :: nfft,nkxc,nspden,n3xccc,paral_kgb,usexcnhat,usepaw 794 real(dp),intent(in) :: ucvol,xc_denpos 795 real(dp),intent(out) :: vxcavg 796 type(electronpositron_type),pointer :: electronpositron 797 !arrays 798 integer,intent(in) :: ngfft(18) 799 real(dp),intent(in) :: gprimd(3,3) 800 real(dp),intent(in) :: nhat(nfft,nspden*usepaw),rhor(nfft,nspden),xccc3d(n3xccc) 801 real(dp),intent(out) :: kxcapn(nfft,nkxc),strsxc(6),vhartr(nfft),vxcapn(nfft,nspden) 802 type(MPI_type),intent(in) :: mpi_enreg 803 804 !Local variables------------------------------- 805 !scalars 806 integer :: cplex,ierr,ifft,ishift,iwarn,iwarnp,nfftot,ngr,ngrad,nspden_ep 807 real(dp) :: exc,excdc,strdiag 808 character(len=500) :: message 809 !arrays 810 real(dp),parameter :: qphon(3)=(/0._dp,0._dp,0._dp/) 811 real(dp) :: vxcavg_tmp(1) 812 real(dp),allocatable :: fxcapn(:),grho2apn(:),rhoe(:,:,:),rhop(:,:),rhotote(:),vxc_ep(:),vxcgr_ep(:) 813 814 ! ************************************************************************* 815 816 if (electronpositron_calctype(electronpositron)/=1) then 817 message = 'Only electronpositron%calctype=1 allowed !' 818 ABI_BUG(message) 819 end if 820 821 if (nkxc>3) then 822 message = 'nkxc>3 (Kxc for GGA) not yet implemented !' 823 ABI_ERROR(message) 824 end if 825 826 !Hartree potential of the positron is zero 827 vhartr=zero 828 829 !Some allocations/inits 830 ngrad=1;if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2 831 ngr=0;if (ngrad==2) ngr=nfft 832 ABI_MALLOC(fxcapn,(nfft)) 833 ABI_MALLOC(grho2apn,(ngr)) 834 nspden_ep=1;cplex=1;ishift=0 835 iwarn=0;iwarnp=1 836 837 !Compute total electronic density 838 ABI_MALLOC(rhotote,(nfft)) 839 rhotote(:)=electronpositron%rhor_ep(:,1) 840 if (n3xccc>0) rhotote(:)=rhotote(:)+xccc3d(:) 841 if (usepaw==1.and.usexcnhat==0) rhotote(:)=rhotote(:)-electronpositron%nhat_ep(:,1) 842 843 !Extra total electron/positron densities; compute gradients for GGA 844 ABI_MALLOC(rhoe,(nfft,nspden_ep,ngrad**2)) 845 ABI_MALLOC(rhop,(nfft,nspden_ep)) 846 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,qphon,rhotote,rhoe) 847 if (ngrad==2) grho2apn(:)=rhoe(:,1,2)**2+rhoe(:,1,3)**2+rhoe(:,1,4)**2 848 rhop(:,1)=rhor(:,1);if (usepaw==1.and.usexcnhat==0) rhop(:,1)=rhop(:,1)-nhat(:,1) 849 ABI_FREE(rhotote) 850 851 !Make the densities positive 852 call mkdenpos(iwarn ,nfft,nspden_ep,1,rhoe(:,1,1),xc_denpos) 853 if (.not.electronpositron%posdensity0_limit) then 854 call mkdenpos(iwarnp,nfft,nspden_ep,1,rhop,xc_denpos) 855 end if 856 857 !Compute electron-positron Vxc_pos, Vxc_el, Fxc, Kxc, ... 858 ABI_MALLOC(vxc_ep,(nfft)) 859 ABI_MALLOC(vxcgr_ep,(ngr)) 860 if (nkxc==0) then 861 call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,& 862 & rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn) 863 else 864 call xcpositron(fxcapn,grho2apn,electronpositron%ixcpositron,ngr,nfft,electronpositron%posdensity0_limit,& 865 & rhoe(:,1,1),rhop(:,1),vxc_ep,vxcgr_ep,vxcapn,dvxce=kxcapn) 866 end if 867 ABI_FREE(rhoe) 868 ABI_FREE(vxc_ep) 869 ABI_FREE(vxcgr_ep) 870 ABI_FREE(grho2apn) 871 872 !Store Vxc and Kxc according to spin components 873 if (nspden>=2) vxcapn(:,2)=vxcapn(:,1) 874 if (nspden==4) vxcapn(:,3:4)=zero 875 if (nkxc==3) then 876 kxcapn(:,1)=two*kxcapn(:,1) 877 kxcapn(:,2)=kxcapn(:,1) 878 kxcapn(:,3)=kxcapn(:,1) 879 end if 880 881 !Compute XC energies and contribution to stress tensor 882 electronpositron%e_xc =zero 883 electronpositron%e_xcdc=zero 884 strdiag=zero 885 nfftot=PRODUCT(ngfft(1:3)) 886 do ifft=1,nfft 887 electronpositron%e_xc =electronpositron%e_xc +fxcapn(ifft) 888 electronpositron%e_xcdc=electronpositron%e_xcdc+vxcapn(ifft,1)*rhor(ifft,1) 889 ! strdiag=strdiag+fxcapn(ifft) ! Already stored in rhotoxc ! 890 strdiag=strdiag-vxcapn(ifft,1)*rhop(ifft,1) 891 end do 892 if (usepaw==1.and.usexcnhat==0) then 893 do ifft=1,nfft 894 electronpositron%e_xcdc=electronpositron%e_xcdc-vxcapn(ifft,1)*nhat(ifft,1) 895 end do 896 end if 897 electronpositron%e_xc =electronpositron%e_xc *ucvol/dble(nfftot) 898 electronpositron%e_xcdc=electronpositron%e_xcdc*ucvol/dble(nfftot) 899 strdiag=strdiag/dble(nfftot) 900 ABI_FREE(fxcapn) 901 ABI_FREE(rhop) 902 903 !Reduction in case of parallelism 904 if(mpi_enreg%paral_kgb==1)then 905 if(paral_kgb/=0)then 906 exc=electronpositron%e_xc;excdc=electronpositron%e_xcdc 907 call xmpi_sum(exc ,mpi_enreg%comm_fft,ierr) 908 call xmpi_sum(excdc,mpi_enreg%comm_fft,ierr) 909 electronpositron%e_xc=exc;electronpositron%e_xcdc=excdc 910 call xmpi_sum(strsxc,mpi_enreg%comm_fft,ierr) 911 end if 912 end if 913 914 !Store stress tensor 915 strsxc(1:3)=strdiag 916 strsxc(4:6)=zero 917 918 !Compute vxcavg 919 call mean_fftr(vxcapn(:,1),vxcavg_tmp,nfft,nfftot,1) 920 vxcavg=vxcavg_tmp(1) 921 922 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
SOURCE
355 subroutine destroy_electronpositron(electronpositron) 356 357 !Arguments ------------------------------------ 358 !scalars 359 type(electronpositron_type),pointer :: electronpositron 360 361 !************************************************************************ 362 363 !@electronpositron_type 364 365 if (associated(electronpositron)) then 366 367 if (allocated(electronpositron%cg_ep)) then 368 ABI_FREE(electronpositron%cg_ep) 369 end if 370 if (allocated(electronpositron%eigen_ep)) then 371 ABI_FREE(electronpositron%eigen_ep) 372 end if 373 if (allocated(electronpositron%occ_ep)) then 374 ABI_FREE(electronpositron%occ_ep) 375 end if 376 if (allocated(electronpositron%rhor_ep)) then 377 ABI_FREE(electronpositron%rhor_ep) 378 end if 379 if (allocated(electronpositron%nhat_ep)) then 380 ABI_FREE(electronpositron%nhat_ep) 381 end if 382 if (allocated(electronpositron%vha_ep)) then 383 ABI_FREE(electronpositron%vha_ep) 384 end if 385 if (allocated(electronpositron%lmselect_ep)) then 386 ABI_FREE(electronpositron%lmselect_ep) 387 end if 388 if (allocated(electronpositron%gred_ep)) then 389 ABI_FREE(electronpositron%gred_ep) 390 end if 391 if (allocated(electronpositron%stress_ep)) then 392 ABI_FREE(electronpositron%stress_ep) 393 end if 394 395 if (electronpositron%has_pawrhoij_ep/=0) then 396 call pawrhoij_free(electronpositron%pawrhoij_ep) 397 end if 398 if (allocated(electronpositron%pawrhoij_ep)) then 399 ABI_FREE(electronpositron%pawrhoij_ep) 400 end if 401 402 if (electronpositron%dimcprj/=0) then 403 call pawcprj_free(electronpositron%cprj_ep) 404 end if 405 if (allocated(electronpositron%cprj_ep)) then 406 ABI_FREE(electronpositron%cprj_ep) 407 end if 408 409 electronpositron%calctype =0 410 electronpositron%particle =-1 411 electronpositron%dimcg =0 412 electronpositron%dimcprj =0 413 electronpositron%dimeigen =0 414 electronpositron%dimocc =0 415 electronpositron%has_pawrhoij_ep=0 416 electronpositron%has_pos_ham =0 417 electronpositron%istep =0 418 electronpositron%istep_scf =0 419 420 electronpositron%posdensity0_limit=.false. 421 electronpositron%scf_converged=.false. 422 423 ABI_FREE(electronpositron) 424 425 end if 426 427 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
SOURCE
730 integer function electronpositron_calctype(electronpositron) 731 732 !Arguments ------------------------------------ 733 !scalars 734 type(electronpositron_type),pointer :: electronpositron 735 736 !************************************************************************ 737 738 if (associated(electronpositron)) then 739 electronpositron_calctype=electronpositron%calctype 740 else 741 electronpositron_calctype=0 742 end if 743 744 745 end function electronpositron_calctype
m_electronpositron/electronpositron_type [ Types ]
[ Top ] [ m_electronpositron ] [ Types ]
NAME
FUNCTION
NOTES
SOURCE
63 type, public :: electronpositron_type 64 65 ! Integer scalars 66 integer :: calctype ! type of electron-positron calculation: 67 ! 0: no calculation 68 ! 1: positron in the electrons potential 69 ! 2: electrons in the positron potential 70 integer :: particle ! current particle stored in electronpositron%xxx_ep arrays 71 ! -1: no particle, 0: electron, 1: positron 72 integer :: dimcg ! Dimension of cg array dimcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 73 integer :: dimcprj ! Dimension of cprj array dimcprj=dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol*usecprj 74 integer :: dimeigen ! Dimension of eigen array dimeigen=dtset%mband*dtset%nkpt*dtset%nsppol 75 integer :: dimocc ! Dimension of occ array dimocc=dtset%mband*dtset%nkpt*dtset%nsppol 76 integer :: has_pawrhoij_ep ! flag for pawrhoij_ep (0: not allocated, 1: allocated, 2: computed) 77 integer :: has_pos_ham ! flag: 1 if current Hamiltonian in memory (vtrial, vpsp, vhartr, vxc, paw_ij%dij) 78 ! is the positronic hamiltonian, 0 is it is the electronic one 79 integer :: ixcpositron ! XC type for electron-positron correlation 80 integer :: istep ! Current index of TC-DFT SCF step 81 integer :: istep_scf ! Current index of DFT SCF step in current electron/positron minimization 82 integer :: lmmax ! Max. number of (l,m) moments over all types of atom 83 integer :: natom ! Number of atoms 84 integer :: nfft ! Number of points in FFT grid 85 integer :: nspden ! Number of spin density components 86 integer :: nstep ! Max. number of steps for the TC-DFT SCF cycle 87 88 ! Logical scalars 89 logical :: posdensity0_limit ! True if we are in the zero positron density limit 90 logical :: scf_converged ! True if the SCF cycle is converged for a positronic/electronic GS calculation 91 92 ! Real(dp) scalars 93 real(dp) :: e_hartree ! Hartree electron-positron interaction energy 94 real(dp) :: e_xc ! XC electron-positron interaction energy 95 real(dp) :: e_xcdc ! Double-counting XC electron-positron interaction energy 96 real(dp) :: e_paw ! PAW electron-positron interaction energy 97 real(dp) :: e_pawdc ! Double-counting PAW electron-positron interaction energy 98 real(dp) :: e0 ! Energy only due to particle(s) currently evolving 99 ! calctype=1, energy due to positron only 100 ! calctype=2, energy due to electrons only 101 real(dp) :: etotal_prev ! Total energy of the previous GS calculation 102 real(dp) :: lambda ! Electron-positron annihilation rate 103 real(dp) :: lifetime ! Positron lifetime 104 real(dp) :: maxfor_prev ! Max. force of the previous GS calculation 105 real(dp) :: posocc ! Occupation number for the positron 106 real(dp) :: postoldfe ! Tolerance on total energy for the TC-DFT SCF cycle 107 real(dp) :: postoldff ! Tolerance on max. force for the TC-DFT SCF cycle 108 109 ! Other scalars 110 type(energies_type) :: energies_ep ! Energies of the previous electronic/positronic SCF step 111 112 ! Logical pointers 113 logical, allocatable :: lmselect_ep(:,:) 114 ! lmselect_ep(lmmax,my_natom) 115 ! flags selecting the non-zero LM-moments of on-site densities 116 117 ! Real(dp) pointers 118 real(dp), allocatable :: cg_ep(:,:) 119 ! cg_ep(2,dimcg) 120 ! if typecalc=1: electronic wavefunctions 121 ! if typecalc=2: positronic wavefunctions 122 123 real(dp), allocatable :: eigen_ep(:) 124 ! eigen(dimeigen) 125 ! if typecalc=1: electronic eigen energies 126 ! if typecalc=2: positronic eigen energies 127 128 real(dp), allocatable :: gred_ep(:,:) 129 ! gred_ep(3,natom) 130 ! if typecalc=1: forces only due to electrons 131 ! if typecalc=2: forces only due to positron 132 133 real(dp), allocatable :: nhat_ep(:,:) 134 ! nhat_ep(nfft,nspden) 135 ! if typecalc=1: electronic compensation charge density in real space 136 ! if typecalc=2: positronic compensation charge density in real space 137 138 real(dp), allocatable :: occ_ep(:) 139 ! occ(dimocc) 140 ! if typecalc=1: electronic occupations 141 ! if typecalc=2: positronic occupations 142 143 real(dp), allocatable :: rhor_ep(:,:) 144 ! rhor_ep(nfft,nspden) 145 ! if typecalc=1: electronic density in real space 146 ! if typecalc=2: positronic density in real space 147 148 real(dp), allocatable :: stress_ep(:) 149 ! stress_ep(6) 150 ! if typecalc=1: stresses only due to electrons 151 ! if typecalc=2: stresses only due to positron 152 153 real(dp), allocatable :: vha_ep(:) 154 ! vha_ep(nfft) 155 ! if typecalc=1: electronic Hartree potential 156 ! if typecalc=2: positronic Hartree potential 157 158 ! Other pointers 159 type(pawrhoij_type), allocatable :: pawrhoij_ep(:) 160 ! pawrhoij_ep(natom) 161 ! Relevant only if PAW 162 ! if typecalc=1: electronic PAW occupation matrix associated with rhor_ep 163 ! if typecalc=2: positronic PAW occupation matrix associated with rhor_ep 164 165 type(pawcprj_type), allocatable :: cprj_ep(:,:) 166 ! cprj_ep(natom,dimcprj) 167 ! Relevant only if PAW 168 ! if typecalc=1: electronic WF projected on nl projectors <p_i|Cnk> 169 ! if typecalc=2: positronic WF projected on nl projectors <p_i|Cnk> 170 171 end type electronpositron_type 172 173 ! public procedures 174 public :: init_electronpositron 175 public :: destroy_electronpositron 176 public :: exchange_electronpositron 177 public :: electronpositron_calctype 178 public :: rhohxcpositron 179 180 CONTAINS 181 182 !===========================================================
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)
gred(3,natom)=gradients wrt nuclear positions 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
SOURCE
470 subroutine exchange_electronpositron(cg,cprj,dtset,eigen,electronpositron,energies,gred,mcg,mcprj,& 471 & mpi_enreg,my_natom,nfft,ngfft,nhat,npwarr,occ,paw_an,pawrhoij,& 472 & rhog,rhor,stress,usecprj,vhartr) 473 474 !Arguments ------------------------------------ 475 !scalars 476 integer,intent(in) :: mcg,mcprj,my_natom,nfft,usecprj 477 type(dataset_type),intent(in) :: dtset 478 type(electronpositron_type),pointer :: electronpositron 479 type(energies_type),intent(inout) :: energies 480 type(MPI_type),intent(in) :: mpi_enreg 481 !arrays 482 integer,intent(in) :: ngfft(18),npwarr(dtset%nkpt) 483 real(dp),intent(inout) :: cg(2,mcg) 484 real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol) 485 real(dp),intent(inout) :: gred(3,dtset%natom),nhat(nfft,dtset%nspden) 486 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 487 real(dp), intent(inout) :: rhog(2,nfft),rhor(nfft,dtset%nspden) 488 real(dp),intent(inout) :: stress(6),vhartr(nfft) 489 type(pawcprj_type) :: cprj(dtset%natom,mcprj*usecprj) 490 type(paw_an_type),intent(inout) :: paw_an(my_natom*dtset%usepaw) 491 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*dtset%usepaw) 492 493 !Local variables------------------------------- 494 !scalars 495 integer :: comm,iatom,ib,ibsp,icg,icgb,ifft,ii,ilm,ikpt 496 integer :: ispden,isppol,ispinor,me,my_nspinor,nband_k,npw_k,sz1,sz2,sz3 497 logical :: ltmp 498 real(dp) :: rtmp 499 type(energies_type) :: energies_tmp 500 !arrays 501 integer,allocatable :: nlmn(:),typ(:) 502 real(dp) :: ctmp(2) 503 type(pawcprj_type),allocatable :: cprj_tmp(:,:) 504 type(pawrhoij_type),allocatable :: pawrhoij_tmp(:) 505 506 !********************************************************************* 507 508 if (associated(electronpositron)) then 509 if (electronpositron%particle/=EP_NOTHING) then 510 511 ! Type of particle stored 512 if (electronpositron%particle==EP_ELECTRON) then 513 electronpositron%particle=EP_POSITRON 514 else if (electronpositron%particle==EP_POSITRON) then 515 electronpositron%particle=EP_ELECTRON 516 end if 517 518 ! Energies 519 ctmp(1)=energies%e_electronpositron 520 ! ctmp(2)=energies%edc_electronpositron 521 call energies_copy(electronpositron%energies_ep,energies_tmp) 522 call energies_copy(energies,electronpositron%energies_ep) 523 call energies_copy(energies_tmp,energies) 524 energies%e_electronpositron=ctmp(1) 525 ! energies%edc_electronpositron=ctmp(2) 526 energies%e0_electronpositron=electronpositron%e0 527 electronpositron%e0=electronpositron%energies_ep%e0_electronpositron 528 529 ! Density and PAW occupation matrix 530 do ispden=1,dtset%nspden 531 do ifft=1,nfft 532 rtmp=rhor(ifft,ispden) 533 rhor(ifft,ispden)=electronpositron%rhor_ep(ifft,ispden) 534 electronpositron%rhor_ep(ifft,ispden)=rtmp 535 end do 536 if (allocated(electronpositron%nhat_ep).and.size(nhat,2)>0) then 537 do ifft=1,nfft 538 rtmp=nhat(ifft,ispden) 539 nhat(ifft,ispden)=electronpositron%nhat_ep(ifft,ispden) 540 electronpositron%nhat_ep(ifft,ispden)=rtmp 541 end do 542 end if 543 end do 544 call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,1,ngfft,0) 545 if (dtset%usepaw==1.and.my_natom>0) then 546 if (electronpositron%has_pawrhoij_ep==1) then 547 ABI_MALLOC(pawrhoij_tmp,(my_natom)) 548 ABI_MALLOC(typ,(my_natom)) 549 ABI_MALLOC(nlmn,(my_natom)) 550 do iatom=1,my_natom 551 typ(iatom)=iatom 552 nlmn(iatom)=pawrhoij(iatom)%lmn_size 553 end do 554 ! Be careful: parallelism over atoms is ignored... 555 call pawrhoij_alloc(pawrhoij_tmp,pawrhoij(1)%cplex_rhoij,pawrhoij(1)%nspden,& 556 & pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,typ, & 557 & lmnsize=nlmn,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,& 558 & qphase=pawrhoij(1)%qphase,use_rhoij_=pawrhoij(1)%use_rhoij_,& 559 & use_rhoijres=pawrhoij(1)%use_rhoijres) 560 ABI_FREE(typ) 561 ABI_FREE(nlmn) 562 call pawrhoij_copy(pawrhoij,pawrhoij_tmp) 563 call pawrhoij_copy(electronpositron%pawrhoij_ep,pawrhoij) 564 call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep) 565 if (pawrhoij_tmp(1)%ngrhoij>0.and.pawrhoij(1)%ngrhoij==0) then 566 do iatom=1,my_natom 567 sz1=pawrhoij_tmp(iatom)%ngrhoij 568 sz2=pawrhoij_tmp(iatom)%cplex_rhoij*pawrhoij_tmp(iatom)%qphase*pawrhoij_tmp(iatom)%lmn2_size 569 sz3=pawrhoij_tmp(iatom)%nspden 570 ABI_MALLOC(pawrhoij(iatom)%grhoij,(sz1,sz2,sz3)) 571 pawrhoij(iatom)%grhoij(:,:,:)=pawrhoij_tmp(iatom)%grhoij(:,:,:) 572 end do 573 end if 574 if (pawrhoij_tmp(1)%use_rhoijres>0.and.pawrhoij(1)%use_rhoijres==0) then 575 do iatom=1,my_natom 576 sz1=pawrhoij_tmp(iatom)%cplex_rhoij*pawrhoij_tmp(iatom)%qphase*pawrhoij_tmp(iatom)%lmn2_size 577 sz2=pawrhoij_tmp(iatom)%nspden 578 ABI_MALLOC(pawrhoij(iatom)%rhoijres,(sz1,sz2)) 579 pawrhoij(iatom)%rhoijres(:,:)=pawrhoij_tmp(iatom)%rhoijres(:,:) 580 end do 581 end if 582 if (pawrhoij_tmp(1)%use_rhoij_>0.and.pawrhoij(1)%use_rhoij_==0) then 583 do iatom=1,my_natom 584 sz1=pawrhoij_tmp(iatom)%cplex_rhoij*pawrhoij_tmp(iatom)%qphase*pawrhoij_tmp(iatom)%lmn2_size 585 sz2=pawrhoij_tmp(iatom)%nspden 586 ABI_MALLOC(pawrhoij(iatom)%rhoij_,(sz1,sz2)) 587 pawrhoij(iatom)%rhoij_(:,:)=pawrhoij_tmp(iatom)%rhoij_(:,:) 588 end do 589 end if 590 if (pawrhoij_tmp(1)%lmnmix_sz>0.and.pawrhoij(1)%lmnmix_sz==0) then 591 do iatom=1,my_natom 592 ABI_MALLOC(pawrhoij(iatom)%kpawmix,(pawrhoij_tmp(iatom)%lmnmix_sz)) 593 pawrhoij(iatom)%kpawmix(:)=pawrhoij_tmp(iatom)%kpawmix(:) 594 end do 595 end if 596 call pawrhoij_free(pawrhoij_tmp) 597 ABI_FREE(pawrhoij_tmp) 598 else 599 do iatom=1,my_natom 600 pawrhoij(iatom)%rhoijp=zero 601 end do 602 end if 603 end if 604 605 ! Hartree potential 606 do ifft=1,nfft 607 rtmp=vhartr(ifft) 608 vhartr(ifft)=electronpositron%vha_ep(ifft) 609 electronpositron%vha_ep(ifft)=rtmp 610 end do 611 612 ! PAW LM-moment selection flags 613 if (dtset%usepaw==1.and.my_natom>0) then 614 if (electronpositron%lmmax>0) then 615 do iatom=1,my_natom 616 do ilm=1,paw_an(iatom)%lm_size 617 ltmp=electronpositron%lmselect_ep(ilm,iatom) 618 electronpositron%lmselect_ep(ilm,iatom)=paw_an(iatom)%lmselect(ilm) 619 paw_an(iatom)%lmselect(ilm)=ltmp 620 end do 621 end do 622 else 623 do iatom=1,my_natom 624 paw_an(iatom)%lmselect(:)=.true. 625 end do 626 end if 627 end if 628 629 ! Wave-functions 630 if (electronpositron%dimcg>0) then 631 do ii=1,electronpositron%dimcg 632 ctmp(1:2)=electronpositron%cg_ep(1:2,ii) 633 electronpositron%cg_ep(1:2,ii)=cg(1:2,ii) 634 cg(1:2,ii)=ctmp(1:2) 635 end do 636 else 637 icg=0 638 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 639 comm=mpi_enreg%comm_cell 640 me=xmpi_comm_rank(comm) 641 do isppol=1,dtset%nsppol 642 do ikpt=1,dtset%nkpt 643 npw_k=npwarr(ikpt);nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 644 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 645 icgb=icg;ibsp=0 646 do ib=1,nband_k 647 cg(:,icgb+1:icgb+my_nspinor*npw_k)=zero 648 do ispinor=1,my_nspinor 649 ibsp=ibsp+1;if (ibsp<my_nspinor*npw_k) cg(1,icgb+ibsp)=one 650 end do 651 icgb=icgb+my_nspinor*npw_k 652 end do 653 if (dtset%mkmem/=0) icg=icg+my_nspinor*npw_k*nband_k 654 end do 655 end do 656 end if 657 if (dtset%usepaw==1) then 658 if(electronpositron%dimcprj>0) then 659 ABI_MALLOC(nlmn,(dtset%natom)) 660 ABI_MALLOC(cprj_tmp,(dtset%natom,electronpositron%dimcprj)) 661 do iatom=1,dtset%natom;nlmn(iatom)=cprj(iatom,1)%nlmn;end do 662 call pawcprj_alloc(cprj_tmp,cprj(1,1)%ncpgr,nlmn) 663 ABI_FREE(nlmn) 664 call pawcprj_copy(electronpositron%cprj_ep,cprj_tmp) 665 call pawcprj_copy(cprj,electronpositron%cprj_ep) 666 call pawcprj_copy(cprj_tmp,cprj) 667 call pawcprj_free(cprj_tmp) 668 ABI_FREE(cprj_tmp) 669 else 670 !TO BE ACTIVATED WHEN cprj IS PRESENT 671 ! call pawcprj_set_zero(cprj) 672 end if 673 end if 674 675 ! Eigenvalues 676 if (electronpositron%dimeigen>0) then 677 do ii=1,electronpositron%dimeigen 678 rtmp=eigen(ii) 679 eigen(ii)=electronpositron%eigen_ep(ii) 680 electronpositron%eigen_ep(ii)=rtmp 681 end do 682 else 683 eigen(:)=9.99999_dp 684 end if 685 686 ! Occupations 687 if (electronpositron%dimocc>0) then 688 do ii=1,electronpositron%dimocc 689 rtmp=occ(ii) 690 occ(ii)=electronpositron%occ_ep(ii) 691 electronpositron%occ_ep(ii)=rtmp 692 end do 693 else 694 occ(:)=9.99999_dp 695 end if 696 697 ! Forces 698 if (allocated(electronpositron%gred_ep)) then 699 do iatom=1,dtset%natom 700 electronpositron%gred_ep(1:3,iatom)=gred(1:3,iatom)-electronpositron%gred_ep(1:3,iatom) 701 end do 702 end if 703 704 ! Stresses 705 if (allocated(electronpositron%stress_ep)) then 706 electronpositron%stress_ep(1:6)=stress(1:6)-electronpositron%stress_ep(1:6) 707 end if 708 709 end if 710 end if 711 712 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
SOURCE
205 subroutine init_electronpositron(ireadwf,dtset,electronpositron,mpi_enreg,nfft,pawrhoij,pawtab) 206 207 !Arguments ------------------------------------ 208 !scalars 209 integer,intent(in) :: ireadwf,nfft 210 type(dataset_type),intent(in) :: dtset 211 type(electronpositron_type),pointer :: electronpositron 212 type(MPI_type),intent(in) :: mpi_enreg 213 !arrays 214 type(pawrhoij_type), intent(in) :: pawrhoij(mpi_enreg%my_natom*dtset%usepaw) 215 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 216 217 !Local variables------------------------------- 218 !scalars 219 integer :: ii,my_nspinor,ncpgr,optfor,optstr 220 logical,parameter :: include_nhat_in_gamma=.false. 221 !arrays 222 integer,allocatable :: nlmn(:) 223 224 !************************************************************************ 225 226 !@electronpositron_type 227 228 if (dtset%positron/=0) then 229 230 ABI_MALLOC(electronpositron,) 231 232 electronpositron%calctype=0 233 electronpositron%particle=-1 234 235 electronpositron%ixcpositron=dtset%ixcpositron 236 electronpositron%natom=dtset%natom 237 electronpositron%nfft=nfft 238 electronpositron%nspden=dtset%nspden 239 electronpositron%istep=0 240 electronpositron%istep_scf=0 241 242 electronpositron%posocc=dtset%posocc 243 electronpositron%nstep=dtset%posnstep 244 electronpositron%postoldfe=dtset%postoldfe 245 electronpositron%postoldff=dtset%postoldff 246 electronpositron%posdensity0_limit=(dtset%ixcpositron/=2) 247 electronpositron%scf_converged=.false. 248 electronpositron%has_pos_ham=0 249 250 call energies_init(electronpositron%energies_ep) 251 252 electronpositron%e_hartree =zero 253 electronpositron%e_xc =zero 254 electronpositron%e_xcdc =zero 255 electronpositron%e_paw =zero 256 electronpositron%e_pawdc =zero 257 electronpositron%e0 =zero 258 electronpositron%etotal_prev=zero 259 electronpositron%maxfor_prev=zero 260 261 electronpositron%lambda=zero 262 electronpositron%lifetime=zero 263 264 ABI_MALLOC(electronpositron%rhor_ep,(nfft,dtset%nspden)) 265 ABI_MALLOC(electronpositron%vha_ep,(nfft)) 266 ABI_MALLOC(electronpositron%pawrhoij_ep,(mpi_enreg%my_natom*dtset%usepaw)) 267 268 if (dtset%usepaw==1) then 269 electronpositron%has_pawrhoij_ep=1 270 if (mpi_enreg%my_natom>0) then 271 call pawrhoij_alloc(electronpositron%pawrhoij_ep,pawrhoij(1)%cplex_rhoij,pawrhoij(1)%nspden,& 272 & pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,dtset%typat,& 273 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,& 274 & pawtab=pawtab,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,& 275 & qphase=pawrhoij(1)%qphase,use_rhoij_=pawrhoij(1)%use_rhoij_,& 276 & use_rhoijres=pawrhoij(1)%use_rhoijres) 277 end if 278 electronpositron%lmmax=0 279 do ii=1,dtset%ntypat 280 electronpositron%lmmax=max(electronpositron%lmmax,pawtab(ii)%lcut_size**2) 281 end do 282 ABI_MALLOC(electronpositron%lmselect_ep,(electronpositron%lmmax,mpi_enreg%my_natom)) 283 if (maxval(pawtab(1:dtset%ntypat)%usexcnhat)==0.or.(.not.include_nhat_in_gamma)) then 284 ABI_MALLOC(electronpositron%nhat_ep,(nfft,dtset%nspden)) 285 end if 286 else 287 electronpositron%has_pawrhoij_ep=0 288 electronpositron%lmmax=0 289 end if 290 291 if (dtset%positron<=-10.or.dtset%posdoppler>0) then 292 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 293 electronpositron%dimcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 294 electronpositron%dimocc=dtset%mband*dtset%nkpt*dtset%nsppol 295 electronpositron%dimeigen=dtset%mband*dtset%nkpt*dtset%nsppol 296 ABI_MALLOC(electronpositron%cg_ep,(2,electronpositron%dimcg)) 297 ABI_MALLOC(electronpositron%eigen_ep,(electronpositron%dimeigen)) 298 ABI_MALLOC(electronpositron%occ_ep,(electronpositron%dimocc)) 299 electronpositron%dimcprj=0 300 ! if (.false.) then !TEMPORARY: will be activated later 301 if (dtset%usepaw==1.and.dtset%pawusecp>0.and.dtset%posdoppler>0) then 302 electronpositron%dimcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 303 if (mpi_enreg%paral_kgb/=0) electronpositron%dimcprj=electronpositron%dimcprj/mpi_enreg%nproc_band 304 ABI_MALLOC(electronpositron%cprj_ep,(dtset%natom,electronpositron%dimcprj)) 305 ABI_MALLOC(nlmn,(dtset%natom)) 306 ncpgr=0 307 do ii=1,dtset%natom;nlmn(ii)=pawtab(dtset%typat(ii))%lmn_size;end do 308 call pawcprj_alloc(electronpositron%cprj_ep,ncpgr,nlmn) 309 ABI_FREE(nlmn) 310 else 311 ABI_MALLOC(electronpositron%cprj_ep,(dtset%natom,electronpositron%dimcprj)) 312 end if 313 else 314 electronpositron%dimcg =0 315 electronpositron%dimcprj =0 316 electronpositron%dimeigen=0 317 electronpositron%dimocc =0 318 end if 319 320 optfor=0;optstr=0 321 if ((dtset%optforces>0.or.dtset%ionmov/=0.or.abs(dtset%toldff)>tiny(0._dp))) optfor=1 322 if (dtset%optstress>0.and.dtset%iscf>0.and.(dtset%nstep>0.or.ireadwf==1)) optstr=1 323 324 if (optfor>0) then 325 ABI_MALLOC(electronpositron%gred_ep,(3,dtset%natom)) 326 electronpositron%gred_ep(:,:)=zero 327 end if 328 329 if (optstr>0) then 330 ABI_MALLOC(electronpositron%stress_ep,(6)) 331 electronpositron%stress_ep(:)=zero 332 end if 333 334 else !dtset%positron==0 335 nullify(electronpositron) 336 end if 337 338 end subroutine init_electronpositron