TABLE OF CONTENTS
- ABINIT/dfpt_etot
- ABINIT/dfpt_newvtr
- ABINIT/dfpt_nselt
- ABINIT/dfpt_nstdy
- ABINIT/dfpt_nsteltwf
- ABINIT/dfpt_rhofermi
- ABINIT/dfpt_scfcv
- ABINIT/dfpt_wfkfermi
- ABINIT/m_dfpt_scfcv
- ABINIT/newfermie1
ABINIT/dfpt_etot [ Functions ]
NAME
dfpt_etot
FUNCTION
Assemble different contributions to the variational part of the 2nd derivative of total energy
INPUTS
berryopt= 4/14: electric field is on; berryopt = 6/7/16/17: electric displacement field is on; eberry=energy associated with Berry phase edocc=correction to 2nd-order total energy coming from changes of occupation ehart1=1st-order Hartree part of 2nd-order total energy eeig0=0th-order eigenenergies part of 2nd-order total energy eew=2nd derivative of Ewald energy (hartree) efrhar=contrib. from frozen-wavefunction, hartree energy, to the 2nd-derivative of total energy efrkin=contrib. from frozen-wavefunction, kinetic energy, to the 2nd-derivative of total energy efrloc=contrib. from frozen-wavefunction, local potential, to the 2nd-derivative of total energy efrnl=contribution from frozen-wavefunction, non-local potential, to the 2nd-derivative of total energy efrx1=contrib. from frozen-wavefunction, xc core correction(1), to the 2nd-derivative of total energy efrx2=contribution from frozen-wavefunction, xc core correction(2), to the second-derivative of total energy. ek0=0th-order kinetic energy part of 2nd-order total energy. ek1=1st-order kinetic energy part of 2nd-order total energy. eii=2nd derivative of pseudopotential core energy (hartree) eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy elpsp1=1st-order local pseudopot. part of 2nd-order total energy. end0=0th-order nuclear dipole energy part of 2nd-order total energy. end1=1st-order nuclear dipole energy part of 2nd-order total energy. enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy. enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy. epaw1=1st-order PAW on-sitew part of 2nd-order total energy. evdw=DFT-D semi-empirical part of 2nd-order total energy evxctau0=0th-order vxctau energy evxctau1=1st-order vxctau energy exc1=1st-order exchange-correlation part of 2nd-order total energy ipert=type of the perturbation natom=number of atoms optene=option for the computation of 2nd-order total energy (-1=no computation; 0=direct scheme; 1=double-counting scheme)
OUTPUT
deltae=change in energy between the previous and present SCF cycle and previous SCF cycle. etotal=2nd-order total energy evar=variational part of the 2nd-order total energy
SIDE EFFECTS
input/output elast=previous value of the 2nd-order total energy, needed to compute deltae, then updated (cannot simply be saved, because set to zero at each new call of dfpt_scfcv).
SOURCE
1604 subroutine dfpt_etot(berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,& 1605 & efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,& 1606 & end0,enl0,enl1,epaw1,etotal,evar,evdw,evxctau0,& 1607 & exc1,ipert,natom,optene) 1608 1609 !Arguments ------------------------------------ 1610 !scalars 1611 integer,intent(in) :: berryopt,ipert,natom,optene 1612 real(dp),intent(in) :: eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1 1613 real(dp),intent(in) :: efrx2,ehart1,eii,ek0,ek1,eloc0,elpsp1,end0,enl0,enl1,epaw1 1614 real(dp),intent(in) :: evdw,evxctau0,exc1 1615 real(dp),intent(inout) :: elast 1616 real(dp),intent(out) :: deltae,etotal,evar 1617 1618 !Local variables------------------------------- 1619 !scalars 1620 ! character(len=500) :: message 1621 1622 ! ********************************************************************* 1623 1624 if (optene==1) then 1625 ABI_BUG('Double-counting scheme not yet allowed!') 1626 end if 1627 1628 if (optene>-1) then 1629 1630 ! Compute 2nd-order variational energy by direct scheme 1631 if (optene==0) then 1632 1633 ! Atomic displ. perturbation 1634 if ( ipert>=1 .and. ipert<=natom ) then 1635 evar=ek0+edocc+eeig0+eloc0+enl0+ehart1+exc1+enl1+epaw1+elpsp1 1636 1637 else if (ipert==natom+1) then 1638 ! evar=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1+end0+end1+evxctau0+evxctau1 1639 ! JWZ debug end1 and evxctau1 already appear in enl1 (?!) 1640 evar=ek0+edocc+eeig0+eloc0+ek1+ehart1+exc1+enl0+enl1+end0+evxctau0 1641 1642 else if (ipert==natom+10 .or. ipert==natom+11) then 1643 evar=ek0+edocc+eeig0+eloc0+enl0 +ek1 ! here ek1 contains a lot of contributions 1644 1645 ! For ipert==natom+2, some contributions vanish, noticeably ek1 1646 else if (ipert==natom+2) then 1647 evar=ek0+edocc+eeig0+eloc0+enl0+ehart1+exc1+enl1+ek1+epaw1 1648 1649 ! All terms enter for strain perturbation 1650 else if ( ipert==natom+3 .or. ipert==natom+4 ) then 1651 evar=ek0+edocc+eeig0+eloc0+enl0+ehart1+exc1+enl1+ek1+epaw1+elpsp1 1652 1653 ! terms for Zeeman perturbation, SPr 2deb 1654 else if ( ipert==natom+5 ) then 1655 evar=ek0+edocc+eeig0+eloc0+enl0+ehart1+exc1+epaw1 1656 end if 1657 end if 1658 1659 ! Compute energy residual 1660 deltae=evar-elast 1661 elast=evar 1662 1663 ! Compute 2nd-order total energy by direct scheme 1664 if (optene==0) then 1665 if (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17) then 1666 if (ipert<=natom) then 1667 etotal=evar+eew+evdw+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2+two*eberry 1668 else if (ipert==natom+2) then 1669 etotal=half*evar+eew+evdw+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2+two*eberry 1670 end if 1671 else 1672 if (ipert/=natom+10 .and. ipert/=natom+11) then 1673 etotal=evar+eew+evdw+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2 1674 else 1675 etotal=evar ! For 2nd order sternheimer equations, the total (4th order) energy is not used (yet) 1676 end if 1677 end if 1678 end if 1679 1680 end if 1681 1682 end subroutine dfpt_etot
ABINIT/dfpt_newvtr [ Functions ]
NAME
dfpt_newvtr
FUNCTION
Compute new first-order trial potential by mixing new and old values. First, compute preconditioned residual first-order potential. Then, call one of the self-consistency drivers, and update vtrial.
INPUTS
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX dielar(7)=input parameters for dielectric matrix: diecut,dielng,diemac,diemix,diegap,dielam,diemixmag. dtset <type(dataset_type)>=all input variables in this dataset | isecur=level of security of the computation | mffmem=governs the number of FFT arrays which are fit in core memory | it is either 1, in which case the array f_fftgr is used, | or 0, in which case the array f_fftgr_disk is used | natom=number of atoms | nspden=number of spin-density components | paral_kgb=option for (kpt,g vectors,bands) parallelism | pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part etotal=the total energy obtained from the input vtrial ffttomix(nfft*(1-nfftmix/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse) initialized= if 0 the initialization of the RF run is not yet finished iscf=( <= 0 =>non-SCF), >0 => SCF) iscf =1 => determination of the largest eigenvalue of the SCF cycle iscf =2 => SCF cycle, simple mixing iscf =3 => SCF cycle, Anderson mixing iscf =4 => SCF cycle, Anderson mixing (order 2) iscf =5 => SCF cycle, CG based on the minimization of the energy iscf =7 => SCF cycle, Pulay mixing ispmix=1 if mixing is done in real space, 2 if mixing is done in reciprocal space istep= number of the step in the SCF cycle mixtofft(nfftmix*(1-nfftmix/nfft))=Index of the points of the FFT grid used for mixing (coarse) on the FFT (fine) grid mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor nfft=(effective) number of FFT grid points (for this processor) nfftmix=dimension of FFT grid used to mix the densities (used in PAW only) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft ngfftmix(18)=contain all needed information about 3D FFT, for the grid corresponding to nfftmix npawmix=-PAW only- number of spherical part elements to be mixed qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2). pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data Use here rhoij residuals (and gradients) rhor(cplex*nfft,nspden)=array for 1st-order electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) usepaw= 0 for non paw calculation; =1 for paw calculation vresid(cplex*nfft,nspden)=array for the residual of the potential xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.
SIDE EFFECTS
vtrial(cplex*nfft,nspden)= at input, it is the "in" trial potential that gave vresid=(v_out-v_in)
at output, it is an updated "mixed" trial potential
==== if usepaw==1
pawrhoij(natom)%nrhoijsel,rhoijselect,rhoijp= several arrays
containing new values of rhoij (augmentation occupancies)
NOTES
In case of PAW calculations:
Computations are done either on the fine FFT grid or the coarse grid (depending on dtset%pawmixdg)
All variables (nfft,ngfft,mgfft) refer to the fine FFT grid.
All arrays (densities/potentials...) are computed on this fine FFT grid.
Developpers have to be careful when introducing others arrays:
they have to be stored on the fine FFT grid.
In case of norm-conserving calculations the FFT grid is the usual FFT grid.
Subtility in PAW and non-collinear magnetism:
Potentials are stored in (up-up,dn-dn,Re[up-dn],Im[up-dn]) format
On-site occupancies (rhoij) are stored in (n,mx,my,mz)
This is compatible provided that the mixing factors for n and m are identical
and that the residual is not a combination of V_res and rhoij_res (pawoptmix=0).
SOURCE
1933 subroutine dfpt_newvtr(cplex,dbl_nnsclo,dielar,dtset,etotal,ffttomix,& 1934 & initialized,iscf,ispmix,istep,mix,mixtofft,& 1935 & mpi_enreg,my_natom,nfft,nfftmix,ngfft,ngfftmix,npawmix,pawrhoij,& 1936 & qphon,rhor,rprimd,usepaw,vresid,vtrial) 1937 1938 !Arguments------------------------------- 1939 !scalars 1940 integer,intent(in) :: cplex,initialized,iscf,ispmix,istep,my_natom,nfft 1941 integer,intent(in) :: nfftmix,npawmix,usepaw 1942 integer,intent(inout) :: dbl_nnsclo !vz_i 1943 real(dp),intent(in) :: etotal 1944 type(MPI_type),intent(in) :: mpi_enreg 1945 type(ab7_mixing_object), intent(inout) :: mix 1946 type(dataset_type),intent(in) :: dtset 1947 !arrays 1948 integer,intent(in) :: ffttomix(nfft*(1-nfftmix/nfft)) 1949 integer,intent(in) :: mixtofft(nfftmix*(1-nfftmix/nfft)),ngfft(18) 1950 integer,intent(in) :: ngfftmix(18) 1951 real(dp),intent(in) :: dielar(7),qphon(3) 1952 real(dp), intent(in), target :: rhor(cplex*nfft,dtset%nspden) 1953 real(dp),intent(in) :: rprimd(3,3) 1954 real(dp),intent(inout) :: vresid(cplex*nfft,dtset%nspden) 1955 real(dp),intent(inout) :: vtrial(cplex*nfft,dtset%nspden) 1956 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw) 1957 1958 !Local variables------------------------------- 1959 !scalars 1960 integer :: cplex_mix,cplex_rhoij,dplex,i_vresid1,i_vrespc1,iatom,ifft,indx,iq,iq0 1961 integer :: irhoij,ispden,jfft,jrhoij,klmn,kklmn,kmix,moved_atm_inside,nfftot,qphase 1962 integer :: mpicomm,errid 1963 logical :: mpi_summarize,reset 1964 real(dp) :: fact,mixfac,mixfac_eff,mixfacmag,ucvol 1965 character(len=500) :: msg 1966 !arrays 1967 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2) 1968 real(dp),allocatable :: rhoijrespc(:),rhoijtmp(:,:) 1969 real(dp),allocatable :: vresid0(:,:),vrespc(:,:),vreswk(:,:) 1970 real(dp), pointer :: vtrial0(:,:),vpaw(:) 1971 real(dp),allocatable :: vtrialg(:,:,:) 1972 1973 ! ************************************************************************* 1974 1975 DBG_ENTER("COLL") 1976 1977 call timab(158,1,tsec) 1978 1979 !Compatibility tests 1980 if(usepaw==1) then 1981 if(dtset%nspden==4.and.dtset%pawoptmix==1) then 1982 ABI_ERROR('pawoptmix=1 is not compatible with nspden=4 !') 1983 end if 1984 if (my_natom>0) then 1985 if (pawrhoij(1)%qphase<cplex) then 1986 ABI_ERROR('pawrhoij()%qphase must be >=cplex !') 1987 end if 1988 end if 1989 end if 1990 1991 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 1992 cplex_mix=max(cplex,ispmix) 1993 if (usepaw==1.and.my_natom>0) then 1994 cplex_rhoij=pawrhoij(1)%cplex_rhoij 1995 qphase=pawrhoij(1)%qphase 1996 end if 1997 1998 !Compute different geometric tensor, as well as ucvol, from rprimd 1999 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2000 moved_atm_inside=0 2001 2002 !Select components of potential to be mixed 2003 ABI_MALLOC(vtrial0,(cplex_mix*nfftmix,dtset%nspden)) 2004 ABI_MALLOC(vresid0,(cplex_mix*nfftmix,dtset%nspden)) 2005 if (ispmix==1.and.nfft==nfftmix) then 2006 vtrial0=vtrial;vresid0=vresid 2007 else if (nfft==nfftmix) then 2008 do ispden=1,dtset%nspden 2009 call fourdp(cplex,vtrial0(:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,1, ngfft, 0) 2010 call fourdp(cplex,vresid0(:,ispden),vresid(:,ispden),-1,mpi_enreg,nfft,1, ngfft, 0) 2011 end do 2012 else 2013 ABI_MALLOC(vtrialg,(2,nfft,dtset%nspden)) 2014 ABI_MALLOC(vreswk,(2,nfft)) 2015 do ispden=1,dtset%nspden 2016 fact=dielar(4);if (ispden>1) fact=dielar(7) 2017 call fourdp(cplex,vtrialg(:,:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,1, ngfft, 0) 2018 call fourdp(cplex,vreswk,vresid(:,ispden),-1,mpi_enreg,nfft,1, ngfft, 0) 2019 do ifft=1,nfft 2020 if (ffttomix(ifft)>0) then 2021 jfft=2*ffttomix(ifft) 2022 vtrial0(jfft-1,ispden)=vtrialg(1,ifft,ispden) 2023 vtrial0(jfft ,ispden)=vtrialg(2,ifft,ispden) 2024 vresid0(jfft-1,ispden)=vreswk(1,ifft) 2025 vresid0(jfft ,ispden)=vreswk(2,ifft) 2026 else 2027 vtrialg(:,ifft,ispden)=vtrialg(:,ifft,ispden)+fact*vreswk(:,ifft) 2028 end if 2029 end do 2030 end do 2031 ABI_FREE(vreswk) 2032 end if 2033 2034 !Precondition the potential residual: 2035 !Use a model dielectric function preconditioning, or simple mixing 2036 ABI_MALLOC(vrespc,(cplex_mix*nfftmix,dtset%nspden)) 2037 call moddiel(cplex_mix,dielar,mpi_enreg,nfftmix,ngfftmix,dtset%nspden,ispmix,0,qphon,rprimd,vresid0,vrespc) 2038 2039 !PAW only : precondition the rhoij quantities (augmentation occupancies) residuals. 2040 !Use a simple preconditionning with the same mixing factor 2041 !as the model dielectric function. 2042 if (usepaw==1.and.my_natom>0) then 2043 ABI_MALLOC(rhoijrespc,(npawmix)) 2044 mixfac=dielar(4);mixfacmag=abs(dielar(7)) 2045 if (cplex_rhoij==1) then 2046 indx=0 2047 do iatom=1,my_natom 2048 do iq=1,qphase 2049 iq0=merge(0,cplex_rhoij*pawrhoij(iatom)%lmn2_size,iq==1) 2050 do ispden=1,pawrhoij(iatom)%nspden 2051 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 2052 do kmix=1,pawrhoij(iatom)%lmnmix_sz 2053 indx=indx+1;klmn=iq0+pawrhoij(iatom)%kpawmix(kmix) 2054 rhoijrespc(indx)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 2055 end do 2056 end do 2057 end do 2058 end do 2059 else 2060 indx=-1 2061 do iatom=1,my_natom 2062 do iq=1,qphase 2063 iq0=merge(0,cplex_rhoij*pawrhoij(iatom)%lmn2_size,iq==1) 2064 do ispden=1,pawrhoij(iatom)%nspden 2065 mixfac_eff=mixfac;if (ispden>1) mixfac_eff=mixfacmag 2066 do kmix=1,pawrhoij(iatom)%lmnmix_sz 2067 indx=indx+2;klmn=iq0+2*pawrhoij(iatom)%kpawmix(kmix)-1 2068 rhoijrespc(indx:indx+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 2069 end do 2070 end do 2071 end do 2072 end do 2073 end if 2074 end if 2075 2076 !------Compute new vtrial 2077 2078 i_vresid1=mix%i_vresid(1) 2079 i_vrespc1=mix%i_vrespc(1) 2080 2081 !Initialise working arrays for the mixing object. 2082 call ab7_mixing_eval_allocate(mix, istep) 2083 2084 !Copy current step arrays. 2085 call ab7_mixing_copy_current_step(mix, vresid0, errid, msg, arr_respc = vrespc) 2086 2087 if (errid /= AB7_NO_ERROR) then 2088 ABI_ERROR(msg) 2089 end if 2090 2091 ABI_FREE(vrespc) 2092 ABI_FREE(vresid0) 2093 2094 !PAW: either use the array f_paw or the array f_paw_disk 2095 ABI_MALLOC(vpaw,(npawmix*usepaw)) 2096 if (usepaw==1.and.my_natom>0) then 2097 dplex=cplex_rhoij-1 ; indx=-dplex 2098 do iatom=1,my_natom 2099 ABI_MALLOC(rhoijtmp,(cplex_rhoij*pawrhoij(iatom)%lmn2_size,1)) 2100 do iq=1,qphase 2101 iq0=merge(0,cplex_rhoij*pawrhoij(iatom)%lmn2_size,iq==1) 2102 do ispden=1,pawrhoij(iatom)%nspden 2103 rhoijtmp=zero ; jrhoij=iq0+1 2104 do irhoij=1,pawrhoij(iatom)%nrhoijsel 2105 klmn=cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 2106 rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 2107 jrhoij=jrhoij+cplex_rhoij 2108 end do 2109 do kmix=1,pawrhoij(iatom)%lmnmix_sz 2110 indx=indx+cplex_rhoij;klmn=cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex ; kklmn=iq0+klmn 2111 vpaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden) 2112 mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden) 2113 mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex) 2114 end do 2115 end do 2116 end do 2117 ABI_FREE(rhoijtmp) 2118 end do 2119 end if 2120 2121 !Unlike for GS, no need to modify the mean of vtrial 2122 2123 mpicomm=0;mpi_summarize=.false. 2124 reset=.false.;if (initialized==0) reset=.true. 2125 call ab7_mixing_eval(mix, vtrial0, istep, nfftot, ucvol, & 2126 & mpicomm, mpi_summarize, errid, msg, & 2127 & reset = reset, isecur = dtset%isecur, & 2128 & pawopt = dtset%pawoptmix, response = 1, pawarr = vpaw, & 2129 & etotal = etotal, potden = rhor, comm_atom=mpi_enreg%comm_atom) 2130 2131 if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then 2132 dbl_nnsclo = 1 2133 else if (errid /= AB7_NO_ERROR) then 2134 ! MG FIXME, Why this? 2135 ! One should propagate the error so that we can handle it 2136 ! in the caller! 2137 ABI_ERROR(msg) 2138 end if 2139 2140 !Do here the mixing of the potential 2141 if(iscf==2 .or. iscf==3 .or. iscf==7)then 2142 ! PAW: restore rhoij from compact storage 2143 if (usepaw==1.and.my_natom>0) then 2144 dplex=cplex_rhoij-1 ; indx=-dplex 2145 do iatom=1,my_natom 2146 ABI_MALLOC(rhoijtmp,(cplex_rhoij*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 2147 rhoijtmp=zero 2148 do iq=1,qphase 2149 iq0=merge(0,cplex_rhoij*pawrhoij(iatom)%lmn2_size,iq==1) 2150 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 2151 do ispden=1,pawrhoij(iatom)%nspden 2152 jrhoij=iq0+1 2153 do irhoij=1,pawrhoij(iatom)%nrhoijsel 2154 klmn=iq0+cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 2155 rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 2156 jrhoij=jrhoij+cplex_rhoij 2157 end do 2158 end do 2159 end if 2160 do ispden=1,pawrhoij(iatom)%nspden 2161 do kmix=1,pawrhoij(iatom)%lmnmix_sz 2162 indx=indx+cplex_rhoij;klmn=iq0+cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex 2163 rhoijtmp(klmn:klmn+dplex,ispden)=vpaw(indx:indx+dplex) 2164 end do 2165 end do 2166 end do 2167 call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,& 2168 & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,& 2169 & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp) 2170 ABI_FREE(rhoijtmp) 2171 end do 2172 end if 2173 2174 else if(iscf==5 .or. iscf==6)then 2175 if(ispmix/=1) then 2176 ABI_ERROR('Mixing on reciprocal space not allowed with iscf=5 or 6.') 2177 end if 2178 ! PAW: apply a simple mixing to rhoij (this is temporary) 2179 if (usepaw==1.and.my_natom>0) then 2180 indx=1-cplex_rhoij 2181 do iatom=1,my_natom 2182 ABI_MALLOC(rhoijtmp,(cplex_rhoij*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden)) 2183 rhoijtmp=zero 2184 do iq=1,qphase 2185 iq0=merge(0,cplex_rhoij*pawrhoij(iatom)%lmn2_size,iq==1) 2186 if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then 2187 do ispden=1,pawrhoij(iatom)%nspden 2188 do kmix=1,pawrhoij(iatom)%lmnmix_sz 2189 indx=indx+cplex_rhoij;klmn=iq0+cplex_rhoij*pawrhoij(iatom)%kpawmix(kmix)-dplex 2190 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) & 2191 & -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden) 2192 end do 2193 end do 2194 end if 2195 do ispden=1,pawrhoij(iatom)%nspden 2196 jrhoij=iq0+1 2197 do irhoij=1,pawrhoij(iatom)%nrhoijsel 2198 klmn=iq0+cplex_rhoij*pawrhoij(iatom)%rhoijselect(irhoij)-dplex 2199 rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) & 2200 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 2201 jrhoij=jrhoij+cplex_rhoij 2202 end do 2203 end do 2204 end do 2205 call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,& 2206 & pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,& 2207 & pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp) 2208 ABI_FREE(rhoijtmp) 2209 end do 2210 end if 2211 end if 2212 2213 ABI_FREE(vpaw) 2214 if (usepaw==1.and.my_natom>0) then 2215 ABI_FREE(rhoijrespc) 2216 end if 2217 2218 !Eventually write the data on disk and deallocate f_fftgr_disk 2219 call ab7_mixing_eval_deallocate(mix) 2220 2221 !Restore potential 2222 if (ispmix==1.and.nfft==nfftmix) then 2223 vtrial=vtrial0 2224 else if (nfft==nfftmix) then 2225 do ispden=1,dtset%nspden 2226 call fourdp(cplex,vtrial0(:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,1, ngfft,0) 2227 end do 2228 else 2229 do ispden=1,dtset%nspden 2230 do ifft=1,nfftmix 2231 jfft=mixtofft(ifft) 2232 vtrialg(1,jfft,ispden)=vtrial0(2*ifft-1,ispden) 2233 vtrialg(2,jfft,ispden)=vtrial0(2*ifft ,ispden) 2234 end do 2235 call fourdp(cplex,vtrialg(:,:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,1,ngfft,0) 2236 end do 2237 ABI_FREE(vtrialg) 2238 end if 2239 ABI_FREE(vtrial0) 2240 2241 call timab(158,2,tsec) 2242 2243 DBG_ENTER("COLL") 2244 2245 end subroutine dfpt_newvtr
ABINIT/dfpt_nselt [ Functions ]
NAME
dfpt_nselt
FUNCTION
This routine compute the non-stationary expression for the second derivative of the total energy, wrt strain for a whole row of mixed strain derivatives.
INPUTS
cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q. cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX ecut=cut-off energy for plane wave basis sphere (Ha) ecutsm=smearing energy for plane wave kinetic energy (Ha) effmass_free=effective mass for electrons (1. in common case) gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations gsqcut=cutoff on (k+G)^2 (bohr^-2) idir=direction of the perturbation ipert=type of the perturbation istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates. kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ kxc(nfft,nkxc)=exchange and correlation kernel mband=maximum number of bands mband_mem=maximum number of bands on this cpu mgfft=maximum size of 1D FFTs mkmem =number of k points treated by this node. mk1mem =number of k points treated by this node (RF data). mpert =maximum number of ipert mpi_enreg=information about MPI parallelization mpsang= 1+maximum angular momentum for nonlocal pseudopotentials mpw=maximum dimensioned size of npw or wfs at k mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs). maximum dimension for q points in grids for nonlocal form factors natom=number of atoms in cell. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin 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 nkpt_rbz=number of k points in the reduced BZ for this perturbation nkxc=second dimension of the kxc array. If /=0, the exchange-correlation kernel must be computed. nloalg(3)=governs the choice of the algorithm for non-local operator. npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point nspden=number of spin-density components nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized nsym1=number of symmetry elements in space group consistent with perturbation ntypat=number of types of atoms in unit cell. occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced Brillouin zone (usually =2) ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information prtbbb=if 1, band-by-band decomposition (also dim of d2bbb) psps <type(pseudopotential_type)>=variables related to pseudopotentials qphon(3)=reduced coordinates for the phonon wavelength rhog(2,nfft)=array for Fourier transform of GS electron density rhor(nfft,nspden)=GS electron density in electrons/bohr**3. rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3. rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) symrc1(3,3,nsym1)=symmetry operations in reciprocal space typat(natom)=type integer for each atom in cell ucvol=unit cell volume in bohr**3. wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang)= real spherical harmonics for each G and k+q point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k point ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k+g point [rfstrs_ref]= if eq 1 the reference energy in vlocalstr is shited to the same valuea as in the FxE routines
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
second order derivatives
d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs
SOURCE
2335 subroutine dfpt_nselt(blkflg,cg,cg1,cplex,& 2336 & d2bbb,d2lo,d2nl,ecut,ecutsm,effmass_free,& 2337 & gmet,gprimd,gsqcut,idir,& 2338 & ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mband,mband_mem,mgfft,& 2339 & mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,mpw1,& 2340 & natom,nband_rbz,nfft,ngfft,& 2341 & nkpt_rbz,nkxc,nloalg,npwarr,npwar1,nspden,nspinor,nsppol,& 2342 & nsym1,ntypat,occ_rbz,& 2343 & ph1d,prtbbb,psps,qphon,rhog,& 2344 & rhor,rhor1,rmet,rprimd,symrc1,typat,ucvol,& 2345 & wtk_rbz,& 2346 & xred,ylm,ylm1,ylmgr,ylmgr1,& 2347 & rfstrs_ref) 2348 2349 !Arguments ------------------------------- 2350 !scalars 2351 integer,intent(in) :: cplex,idir,ipert,mband,mgfft,mk1mem 2352 integer,intent(in) :: mband_mem 2353 integer,intent(in) :: mkmem,mpert,mpsang,mpw,mpw1,natom,nfft,nkpt_rbz 2354 integer,intent(in) :: nkxc,nspden,nspinor,nsppol,nsym1,ntypat 2355 integer,intent(in) :: prtbbb 2356 integer,intent(in),optional :: rfstrs_ref 2357 real(dp),intent(in) :: ecut,ecutsm,effmass_free,gsqcut,ucvol 2358 type(MPI_type),intent(in) :: mpi_enreg 2359 type(pseudopotential_type),intent(in) :: psps 2360 !arrays 2361 integer,intent(in) :: istwfk_rbz(nkpt_rbz) 2362 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 2363 integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18) 2364 integer,intent(in) :: nloalg(3),npwar1(nkpt_rbz),npwarr(nkpt_rbz) 2365 integer,intent(in) :: symrc1(3,3,nsym1),typat(natom) 2366 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 2367 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol) 2368 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol) 2369 real(dp),intent(in) :: gmet(3,3) 2370 real(dp),intent(in) :: gprimd(3,3),kpt_rbz(3,nkpt_rbz),kxc(nfft,nkxc) 2371 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol) 2372 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qphon(3),rhog(2,nfft) 2373 real(dp),intent(in) :: rhor(nfft,nspden) 2374 real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3) 2375 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom) 2376 real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm) 2377 real(dp),intent(in) :: ylm1(mpw1*mk1mem,mpsang*mpsang*psps%useylm) 2378 real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm) 2379 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*psps%useylm) 2380 real(dp),intent(out) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb) 2381 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) 2382 real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert) 2383 2384 !Local variables------------------------------- 2385 !scalars 2386 integer :: ban2tot,bantot,bd2tot_index,bdtot_index,g0term 2387 integer :: icg,icg1,idir1,ifft,ii,ikg,ikg1,ikpt,comm 2388 integer :: ilm,ipert1,ispden,isppol,istr1,istwf_k 2389 integer :: mbd2kpsp,mbdkpsp,me,n1,n2,n3,n3xccc,n4,n5,n6 2390 integer :: nband_k,nfftot,npw1_k,npw_k,option,rfstrs_ref_ 2391 logical :: nmxc=.false. 2392 real(dp) :: doti,dotr 2393 real(dp) :: wtk_k 2394 character(len=500) :: message 2395 type(gs_hamiltonian_type) :: gs_hamk 2396 !arrays 2397 integer :: ikpt_fbz(3) 2398 integer,allocatable :: kg1_k(:,:),kg_k(:,:) 2399 real(dp) :: kpoint(3),restr(6),dummy(0,0) 2400 real(dp),allocatable :: d2bbb_k(:,:,:,:),d2nl_k(:,:,:) 2401 real(dp),allocatable :: occ_k(:) 2402 real(dp),allocatable :: vhartr01(:),vpsp1(:),vxc1(:,:),xccc3d1(:),ylm1_k(:,:) 2403 real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:),ylmgr_k(:,:,:) 2404 type(pawtab_type) :: pawtab_dum(0) 2405 2406 2407 ! ********************************************************************* 2408 2409 !Init me 2410 comm = mpi_enreg%comm_cell 2411 me = mpi_enreg%me_kpt 2412 2413 !Unit numbers 2414 2415 !Zero only portion of nonlocal matrix to be computed here 2416 d2nl(:,:,natom+3:natom+4,idir,ipert)=zero 2417 bdtot_index=0 2418 bd2tot_index=0 2419 icg=0 2420 icg1=0 2421 mbdkpsp=mband*nkpt_rbz*nsppol 2422 mbd2kpsp=2*mband**2*nkpt_rbz*nsppol 2423 2424 !Update list of computed matrix elements 2425 if((ipert==natom+3) .or. (ipert==natom+4)) then 2426 ! Eventually expand when strain coupling to other perturbations is implemented 2427 do ipert1=natom+3,natom+4 2428 do idir1=1,3 2429 blkflg(idir1,ipert1,idir,ipert)=1 2430 end do 2431 end do 2432 end if 2433 2434 ABI_MALLOC(d2bbb_k,(2,3,mband,mband*prtbbb)) 2435 ABI_MALLOC(d2nl_k,(2,3,mpert)) 2436 2437 ABI_MALLOC(kg_k,(3,mpw)) 2438 ABI_MALLOC(kg1_k,(3,mpw1)) 2439 2440 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 2441 n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6) 2442 nfftot=n1*n2*n3 2443 2444 !Initialize Hamiltonian (k-independent terms) - NCPP only 2445 call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,& 2446 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,ph1d=ph1d) 2447 2448 bantot = 0 2449 ban2tot = 0 2450 2451 !LOOP OVER SPINS 2452 do isppol=1,nsppol 2453 2454 if (nsppol/=1) then 2455 write(message,*)' **** In dfpt_nselt for isppol=',isppol 2456 call wrtout(std_out,message,'COLL') 2457 end if 2458 2459 ikg=0 2460 ikg1=0 2461 2462 ikpt_fbz(1:3)=0 2463 2464 ! BIG FAT k POINT LOOP 2465 do ikpt=1,nkpt_rbz 2466 2467 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 2468 istwf_k=istwfk_rbz(ikpt) 2469 npw_k=npwarr(ikpt) 2470 npw1_k=npwar1(ikpt) 2471 kpoint(:)=kpt_rbz(:,ikpt) 2472 2473 bantot = bantot + nband_k 2474 ban2tot = ban2tot + 2*nband_k**2 2475 2476 ! asserts at least 1 band of the current k and spin is on present processor 2477 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 2478 bdtot_index=bdtot_index+nband_k 2479 bd2tot_index=bd2tot_index+2*nband_k**2 2480 ! Skip the rest of the k-point loop 2481 cycle 2482 end if 2483 2484 ABI_MALLOC(occ_k,(nband_k)) 2485 2486 ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang)) 2487 ABI_MALLOC(ylm1_k,(npw1_k,mpsang*mpsang)) 2488 if (ipert==natom+3.or.ipert==natom+4) then 2489 ABI_MALLOC(ylmgr_k,(npw_k,3,mpsang*mpsang)) 2490 ABI_MALLOC(ylmgr1_k,(npw1_k,3,mpsang*mpsang)) 2491 end if 2492 2493 ! enl1_k(:)=zero 2494 d2nl_k(:,:,:)=zero 2495 if(prtbbb==1)d2bbb_k(:,:,:,:)=zero 2496 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 2497 2498 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 2499 if (psps%useylm==1) then 2500 do ilm=1,mpsang*mpsang 2501 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 2502 end do 2503 if (ipert==natom+3.or.ipert==natom+4) then 2504 do ilm=1,mpsang*mpsang 2505 do ii=1,3 2506 ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm) 2507 end do 2508 end do 2509 end if 2510 end if 2511 2512 wtk_k=wtk_rbz(ikpt) 2513 2514 kg1_k(:,:) = 0 2515 2516 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 2517 if (psps%useylm==1) then 2518 do ilm=1,mpsang*mpsang 2519 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 2520 end do 2521 if (ipert==natom+3.or.ipert==natom+4) then 2522 do ilm=1,mpsang*mpsang 2523 do ii=1,3 2524 ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm) 2525 end do 2526 end do 2527 end if 2528 end if 2529 2530 ! Compute the eigenvalues, wavefunction, 2531 ! contributions to kinetic energy, nonlocal energy, forces, 2532 ! and update of rhor1 to this k-point and this spin polarization. 2533 2534 ! Note that dfpt_nsteltwf is called with kpoint, while kpt is used inside dfpt_vtowfk 2535 call dfpt_nsteltwf(cg,cg1,d2nl_k,ecut,ecutsm,effmass_free,gs_hamk,icg,icg1,ikpt,isppol,& 2536 & istwf_k,kg_k,kg1_k,kpoint,mband,mband_mem,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,natom,nband_k,& 2537 & npw_k,npw1_k,nspinor,nsppol,ntypat,occ_k,psps,rmet,wtk_k,ylm_k,ylmgr_k) 2538 d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:) 2539 if(prtbbb==1)then 2540 d2bbb(:,:,idir,ipert,:,:) = d2bbb(:,:,idir,ipert,:,:) + & 2541 & d2bbb_k(:,:,:,:) 2542 end if 2543 2544 ABI_FREE(occ_k) 2545 2546 ! Keep track of total number of bands (all k points so far, even for 2547 ! k points not treated by me) 2548 bdtot_index=bdtot_index+nband_k 2549 bd2tot_index=bd2tot_index+2*nband_k**2 2550 2551 ! Shift array memory 2552 if (mkmem/=0) then 2553 icg=icg+npw_k*nspinor*proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 2554 ikg=ikg+npw_k 2555 end if 2556 if (mk1mem/=0) then 2557 icg1=icg1+npw1_k*nspinor*proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 2558 ikg1=ikg1+npw1_k 2559 end if 2560 ABI_FREE(ylm_k) 2561 ABI_FREE(ylm1_k) 2562 if (ipert==natom+3.or.ipert==natom+4) then 2563 ABI_FREE(ylmgr_k) 2564 ABI_FREE(ylmgr1_k) 2565 end if 2566 2567 end do ! End big k point loop 2568 end do ! End loop over spins 2569 2570 if(xmpi_paral==1)then 2571 call xmpi_barrier(comm) 2572 call wrtout(std_out,' dfpt_nselt: loop on k-points and spins done in parallel','COLL') 2573 end if 2574 2575 !Treat now varying occupation numbers 2576 !if(occopt>=3 .and. occopt <=8) then 2577 !SUPPRESSED metallic coding of vtorho 2578 2579 !Treat fixed occupation numbers 2580 !else 2581 2582 !Accumulation over parallel processed now carried out for all terms 2583 !in dfpt_nstdy.f 2584 2585 !End of test on varying or fixed occupation numbers 2586 !end if 2587 2588 !The imaginary part of d2nl will be must be set to zero here since 2589 !time-reversal symmetry will always be true for the strain peturbation. 2590 !The symmetry-reduced kpt set will leave a non-zero imaginary part. 2591 2592 d2nl(2,:,natom+3:natom+4,idir,ipert)=zero 2593 2594 !Symmetrize the non-local contributions, 2595 !as was needed for the stresses in a ground-state calculation 2596 2597 if (nsym1>1) then 2598 ! Pack like symmetric-storage cartesian stress tensor 2599 ii=0 2600 do ipert1=natom+3,natom+4 2601 do idir1=1,3 2602 ii=ii+1 2603 restr(ii)=d2nl(1,idir1,ipert1,idir,ipert) 2604 end do 2605 end do 2606 ! Do the symmetrization using the ground state routine 2607 call stresssym(gprimd,nsym1,restr,symrc1) 2608 ! Unpack symmetrized stress tensor 2609 ii=0 2610 do ipert1=natom+3,natom+4 2611 do idir1=1,3 2612 ii=ii+1 2613 d2nl(1,idir1,ipert1,idir,ipert)=restr(ii) 2614 end do 2615 end do 2616 end if !nsym>1 2617 2618 !---------------------------------------------------------------------------- 2619 !Now, treat the local contribution 2620 2621 ABI_MALLOC(vpsp1,(cplex*nfft)) 2622 n3xccc=0 2623 if(psps%n1xccc/=0)n3xccc=nfft 2624 ABI_MALLOC(xccc3d1,(cplex*n3xccc)) 2625 ABI_MALLOC(vxc1,(cplex*nfft,nspden)) 2626 ABI_MALLOC(vhartr01,(nfft)) 2627 xccc3d1(:)=zero 2628 2629 !To compute Absolute Deformation Potentials toghether with FxE tensor 2630 !the reference has to be the same as in the FxE routines 2631 rfstrs_ref_=0; if (present(rfstrs_ref)) rfstrs_ref_=rfstrs_ref 2632 g0term=0; if (rfstrs_ref_==1) g0term=1 2633 2634 !Double loop over strain perturbations 2635 do ipert1=natom+3,natom+4 2636 do idir1=1,3 2637 if(ipert1==natom+3) then 2638 istr1=idir1 2639 else 2640 istr1=idir1+3 2641 end if 2642 2643 ! Get first-order local potential. 2644 call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfft,mpi_enreg,& 2645 & psps%mqgrid_vl,natom,gs_hamk%nattyp,nfft,ngfft,ntypat,ph1d,psps%qgrid_vl,& 2646 & ucvol,psps%vlspl,vpsp1,g0term=g0term) 2647 2648 ! Get first-order hartree potential. 2649 call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,natom,nfft,ngfft,& 2650 & rhog,rprimd,vhartr01) 2651 2652 ! Get first-order exchange-correlation potential 2653 if(psps%n1xccc/=0)then 2654 call dfpt_mkcore(cplex,idir1,ipert1,natom,ntypat,n1,psps%n1xccc,& 2655 & n2,n3,qphon,rprimd,typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 2656 end if ! psps%n1xccc/=0 2657 2658 option=0 2659 call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,natom,nfft,ngfft,& 2660 & dummy,dummy,nkxc,nmxc,nspden,n3xccc,option,qphon,rhor,rhor1,rprimd,& 2661 & 0,0,vxc1,xccc3d1) 2662 2663 ! Combines density j2 with local potential j1 2664 do ispden=1,min(nspden,2) 2665 do ifft=1,cplex*nfft 2666 vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft)+vhartr01(ifft) 2667 end do 2668 end do 2669 call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol) 2670 write(std_out,*) 2671 d2lo(1,idir1,ipert1,idir,ipert)=dotr 2672 d2lo(2,idir1,ipert1,idir,ipert)=doti 2673 end do ! istr1 2674 end do ! ipert1 2675 2676 call gs_hamk%free() 2677 2678 ABI_FREE(vxc1) 2679 ABI_FREE(xccc3d1) 2680 ABI_FREE(vhartr01) 2681 2682 ABI_FREE(d2bbb_k) 2683 ABI_FREE(d2nl_k) 2684 ABI_FREE(kg_k) 2685 ABI_FREE(kg1_k) 2686 ABI_FREE(vpsp1) 2687 2688 end subroutine dfpt_nselt
ABINIT/dfpt_nstdy [ Functions ]
NAME
dfpt_nstdy
FUNCTION
This routine compute the non-stationary expression for the second derivative of the total energy, for a whole row of mixed derivatives. Only for norm-conserving pseudopotentials (no PAW)
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions at k cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q. cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues gmet(3,3)=reciprocal space metric tensor in bohr**-2. gsqcut=cutoff on (k+G)^2 (bohr^-2) idir=direction of the perturbation indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points indsy1(4,nsym1,natom)=indirect indexing array for atom labels ipert=type of the perturbation istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates. kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ kxc(nfft,nkxc)=exchange and correlation kernel mkmem =number of k points treated by this node (GS data) mk1mem =number of k points treated by this node (RF data) mpert =maximum number of ipert mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs). nattyp(ntypat)= # atoms of each type. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt=number of k points in the full BZ nkpt_rbz=number of k points in the reduced BZ for this perturbation nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym1=number of symmetry elements in space group consistent with i perturbation occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced Brillouin zone (usually =2) ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information psps <type(pseudopotential_type)>=variables related to pseudopotentials rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3. rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) symrc1(3,3,nsym1)=symmetry operations in reciprocal space ucvol=unit cell volume in bohr**3. wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
second order derivatives
d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs
NOTES
Note that the ddk perturbation should not be treated here.
SOURCE
2996 subroutine dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,& 2997 & gmet,gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mkmem,mk1mem,& 2998 & mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,nfft,ngfft,nkpt,nkpt_rbz,nkxc,& 2999 & npwarr,npwar1,nspden,nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,& 3000 & symrc1,ucvol,wtk_rbz,xred,ylm,ylm1,rhor,vxc) 3001 3002 !Arguments ------------------------------- 3003 !scalars 3004 integer,intent(in) :: cplex,idir,ipert,mk1mem,mkmem,mpert,mpw,mpw1,nfft,nkpt,nkpt_rbz,nkxc,nspden,nsppol,nsym1 3005 integer,intent(in) :: mband_mem_rbz 3006 real(dp),intent(in) :: gsqcut,ucvol 3007 type(MPI_type),intent(in) :: mpi_enreg 3008 type(datafiles_type),intent(in) :: dtfil 3009 type(dataset_type),intent(in) :: dtset 3010 type(pseudopotential_type),intent(in) :: psps 3011 !arrays 3012 integer,intent(in) :: atindx(dtset%natom),indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom) 3013 integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 3014 integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol),ngfft(18) 3015 integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz),symrc1(3,3,nsym1) 3016 integer,intent(inout) :: blkflg(3,mpert,3,mpert) !vz_i 3017 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband_mem_rbz*mkmem*nsppol) 3018 real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*nsppol) 3019 real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*nsppol) 3020 real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol) 3021 real(dp),intent(in) :: gmet(3,3),kpt_rbz(3,nkpt_rbz) 3022 real(dp),intent(in) :: kxc(nfft,nkxc),occ_rbz(dtset%mband*nkpt_rbz*nsppol) 3023 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 3024 real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3) 3025 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom) 3026 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 3027 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 3028 real(dp),intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*dtset%prtbbb)!vz_i 3029 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i 3030 ! optional 3031 real(dp),optional,intent(in) :: rhor(nfft,nspden) 3032 real(dp),optional,intent(in) :: vxc(cplex*nfft,nspden) 3033 3034 !Local variables------------------------------- 3035 !scalars 3036 integer,parameter :: formeig1=1 3037 integer :: ban2tot,bantot,bdtot_index,ddkcase,iband,icg,icg1,idir1 3038 integer :: ierr,ifft,ii,ikg,ikg1,ikpt,ilm,ipert1,ispden,isppol 3039 integer :: istwf_k,isym,jj,master,me,n1,n2,n3,n3xccc,n4,n5,n6 3040 integer :: nband_k,nfftot,npw1_k,npw_k,nspinor_,option,spaceworld,optnc 3041 real(dp) :: doti,dotr,wtk_k 3042 logical :: nmxc=.false.,t_exist 3043 character(len=500) :: msg 3044 character(len=fnlen) :: fiwfddk 3045 type(gs_hamiltonian_type) :: gs_hamkq 3046 !arrays 3047 integer :: ddkfil(3) 3048 integer,allocatable :: kg1_k(:,:),kg_k(:,:),symrl1(:,:,:) 3049 real(dp) :: d2nl_elfd(2,3),d2nl_mgfd(2,3),kpoint(3),kpq(3),sumelfd(2),summgfd(2),tsec(2) 3050 real(dp),allocatable :: buffer1(:),buffer2(:),d2bbb_k(:,:,:,:),d2nl_k(:,:,:) 3051 real(dp),allocatable :: eig1_k(:),eig_k(:),occ_k(:) 3052 real(dp) :: rhodummy(0,0) 3053 real(dp),allocatable :: vpsp1(:),vxc1(:,:),work1(:,:,:),xccc3d1(:),ylm1_k(:,:),ylm_k(:,:) 3054 type(pawtab_type) :: pawtab(dtset%ntypat*psps%usepaw) 3055 type(wfk_t) :: ddks(3) 3056 3057 3058 ! ********************************************************************* 3059 3060 ABI_UNUSED(nkpt) 3061 3062 DBG_ENTER("COLL") 3063 3064 !Not valid for PAW 3065 if (psps%usepaw==1) then 3066 msg='This routine cannot be used for PAW (use dfpt_nstpaw instead) !' 3067 ABI_BUG(msg) 3068 end if 3069 3070 3071 !Keep track of total time spent in dfpt_nstdy 3072 call timab(111,1,tsec) 3073 3074 !Init parallelism 3075 spaceworld=mpi_enreg%comm_cell 3076 me=mpi_enreg%me_kpt 3077 3078 master =0 3079 3080 !Zero only portion of nonlocal matrix to be computed here 3081 d2nl(:,:,1:dtset%natom+2,idir,ipert)=zero 3082 3083 ABI_MALLOC(d2bbb_k,(2,3,dtset%mband,dtset%mband*dtset%prtbbb)) 3084 ABI_MALLOC(d2nl_k,(2,3,mpert)) 3085 ABI_MALLOC(eig_k,(nsppol*dtset%mband)) 3086 ABI_MALLOC(eig1_k,(2*nsppol*dtset%mband**2)) 3087 ABI_MALLOC(kg_k,(3,mpw)) 3088 ABI_MALLOC(kg1_k,(3,mpw1)) 3089 3090 !Do not try to open electric field file 3091 ddkfil(:)=0 3092 !The treatment of homogeneous electric field potential need the existence of d/dk files. 3093 do idir1=1,3 3094 ddkcase=idir1+dtset%natom*3 3095 call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk) 3096 3097 ! Check that ddk file exists 3098 t_exist = file_exists(fiwfddk) 3099 if (.not. t_exist) then 3100 ! Try netcdf file. 3101 t_exist = file_exists(nctk_ncify(fiwfddk)) 3102 if (t_exist) then 3103 fiwfddk = nctk_ncify(fiwfddk) 3104 write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name." 3105 call wrtout(std_out,msg,'COLL') 3106 end if 3107 end if 3108 3109 if (t_exist) then 3110 ! Note the use of unit numbers 21, 22 and 23 3111 ! Open files in sequential mode 3112 ddkfil(idir1)=20+idir1 3113 write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk) 3114 call wrtout([std_out, ab_out], msg) 3115 call wfk_open_read(ddks(idir1),fiwfddk,formeig1,dtset%iomode,ddkfil(idir1), xmpi_comm_self) 3116 end if 3117 end do 3118 3119 !Update list of computed matrix elements 3120 if (ipert /= dtset%natom + 1) then 3121 do ipert1=1,mpert 3122 do idir1=1,3 3123 if(ipert1 <= dtset%natom .or. ipert1==dtset%natom+2 .and. ddkfil(idir1)/=0) then 3124 blkflg(idir1,ipert1,idir,ipert)=1 3125 end if 3126 end do 3127 end do 3128 else 3129 ipert1 = dtset%natom + 1 3130 do idir1=1,3 3131 ! If was already computed in another run or dataset, or if is to be computed in the present one 3132 if ((ddkfil(idir1) /= 0).or. (dtset%rfdir(idir1)/=0.and. idir1<=idir) ) then 3133 ! if ((ddkfil(idir1) /= 0).or. (idir1==idir) ) then 3134 blkflg(idir1,ipert1,idir,ipert)=1 3135 end if 3136 end do 3137 end if 3138 3139 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 3140 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 3141 nspinor_=dtset%nspinor 3142 3143 bantot = 0 3144 ban2tot = 0 3145 3146 !==== Initialize most of the Hamiltonian ==== 3147 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 3148 !2) Perform the setup needed for the non-local factors: 3149 !3) Constant kleimann-Bylander energies are copied from psps to gs_hamk. 3150 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,dtset%natom,& 3151 & dtset%typat,xred,nfft,dtset%mgfft,ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 3152 & gpu_option=dtset%gpu_option) 3153 3154 !LOOP OVER SPINS 3155 bdtot_index=0 3156 icg=0;icg1=0 3157 do isppol=1,nsppol 3158 3159 ikg=0;ikg1=0 3160 3161 ! Continue to initialize the Hamiltonian 3162 call gs_hamkq%load_spin(isppol,with_nonlocal=.true.) 3163 3164 ! BIG FAT k POINT LOOP 3165 do ikpt=1,nkpt_rbz 3166 3167 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 3168 istwf_k=istwfk_rbz(ikpt) 3169 npw_k=npwarr(ikpt) 3170 npw1_k=npwar1(ikpt) 3171 3172 eig_k(1:nband_k) = eigen0(1+bantot:nband_k+bantot) 3173 eig1_k(1:2*nband_k**2) = eigen1(1+ban2tot:2*nband_k**2+ban2tot) 3174 bantot = bantot + nband_k 3175 ban2tot = ban2tot + 2*nband_k**2 3176 3177 3178 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 3179 bdtot_index=bdtot_index+nband_k 3180 ! The wavefunction blocks for ddk file is skipped elsewhere in the loop 3181 ! Skip the rest of the k-point loop 3182 cycle 3183 end if 3184 3185 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 3186 ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)) 3187 3188 ! In case of electric field pert1, read ddk wfs file 3189 ! Note that the symmetries are not used for ddk, so read each k point 3190 ! Also take into account implicitly the parallelism over k points 3191 3192 do idir1=1,3 3193 if (ddkfil(idir1)/=0) then 3194 ii = ddks(idir1)%findk(kpt_rbz(:, ikpt)) 3195 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1") 3196 end if 3197 end do 3198 3199 ABI_MALLOC(occ_k,(nband_k)) 3200 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 3201 kpoint(:)=kpt_rbz(:,ikpt) 3202 kpq(:)=kpoint(:)+dtset%qptn(:) 3203 wtk_k=wtk_rbz(ikpt) 3204 d2nl_k(:,:,:)=zero 3205 if(dtset%prtbbb==1)d2bbb_k(:,:,:,:)=zero 3206 3207 ! Get plane-wave vectors and related data at k 3208 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 3209 if (psps%useylm==1) then 3210 do ilm=1,psps%mpsang*psps%mpsang 3211 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 3212 end do 3213 end if 3214 3215 ! Get plane-wave vectors and related data at k+q 3216 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 3217 if (psps%useylm==1) then 3218 do ilm=1,psps%mpsang*psps%mpsang 3219 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 3220 end do 3221 end if 3222 3223 ! Compute the eigenvalues, wavefunction, 3224 ! contributions to kinetic energy, nonlocal energy, forces, 3225 ! and update of rhor1 to this k-point and this spin polarization. 3226 ! Note that dfpt_nstwf is called with kpoint, while kpt is used inside dfpt_vtowfk 3227 call dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,& 3228 & icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpoint,kpq,mband_mem_rbz,mkmem,mk1mem,mpert,& 3229 & mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,& 3230 & occ_k,psps,rmet,ddks,wtk_k,ylm_k,ylm1_k) 3231 3232 d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:) 3233 if(dtset%prtbbb==1)d2bbb(:,:,idir,ipert,:,:)=d2bbb(:,:,idir,ipert,:,:)+d2bbb_k(:,:,:,:) 3234 3235 ! Keep track of total number of bands 3236 bdtot_index=bdtot_index+nband_k 3237 3238 ! Shift arrays memory 3239 if (mkmem/=0) then 3240 icg=icg+npw_k*dtset%nspinor*proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 3241 ikg=ikg+npw_k 3242 end if 3243 if (mk1mem/=0) then 3244 icg1=icg1+npw1_k*dtset%nspinor*proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 3245 ikg1=ikg1+npw1_k 3246 end if 3247 3248 ABI_FREE(occ_k) 3249 ABI_FREE(ylm_k) 3250 ABI_FREE(ylm1_k) 3251 3252 ! End big k point loop 3253 end do 3254 3255 ! End loop over spins 3256 end do 3257 3258 call gs_hamkq%free() 3259 3260 !Treat fixed occupation numbers (as in vtorho) 3261 if(xmpi_paral==1)then 3262 ABI_MALLOC(buffer1,(2*3*mpert)) 3263 ABI_MALLOC(buffer2,(2*3*mpert)) 3264 ! Pack d2nl 3265 buffer1(1:2*3*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/2*3*mpert/)) 3266 ! Build sum of everything 3267 call timab(48,1,tsec) 3268 call xmpi_sum(buffer1,buffer2,2*3*mpert,spaceworld,ierr) 3269 call timab(48,2,tsec) 3270 ! Unpack the final result 3271 d2nl(:,:,:,idir,ipert)=reshape(buffer2(:),(/2,3,mpert/)) 3272 ABI_FREE(buffer1) 3273 ABI_FREE(buffer2) 3274 3275 if(dtset%prtbbb==1)then 3276 ABI_MALLOC(buffer1,(2*3*dtset%mband*dtset%mband)) 3277 ABI_MALLOC(buffer2,(2*3*dtset%mband*dtset%mband)) 3278 ! Pack d2bbb 3279 buffer1(1:2*3*dtset%mband*dtset%mband)=reshape(d2bbb(:,:,idir,ipert,:,:),(/2*3*dtset%mband*dtset%mband/)) 3280 ! Build sum of everything 3281 call timab(48,1,tsec) 3282 call xmpi_sum(buffer1,buffer2,2*3*dtset%mband*dtset%mband,spaceworld,ierr) 3283 call timab(48,2,tsec) 3284 ! Unpack the final result 3285 d2bbb(:,:,idir,ipert,:,:)=reshape(buffer2(:),(/2,3,dtset%mband,dtset%mband/)) 3286 ABI_FREE(buffer1) 3287 ABI_FREE(buffer2) 3288 end if 3289 end if ! xmpi_paral==1 3290 3291 !In the case of the strain perturbation time-reversal symmetry will always 3292 !be true so imaginary part of d2nl will be must be set to zero here since 3293 !the symmetry-reduced kpt set will leave a non-zero imaginary part. 3294 if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) d2nl(2,:,:,idir,ipert)=zero 3295 3296 !In case of electric field ipert1, close the ddk wf files 3297 do idir1=1,3 3298 if (ddkfil(idir1)/=0) call ddks(idir1)%close() 3299 end do 3300 3301 !Symmetrize the non-local contributions, 3302 !as was needed for the forces in a ground-state calculation 3303 !However, here the quantity is complex, and there are phases ! 3304 3305 !Do the transform 3306 ABI_MALLOC(work1,(2,3,dtset%natom)) 3307 do ipert1=1,dtset%natom 3308 do idir1=1,3 3309 work1(1,idir1,ipert1)=d2nl(1,idir1,ipert1,idir,ipert) 3310 work1(2,idir1,ipert1)=d2nl(2,idir1,ipert1,idir,ipert) 3311 end do 3312 end do 3313 call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work1,indsy1,ipert,nsym1,dtset%qptn,symrc1) 3314 ABI_FREE(work1) 3315 3316 !Must also symmetrize the electric/magnetic field perturbation response ! 3317 !(XG 000803 This was not implemented until now) 3318 if(sum(ddkfil(:))/=0)then 3319 ! Get the symmetry matrices in terms of real space basis 3320 ABI_MALLOC(symrl1,(3,3,nsym1)) 3321 do isym=1,nsym1 3322 call mati3inv(symrc1(:,:,isym),symrl1(:,:,isym)) 3323 end do 3324 ! There should not be any imaginary part, but stay general (for debugging) 3325 d2nl_elfd(:,:)=d2nl(:,:,dtset%natom+2,idir,ipert) 3326 do ii=1,3 3327 sumelfd(:)=zero 3328 summgfd(:)=zero 3329 do isym=1,nsym1 3330 do jj=1,3 3331 if(symrl1(ii,jj,isym)/=0)then 3332 if(ddkfil(jj)==0)then 3333 blkflg(ii,dtset%natom+2,idir,ipert)=0 3334 end if 3335 end if 3336 end do 3337 sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+& 3338 & dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+& 3339 & dble(symrl1(ii,3,isym))*d2nl_elfd(:,3) 3340 summgfd(:)=summgfd(:)+dble(symrl1(ii,1,isym))*d2nl_mgfd(:,1)+& 3341 & dble(symrl1(ii,2,isym))*d2nl_mgfd(:,2)+& 3342 & dble(symrl1(ii,3,isym))*d2nl_mgfd(:,3) 3343 end do 3344 d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1) 3345 end do 3346 3347 if ((dtset%prtbbb==1).and.(ipert<=dtset%natom)) then 3348 do iband = 1,dtset%mband 3349 d2nl_elfd(:,:)=d2bbb(:,:,idir,ipert,iband,iband) 3350 do ii=1,3 3351 sumelfd(:)=zero 3352 do isym=1,nsym1 3353 sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+& 3354 & dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+& 3355 & dble(symrl1(ii,3,isym))*d2nl_elfd(:,3) 3356 end do 3357 d2bbb(:,ii,idir,ipert,iband,iband)=sumelfd(:)/dble(nsym1) 3358 end do 3359 end do !iband 3360 end if 3361 3362 ABI_FREE(symrl1) 3363 end if 3364 3365 !---------------------------------------------------------------------------- 3366 !Now, treat the local contribution 3367 3368 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 3369 ABI_MALLOC(vpsp1,(cplex*nfft)) 3370 if (ipert /= dtset%natom + 1) then 3371 n3xccc=0;if(psps%n1xccc/=0) n3xccc=nfft 3372 ABI_MALLOC(xccc3d1,(cplex*n3xccc)) 3373 ABI_MALLOC(vxc1,(cplex*nfft,nspden)) 3374 3375 do ipert1=1,mpert 3376 do idir1=1,3 3377 if(ipert1 <= dtset%natom)then 3378 3379 ! Get first-order local potential and first-order pseudo core density 3380 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir1,ipert1,mpi_enreg,psps%mqgrid_ff,dtset%natom,& 3381 & nattyp,nfft,ngfft,dtset%ntypat,n1,n2,n3,ph1d,psps%qgrid_ff,& 3382 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 3383 if(psps%n1xccc/=0)then 3384 call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,n1,psps%n1xccc,& 3385 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 3386 end if 3387 3388 ! Get first-order exchange-correlation potential (core-correction contribution only !) 3389 if(psps%n1xccc/=0)then 3390 option=0 3391 !FR SPr EB non-collinear magnetism 3392 if (nspden==4.and.present(rhor).and.present(vxc)) then 3393 optnc=1 3394 call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,rhodummy,0,& 3395 & nkxc,nmxc,nspden,n3xccc,optnc,option,dtset%qptn,rhor,rhor1,& 3396 & rprimd,0,vxc,vxc1,xccc3d1) 3397 else 3398 call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,& 3399 & nkxc,nmxc,nspden,n3xccc,option,dtset%qptn,rhodummy,& 3400 & rprimd,0,vxc1,xccc3d1) 3401 end if 3402 else 3403 vxc1(:,:)=zero 3404 end if 3405 3406 ! Norm-conserving pseudpopotential case: 3407 ! Combines density j2 with local potential j1 (vpsp1 and vxc1) 3408 ! XG030514 : this is a first possible coding, however, each dotprod contains 3409 ! a parallel section (reduction), so it is better to use only one dotprod ... 3410 ! call dotprod_vn(cplex,rhor1,dr_psp1,di_psp1,mpi_enreg,nfft,nfftot,1,2,vpsp1,ucvol) 3411 ! call dotprod_vn(cplex,rhor1,dr_xc1,di_xc1,mpi_enreg,nfft,nfftot,nspden,2,vxc1,ucvol) 3412 ! dotr=dr_psp1+dr_xc1;doti=di_psp1+di_xc1... but then, one needs to overload vxc1 3413 do ispden=1,min(nspden,2) 3414 do ifft=1,cplex*nfft 3415 vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft) 3416 end do 3417 end do 3418 call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol) 3419 3420 ! MVeithen 021212 : in case ipert = 2, these lines compute the local part 3421 ! of the Born effective charges from phonon and electric 3422 ! field type perturbations, see eq. 43 of 3423 ! X. Gonze and C. Lee, PRB 55, 10355 (1997) [[cite:Gonze1997a]] 3424 ! The minus sign is due to the fact that the effective charges 3425 ! are minus the second derivatives of the energy 3426 if (ipert == dtset%natom+2) then 3427 d2lo(1,idir1,ipert1,idir,ipert)=-dotr 3428 d2lo(2,idir1,ipert1,idir,ipert)=-doti 3429 else 3430 d2lo(1,idir1,ipert1,idir,ipert)=dotr 3431 d2lo(2,idir1,ipert1,idir,ipert)=doti 3432 end if 3433 ! Endif ipert1<=natom 3434 end if 3435 end do 3436 end do 3437 3438 ABI_FREE(vxc1) 3439 ABI_FREE(xccc3d1) 3440 3441 end if ! ipert /= natom +1 3442 3443 ABI_FREE(d2bbb_k) 3444 ABI_FREE(d2nl_k) 3445 ABI_FREE(kg_k) 3446 ABI_FREE(kg1_k) 3447 ABI_FREE(vpsp1) 3448 ABI_FREE(eig_k) 3449 ABI_FREE(eig1_k) 3450 3451 call timab(111,2,tsec) 3452 3453 DBG_EXIT("COLL") 3454 3455 end subroutine dfpt_nstdy
ABINIT/dfpt_nsteltwf [ Functions ]
NAME
dfpt_nsteltwf
FUNCTION
This routine computes the non-local and kinetic contribution to the 2DTE matrix elements, in the non-stationary formulation
INPUTS
cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions
cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)=pw coefficients of RF
wavefunctions at k,q.
ecut=cut-off energy for plane wave basis sphere (Ha)
ecutsm=smearing energy for plane wave kinetic energy (Ha) (NOT NEEDED !)
effmass_free=effective mass for electrons (1. in common case)
gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
icg=shift to be applied on the location of data in the array cg
icg1=shift to be applied on the location of data in the array cg1
ikpt=number of the k-point
isppol=1 for unpolarized, 2 for spin-polarized
istwf_k=flag controlling the storage of WFs
kg_k(3,npw_k)=reduced planewave coordinates.
kg1_k(3,npw1_k)=reduced planewave coordinates at k+q, with RF k points
kpoint(3)=k-point in reduced coordinates
mband=maximum number of bands
mband_mem=maximum number of bands on this cpu
mkmem =number of k points treated by this node.
mk1mem =number of k points treated by this node (RF data).
mpert =maximum number of ipert
mpi_enreg=information about MPI parallelization
mpw=maximum dimensioned size of npw or wfs at k
mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
natom=number of atoms in cell.
nband_k=number of bands at this k point for that spin polarization
npw_k=number of plane waves at this k point
npw1_k=number of plane waves at this k+q point
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
ntypat=number of types of atoms in unit cell.
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
psps <type(pseudopotential_type)>=variables related to pseudopotentials
rmet(3,3)=real space metric (bohr**2)
wtk_k=weight assigned to the k point.
ylm(npw_k,mpsang*mpsang)= real spherical harmonics for each G and k point
ylmgr(npw_k,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k point
OUTPUT
d2nl_k(2,3,mpert)=non-local contributions to non-stationary 2DTE, for the present k point, and perturbation idir, ipert
SOURCE
2743 subroutine dfpt_nsteltwf(cg,cg1,d2nl_k,ecut,ecutsm,effmass_free,gs_hamk,icg,icg1,ikpt,isppol,& 2744 & istwf_k,kg_k,kg1_k,kpoint,mband,mband_mem,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,natom,nband_k,& 2745 & npw_k,npw1_k,nspinor,nsppol,ntypat,occ_k,psps,rmet,wtk_k,ylm,ylmgr) 2746 2747 2748 2749 !Arguments ------------------------------------ 2750 !scalars 2751 integer,intent(in) :: icg,icg1,ikpt,isppol,istwf_k,mband,mk1mem,mkmem,mpert,mpw,mpw1,natom 2752 integer,intent(in) :: mband_mem 2753 integer,intent(in) :: nspinor,nsppol,ntypat 2754 integer,intent(inout) :: nband_k,npw1_k,npw_k 2755 real(dp),intent(in) :: ecut,ecutsm,effmass_free,wtk_k 2756 type(MPI_type),intent(in) :: mpi_enreg 2757 type(pseudopotential_type),intent(in) :: psps 2758 !arrays 2759 integer,intent(in) :: kg1_k(3,npw1_k),kg_k(3,npw_k) 2760 real(dp),intent(in) :: kpoint(3) 2761 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol) 2762 real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol) 2763 real(dp),intent(in) :: occ_k(nband_k),rmet(3,3) 2764 real(dp),intent(in) :: ylm(npw_k,psps%mpsang*psps%mpsang) 2765 real(dp),intent(in) :: ylmgr(npw_k,3,psps%mpsang*psps%mpsang) 2766 real(dp),intent(inout) :: d2nl_k(2,3,mpert) 2767 2768 !Local variables------------------------------- 2769 !scalars 2770 integer :: choice,cpopt,dimffnl,dimffnl2,iband 2771 integer :: iband_me 2772 integer :: ider,idir0,idir1,ilmn,ipert1,ipw,ipws,ispinor,istr1,itypat 2773 integer :: nkpg,nnlout,paw_opt,signs,tim_nonlop 2774 real(dp) :: doti,dotr 2775 type(gs_hamiltonian_type) :: gs_hamk 2776 !arrays 2777 real(dp) :: enlout(6),dum_svectout(1,1),dum(1),kpg_dum(0,0) 2778 real(dp),allocatable :: cwave0(:,:),cwavef(:,:),dkinpw(:),eig2_k(:) 2779 real(dp),allocatable :: ffnl(:,:,:,:),ffnl_ylm(:,:,:,:),ghc(:,:) 2780 real(dp),allocatable :: gvnlx1(:,:),gvnlxc(:,:),kinpw1(:),ph3d(:,:,:) 2781 type(pawcprj_type) :: cprj_dum(0,0) 2782 2783 ! ********************************************************************* 2784 2785 !Init me 2786 ABI_MALLOC(ghc,(2,npw1_k*nspinor)) 2787 ABI_MALLOC(gvnlxc,(2,npw1_k*nspinor)) 2788 ABI_MALLOC(gvnlx1,(2,npw1_k*nspinor)) 2789 ABI_MALLOC(eig2_k,(2*nsppol*mband**2)) 2790 ABI_MALLOC(kinpw1,(npw1_k)) 2791 ABI_MALLOC(dkinpw,(npw_k)) 2792 nkpg=0 2793 2794 !Compute nonlocal form factors ffnl at (k+G), for all atoms 2795 dimffnl=2 2796 ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat)) 2797 if (psps%useylm==0) then 2798 ider=1;idir0=0 2799 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,ider,idir0,& 2800 & psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 2801 & npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,ylmgr) 2802 else 2803 ider=1;idir0=-7;dimffnl2=7 2804 ABI_MALLOC(ffnl_ylm,(npw_k,dimffnl2,psps%lmnmax,ntypat)) 2805 call mkffnl(psps%dimekb,dimffnl2,psps%ekb,ffnl_ylm,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,& 2806 & ider,idir0,psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,& 2807 & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,ylmgr) 2808 do itypat=1,ntypat 2809 do ilmn=1,psps%lmnmax 2810 ffnl(:,1,ilmn,itypat)=ffnl_ylm(:,1,ilmn,itypat) 2811 end do 2812 end do 2813 end if 2814 2815 !Compute kinetic contributions (1/2) (2 Pi)**2 (k+G)**2: 2816 ! call mkkin(ecut,ecutsm,effmass_free,gs_hamk%gmet,kg1_k,kinpw1,kpoint,npw1_k) 2817 call mkkin(ecut,ecutsm,effmass_free,gs_hamk%gmet,kg1_k,kinpw1,kpoint,npw1_k,0,0) 2818 2819 !Load k/k+q-dependent part in the Hamiltonian datastructure 2820 ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk)) 2821 call gs_hamk%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,ffnl_k=ffnl,& 2822 & ph3d_k=ph3d,compute_ph3d=.true.) 2823 2824 ABI_MALLOC(cwave0,(2,npw_k*nspinor)) 2825 ABI_MALLOC(cwavef,(2,npw1_k*nspinor)) 2826 2827 !Loop over bands 2828 iband_me = 0 2829 do iband=1,nband_k 2830 2831 if(mpi_enreg%proc_distrb(ikpt, iband, isppol) /= mpi_enreg%me_kpt) then 2832 ! Skip the eigenvalue and the gvnl records of this band 2833 cycle 2834 end if 2835 iband_me = iband_me + 1 2836 2837 ! Get ground-state and first-order wavefunctions 2838 cwave0(:,:)=cg(:,1+(iband_me-1)*npw_k*nspinor+icg:iband*npw_k*nspinor+icg) 2839 cwavef(:,:)=cg1(:,1+(iband_me-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1) 2840 2841 ! Double loop over strain perturbations 2842 do ipert1=natom+3,natom+4 2843 do idir1=1,3 2844 if (ipert1==natom+3) istr1=idir1 2845 if (ipert1==natom+4) istr1=idir1+3 2846 2847 ! Compute the derivative of the kinetic operator vs strain in dkinpw 2848 call kpgstr(dkinpw,ecut,ecutsm,effmass_free,gs_hamk%gmet,gs_hamk%gprimd,istr1,& 2849 & kg1_k,kpoint,npw1_k) 2850 2851 ! Get |vnon-locj1|u0> : 2852 ! first-order non-local, applied to zero-order wavefunction 2853 ! (??) this routine gives MINUS the non-local contribution 2854 2855 ! When using Ylms, load the correct ffnl derivative 2856 if (psps%useylm==1) then 2857 do itypat=1,ntypat 2858 do ilmn=1,psps%lmnmax 2859 ffnl(:,2,ilmn,itypat)=ffnl_ylm(:,1+istr1,ilmn,itypat) 2860 end do 2861 end do 2862 end if 2863 2864 signs=2 ; choice=3 ; nnlout=6 ; paw_opt=0 ; cpopt=-1 ; tim_nonlop=5 2865 call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,istr1,dum,mpi_enreg,1,nnlout,paw_opt,& 2866 & signs,dum_svectout,tim_nonlop,cwave0,gvnlx1) 2867 ! <G|Vnl1|Cnk> is contained in gvnlx1 2868 2869 ! Kinetic contribution 2870 do ispinor=1,nspinor 2871 do ipw=1,npw1_k 2872 ipws=ipw+npw1_k*(ispinor-1) 2873 if(kinpw1(ipw)<huge(zero)*1.d-11)then 2874 gvnlx1(1,ipws)=gvnlx1(1,ipws)+dkinpw(ipw)*cwave0(1,ipws) 2875 gvnlx1(2,ipws)=gvnlx1(2,ipws)+dkinpw(ipw)*cwave0(2,ipws) 2876 else 2877 gvnlx1(1,ipws)=0.0_dp 2878 gvnlx1(2,ipws)=0.0_dp 2879 end if 2880 end do 2881 end do 2882 2883 ! construct the matrix element (<uj2|vj1|u0>)complex conjug. 2884 ! and add it to the 2nd-order matrix 2885 ! imaginary term should be zero for strain-strain 2nd derivatives, 2886 ! but keep it as a test for now 2887 call dotprod_g(dotr,doti,gs_hamk%istwf_k,npw1_k*nspinor,2,cwavef,gvnlx1,& 2888 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2889 2890 d2nl_k(1,idir1,ipert1)= d2nl_k(1,idir1,ipert1)+wtk_k*occ_k(iband)*2.0_dp*dotr 2891 d2nl_k(2,idir1,ipert1)= d2nl_k(2,idir1,ipert1)-wtk_k*occ_k(iband)*2.0_dp*doti 2892 2893 end do !idir1 2894 end do !ipert1 2895 2896 ! UNTIL NOW, DO NOT TAKE INTO ACCOUNT istwf_k 2897 2898 ! End loop over bands 2899 end do 2900 2901 ABI_FREE(cwave0) 2902 ABI_FREE(cwavef) 2903 2904 !################################################################### 2905 2906 ABI_FREE(eig2_k) 2907 ABI_FREE(ghc) 2908 ABI_FREE(gvnlxc) 2909 ABI_FREE(gvnlx1) 2910 ABI_FREE(kinpw1) 2911 ABI_FREE(dkinpw) 2912 ABI_FREE(ffnl) 2913 ABI_FREE(ph3d) 2914 if (psps%useylm==1) then 2915 ABI_FREE(ffnl_ylm) 2916 end if 2917 2918 end subroutine dfpt_nsteltwf
ABINIT/dfpt_rhofermi [ Functions ]
NAME
dfpt_rhofermi
FUNCTION
This routine computes the fixed contribution to the first-order Fermi energy for metallic occupation and Q=0, as well as the Fermi level charge density needed to compute the remainder of the first-order Fermi energy from the self-consistent local potential at each step in the iteration process.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions. cgq(2,mpw1*nspinor*mband_mem*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q. cplex: if 1, real space 1-order functions on FFT grid are REAL; if 2, COMPLEX cprj(natom,nspinor*mband_mem*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors: cprj=<p_i|Cnk> cprjq(natom,nspinor*mband_mem*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors: cprjq=<p_i|Cnk+q> doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree) eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) gmet(3,3)=reciprocal space metric tensor in bohr**-2. gprimd(3,3)=dimensional reciprocal space primitive translations idir=direction of the perturbation indsy1(4,nsym1,natom)=indirect indexing array for atom labels ipert=type of the perturbation irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates. kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points. mband=maximum number of bands mband_mem=maximum number of bands on this cpu mkmem =number of k points treated by this node (GS data) mkqmem =number of k+q points treatede by this node (GS data) mk1mem =number of k points treated by this node. mpw=maximum dimensioned size of npw or wfs at k mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs). maximum dimension for q points in grids for nonlocal form factors natom=number of atoms in cell. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>) nfftf= -PAW ONLY- number of FFT grid points for the fine grid (nfftf=nfft for norm-conserving potential runs - see comment in respfn.F90) ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid nhatfermi(nfft,nspden)=array for fermi-level compensation charge density (PAW only) nkpt_rbz=number of k points in the IBZ for this perturbation mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym1=number of symmetry elements in space group consistent with perturbation occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2) at each k+q point of the reduced Brillouin zone. occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k (usually 2) paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawang <type(pawang_type)>=paw angular mesh and related data pawang1 <type(pawang_type)>=pawang datastr. containing only symmetries preserving the perturbation pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data phnons1(2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information prtvol=control print volume and debugging output psps <type(pseudopotential_type)>=variables related to pseudopotentials rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional real space primitive translations symaf1(nsym1)=(anti)ferromagnetic part of symmetry operations symrc1(3,3,nsym1)=symmetry operations in reciprocal space symrl1(3,3,nsym1)=3x3 matrices of the group symmetries tnons1(3,nsym1)=non-symmorphic translations ucvol=volume of the unit cell usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory useylmgr1= 1 if ylmgr1 array is allocated vtrial(nfftf,nspden)=GS potential (Hartree). vxc(nfftf,nspden)=XC potential (Hartree). wtk_rbz(nkpt_rbz)=weight assigned to each k point. xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= spherical harmonics for each G and k+g point ylmgr1(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k+q
OUTPUT
eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues (hartree) - only digonal elements computed here fe1fixed=fixed contribution to the first-order Fermi energy (nonlocal and kinetic in the case of strain) nhatfermi(cplex*nfftf,nspden)=fermi-level compensation charge density (PAW only) rhorfermi(cplex*nfftf,nspden)=fermi-level electronic density
NOTES
This routine will NOT work with nspden==4:
at least the use of fftpac should be modified.
SOURCE
3564 subroutine dfpt_rhofermi(cg,cgq,cplex,cprj,cprjq,& 3565 & doccde_rbz,docckqde,dtfil,dtset,eigenq,eigen0,eigen1,fe1fixed,gmet,gprimd,idir,& 3566 & indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,mband,mband_mem,mkmem,mkqmem,mk1mem,mpi_enreg,& 3567 & mpw,mpw1,my_natom,natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1,nspden,& 3568 & nsppol,nsym1,occkq,occ_rbz,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,& 3569 & phnons1,ph1d,prtvol,psps,rhorfermi,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,& 3570 & ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1) 3571 3572 !Arguments ------------------------------- 3573 !scalars 3574 integer,intent(in) :: cplex,idir,ipert,mband,mk1mem,mkmem,mkqmem 3575 integer,intent(in) :: mband_mem 3576 integer,intent(in) :: mpw,mpw1,my_natom,natom,ncpgr,nfftf,nkpt_rbz,nspden,nsppol,nsym1 3577 integer,intent(in) :: prtvol,usecprj,useylmgr1 3578 real(dp),intent(in) :: ucvol 3579 real(dp),intent(out) :: fe1fixed 3580 type(MPI_type),intent(in) :: mpi_enreg 3581 type(datafiles_type),intent(in) :: dtfil 3582 type(dataset_type),intent(in) :: dtset 3583 type(pawang_type),intent(in) :: pawang,pawang1 3584 type(pawfgr_type),intent(in) :: pawfgr 3585 type(pseudopotential_type),intent(in) :: psps 3586 !arrays 3587 integer,intent(in) :: indsy1(4,nsym1,natom) 3588 integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4)) 3589 integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 3590 integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfftf(18) 3591 integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz),symaf1(nsym1) 3592 integer,intent(in) :: symrc1(3,3,nsym1),symrl1(3,3,nsym1) 3593 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband_mem*mkmem*nsppol) 3594 real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*mband_mem*mkqmem*nsppol) 3595 real(dp),intent(in) :: doccde_rbz(mband*nkpt_rbz*nsppol) 3596 real(dp),intent(in) :: docckqde(mband*nkpt_rbz*nsppol) 3597 real(dp),intent(in) :: eigen0(mband*nkpt_rbz*nsppol) 3598 real(dp),intent(in) :: eigenq(mband*nkpt_rbz*nsppol),gmet(3,3),gprimd(3,3) 3599 real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz) 3600 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol),occkq(mband*nkpt_rbz*nsppol) 3601 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*natom) 3602 real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4)) 3603 real(dp),intent(in) :: rmet(3,3),rprimd(3,3),tnons1(3,nsym1) 3604 real(dp),intent(in) :: vtrial(nfftf,nspden),vxc(nfftf,nspden),wtk_rbz(nkpt_rbz) 3605 real(dp),intent(in) :: xred(3,natom),ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 3606 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 3607 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 3608 real(dp),intent(out) :: eigen1(2*mband*mband*nkpt_rbz*nsppol) 3609 real(dp),intent(out) :: nhatfermi(:,:) 3610 real(dp),intent(out) :: rhorfermi(cplex*nfftf,nspden) 3611 type(pawcprj_type),intent(in) :: cprj (natom,dtset%nspinor*mband_mem*mkmem *nsppol*usecprj) 3612 type(pawcprj_type),intent(in) :: cprjq(natom,dtset%nspinor*mband_mem*mkqmem*nsppol*usecprj) 3613 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 3614 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 3615 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw) 3616 type(pawrhoij_type),target,intent(inout)::pawrhoijfermi(my_natom*psps%usepaw)!vz_i 3617 type(pawtab_type), intent(in) :: pawtab(dtset%ntypat*psps%usepaw) 3618 3619 !Local variables------------------------------- 3620 !scalars 3621 integer,parameter :: level=17 3622 integer :: bd2tot_index,bdtot_index,buffer_size,cplex_rhoij 3623 integer :: dimffnl1,dimffnlk,iatom,iband,ibg,ibgq 3624 integer :: icg,icgq,ider,idir0,ierr,ii,ikg,ikg1,ikpt,ilm,ilmn,indx 3625 integer :: ispden,isppol,istr,istwf_k 3626 integer :: mbd2kpsp,mcgq,mcgq_disk,mcprjq,mcprjq_disk 3627 integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nkpg1,npw1_k,npw_k,nspden_rhoij 3628 integer :: optfr,qphase_rhoij,spaceworld 3629 integer :: nband_me 3630 logical :: paral_atom,qne0 3631 real(dp) :: arg,fe1norm,invfe1norm,wtk_k 3632 type(gs_hamiltonian_type) :: gs_hamkq 3633 type(rf_hamiltonian_type) :: rf_hamkq 3634 !arrays 3635 integer,allocatable :: kg1_k(:,:),kg_k(:,:) 3636 real(dp) :: kpoint(3),kpq(3),tsec(2) 3637 real(dp) :: ylmgr_dum(1,1,1) 3638 real(dp),allocatable :: buffer1(:),dkinpw(:),doccde_k(:) 3639 real(dp),allocatable :: doccde_kq(:),eig0_k(:),eig0_kq(:),eig1_k(:) 3640 real(dp),allocatable :: fe1fixed_k(:),fe1norm_k(:) 3641 real(dp),allocatable :: ffnl1(:,:,:,:),ffnlk(:,:,:,:) 3642 real(dp),allocatable :: kinpw1(:),kpg1_k(:,:),kpg_k(:,:) 3643 real(dp),allocatable :: occ_k(:),occ_kq(:),ph3d(:,:,:),ph3d1(:,:,:) 3644 real(dp),allocatable :: rhoaug(:,:,:),rhogfermi(:,:),rhowfr(:,:) 3645 real(dp),allocatable :: rhoaug4(:,:,:,:) 3646 real(dp),allocatable :: rocceig(:,:),ylm1_k(:,:),ylm_k(:,:),ylmgr1_k(:,:,:) 3647 type(paw_ij_type),allocatable :: paw_ij1fr(:) 3648 type(pawrhoij_type),pointer :: pawrhoijfermi_unsym(:) 3649 ! real(dp),allocatable :: vlocal1(:,:,:,:),vlocal_tmp(:,:,:,:) 3650 ! real(dp),allocatable :: v1zeeman(:,:),vtrial_tmp(:,:) 3651 3652 ! ********************************************************************* 3653 3654 DBG_ENTER('COLL') 3655 3656 !Check arguments validity 3657 if (ipert>natom.and.ipert/=natom+3.and.ipert/=natom+4.and.ipert/=natom+5) then 3658 ABI_BUG('wrong ipert argument!') 3659 end if 3660 if (cplex/=1) then 3661 ABI_BUG('wrong cplex/=1 argument !') 3662 end if 3663 3664 !Keep track of total time spent in this routine 3665 call timab(121,1,tsec) 3666 call timab(124,1,tsec) 3667 3668 !Retrieve parallelism data 3669 spaceworld=mpi_enreg%comm_cell 3670 me=mpi_enreg%me_kpt 3671 paral_atom=(my_natom/=dtset%natom) 3672 3673 !Initialize output variables 3674 fe1fixed=zero 3675 if (psps%usepaw==0) rhorfermi(:,:)=zero 3676 3677 !Initialisations/allocation of temporary variables 3678 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 3679 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 3680 bdtot_index=0 ; bd2tot_index=0 ; ibg=0 ; ibgq=0 ; icg=0 ; icgq=0 3681 qne0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2>=tol14) 3682 mbd2kpsp=2*mband**2*nkpt_rbz*nsppol 3683 fe1norm=zero 3684 if (nspden/=4) then 3685 ABI_MALLOC(rhoaug,(cplex*n4,n5,n6)) 3686 else 3687 ABI_MALLOC(rhoaug4,(cplex*n4,n5,n6,nspden)) 3688 end if 3689 ABI_MALLOC(kg_k,(3,mpw)) 3690 ABI_MALLOC(kg1_k,(3,mpw1)) 3691 if (psps%usepaw==1) then 3692 ABI_MALLOC(rhowfr,(cplex*dtset%nfft,dtset%nspden)) 3693 rhowfr(:,:)=zero 3694 end if 3695 3696 mcgq=mpw1*dtset%nspinor*mband_mem*mkqmem*nsppol;mcgq_disk=0 3697 3698 !Prepare RF PAW files for reading and writing if mkmem, mkqmem or mk1mem==0 3699 if (psps%usepaw==1) then 3700 mcprjq=dtset%nspinor*mband_mem*mkqmem*nsppol*usecprj;mcprjq_disk=0 3701 else 3702 mcprjq=0;mcprjq_disk=0 3703 end if 3704 3705 !PAW:has to compute frozen part of Dij^(1) (without Vpsp(1) contribution) 3706 if (psps%usepaw==1) then 3707 ABI_MALLOC(paw_ij1fr,(my_natom)) 3708 call paw_ij_nullify(paw_ij1fr) 3709 call paw_ij_init(paw_ij1fr,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,& 3710 & dtset%natom,dtset%ntypat,dtset%typat,pawtab,has_dijfr=1,& 3711 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom ) 3712 optfr=1 3713 ABI_MALLOC(buffer1,(0)) 3714 call pawdijfr(gprimd,idir,ipert,my_natom,natom,nfftf,ngfftf,dtset%nspden,dtset%nsppol,& 3715 & dtset%ntypat,optfr,paw_ij1fr,pawang,pawfgrtab,pawrad,pawtab,& 3716 & cplex,dtset%qptn,rprimd,ucvol,buffer1,vtrial,vxc,xred,& 3717 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 3718 ABI_FREE(buffer1) 3719 end if 3720 3721 !PAW:allocate memory for non-symetrized occupancies matrix at EFermi (pawrhoijfermi) 3722 pawrhoijfermi_unsym => pawrhoijfermi 3723 if (psps%usepaw==1) then 3724 if (paral_atom) then 3725 ABI_MALLOC(pawrhoijfermi_unsym,(natom)) 3726 !Q phase should be 1 because Q=0 3727 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 3728 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 3729 call pawrhoij_alloc(pawrhoijfermi_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 3730 & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1) 3731 else 3732 call pawrhoij_init_unpacked(pawrhoijfermi_unsym) 3733 end if 3734 end if 3735 3736 !Initialize most of the Hamiltonian (arrays and quantities that do not depend on k + nl form factors) 3737 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,natom,& 3738 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,& 3739 & paw_ij=paw_ij,usecprj=usecprj,ph1d=ph1d,gpu_option=dtset%gpu_option,& 3740 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_spintab=mpi_enreg%my_isppoltab) 3741 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,paw_ij1=paw_ij1fr,& 3742 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,mpi_spintab=mpi_enreg%my_isppoltab) 3743 3744 3745 3746 !LOOP OVER SPINS 3747 do isppol=1,nsppol 3748 ikg=0;ikg1=0 3749 ! Continue to initialize the Hamiltonian at k+q 3750 call gs_hamkq%load_spin(isppol,with_nonlocal=.true.) 3751 3752 call rf_hamkq%load_spin(isppol,with_nonlocal=.true.) 3753 3754 3755 ! Nullify contribution to density at EFermi from this k-point 3756 if (nspden/=4) then 3757 rhoaug(:,:,:)=zero 3758 else 3759 rhoaug4(:,:,:,:)=zero 3760 end if 3761 call timab(125,1,tsec) 3762 3763 ! BIG FAT k POINT LOOP 3764 do ikpt=1,nkpt_rbz 3765 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 3766 istwf_k=istwfk_rbz(ikpt) 3767 npw_k=npwarr(ikpt) 3768 npw1_k=npwar1(ikpt) 3769 wtk_k=wtk_rbz(ikpt) 3770 3771 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 3772 eigen1(1+bd2tot_index : 2*nband_k**2+bd2tot_index) = zero 3773 bdtot_index=bdtot_index+nband_k 3774 bd2tot_index=bd2tot_index+2*nband_k**2 3775 ! Skip the rest of the k-point loop 3776 cycle 3777 end if 3778 3779 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 3780 ABI_MALLOC(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)) 3781 ABI_MALLOC(ylmgr1_k,(npw1_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 3782 3783 ! Continue to initialize the Hamiltonian at k+q 3784 kpoint(:)=kpt_rbz(:,ikpt) 3785 kpq(:)=kpoint(:)+dtset%qptn(1:3) 3786 3787 ABI_MALLOC(doccde_k,(nband_k)) 3788 ABI_MALLOC(doccde_kq,(nband_k)) 3789 ABI_MALLOC(eig0_k,(nband_k)) 3790 ABI_MALLOC(eig0_kq,(nband_k)) 3791 ABI_MALLOC(eig1_k,(2*nband_k**2)) 3792 ABI_MALLOC(fe1fixed_k,(nband_k)) 3793 ABI_MALLOC(fe1norm_k,(nband_k)) 3794 ABI_MALLOC(occ_k,(nband_k)) 3795 ABI_MALLOC(occ_kq,(nband_k)) 3796 ABI_MALLOC(rocceig,(nband_k,nband_k)) 3797 3798 eig1_k(:)=zero 3799 eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 3800 eig0_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index) 3801 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 3802 occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index) 3803 doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index) 3804 doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index) 3805 3806 ! For each pair of active bands (m,n), generates the ratios 3807 ! rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)) 3808 ! and decide to which band to attribute it. 3809 call occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,& 3810 & dtset%occopt,occ_k,occ_kq,rocceig) 3811 3812 ! Get plane-wave coeffs and related data at k 3813 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 3814 if (psps%useylm==1) then 3815 do ilm=1,psps%mpsang*psps%mpsang 3816 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 3817 end do 3818 end if 3819 3820 ! Get plane-wave coeffs and related data at k+q 3821 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 3822 if (psps%useylm==1) then 3823 do ilm=1,psps%mpsang*psps%mpsang 3824 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 3825 end do 3826 if (useylmgr1==1) then 3827 do ilm=1,psps%mpsang*psps%mpsang 3828 do ii=1,3 3829 ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm) 3830 end do 3831 end do 3832 end if 3833 end if 3834 3835 ! Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 3836 3837 ! Compute (k+G) vectors 3838 nkpg=0;if(ipert>=1.and.ipert<=natom) nkpg=3*dtset%nloalg(3) 3839 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 3840 if (nkpg>0) then 3841 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 3842 end if 3843 3844 ! Compute (k+q+G) vectors 3845 nkpg1=0;if(ipert>=1.and.ipert<=natom) nkpg1=3*dtset%nloalg(3) 3846 ABI_MALLOC(kpg1_k,(npw1_k,nkpg1)) 3847 if (nkpg1>0) then 3848 call mkkpg(kg1_k,kpg1_k,kpq,nkpg1,npw1_k) 3849 end if 3850 3851 ! ===== Preparation of non-local contributions 3852 3853 dimffnlk=0;if (ipert<=natom) dimffnlk=1 3854 ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,dtset%ntypat)) 3855 3856 ! Compute nonlocal form factors ffnlk at (k+G) 3857 if (ipert<=natom ) then 3858 ider=0;idir0=0 3859 call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,& 3860 & gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,& 3861 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,dtset%ntypat,& 3862 & psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum) 3863 end if 3864 3865 ! Compute nonlocal form factors ffnl1 at (k+q+G) 3866 !-- Atomic displacement perturbation 3867 if (ipert<=natom) then 3868 ider=0;idir0=0 3869 !-- Strain perturbation 3870 else if (ipert==natom+3.or.ipert==natom+4) then 3871 if (ipert==natom+3) istr=idir 3872 if (ipert==natom+4) istr=idir+3 3873 ider=1;idir0=-istr 3874 else if (ipert==natom+5) then !SPr deb rfmagn 3875 ider=0;idir0=0 3876 end if 3877 dimffnl1=1+ider;if (ider==1.and.idir0==0) dimffnl1=dimffnl1+2*psps%useylm 3878 ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,dtset%ntypat)) 3879 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,& 3880 & psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,& 3881 & npw1_k,dtset%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k) 3882 3883 ! ===== Preparation of kinetic contributions 3884 3885 ABI_MALLOC(dkinpw,(npw_k)) 3886 ABI_MALLOC(kinpw1,(npw1_k)) 3887 3888 ! Compute the derivative of the kinetic operator vs strain in dkinpw 3889 if (ipert==natom+3.or.ipert==natom+4) then 3890 if (ipert==natom+3) istr=idir 3891 if (ipert==natom+4) istr=idir+3 3892 call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr,& 3893 & kg_k,kpoint,npw_k) 3894 end if 3895 3896 ! Compute (1/2) (2 Pi)**2 (k+q+G)**2: 3897 ! call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k) 3898 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0) 3899 3900 ! ===== Load the k/k+q dependent parts of the Hamiltonian 3901 3902 ! Load k-dependent part in the Hamiltonian datastructure 3903 ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk)) 3904 call gs_hamkq%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,& 3905 & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.) 3906 if (size(ffnlk)>0) then 3907 call gs_hamkq%load_k(ffnl_k=ffnlk) 3908 else 3909 call gs_hamkq%load_k(ffnl_k=ffnl1) 3910 end if 3911 3912 ! Load k+q-dependent part in the Hamiltonian datastructure 3913 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 3914 call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 3915 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,& 3916 & compute_gbound=.true.) 3917 if (qne0) then 3918 ABI_MALLOC(ph3d1,(2,npw1_k,gs_hamkq%matblk)) 3919 call gs_hamkq%load_kprime(ph3d_kp=ph3d1,compute_ph3d=.true.) 3920 end if 3921 3922 ! Load k-dependent part in the 1st-order Hamiltonian datastructure 3923 call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw) 3924 3925 ! Compute fixed contributions to 1st-order Fermi energy 3926 ! and Fermi level charge density 3927 fe1fixed_k(:)=zero ; fe1norm_k(:)=zero 3928 3929 ! Note that dfpt_wfkfermi is called with kpoint, while kpt is used inside dfpt_wfkfermi 3930 if (nspden/=4) then 3931 call dfpt_wfkfermi(cg,cgq,cplex,cprj,cprjq,dtfil,eig0_k,eig1_k,fe1fixed_k,& 3932 & fe1norm_k,gs_hamkq,ibg,ibgq,icg,icgq,idir,ikpt,ipert,isppol,dtset%kptopt,mband_mem,& 3933 & mcgq,mcprjq,mkmem,mpi_enreg,mpw,nband_k,ncpgr,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k,& 3934 & pawrhoijfermi_unsym,prtvol,rf_hamkq,rhoaug,rocceig,wtk_k) 3935 else 3936 call dfpt_wfkfermi(cg,cgq,cplex,cprj,cprjq,dtfil,eig0_k,eig1_k,fe1fixed_k,& 3937 & fe1norm_k,gs_hamkq,ibg,ibgq,icg,icgq,idir,ikpt,ipert,isppol,dtset%kptopt,mband_mem,& 3938 & mcgq,mcprjq,mkmem,mpi_enreg,mpw,nband_k,ncpgr,npw_k,npw1_k,dtset%nspinor,nsppol,occ_k,& 3939 & pawrhoijfermi_unsym,prtvol,rf_hamkq,rhoaug4,rocceig,wtk_k) 3940 end if 3941 ! Free temporary storage 3942 ABI_FREE(kpg_k) 3943 ABI_FREE(kpg1_k) 3944 ABI_FREE(dkinpw) 3945 ABI_FREE(ffnlk) 3946 ABI_FREE(ffnl1) 3947 ABI_FREE(kinpw1) 3948 ABI_FREE(doccde_k) 3949 ABI_FREE(doccde_kq) 3950 ABI_FREE(eig0_k) 3951 ABI_FREE(eig0_kq) 3952 ABI_FREE(occ_kq) 3953 ABI_FREE(rocceig) 3954 ABI_FREE(ylm_k) 3955 ABI_FREE(ylm1_k) 3956 ABI_FREE(ylmgr1_k) 3957 ABI_FREE(ph3d) 3958 if (allocated(ph3d1)) then 3959 ABI_FREE(ph3d1) 3960 end if 3961 3962 ! Save eigenvalues (hartree) 3963 eigen1 (1+bd2tot_index : 2*nband_k**2+bd2tot_index) = eig1_k(:) 3964 3965 ! Accumulate sum over k points for 1st-order Fermi energy components 3966 do iband=1,nband_k 3967 fe1fixed=fe1fixed+wtk_k*occ_k(iband)*fe1fixed_k(iband) 3968 fe1norm=fe1norm+wtk_k*occ_k(iband)*fe1norm_k(iband) 3969 end do 3970 3971 ABI_FREE(eig1_k) 3972 ABI_FREE(occ_k) 3973 ABI_FREE(fe1fixed_k) 3974 ABI_FREE(fe1norm_k) 3975 3976 ! Keep track of total number of bands 3977 ! (all k points so far, even for k points not treated by me) 3978 bdtot_index=bdtot_index+nband_k 3979 bd2tot_index=bd2tot_index+2*nband_k**2 3980 3981 nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 3982 ! Shift array memory 3983 if (mkmem/=0) then 3984 ibg=ibg+nband_me 3985 icg=icg+npw_k*dtset%nspinor*nband_me 3986 ikg=ikg+npw_k 3987 end if 3988 if (mkqmem/=0) then 3989 ibgq=ibgq+dtset%nspinor*nband_me 3990 icgq=icgq+npw1_k*dtset%nspinor*nband_me 3991 end if 3992 if (mk1mem/=0) then 3993 ikg1=ikg1+npw1_k 3994 end if 3995 3996 ! End big k point loop 3997 end do 3998 3999 call timab(125,2,tsec) 4000 4001 ! Transfer density on augmented fft grid to normal fft grid in real space 4002 ! Also take into account the spin. 4003 if (nspden/=4) then 4004 if (psps%usepaw==0) then 4005 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhorfermi,rhoaug,1) 4006 else 4007 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhowfr ,rhoaug,1) 4008 end if 4009 else 4010 if (psps%usepaw==0) then 4011 do ispden=1,4 4012 call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhorfermi,rhoaug4(:,:,:,ispden),1) 4013 end do 4014 end if 4015 end if 4016 4017 end do ! End loop over spins 4018 4019 !More memory cleaning 4020 call gs_hamkq%free() 4021 call rf_hamkq%free() 4022 if(psps%usepaw==1) then 4023 call paw_ij_free(paw_ij1fr) 4024 ABI_FREE(paw_ij1fr) 4025 end if 4026 if (nspden/=4) then 4027 ABI_FREE(rhoaug) 4028 else 4029 ABI_FREE(rhoaug4) 4030 end if 4031 ABI_FREE(kg_k) 4032 ABI_FREE(kg1_k) 4033 4034 call timab(124,2,tsec) 4035 4036 4037 !=== MPI communications ================== 4038 if(xmpi_paral==1)then 4039 call timab(129,1,tsec) 4040 4041 ! Identify MPI buffer size 4042 buffer_size=cplex*dtset%nfft*nspden+2+mbd2kpsp 4043 ABI_MALLOC(buffer1,(buffer_size)) 4044 4045 ! Pack rhorfermi, fe1fixed, fe1norm 4046 indx=cplex*dtset%nfft*nspden 4047 if (psps%usepaw==0) then 4048 buffer1(1:indx)=reshape(rhorfermi,(/indx/)) 4049 else 4050 buffer1(1:indx)=reshape(rhowfr,(/indx/)) 4051 end if 4052 buffer1(indx+1)=fe1fixed ; buffer1(indx+2)=fe1norm 4053 indx=indx+2 ; bd2tot_index=0 4054 do isppol=1,nsppol 4055 do ikpt=1,nkpt_rbz 4056 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 4057 buffer1(indx+1:indx+2*nband_k**2)=eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2) 4058 bd2tot_index=bd2tot_index+2*nband_k**2 4059 indx=indx+2*nband_k**2 4060 end do 4061 end do 4062 if(indx<buffer_size)buffer1(indx+1:buffer_size)=zero 4063 4064 ! Build sum of everything 4065 call timab(48,1,tsec) 4066 call xmpi_sum(buffer1,buffer_size,spaceworld,ierr) 4067 call timab(48,2,tsec) 4068 4069 ! Unpack the final result 4070 indx=cplex*dtset%nfft*nspden 4071 if (psps%usepaw==0) then 4072 rhorfermi(:,:)=reshape(buffer1(1:indx),(/cplex*dtset%nfft,nspden/)) 4073 else 4074 rhowfr(:,:)=reshape(buffer1(1:indx),(/cplex*dtset%nfft,nspden/)) 4075 end if 4076 fe1fixed=buffer1(indx+1) ; fe1norm =buffer1(indx+2) 4077 indx=indx+2 ; bd2tot_index=0 4078 do isppol=1,nsppol 4079 do ikpt=1,nkpt_rbz 4080 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 4081 eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2)=buffer1(indx+1:indx+2*nband_k**2) 4082 bd2tot_index=bd2tot_index+2*nband_k**2 4083 indx=indx+2*nband_k**2 4084 end do 4085 end do 4086 ABI_FREE(buffer1) 4087 4088 ! Accumulate PAW occupancies 4089 if (psps%usepaw==1) then 4090 call pawrhoij_mpisum_unpacked(pawrhoijfermi_unsym,spaceworld) 4091 end if 4092 4093 call timab(129,2,tsec) 4094 end if ! if kpt parallel 4095 !=== End communications ================== 4096 4097 call timab(127,1,tsec) 4098 4099 !Normalize the fixed part of fermie1 4100 invfe1norm = zero ; if (abs(fe1norm) > tol10) invfe1norm=one/fe1norm 4101 fe1fixed=fe1fixed*invfe1norm 4102 4103 4104 if(nspden==4) then 4105 ! FR SPr symrhg will manage correctly this rearrangement 4106 rhorfermi(:,2)=rhorfermi(:,2)+(rhorfermi(:,1)+rhorfermi(:,4)) !(n+mx) 4107 rhorfermi(:,3)=rhorfermi(:,3)+(rhorfermi(:,1)+rhorfermi(:,4)) !(n+my) 4108 call timab(17,2,tsec) 4109 end if 4110 4111 !Symmetrize the density 4112 !In order to have the symrhg working in parallel on FFT coefficients, the size 4113 !of irzzon1 and phnons1 should be set to nfftot. Therefore, nsym\=1 does not work. 4114 !We also have the spin-up density, symmetrized, in rhorfermi(:,2). 4115 ABI_MALLOC(rhogfermi,(2,dtset%nfft)) 4116 if (psps%usepaw==0) then 4117 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,nspden,& 4118 & nsppol,nsym1,phnons1,rhogfermi,rhorfermi,rprimd,symaf1,symrl1,tnons1) 4119 else 4120 call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,nspden,& 4121 & nsppol,nsym1,phnons1,rhogfermi,rhowfr,rprimd,symaf1,symrl1,tnons1) 4122 end if 4123 4124 !PAW: Build new rhoij quantities then symetrize them 4125 !Compute and add the compensation density to rhowfr to get the total density 4126 if (psps%usepaw == 1) then 4127 if (size(nhatfermi)>0) then 4128 call pawmkrho(1,arg,cplex,gprimd,0,indsy1,0,mpi_enreg,& 4129 & my_natom,natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,pawfgr,& 4130 & pawfgrtab,-10001,pawrhoijfermi,pawrhoijfermi_unsym,pawtab,dtset%qptn,& 4131 & rhogfermi,rhowfr,rhorfermi,rprimd,symaf1,symrc1,dtset%typat,ucvol,& 4132 & dtset%usewvl,xred,pawang_sym=pawang1,pawnhat=nhatfermi) 4133 else 4134 call pawmkrho(1,arg,cplex,gprimd,0,indsy1,0,mpi_enreg,& 4135 & my_natom,natom,nspden,nsym1,dtset%ntypat,dtset%paral_kgb,pawang,pawfgr,& 4136 & pawfgrtab,-10001,pawrhoijfermi,pawrhoijfermi_unsym,pawtab,dtset%qptn,& 4137 & rhogfermi,rhowfr,rhorfermi,rprimd,symaf1,symrc1,dtset%typat,ucvol,& 4138 & dtset%usewvl,xred,pawang_sym=pawang1) 4139 end if 4140 ABI_FREE(rhowfr) 4141 call pawrhoij_free_unpacked(pawrhoijfermi_unsym) 4142 if (paral_atom) then 4143 call pawrhoij_free(pawrhoijfermi_unsym) 4144 ABI_FREE(pawrhoijfermi_unsym) 4145 end if 4146 end if 4147 ABI_FREE(rhogfermi) 4148 4149 !Normalize the Fermi level charge density (and associated PAW occupancies) 4150 rhorfermi(:,:)=invfe1norm*rhorfermi(:,:) 4151 if (psps%usepaw==1) then 4152 if (size(nhatfermi)>0) nhatfermi(:,:)=invfe1norm*nhatfermi(:,:) 4153 do iatom=1,my_natom 4154 do ispden=1,nspden 4155 do ilmn=1,pawrhoijfermi(iatom)%nrhoijsel 4156 pawrhoijfermi(iatom)%rhoijp(ilmn,ispden)=& 4157 & pawrhoijfermi(iatom)%rhoijp(ilmn,ispden)*invfe1norm 4158 end do 4159 end do 4160 end do 4161 end if 4162 4163 call timab(127,2,tsec) 4164 call timab(121,2,tsec) 4165 4166 DBG_EXIT('COLL') 4167 4168 end subroutine dfpt_rhofermi
ABINIT/dfpt_scfcv [ Functions ]
NAME
dfpt_scfcv
FUNCTION
Conducts set of passes or overall iterations of preconditioned conjugate gradient algorithm to converge wavefunctions to optimum and optionally to compute mixed derivatives of energy.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol)=pw coefficients of GS wavefunctions at k. cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q. cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k projected with non-local projectors: cprj=<p_i|Cnk> cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q projected with non-local projectors: cprjq=<p_i|Cnk+q> cpus= cpu time limit in seconds doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eew=2nd derivative of Ewald energy (hartree) efrhar=Contribution from frozen-wavefunction, hartree energy, to the second-derivative of total energy. efrkin=Contribution from frozen-wavefunction, kinetic energy, to the second-derivative of total energy. efrloc=Contribution from frozen-wavefunction, local potential, to the second-derivative of total energy. efrnl=Contribution from frozen-wavefunction, non-local potential, to the second-derivative of total energy. efrx1=Contribution from frozen-wavefunction, xc core correction(1), to the second-derivative of total energy. efrx2=Contribution from frozen-wavefunction, xc core correction(2), to the second-derivative of total energy. eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree) eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) eii=2nd derivative of pseudopotential core energy (hartree) evdw=DFT-D semi-empirical part of 2nd-order total energy fermie=fermi energy (Hartree) hdr <type(hdr_type)>=the header of wf, den and pot files idir=direction of the current perturbation indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points indsy1(4,nsym1,natom)=indirect indexing array for atom labels ipert=type of the perturbation irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for RF symmetries istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates at k kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points. kxc(nfftf,nkxc)=exchange and correlation kernel (see rhotoxc.f) mband_mem_rbz=maximum number of bands per processor in memory for cg mgfftf=maximum size of 1D FFTs for the "fine" grid (see NOTES in respfn.F90) mkmem =number of k points treated by this node (GS data) mkqmem =number of k+q points which can fit in memory (GS data); 0 if use disk mk1mem =number of k points which can fit in memory (RF data); 0 if use disk mpert=maximum number of ipert mpw=maximum dimensioned size of npw for wfs at k. mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs). nattyp(ntypat)= # atoms of each type. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point, for each polarization ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>) nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) ngfftf(1:18)=integer array with FFT box dimensions and other for the "fine" grid (see NOTES in respfn.F90) nkpt=number of k points in the full BZ nkpt_rbz=number of k points in the reduced BZ for this perturbation nkxc=second dimension of the kxc array. mpi_enreg=information about MPI parallelization my_natom=number of atoms treated by current processor npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point nspden=number of spin-density components nspinor=number of spinorial components of the wavefunctions nsym1=number of symmetry elements in space group consistent with perturbation n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used otherwise, nfftf occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2) at each k+q point of the reduced Brillouin zone. occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2) at each k point of the reduced Brillouin zone. paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawang <type(pawang_type)>=paw angular mesh and related data pawang1 <type(pawang_type)>=pawang datastr. containing only symmetries preserving the perturbation pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data pertcase=fuill index of the perturbation phnons1(2,nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic transl. phases, for RF symmetries ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information ph1df(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information for the "fine" grid prtbbb=if 1, band-by-band decomposition (also dim of d2bbb) psps <type(pseudopotential_type)>=variables related to pseudopotentials qphon(3)=reduced coordinates for the phonon wavelength rhog(2,nfftf)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. rprimd(3,3)=dimensional primitive translations in real space (bohr) symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations symrc1(3,3,nsym1)=symmetry operations in reciprocal space symrl1(3,3,nsym1)=symmetry operations in real space in terms of primitive translations tnons1(3,nsym1)=non-symmorphic translations usecprj= 1 if cprj, cprjq arrays are stored in memory useylmgr = 1 if ylmgr array is allocated useylmgr1= 1 if ylmgr1 array is allocated ddk<wfk_t>=ddk file vpsp1(cplex*nfftf)=first-order derivative of the ionic potential vtrial(nfftf,nspden)=GS potential (Hartree). vxc(nfftf,nspden)=Exchange-Correlation GS potential (Hartree) vxctau(nfftf,nspden,4*usekden)=derivative of e_xc with respect to kinetic energy density, for mGGA wtk_rbz(nkpt_rbz)=weight for each k point in the reduced Brillouin zone xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm*useylmgr)= gradients of real spherical harmonics at k ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm*useylmgr1)= gradients of real spherical harmonics at k+q
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF
wavefunctions at k,q. They are orthogonalized to the active.
d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
second order derivatives
d2lo(2,mpert,3,mpert)=local contributions to the 2DTEs
d2nl(2,mpert,3,mpert)=non-local contributions to the 2DTEs
d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs
eberry=energy associated with Berry phase
edocc=correction to 2nd-order total energy coming from changes of occupation
eeig0=0th-order eigenenergies part of 2nd-order total energy
ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
for strain perturbation only (zero otherwise, and not used)
ehart1=1st-order Hartree part of 2nd-order total energy
eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues (hartree)
ek0=0th-order kinetic energy part of 2nd-order total energy.
ek1=1st-order kinetic energy part of 2nd-order total energy.
eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
end0=0th-order nuclear dipole part of 2nd-order total energy
end1=1st-order nuclear dipole part of 2nd-order total energy
enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
eovl1=1st-order change of wave-functions overlap, part of 2nd-order energy
PAW only - Eq(79) and Eq(80) of PRB 78, 035105 (2008) [[cite:Audouze2008]]
epaw1=1st-order PAW on-site part of 2nd-order total energy.
etotal=total energy (sum of 7 contributions) (hartree)
evxctau0=0th-order energy from vxctau
evxctau1=1st-order energy from vxctau
exc1=1st-order exchange-correlation part of 2nd-order total energy.
gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}|\Psi^{(1)}>
The wavefunction is orthogonal to the active space (for metals). It is not
coherent with cg1.
resid(mband*nkpt_rbz*nsppol)=residuals for each band over all k points
of the reduced Brillouin zone, and spins
residm=maximum value from resid array (except for nbdbuf highest bands)
conv_retcode=return code, 0 if convergence was achieved.
SIDE EFFECTS
cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=updated wavefunctions (ortho. to occ. states);
initialized= if 0 the initialization of the RF run is not yet finished
mpi_enreg=information about MPI parallelization
rhog1(2,nfftf)=array for Fourier transform of RF electron density
rhor1(cplex*nfftf,nspden)=array for RF electron density in electrons/bohr**3.
=== if psps%usepaw==1
pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
SOURCE
272 subroutine dfpt_scfcv(atindx,blkflg,cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cpus,& 273 & dielt,dim_eig2rf,doccde_rbz,docckqde,dtfil,dtset,& 274 & d2bbb,d2lo,d2nl,d2ovl,eberry,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,& 275 & ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek0,ek1,eloc0,elpsp1,& 276 & end0,end1,enl0,enl1,eovl1,epaw1,etotal,evxctau0,evxctau1,evdw,exc1,& 277 & fermie,gh0c1_set,gh1c_set,hdr,idir,indkpt1,& 278 & indsy1,initialized,ipert,irrzon1,istwfk_rbz,& 279 & kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mgfftf,mkmem,mkqmem,mk1mem,& 280 & mpert,mpi_enreg,mpw,mpw1,mpw1_mq,my_natom,nattyp,nband_rbz,ncpgr,& 281 & nfftf,ngfftf,nhat,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,& 282 & nsym1,n3xccc,occkq,occ_rbz,& 283 & paw_an,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawrhoij1,pawtab,& 284 & pertcase,phnons1,ph1d,ph1df,& 285 & prtbbb,psps,qphon,resid,residm,rhog,rhog1,& 286 & rhor,rhor1,rprimd,symaf1,symrc1,symrl1,tnons1,& 287 & usecprj,useylmgr,useylmgr1,ddk_f,vpsp1,vtrial,vxc,vxctau,& 288 & wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,zeff,conv_retcode,& 289 & kramers_deg,& 290 & cg_mq,cg1_mq,cg1_active_mq,docckde_mq,eigen_mq,eigen1_mq,gh0c1_set_mq,gh1c_set_mq,& 291 & kg1_mq,npwar1_mq,occk_mq,resid_mq,residm_mq,rhog1_pq,rhog1_mq,rhor1_pq,rhor1_mq) 292 293 !Arguments ------------------------------------ 294 type(dataset_type),intent(in) :: dtset 295 type(pseudopotential_type),intent(in) :: psps 296 integer,intent(in) :: cplex,dim_eig2rf,idir,ipert,mgfftf,mk1mem,mkmem,mkqmem 297 integer,intent(in) :: mpert,mpw,mpw1,my_natom,n3xccc,ncpgr,nfftf 298 integer,intent(in) :: mband_mem_rbz 299 integer,intent(in) :: mpw1_mq !-q duplicate 300 integer,intent(in) :: nkpt,nkpt_rbz,nkxc,nspden 301 integer,intent(in) :: nsym1,pertcase,prtbbb,usecprj,useylmgr,useylmgr1 302 logical,intent(in) :: kramers_deg 303 integer,intent(inout) :: initialized 304 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise 305 integer,intent(in) :: atindx(dtset%natom) 306 integer,intent(out) :: blkflg(3,mpert,3,mpert) 307 integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom) 308 integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 309 integer,intent(in) :: istwfk_rbz(nkpt_rbz) 310 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem),nattyp(psps%ntypat) 311 integer,intent(in) :: nband_rbz(nkpt_rbz*dtset%nsppol) 312 integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz) 313 integer,optional,intent(in) :: npwar1_mq(nkpt_rbz) !-q duplicate 314 integer,optional,intent(in) :: kg1_mq(3,mpw1_mq*mk1mem)! 315 integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1) 316 integer,intent(out) :: conv_retcode 317 real(dp),intent(in) :: cpus,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,eii 318 real(dp),intent(out) :: eberry,edocc,eeig0,ehart01,ehart1,ek0,ek1,eloc0,elpsp1,end0,end1 319 real(dp),intent(out) :: enl0,enl1,eovl1,epaw1,etotal,evdw,evxctau0,evxctau1,exc1,residm 320 real(dp),optional,intent(out) :: residm_mq !-q duplicate 321 real(dp),intent(inout) :: fermie 322 real(dp),intent(in) :: qphon(3) 323 ! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise 324 integer,intent(in) :: ngfftf(18) 325 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband_mem_rbz*mkmem*dtset%nsppol) 326 real(dp),intent(inout) :: cg1(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol) 327 real(dp),intent(out) :: cg1_active(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol*dim_eig2rf) 328 real(dp),intent(out) :: gh1c_set(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol*dim_eig2rf) 329 real(dp),intent(out) :: gh0c1_set(2,mpw1*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol*dim_eig2rf) 330 real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*mband_mem_rbz*mkqmem*dtset%nsppol) 331 real(dp),optional,intent(inout) :: cg1_mq(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol) !start -q duplicates 332 real(dp),optional,intent(out) :: cg1_active_mq(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol*dim_eig2rf)! 333 real(dp),optional,intent(out) :: gh1c_set_mq(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol*dim_eig2rf) ! 334 real(dp),optional,intent(out) :: gh0c1_set_mq(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol*dim_eig2rf) ! 335 real(dp),optional,intent(in) :: cg_mq(2,mpw1_mq*dtset%nspinor*mband_mem_rbz*mkqmem*dtset%nsppol) ! 336 real(dp),optional,intent(in) :: eigen_mq(dtset%mband*nkpt_rbz*dtset%nsppol) ! 337 real(dp),optional,intent(in) :: docckde_mq(dtset%mband*nkpt_rbz*dtset%nsppol) ! 338 real(dp),optional,intent(out) :: eigen1_mq(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol) ! 339 real(dp),optional,intent(in) :: occk_mq(dtset%mband*nkpt_rbz*dtset%nsppol) ! 340 real(dp),optional,intent(out) :: resid_mq(dtset%mband*nkpt_rbz*nspden) !end 341 real(dp),intent(out) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb) 342 real(dp),intent(out) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) 343 real(dp),intent(out) :: d2ovl(2,3,mpert,3,mpert*psps%usepaw) 344 real(dp),intent(in) :: dielt(3,3) 345 real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*dtset%nsppol) 346 real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*dtset%nsppol) 347 real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*dtset%nsppol) 348 real(dp),intent(out) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol) 349 real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*dtset%nsppol) 350 real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),kxc(nfftf,nkxc) 351 real(dp),intent(in) :: nhat(nfftf,dtset%nspden) 352 real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*dtset%nsppol) 353 real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol) 354 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom),ph1df(2,3*(2*mgfftf+1)*dtset%natom) 355 real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 356 real(dp),intent(out) :: resid(dtset%mband*nkpt_rbz*nspden) 357 real(dp),intent(in) :: rhog(2,nfftf),rhor(nfftf,nspden),rprimd(3,3) 358 real(dp),intent(inout) :: rhog1(2,nfftf),rhor1(cplex*nfftf,nspden),xred(3,dtset%natom) 359 real(dp),optional,intent(inout) :: rhog1_pq(2,nfftf),rhor1_pq(cplex*nfftf,nspden) !+q/-q duplicates 360 real(dp),optional,intent(inout) :: rhog1_mq(2,nfftf),rhor1_mq(cplex*nfftf,nspden) ! 361 real(dp),intent(in) :: tnons1(3,nsym1) 362 real(dp),target,intent(in) :: vtrial(nfftf,nspden) 363 real(dp),intent(in) :: vpsp1(cplex*nfftf),vxc(nfftf,nspden) 364 real(dp),intent(inout) :: vxctau(nfftf,dtset%nspden,4*dtset%usekden) 365 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xccc3d1(cplex*n3xccc) 366 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 367 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 368 real(dp),intent(in) :: ylmgr(mpw*mkmem,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 369 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3+6*((ipert-dtset%natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 370 real(dp),intent(in) :: zeff(3,3,dtset%natom) 371 type(pawcprj_type),intent(in) :: cprj(dtset%natom,dtset%nspinor*mband_mem_rbz*mkmem*dtset%nsppol*usecprj) 372 type(pawcprj_type),intent(in) :: cprjq(dtset%natom,dtset%nspinor*mband_mem_rbz*mkqmem*dtset%nsppol*usecprj) 373 type(datafiles_type),intent(in) :: dtfil 374 type(hdr_type),intent(inout) :: hdr 375 type(pawang_type),intent(in) :: pawang,pawang1 376 type(pawfgr_type),intent(in) :: pawfgr 377 type(paw_an_type),intent(in) :: paw_an(my_natom*psps%usepaw) 378 type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw) 379 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw) 380 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 381 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*psps%usepaw) 382 type(pawrhoij_type),intent(inout) :: pawrhoij1(my_natom*psps%usepaw) 383 type(pawtab_type), intent(in) :: pawtab(psps%ntypat*psps%usepaw) 384 type(MPI_type),intent(inout) :: mpi_enreg 385 type(wfk_t),intent(inout) :: ddk_f(4) 386 387 !Local variables------------------------------- 388 !scalars 389 integer,parameter :: level=12,response=1 390 integer :: afford,bantot_rbz,choice,cplex_rhoij,dbl_nnsclo 391 integer :: has_dijfr,has_diju,iatom,ider,idir_dum,idir_paw1,ierr,errid,denpot 392 integer :: iprcel,iscf10_mod,iscf_mod,ispden,ispmix 393 integer :: istep,istep_fock_outer,istep_mix,itypat,izero,me,mgfftdiel,mvdum !lmn2_size, 394 integer :: nfftdiel,nfftmix,nfftotf,nhat1grdim,npawmix,npwdiel,nspden_rhoij,nstep,nzlmopt 395 integer :: optene,optfr,option,optres,prtfor,qphase_rhoij,quit,quit_sum,qzero 396 integer :: my_quit,quitsum_request,timelimit_exit,varid,ncerr,ncid 397 integer ABI_ASYNC :: quitsum_async 398 integer :: rdwrpaw,spaceComm,sz1,sz2,usexcnhat,with_vectornd,Z_kappa 399 integer :: dbl_nnsclo_mq,ifft !-q duplicate for dbl_nnsclo 400 !integer :: pqmq ! pqmq = indicator for potential mixing 401 logical :: need_fermie1,nmxc,paral_atom,use_nhat_gga 402 real(dp) :: wtime_step,now,prev 403 real(dp) :: born,born_bar,boxcut,deltae,diffor,diel_q,dum,ecut,ecutf,elast 404 real(dp) :: epawdc1_dum,evar,fe1fixed,fermie1,gsqcut,qphon_norm,maxfor,renorm,res2,res3,residm2 405 real(dp) :: ucvol,vxcavg,elmag1 406 real(dp) :: res2_mq,fe1fixed_mq,elast_mq 407 real(dp) :: eberry_mq,edocc_mq,eeig0_mq,ehart01_mq,ehart1_mq,ek0_mq,ek1_mq,eloc0_mq,elpsp1_mq 408 real(dp) :: end0_mq,end1_mq,enl0_mq,enl1_mq,eovl1_mq,epaw1_mq,exc1_mq,fermie1_mq,deltae_mq,elmag1_mq 409 real(dp) :: evxctau0_mq,evxctau1_mq 410 character(len=500) :: msg 411 character(len=500),parameter :: MY_NAME="dfpt_scfcv" 412 character(len=fnlen) :: fi1o 413 !character(len=fnlen) :: fi1o_vtk 414 integer :: prtopt 415 type(ab7_mixing_object) :: mix 416 type(efield_type) :: dtefield 417 !arrays 418 integer :: ngfftmix(18) 419 integer,allocatable :: dimcprj(:),pwindall(:,:,:) 420 integer,pointer :: my_atmtab(:) 421 real(dp) :: dielar(7) 422 real(dp) :: favg(3),gmet(3,3),gprimd(3,3),q_cart(3),qphon2(3),qred2cart(3,3) 423 real(dp) :: rhomag(2,nspden),rmet(3,3),tollist(12),tsec(2) 424 real(dp) :: zeff_red(3),zeff_bar(3,3) 425 real(dp) :: intgden(dtset%nspden,dtset%natom),dentot(dtset%nspden) 426 !real(dp) :: zdmc_red(3),zdmc_bar(3,3),mean_rhor1(1) !dynamic magnetic charges and mean density 427 real(dp),allocatable :: dielinv(:,:,:,:,:) 428 real(dp),allocatable :: fcart(:,:),nhat1(:,:),nhat1gr(:,:,:),nhatfermi(:,:),nvresid1(:,:),nvresid2(:,:) 429 real(dp),allocatable :: qmat(:,:,:,:,:,:),resid2(:),rhog2(:,:),rhor2(:,:),rhorfermi(:,:) 430 real(dp),allocatable :: susmat(:,:,:,:,:),vectornd(:,:,:),vhartr1(:),vxc1(:,:) 431 real(dp),allocatable :: vhartr1_tmp(:,:) 432 real(dp),allocatable,target :: vtrial1(:,:),vtrial2(:,:) 433 real(dp),allocatable :: vtrial1_pq(:,:),vtrial1_mq(:,:),rhorfermi_mq(:,:) 434 real(dp),allocatable :: nvresid1_mq(:,:),vxc1_mq(:,:),vhartr1_mq(:) 435 real(dp),pointer :: vtrial1_tmp(:,:) 436 type(pawcprj_type),allocatable :: cprj1(:,:) 437 type(paw_an_type),allocatable :: paw_an1(:) 438 type(paw_ij_type),allocatable :: paw_ij1(:) 439 type(pawrhoij_type),allocatable :: pawrhoijfermi(:) 440 441 ! ********************************************************************* 442 443 DBG_ENTER("COLL") 444 445 if (dtset%occopt == 9) then 446 write(msg,'(a)') "Cannot perform dfpt with occopt = 9: not yet implemented" 447 ABI_ERROR(msg) 448 end if 449 450 call timab(120,1,tsec) 451 call timab(154,1,tsec) 452 453 ! intel 18 really needs this to be initialized 454 maxfor = zero 455 456 ! enable time limit handler if not done in callers. 457 if (enable_timelimit_in(MY_NAME) == MY_NAME) then 458 write(std_out,*)"Enabling timelimit check in function: ",trim(MY_NAME)," with timelimit: ",trim(sec2str(get_timelimit())) 459 end if 460 461 !Parallelism data 462 spaceComm=mpi_enreg%comm_cell 463 me=mpi_enreg%me_kpt 464 paral_atom=(my_natom/=dtset%natom) 465 my_atmtab=>mpi_enreg%my_atmtab 466 467 !Save some variables from dataset definition 468 ecut=dtset%ecut 469 ecutf=ecut;if (psps%usepaw==1) ecutf=dtset%pawecutdg 470 iprcel=dtset%iprcel 471 tollist(1)=dtset%tolmxf;tollist(2)=dtset%tolwfr 472 tollist(3)=dtset%toldff;tollist(4)=dtset%toldfe 473 tollist(6)=dtset%tolvrs;tollist(7)=dtset%tolrff 474 nfftotf=product(ngfftf(1:3)) 475 nstep=dtset%nstep 476 iscf_mod=dtset%iscf 477 iscf10_mod=mod(iscf_mod,10) 478 479 qzero=0; if(qphon(1)**2+qphon(2)**2+qphon(3)**2 < tol14) qzero=1 480 481 need_fermie1=((qzero==1.and.dtset%frzfermi==0.and.nstep>0).and.& 482 & (dtset%occopt>=3.and.dtset%occopt<=8).and. & 483 & (ipert<=dtset%natom.or.ipert==dtset%natom+3.or.ipert==dtset%natom+4.or.ipert==dtset%natom+5)) 484 485 !The value of iscf must be modified if ddk perturbation, see dfpt_looppert.f 486 if (ipert==dtset%natom+1.or.ipert==dtset%natom+10.or.ipert==dtset%natom+11) iscf_mod=-3 487 488 !Compute different geometric tensor, as well as ucvol, from rprimd 489 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 490 491 !Compute large sphere cut-off gsqcut 492 qphon2(:)=zero;if (psps%usepaw==1) qphon2(:)=qphon(:) 493 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,qphon2,ngfftf) 494 495 !Some variables need to be initialized/nullify at start 496 quit=0 ; dbl_nnsclo=0 ; elast=zero; conv_retcode = -1 497 optres=merge(0,1,abs(iscf_mod)<10) 498 nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 499 usexcnhat=0 500 !This might be taken away later 501 edocc=zero ; eeig0=zero ; ehart01=zero ; ehart1=zero ; ek0=zero ; ek1=zero 502 eloc0=zero ; elpsp1=zero ; end0=zero; end1=zero; 503 enl0=zero ; enl1=zero ; eovl1=zero; evxctau0=zero; evxctau1=zero; exc1=zero 504 deltae=zero ; fermie1=zero ; epaw1=zero ; eberry=zero ; elmag1=zero 505 elast_mq=zero 506 dbl_nnsclo_mq=0 507 !This might be taken away later 508 edocc_mq=zero ; eeig0_mq=zero ; ehart01_mq=zero ; ehart1_mq=zero ; ek0_mq=zero ; ek1_mq=zero 509 eloc0_mq=zero ; elpsp1_mq=zero ; enl0_mq=zero ; enl1_mq=zero ; 510 end0_mq=zero; end1_mq=zero; eovl1_mq=zero; evxctau0_mq=zero; evxctau1_mq=zero; exc1_mq=zero 511 deltae_mq=zero ; fermie1_mq=zero ; epaw1_mq=zero ; eberry_mq=zero ; elmag1_mq=zero 512 res2_mq=zero 513 514 !Examine tolerance criteria, and eventually print a line to the output 515 !file (with choice=1, the only non-dummy arguments of scprqt are 516 !nstep, tollist and iscf - still, diffor,res2,prtfor,fcart are here initialized to 0) 517 choice=1 ; prtfor=0 ; diffor=zero ; res2=zero 518 ABI_MALLOC(fcart,(3,dtset%natom)) 519 520 !At present, no double loop 521 istep_mix=1 ; istep_fock_outer=1 522 523 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 524 & etotal,favg,fcart,fermie,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 525 & 1,iscf_mod,istep,istep_fock_outer,istep_mix,kpt_rbz,maxfor,& 526 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 527 & nstep,occ_rbz,0,prtfor,0,& 528 & quit,res2,resid,residm,response,& 529 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 530 531 !Allocations/initializations for PAW only 532 if(psps%usepaw==1) then 533 usexcnhat=maxval(pawtab(:)%usexcnhat) 534 use_nhat_gga=(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) 535 ! 1st-order compensation density 536 ABI_MALLOC(nhat1,(cplex*nfftf,dtset%nspden)) 537 nhat1=zero 538 ! Projections of 1-st order WF on nl projectors 539 ABI_MALLOC(cprj1,(dtset%natom,dtset%nspinor*mband_mem_rbz*mk1mem*dtset%nsppol*usecprj)) 540 if (usecprj==1.and.mk1mem/=0) then 541 !cprj ordered by atom-type 542 ABI_MALLOC(dimcprj,(dtset%natom)) 543 call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 544 call pawcprj_alloc(cprj1,0,dimcprj) 545 ABI_FREE(dimcprj) 546 end if 547 ! 1st-order arrays/variables related to the PAW spheres 548 ABI_MALLOC(paw_an1,(my_natom)) 549 ABI_MALLOC(paw_ij1,(my_natom)) 550 call paw_an_nullify(paw_an1) 551 call paw_ij_nullify(paw_ij1) 552 553 has_dijfr=0;if (ipert/=dtset%natom+1.and.ipert/=dtset%natom+10) has_dijfr=1 554 has_diju=merge(0,1,dtset%usepawu==0) 555 call paw_an_init(paw_an1,dtset%natom,dtset%ntypat,0,0,dtset%nspden,& 556 & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,has_vxctau=dtset%usekden,& 557 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 558 call paw_ij_init(paw_ij1,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,0,dtset%natom,& 559 & dtset%ntypat,dtset%typat,pawtab,& 560 & has_dij=1,has_dijhartree=1,has_dijfr=has_dijfr,has_dijU=has_diju,& 561 & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom) 562 else 563 ABI_MALLOC(nhat1,(0,0)) 564 ABI_MALLOC(cprj1,(0,0)) 565 ABI_MALLOC(paw_an1,(0)) 566 ABI_MALLOC(paw_ij1,(0)) 567 end if ! PAW 568 569 !Various allocations (potentials) 570 ABI_MALLOC(vhartr1,(cplex*nfftf)) 571 ABI_MALLOC(vtrial1,(cplex*nfftf,nspden)) 572 if(.not.kramers_deg) then 573 ABI_MALLOC(vhartr1_mq,(cplex*nfftf)) 574 ABI_MALLOC(vtrial1_pq,(cplex*nfftf,nspden)) 575 ABI_MALLOC(vtrial1_mq,(cplex*nfftf,nspden)) 576 end if 577 ! TODO: for non collinear case this should always be nspden, in NCPP case as well!!! 578 ABI_MALLOC(vxc1,(cplex*nfftf,nspden*(1-usexcnhat))) ! Not always needed 579 vtrial1_tmp => vtrial1 ! this is to avoid errors when vtrial1_tmp is unused 580 581 if (.not.kramers_deg) then 582 ABI_MALLOC(vxc1_mq,(cplex*nfftf,nspden*(1-usexcnhat))) 583 end if 584 585 !Several parameters and arrays for the SCF mixing: 586 !These arrays are needed only in the self-consistent case 587 if (iscf_mod>0.or.iscf_mod==-3) then 588 ABI_MALLOC(nvresid1,(cplex*nfftf,dtset%nspden)) 589 if (nstep==0) nvresid1=zero 590 if ((dtset%getddb .ne. 0 .or. dtset%irdddb .ne.0) .and. qzero .ne. 1) then 591 ABI_MALLOC(nvresid2,(cplex*nfftf,dtset%nspden)) 592 if (nstep==0) nvresid2=zero 593 end if 594 if (.not.kramers_deg) then 595 ABI_MALLOC(nvresid1_mq,(cplex*nfftf,dtset%nspden)) 596 if (nstep==0) nvresid1_mq=zero 597 end if 598 else 599 ABI_MALLOC(nvresid1,(0,0)) 600 if(.not.kramers_deg) then 601 ABI_MALLOC(nvresid1_mq,(0,0)) 602 end if 603 end if 604 if(nstep>0 .and. iscf_mod>0) then 605 dielar(1)=dtset%diecut;dielar(2)=dtset%dielng 606 dielar(3)=dtset%diemac;dielar(4)=dtset%diemix 607 dielar(5)=dtset%diegap;dielar(6)=dtset%dielam 608 dielar(7)=dtset%diemix;if (dtset%iscf>=10) dielar(7)=dtset%diemixmag 609 ! Additional allocation for mixing within PAW 610 npawmix=0 611 if(psps%usepaw==1) then 612 do iatom=1,my_natom 613 itypat=pawrhoij1(iatom)%itypat 614 pawrhoij1(iatom)%use_rhoijres=1 615 sz1=pawrhoij1(iatom)%cplex_rhoij*pawrhoij1(iatom)%qphase*pawrhoij1(iatom)%lmn2_size 616 sz2=pawrhoij1(iatom)%nspden 617 ABI_MALLOC(pawrhoij1(iatom)%rhoijres,(sz1,sz2)) 618 do ispden=1,pawrhoij1(iatom)%nspden 619 pawrhoij1(iatom)%rhoijres(:,ispden)=zero 620 end do 621 ABI_MALLOC(pawrhoij1(iatom)%kpawmix,(pawtab(itypat)%lmnmix_sz)) 622 pawrhoij1(iatom)%lmnmix_sz=pawtab(itypat)%lmnmix_sz 623 pawrhoij1(iatom)%kpawmix=pawtab(itypat)%kmix 624 npawmix=npawmix+pawrhoij1(iatom)%nspden*pawtab(itypat)%lmnmix_sz & 625 & *pawrhoij1(iatom)%cplex_rhoij*pawrhoij1(iatom)%qphase 626 end do 627 end if 628 denpot = AB7_MIXING_POTENTIAL 629 if (dtset%iscf > 10) denpot = AB7_MIXING_DENSITY 630 if (psps%usepaw==1.and.dtset%pawmixdg==0) then 631 ispmix=AB7_MIXING_FOURRIER_SPACE;nfftmix=dtset%nfft;ngfftmix(:)=dtset%ngfft(:) 632 else 633 ispmix=AB7_MIXING_REAL_SPACE;nfftmix=nfftf;ngfftmix(:)=ngfftf(:) 634 end if 635 if (iscf10_mod == 5 .or. iscf10_mod == 6) then 636 call ab7_mixing_new(mix, iscf10_mod, denpot, cplex, & 637 & nfftf, dtset%nspden, npawmix, errid, msg, dtset%npulayit) 638 else 639 call ab7_mixing_new(mix, iscf10_mod, denpot, max(cplex, ispmix), & 640 & nfftmix, dtset%nspden, npawmix, errid, msg, dtset%npulayit) 641 end if 642 if (errid /= AB7_NO_ERROR) then 643 ABI_ERROR(msg) 644 end if 645 if (dtset%mffmem == 0) then 646 call ab7_mixing_use_disk_cache(mix, dtfil%fnametmp_fft) 647 end if 648 end if ! iscf, nstep 649 650 !Here, allocate arrays for computation of susceptibility and dielectric matrix or for TDDFT 651 if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then 652 ! Here, for TDDFT, artificially set iprcel . Also set a variable to reduce the memory needs. 653 afford=1 654 if(iscf_mod==-1) then 655 iprcel=21 656 afford=0 657 end if 658 npwdiel=1 659 mgfftdiel=1 660 nfftdiel=1 661 ! Now, performs allocation 662 ! CAUTION : the dimensions are still those of GS, except for phnonsdiel 663 ABI_MALLOC(dielinv,(2,npwdiel*afford,nspden,npwdiel,nspden)) 664 ABI_MALLOC(susmat,(2,npwdiel*afford,nspden,npwdiel,nspden)) 665 end if 666 667 !Initialize Berry-phase related stuffs 668 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 669 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 670 ABI_MALLOC(pwindall,(max(mpw,mpw1)*mkmem,8,3)) 671 call dfptff_initberry(dtefield,dtset,gmet,kg,kg1,dtset%mband,mkmem,mpi_enreg,& 672 & mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,occ_rbz,pwindall,rprimd) 673 ! calculate inverse of the overlap matrix 674 ABI_MALLOC(qmat,(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)) 675 call qmatrix(cg,dtefield,qmat,mpi_enreg,mpw,mpw1,mkmem,dtset%mband,mband_mem_rbz,& 676 & npwarr,nkpt,dtset%nspinor,dtset%nsppol,pwindall) 677 else 678 ABI_MALLOC(pwindall,(0,0,0)) 679 ABI_MALLOC(qmat,(0,0,0,0,0,0)) 680 end if 681 682 ! if any nuclear dipoles are nonzero, compute the vector potential in real space 683 with_vectornd = 0 684 ! nuclear dipoles only work with the DDK response function 685 if ( (ANY(ABS(dtset%nucdipmom(:,:))>tol8)) .AND. (ipert.EQ.dtset%natom+1) ) with_vectornd = 1 686 if(allocated(vectornd)) then 687 ABI_FREE(vectornd) 688 end if 689 ABI_MALLOC(vectornd,(with_vectornd*nfftf,dtset%nspden,3)) 690 if(with_vectornd .EQ. 1) then 691 call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,& 692 & ngfftf,dtset%nspden,dtset%nucdipmom,rprimd,vectornd,xred) 693 endif 694 695 call timab(154,2,tsec) 696 697 !###################################################################### 698 !PERFORM ELECTRONIC ITERATIONS 699 !###################################################################### 700 701 !Offer option of computing 2nd-order total energy with existing 702 !wavefunctions when nstep<=0, else do nstep iterations 703 !Note that for non-self-consistent calculations, this loop will be exited 704 !after the first call to dfpt_vtorho 705 706 !Pass through the first routines even when nstep==0 707 !write(std_out,*) 'dfpt_scfcv, nstep=', max(1,nstep) 708 709 quitsum_request = xmpi_request_null; timelimit_exit = 0 710 711 do istep=1,max(1,nstep) 712 713 ! Handle time limit condition. 714 if (istep == 1) prev = abi_wtime() 715 if (istep > 1) then 716 now = abi_wtime() 717 wtime_step = now - prev 718 prev = now 719 call wrtout(std_out,sjoin(" dfpt_scfcv: previous iteration took ",sec2str(wtime_step))) 720 721 if (have_timelimit_in(MY_NAME)) then 722 if (istep > 2) then 723 call xmpi_wait(quitsum_request,ierr) 724 if (quitsum_async > 0) then 725 write(msg,"(3a)")" Approaching time limit ",trim(sec2str(get_timelimit())),". Will exit istep loop in dfpt_scfcv." 726 ABI_COMMENT(msg) 727 call wrtout(ab_out, msg, "COLL") 728 timelimit_exit = 1 729 exit 730 end if 731 end if 732 733 my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1 734 call xmpi_isum(my_quit,quitsum_async,spacecomm,quitsum_request,ierr) 735 end if 736 end if 737 738 ! ###################################################################### 739 ! The following steps are done once 740 ! ---------------------------------------------------------------------- 741 if (istep==1)then 742 743 ! PAW only: compute frozen part of 1st-order compensation density 744 ! and frozen part of psp strengths Dij 745 ! ---------------------------------------------------------------------- 746 if (psps%usepaw==1) then 747 optfr=0 748 idir_paw1 = idir 749 if (ipert==dtset%natom+11) then 750 call rf2_getidirs(idir,idir_dum,idir_paw1) 751 end if 752 call pawdijfr(gprimd,idir_paw1,ipert,my_natom,dtset%natom,nfftf,ngfftf,nspden,dtset%nsppol,& 753 & psps%ntypat,optfr,paw_ij1,pawang,pawfgrtab,pawrad,pawtab,cplex,qphon,& 754 & rprimd,ucvol,vpsp1,vtrial,vxc,xred,& 755 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 756 757 if ((iscf_mod>=0.or.usexcnhat==0).and.(dtset%pawstgylm/=0)) then 758 ider=0;if ((ipert<=dtset%natom).and.(use_nhat_gga)) ider=1 759 call pawnhatfr(ider,idir_paw1,ipert,my_natom,dtset%natom,nspden,psps%ntypat,& 760 & pawang,pawfgrtab,pawrhoij,pawtab,rprimd,& 761 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 762 end if 763 end if 764 ! PAW only: we sometimes have to compute 1st-order compensation density 765 ! and eventually add it to density from 1st-order WFs 766 ! ---------------------------------------------------------------------- 767 nhat1grdim=0 768 ABI_MALLOC(nhat1gr,(0,0,0)) 769 if (psps%usepaw==1.and.ipert/=dtset%natom+1.and.ipert/=dtset%natom+10) then 770 call timab(564,1,tsec) 771 nhat1grdim=0;if (dtset%xclevel==2) nhat1grdim=usexcnhat*dtset%pawnhatxc 772 ider=2*nhat1grdim;izero=0 773 if (nhat1grdim>0) then 774 ABI_FREE(nhat1gr) 775 ABI_MALLOC(nhat1gr,(cplex*nfftf,dtset%nspden,3*nhat1grdim)) 776 end if 777 call pawmknhat(dum,cplex,ider,idir_paw1,ipert,izero,gprimd,my_natom,dtset%natom,& 778 & nfftf,ngfftf,nhat1grdim,nspden,psps%ntypat,pawang,pawfgrtab,nhat1gr,nhat1,& 779 & pawrhoij1,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred,& 780 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 781 if (dtfil%ireadwf/=0.and.dtset%get1den==0.and.dtset%ird1den==0.and.initialized==0) then 782 rhor1(:,:)=rhor1(:,:)+nhat1(:,:) 783 call fourdp(cplex,rhog1,rhor1(:,1),-1,mpi_enreg,nfftf,1, ngfftf,0) 784 end if 785 call timab(564,2,tsec) 786 end if 787 ! Set initial guess for 1st-order potential 788 ! ---------------------------------------------------------------------- 789 option=1;optene=0;if (iscf_mod==-2) optene=1 790 call dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,& 791 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,nhat1,nhat1gr,nhat1grdim,& 792 & nkxc,nspden,n3xccc,nmxc,optene,option,dtset%qptn,& 793 & rhog,rhog1,rhor,rhor1,rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,& 794 & nvresid1,res2,vtrial1,vxc,vxc1,xccc3d1,dtset%ixcrot) 795 796 if(.not.kramers_deg) then 797 vtrial1_pq=vtrial1 !save trial potential at +q 798 !rhor1_mq=rhor1 799 !rhog1_mq=rhog1 800 !get initial guess for vtrial1 at -q 801 do ifft=1,nfftf 802 vtrial1_mq(2*ifft-1,1)=+vtrial1(2*ifft-1,1) 803 vtrial1_mq(2*ifft-1,2)=+vtrial1(2*ifft-1,2) 804 vtrial1_mq(2*ifft ,1)=-vtrial1(2*ifft ,1) 805 vtrial1_mq(2*ifft ,2)=-vtrial1(2*ifft ,2) 806 vtrial1_mq(2*ifft-1,3)= vtrial1(2*ifft ,4) !Re[V^12] 807 vtrial1_mq(2*ifft ,3)= vtrial1(2*ifft-1,4) !Im[V^12],see definition of v(:,4) cplex=2 case 808 vtrial1_mq(2*ifft ,4)= vtrial1(2*ifft-1,3) !Re[V^21]=Re[V^12] 809 vtrial1_mq(2*ifft-1,4)= vtrial1(2*ifft ,3) !Re[V^21]=Re[V^12] 810 end do 811 end if 812 813 ! For Q=0 and metallic occupation, initialize quantities needed to 814 ! compute the first-order Fermi energy 815 ! ---------------------------------------------------------------------- 816 if (need_fermie1) then 817 ABI_MALLOC(rhorfermi,(cplex*nfftf,nspden)) 818 if(.not.kramers_deg) then 819 ABI_MALLOC(rhorfermi_mq,(cplex*nfftf,nspden)) 820 end if 821 if (psps%usepaw==1.and.usexcnhat==0) then 822 ABI_MALLOC(nhatfermi,(cplex*nfftf,nspden)) 823 else 824 ABI_MALLOC(nhatfermi,(0,0)) 825 end if 826 ABI_MALLOC(pawrhoijfermi,(my_natom*psps%usepaw)) 827 if (psps%usepaw==1) then 828 !Q phase should be 1 because Q=0 829 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 830 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 831 call pawrhoij_alloc(pawrhoijfermi,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 832 & dtset%nsppol,dtset%typat,pawtab=pawtab,mpi_atmtab=mpi_enreg%my_atmtab,& 833 & comm_atom=mpi_enreg%comm_atom) 834 end if 835 836 call dfpt_rhofermi(cg,cgq,cplex,cprj,cprjq,& 837 & doccde_rbz,docckqde,dtfil,dtset,eigenq,eigen0,eigen1,fe1fixed,gmet,gprimd,idir,& 838 & indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,mband_mem_rbz,mkmem,mkqmem,mk1mem,mpi_enreg,& 839 & mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1,& 840 & nspden,dtset%nsppol,nsym1,occkq,occ_rbz,& 841 & paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,& 842 & phnons1,ph1d,dtset%prtvol,psps,rhorfermi,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,& 843 & ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1) 844 if (.not.kramers_deg) then 845 call dfpt_rhofermi(cg,cg_mq,cplex,cprj,cprjq,& 846 & doccde_rbz,docckde_mq,dtfil,dtset,eigen_mq,eigen0,eigen1_mq,fe1fixed_mq,gmet,gprimd,idir,& 847 & indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1_mq,kpt_rbz,dtset%mband,mband_mem_rbz,mkmem,mkqmem,mk1mem,mpi_enreg,& 848 & mpw,mpw1_mq,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,ngfftf,nhatfermi,nkpt_rbz,npwarr,npwar1_mq,& 849 & nspden,dtset%nsppol,nsym1,occk_mq,occ_rbz,& 850 & paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoijfermi,pawtab,& 851 & phnons1,ph1d,dtset%prtvol,psps,rhorfermi_mq,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,& 852 & ucvol,usecprj,useylmgr1,vtrial,vxc,wtk_rbz,xred,ylm,ylm1,ylmgr1) 853 end if 854 855 end if 856 857 end if ! End the condition of istep==1 858 859 ! ###################################################################### 860 ! The following steps are done at every iteration 861 ! ---------------------------------------------------------------------- 862 863 if (psps%usepaw==1)then 864 ! Computation of "on-site" 2nd-order energy, first-order potentials, first-order densities 865 nzlmopt=0;if (istep==2.and.dtset%pawnzlm>0) nzlmopt=-1 866 if (istep>2) nzlmopt=dtset%pawnzlm 867 call paw_an_reset_flags(paw_an1) ! Force the recomputation of on-site potentials 868 call paw_ij_reset_flags(paw_ij1,self_consistent=.true.) ! Force the recomputation of Dij 869 option=0;if (dtset%iscf>0.and.dtset%iscf<10.and.nstep>0) option=1 870 call pawdenpot(dum,epaw1,epawdc1_dum,ipert,dtset%ixc,my_natom,dtset%natom,& 871 & dtset%nspden,psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an1,paw_an,paw_ij1,pawang,& 872 & dtset%pawprtvol,pawrad,pawrhoij1,dtset%pawspnorb,pawtab,dtset%pawxcdev,& 873 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp, & 874 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 875 876 ! First-order Dij computation 877 call timab(561,1,tsec) 878 if (has_dijfr>0) then 879 !vpsp1 contribution to Dij already stored in frozen part of Dij 880 ABI_MALLOC(vtrial1_tmp,(cplex*nfftf,nspden)) 881 vtrial1_tmp=vtrial1 882 do ispden=1,min(dtset%nspden,2) 883 vtrial1_tmp(:,ispden)=vtrial1_tmp(:,ispden)-vpsp1(:) 884 end do 885 else 886 vtrial1_tmp => vtrial1 887 end if 888 call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,dtset%natom,& 889 & nfftf,nfftotf,dtset%nspden,psps%ntypat,paw_an1,paw_ij1,pawang,& 890 & pawfgrtab,dtset%pawprtvol,pawrad,pawrhoij1,dtset%pawspnorb,pawtab,& 891 & dtset%pawxcdev,qphon,dtset%spnorbscl,ucvol,dtset%cellcharge(1),vtrial1_tmp,vxc1,xred,& 892 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 893 if (has_dijfr>0) then 894 ABI_FREE(vtrial1_tmp) 895 end if 896 897 call symdij(gprimd,indsy1,ipert,my_natom,dtset%natom,nsym1,psps%ntypat,0,& 898 & paw_ij1,pawang1,dtset%pawprtvol,pawtab,rprimd,symaf1,symrc1, & 899 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom,& 900 & qphon=qphon) 901 call timab(561,2,tsec) 902 end if ! end usepaw section 903 904 ! ###################################################################### 905 ! The following steps are done only when nstep>0 906 ! ---------------------------------------------------------------------- 907 908 if(iscf_mod>0.and.nstep>0)then 909 write(msg, '(a,a,i4)' )ch10,' ITER STEP NUMBER ',istep 910 call wrtout(std_out,msg,'COLL') 911 end if 912 913 ! For Q=0 and metallic occupation, calculate the first-order Fermi energy 914 if (need_fermie1) then 915 call newfermie1(cplex,fermie1,fe1fixed,ipert,istep,dtset%ixc,my_natom,dtset%natom,& 916 & nfftf,nfftotf,nhatfermi,nspden,dtset%ntypat,dtset%occopt,paw_an,paw_an1,paw_ij1,pawang,& 917 & dtset%pawnzlm,pawrad,pawrhoij1,pawrhoijfermi,pawtab,dtset%pawxcdev,& 918 & dtset%prtvol,rhorfermi,ucvol,psps%usepaw,usexcnhat,vtrial1,vxc1,dtset%xclevel,& 919 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 920 if (.not.kramers_deg) then 921 !fermie1_mq is updated as well at "-q" 922 call newfermie1(cplex,fermie1_mq,fe1fixed_mq,ipert,istep,dtset%ixc,my_natom,dtset%natom,& 923 & nfftf,nfftotf,nhatfermi,nspden,dtset%ntypat,dtset%occopt,paw_an,paw_an1,paw_ij1,pawang,& 924 & dtset%pawnzlm,pawrad,pawrhoij1,pawrhoijfermi,pawtab,dtset%pawxcdev,& 925 & dtset%prtvol,rhorfermi_mq,ucvol,psps%usepaw,usexcnhat,vtrial1_mq,vxc1_mq,dtset%xclevel,& 926 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 927 end if 928 end if 929 930 ! No need to continue and call dfpt_vtorho, when nstep==0 931 if(nstep==0) exit 932 933 ! #######################e1magh############################################### 934 ! Compute the 1st-order density rho1 from the 1st-order trial potential 935 ! ---------------------------------------------------------------------- 936 937 call dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,& 938 & dbl_nnsclo,dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,dtset%qptn,edocc,& 939 & eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,end0,end1,enl0,enl1,evxctau0,evxctau1,& 940 & fermie1,gh0c1_set,gh1c_set,& 941 & gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,mband_mem_rbz,& 942 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,& 943 & nhat1,nkpt_rbz,npwarr,npwar1,res2,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid1,& 944 & occkq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,& 945 & pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid,residm,rhog1,& 946 & rhor1,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,ucvol,usecprj,useylmgr1,ddk_f,& 947 & vectornd,vtrial,vtrial1,vxctau,with_vectornd,wtk_rbz,xred,ylm,ylm1,ylmgr1) 948 if (.not.kramers_deg) then 949 rhor1_pq=rhor1 !at this stage rhor1_pq contains only one term of the 1st order density at +q 950 rhog1_pq=rhog1 !same for rhog1_pq 951 !get the second term related to 1st order wf at -q 952 953 call dfpt_vtorho(cg,cg_mq,cg1_mq,cg1_active_mq,cplex,cprj,cprjq,cprj1,& 954 & dbl_nnsclo_mq,dim_eig2rf,doccde_rbz,docckde_mq,dtefield,dtfil,dtset,-dtset%qptn,edocc_mq,& 955 & eeig0_mq,eigen_mq,eigen0,eigen1_mq,ek0_mq,ek1_mq,eloc0_mq,end0_mq,end1_mq,& 956 & enl0_mq,enl1_mq,evxctau0_mq,evxctau1_mq,fermie1_mq,gh0c1_set_mq,gh1c_set_mq,& 957 & gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1_mq,kpt_rbz,dtset%mband,mband_mem_rbz,& 958 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1_mq,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,& 959 & nhat1,nkpt_rbz,npwarr,npwar1_mq,res2_mq,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid1_mq,& 960 & occk_mq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,& 961 & pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid_mq,residm_mq,rhog1_mq,& 962 & rhor1_mq,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,ucvol,usecprj,useylmgr1,ddk_f,& 963 & vectornd,vtrial,vtrial1_mq,vxctau,with_vectornd,wtk_rbz,xred,ylm,ylm1,ylmgr1) 964 !reconstruct the +q and -q densities, this might bug if fft parallelization is used, todo... 965 do ifft=1,nfftf 966 rhor1_pq(2*ifft-1,:) = half*(rhor1(2*ifft-1,:)+rhor1_mq(2*ifft-1,:)) 967 rhor1_pq(2*ifft ,:) = half*(rhor1(2*ifft ,:)-rhor1_mq(2*ifft ,:)) 968 rhor1_mq(2*ifft-1,:) = rhor1_pq(2*ifft-1,:) 969 rhor1_mq(2*ifft ,:) =-rhor1_pq(2*ifft ,:) 970 end do 971 rhor1=rhor1_pq 972 call fourdp(cplex,rhog1,rhor1(:,1),-1,mpi_enreg,nfftf,1, ngfftf, 0) 973 call fourdp(cplex,rhog1_mq,rhor1_mq(:,1),-1,mpi_enreg,nfftf,1, ngfftf, 0) 974 end if 975 976 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 977 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 978 979 ! calculate \Omega E \cdot P term 980 if (ipert<=dtset%natom) then 981 ! phonon perturbation 982 call dfptff_ebp(cg,cg1,dtefield,eberry,dtset%mband,mband_mem_rbz,mkmem,& 983 & mpi_enreg,mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat) 984 else if (ipert==dtset%natom+2) then 985 ! electric field perturbation 986 call dfptff_edie(cg,cg1,dtefield,eberry,idir,dtset%mband,mband_mem_rbz,mkmem,& 987 & mpi_enreg,mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd) 988 end if 989 end if 990 991 ! SPr: don't remove the following comments for debugging 992 ! call calcdenmagsph(mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,& 993 !& dtset%ntypat,dtset%ratsm,dtset%ratsph,rhor1,rprimd,dtset%typat,xred,& 994 !& idir+1,cplex,intgden=intgden,rhomag=rhomag) 995 ! call prtdenmagsph(cplex,intgden,dtset%natom,nspden,dtset%ntypat,ab_out,idir+1,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat) 996 997 ! write(*,*) ' n ( 1,2)',intgden(1,1),' ',intgden(1,2) 998 ! write(*,*) ' mx( 1,2)',intgden(2,1),' ',intgden(2,2) 999 ! write(*,*) ' my( 1,2)',intgden(3,1),' ',intgden(3,2) 1000 ! write(*,*) ' mz( 1,2)',intgden(4,1),' ',intgden(4,2) 1001 ! call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,& 1002 !& efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,& 1003 !& enl0,enl1,epaw1,etotal,evar,evdw,exc1,elmag1,ipert,dtset%natom,optene) 1004 ! write(*,*) 'SPr: ek1=',ek1,' exc1=',exc1,' elmag1=',elmag 1005 ! if (ipert==dtset%natom+5) then 1006 ! !calculate 1st order magnetic potential contribution to the energy 1007 ! call dfpt_e1mag(e1mag,rhor1,rhog1); 1008 ! endif 1009 1010 ! ###################################################################### 1011 ! Skip out of step loop if non-SCF (completed) 1012 ! ---------------------------------------------------------------------- 1013 1014 ! Indeed, nstep loops have been done inside dfpt_vtorho 1015 if (iscf_mod<=0 .and. iscf_mod/=-3) exit 1016 1017 ! ###################################################################### 1018 ! In case of density mixing , compute the total 2nd-order energy, 1019 ! check the exit criterion, then mix the 1st-order density 1020 ! ---------------------------------------------------------------------- 1021 1022 if (iscf_mod>=10) then 1023 optene = 1 ! use double counting scheme 1024 call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,& 1025 & efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,& 1026 & end0,enl0,enl1,epaw1,etotal,evar,evdw,evxctau0,exc1,ipert,dtset%natom,optene) 1027 call timab(152,1,tsec) 1028 choice=2 1029 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 1030 & etotal,favg,fcart,fermie,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 1031 & 1,iscf_mod,istep,istep_fock_outer,istep_mix,kpt_rbz,maxfor,& 1032 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 1033 & nstep,occ_rbz,0,prtfor,0,& 1034 & quit,res2,resid,residm,response,& 1035 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 1036 call timab(152,2,tsec) 1037 1038 if (istep==nstep) quit=1 1039 ! If criteria in scprqt say to quit, then exit the loop over istep 1040 quit_sum=quit 1041 call xmpi_sum(quit_sum,spaceComm,ierr) 1042 1043 if (quit_sum>0) exit 1044 ! INSERT HERE CALL TO NEWRHO3 : to be implemented 1045 if (psps%usepaw==1) then 1046 ABI_BUG("newrho3 not implemented: use potential mixing!") 1047 end if 1048 initialized=1 1049 end if 1050 1051 ! ###################################################################### 1052 ! Compute the new 1st-order potential from the 1st-order density 1053 ! ---------------------------------------------------------------------- 1054 1055 if (ipert<dtset%natom+10) then 1056 optene=1 1057 call dfpt_rhotov(cplex,ehart01,ehart1,elpsp1,exc1,elmag1,gsqcut,idir,ipert,& 1058 & dtset%ixc,kxc,mpi_enreg,dtset%natom,nfftf,ngfftf,nhat,nhat1,nhat1gr,nhat1grdim,nkxc,& 1059 & nspden,n3xccc,nmxc,optene,optres,dtset%qptn,rhog,rhog1,rhor,rhor1,& 1060 & rprimd,ucvol,psps%usepaw,usexcnhat,vhartr1,vpsp1,nvresid1,res2,vtrial1,vxc,vxc1,xccc3d1,dtset%ixcrot) 1061 end if 1062 1063 ! ###################################################################### 1064 ! In case of potential mixing , compute the total 2nd-order energy, 1065 ! check the exit criterion, then mix the 1st-order potential 1066 ! ---------------------------------------------------------------------- 1067 1068 if (iscf_mod<10) then 1069 1070 ! PAW: has to compute here the "on-site" 2nd-order energy 1071 if (psps%usepaw==1) then 1072 nzlmopt=0;if (istep==1.and.dtset%pawnzlm>0) nzlmopt=-1 1073 if (istep>1) nzlmopt=dtset%pawnzlm 1074 call paw_an_reset_flags(paw_an1) ! Force the recomputation of on-site potentials 1075 option=2 1076 call pawdenpot(dum,epaw1,epawdc1_dum,ipert,dtset%ixc,my_natom,dtset%natom,dtset%nspden,& 1077 & psps%ntypat,dtset%nucdipmom,nzlmopt,option,paw_an1,paw_an,paw_ij1,pawang,dtset%pawprtvol,& 1078 & pawrad,pawrhoij1,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%spnorbscl,& 1079 & dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp,& 1080 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1081 end if 1082 1083 optene = 0 ! use direct scheme 1084 call dfpt_etot(dtset%berryopt,deltae,eberry,edocc,eeig0,eew,efrhar,efrkin,& 1085 & efrloc,efrnl,efrx1,efrx2,ehart1,ek0,ek1,eii,elast,eloc0,elpsp1,& 1086 !& enl0,enl1,epaw1,etotal,evar,evdw,exc1,elmag1,ipert,dtset%natom,optene) 1087 ! !debug: compute the d2E/d-qd+q energy, should be equal to the one from previous line 1088 ! if(.not.kramers_deg) then 1089 ! call dfpt_etot(dtset%berryopt,deltae_mq,eberry_mq,edocc_mq,eeig0_mq,eew,efrhar,efrkin,& 1090 !& efrloc,efrnl,efrx1,efrx2,ehart1_mq,ek0_mq,ek1_mq,eii,elast_mq,eloc0_mq,elpsp1_mq,& 1091 !& enl0_mq,enl1_mq,epaw1_mq,etotal_mq,evar_mq,evdw,exc1_mq,elmag1_mq,ipert,dtset%natom,optene) 1092 ! end if 1093 & end0,enl0,enl1,epaw1,etotal,evar,evdw,evxctau0,exc1,ipert,dtset%natom,optene) 1094 1095 call timab(152,1,tsec) 1096 choice=2 1097 ! To take into account new definition of hdr_update; 1098 ! test to avoid dfpt and occopt 9 was already done 1099 ! so we can just set fermih = fermie 1100 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 1101 & etotal,favg,fcart,fermie,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 1102 & 1,iscf_mod,istep,istep_fock_outer,istep_mix,kpt_rbz,maxfor,& 1103 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 1104 & nstep,occ_rbz,0,prtfor,0,& 1105 & quit,res2,resid,residm,response,& 1106 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 1107 ! endif 1108 call timab(152,2,tsec) 1109 1110 ! If criteria in scprqt say to quit, then exit the loop over istep 1111 quit_sum=quit 1112 call xmpi_sum(quit_sum,spaceComm,ierr) 1113 if (quit_sum>0) exit 1114 1115 ! TODO 1116 ! Better error handling is the SCF cycle goes bananas: 1117 ! Write a BIG warning in the output file and save the wavefunctions. 1118 ! so that we can restart. 1119 if(iscf_mod/=-3)then 1120 ! Note that nvresid1 and vtrial1 are called vresid and vtrial inside this routine 1121 call dfpt_newvtr(cplex,dbl_nnsclo,dielar,dtset,etotal,pawfgr%fintocoa,& 1122 & initialized,iscf_mod,ispmix,istep,mix,pawfgr%coatofin,& 1123 & mpi_enreg,my_natom,nfftf,nfftmix,ngfftf,ngfftmix,npawmix,pawrhoij1,& 1124 & qphon,rhor1,rprimd,psps%usepaw,nvresid1,vtrial1) 1125 if (.not.kramers_deg) then 1126 !same problem as with density reconstruction, TODO proper fft parallelization... 1127 do ifft=1,nfftf 1128 vtrial1_mq(2*ifft-1,1)=+vtrial1(2*ifft-1,1) 1129 vtrial1_mq(2*ifft-1,2)=+vtrial1(2*ifft-1,2) 1130 vtrial1_mq(2*ifft ,1)=-vtrial1(2*ifft ,1) 1131 vtrial1_mq(2*ifft ,2)=-vtrial1(2*ifft ,2) 1132 vtrial1_mq(2*ifft-1,3)= vtrial1(2*ifft ,4) !Re[V^12] 1133 vtrial1_mq(2*ifft ,3)= vtrial1(2*ifft-1,4) !Im[V^12],see definition of v(:,4) cplex=2 case 1134 vtrial1_mq(2*ifft ,4)= vtrial1(2*ifft-1,3) !Re[V^21]=Re[V^12] 1135 vtrial1_mq(2*ifft-1,4)= vtrial1(2*ifft ,3) !Re[V^21]=Re[V^12] 1136 end do 1137 end if 1138 initialized=1 1139 end if 1140 end if 1141 1142 ! ###################################################################### 1143 ! END MINIMIZATION ITERATIONS 1144 ! Note that there are different "exit" instructions within the loop 1145 ! ###################################################################### 1146 end do ! istep 1147 1148 ! Avoid pending requests if itime == ntime. 1149 call xmpi_wait(quitsum_request,ierr) 1150 if (timelimit_exit == 1) istep = istep - 1 1151 1152 !SP : Here read the _DDB file and extract the Born effective charge and 1153 ! dielectric constant. 1154 ! The idea is to supress the divergence due to a residual Born effective charge 1155 ! by renormalizing the v_hart1. For this, the difference between the ionic 1156 ! Z_kappa and the Born effective charge divided by the dielectric constant is used. 1157 ! --------------------------------------------------------------------------------- 1158 if ((dtset%getddb .ne. 0 .or. dtset%irdddb .ne.0) .and. qzero .ne. 1) then 1159 ABI_MALLOC(rhor2,(cplex*nfftf,nspden)) 1160 ABI_MALLOC(resid2,(dtset%mband*nkpt_rbz*nspden)) 1161 ABI_MALLOC(rhog2,(2,nfftf)) 1162 ABI_MALLOC(vtrial2,(cplex*nfftf,nspden)) 1163 1164 Z_kappa = nint(psps%ziontypat(dtset%typat(ipert))) ! Charge ionic from the psp 1165 qred2cart = two_pi*gprimd 1166 q_cart = MATMUL(qred2cart,qphon) 1167 q_cart = q_cart/SQRT(dot_product(q_cart,q_cart)) 1168 diel_q = dot_product(MATMUL(dielt,q_cart),q_cart) 1169 zeff_bar = SUM(zeff(:,:,:),DIM=3)/dtset%natom 1170 zeff_red = MATMUL(zeff_bar(:,:),rprimd(:,idir))/two_pi 1171 qphon_norm = SQRT(dot_product(qphon,qphon)) 1172 q_cart = MATMUL(qred2cart,qphon) 1173 born_bar = dot_product(q_cart,zeff_red(:)) 1174 zeff_red = MATMUL(zeff(:,:,ipert),rprimd(:,idir))/two_pi 1175 born = dot_product(q_cart,zeff_red(:)) 1176 1177 ! To avoid problem of divergence (0/0) we add a small value to qphon 1178 qphon2 = qphon + tol6 1179 renorm = (1-(qphon2(idir)*Z_kappa-(born-born_bar)/diel_q)/(qphon2(idir)*Z_kappa-born/diel_q)) 1180 1181 vtrial2(:,1) = vtrial1(:,1) -renorm*vhartr1 1182 1183 call dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,& 1184 & dbl_nnsclo,dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,dtset%qptn,edocc,& 1185 & eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,end0,end1,enl0,enl1,evxctau0,evxctau1,& 1186 & fermie1,gh0c1_set,gh1c_set,& 1187 & gmet,gprimd,idir,indsy1,ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,mband_mem_rbz,& 1188 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,dtset%natom,nband_rbz,ncpgr,nfftf,& 1189 & nhat1,nkpt_rbz,npwarr,npwar1,res3,nspden,dtset%nsppol,nsym1,dtset%ntypat,nvresid2,& 1190 & occkq,occ_rbz,optres,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,& 1191 & pawrhoij1,pawtab,phnons1,ph1d,dtset%prtvol,psps,pwindall,qmat,resid2,residm2,rhog2,& 1192 & rhor2,rmet,rprimd,symaf1,symrc1,symrl1,tnons1,ucvol,usecprj,useylmgr1,ddk_f,& 1193 & vectornd,vtrial,vtrial2,vxctau,with_vectornd,wtk_rbz,xred,ylm,ylm1,ylmgr1,1) 1194 1195 write(msg,'(a)') ' '//char(10)//& 1196 ' ---------------------------------' 1197 call wrtout(ab_out,msg,'COLL') 1198 write(msg,'(a,a)')' The charge sum rule is activated'//char(10)//& 1199 ' ---------------------------------' 1200 call wrtout(ab_out,msg,'COLL') 1201 write(msg,'(a,i4)') ' Z_ion (psp):',Z_kappa 1202 call wrtout(ab_out,msg,'COLL') 1203 write(msg,'(a,f12.8)') ' Residual Born effective charge: ',born 1204 call wrtout(ab_out,msg,'COLL') 1205 write(msg,'(a,f12.8)') ' Renormalisation: ',renorm 1206 call wrtout(ab_out,msg,'COLL') 1207 if (renorm > 0.01 ) then 1208 write(msg,'(a,a)')' WARNING: The renormalisation seems large (> 0.01).'//char(10)//& 1209 ' You might consider increasing the k-point grid.' 1210 ABI_WARNING(msg) 1211 call wrtout(ab_out,msg,'COLL') 1212 end if 1213 write(msg,'(a)') ' ' 1214 call wrtout(ab_out,msg,'COLL') 1215 1216 ABI_FREE(nvresid2) 1217 ABI_FREE(rhor2) 1218 ABI_FREE(resid2) 1219 ABI_FREE(rhog2) 1220 ABI_FREE(vtrial2) 1221 end if 1222 1223 if (iscf_mod>0.or.iscf_mod==-3) then 1224 ABI_FREE(nvresid1) 1225 if (.not.kramers_deg) then 1226 ABI_FREE(nvresid1_mq) 1227 end if 1228 end if 1229 1230 !###################################################################### 1231 !Additional steps after SC iterations 1232 !---------------------------------------------------------------------- 1233 1234 call timab(160,1,tsec) 1235 1236 !Compute Dynamic magnetic charges (dmc) in case of rfphon, 1237 !and magnetic susceptibility in case of rfmagn from first order density 1238 !(results to be comapred to dmc from d2e) 1239 !SPr deb 1240 !if (ipert<=dtset%natom.and.dtset%nspden>=2) then 1241 ! 1242 ! mpi_comm_sphgrid=mpi_enreg%comm_fft 1243 ! call mean_fftr(rhor1(:,1),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1244 ! write(*,*) ' Mean 1st order density: ', mean_rhor1 1245 ! call mean_fftr(rhor1(:,2),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1246 ! if (dtset%nspden==2) then 1247 ! write(*,*) ' 1st order m_z : ', mean_rhor1 1248 ! else !nspden==4 1249 ! write(*,*) ' 1st order m_x : ', mean_rhor1 1250 ! call mean_fftr(rhor1(:,3),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1251 ! write(*,*) ' 1st order m_y : ', mean_rhor1 1252 ! call mean_fftr(rhor1(:,4),mean_rhor1,nfftf,nfftotf,1,mpi_comm_sphgrid) 1253 ! write(*,*) ' 1st order m_z : ', mean_rhor1 1254 ! endif 1255 ! 1256 !endif 1257 1258 1259 !Eventually close the dot file, before calling dfpt_nstdy 1260 if ((ipert==dtset%natom+2.and.sum((dtset%qptn(1:3))**2)<=1.0d-7.and.& 1261 & (dtset%berryopt/=4 .and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.& 1262 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17)).or.& 1263 & ipert==dtset%natom+10.or.ipert==dtset%natom+11) then 1264 call ddk_f(1)%close() 1265 end if 1266 if ((ipert==dtset%natom+10 .and. idir>3) .or. ipert==dtset%natom+11) then 1267 call ddk_f(2)%close() 1268 end if 1269 if (ipert==dtset%natom+11) then 1270 call ddk_f(3)%close() 1271 if(idir>3) call ddk_f(4)%close() 1272 end if 1273 1274 !Deallocate the no more needed arrays 1275 if (iscf_mod>0.and.nstep>0) then 1276 call ab7_mixing_deallocate(mix) 1277 end if 1278 if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then 1279 ABI_FREE(dielinv) 1280 ABI_FREE(susmat) 1281 end if 1282 if(allocated(rhorfermi)) then 1283 ABI_FREE(rhorfermi) 1284 end if 1285 if(allocated(rhorfermi_mq)) then 1286 ABI_FREE(rhorfermi_mq) 1287 end if 1288 if(allocated(nhatfermi)) then 1289 ABI_FREE(nhatfermi) 1290 end if 1291 if(allocated(pawrhoijfermi)) then 1292 call pawrhoij_free(pawrhoijfermi) 1293 ABI_FREE(pawrhoijfermi) 1294 end if 1295 if(psps%usepaw==1) then 1296 if (mk1mem/=0.and.usecprj==1) then 1297 call pawcprj_free(cprj1) 1298 end if 1299 do iatom=1,my_natom 1300 if (pawfgrtab(iatom)%nhatfr_allocated>0) then 1301 ABI_FREE(pawfgrtab(iatom)%nhatfr) 1302 end if 1303 pawfgrtab(iatom)%nhatfr_allocated=0 1304 end do 1305 if (nstep>0.and.iscf_mod>0) then 1306 do iatom=1,my_natom 1307 pawrhoij1(iatom)%lmnmix_sz=0 1308 pawrhoij1(iatom)%use_rhoijres=0 1309 ABI_FREE(pawrhoij1(iatom)%kpawmix) 1310 ABI_FREE(pawrhoij1(iatom)%rhoijres) 1311 end do 1312 end if 1313 end if ! PAW 1314 ABI_FREE(cprj1) 1315 ABI_FREE(nhat1gr) 1316 1317 call timab(160,2,tsec) 1318 call timab(150,1,tsec) 1319 1320 if (psps%usepaw==0.and.dtset%userie/=919.and. & 1321 & (ipert==dtset%natom+3.or.ipert==dtset%natom+4)) then 1322 call dfpt_nselt(blkflg,cg,cg1,cplex,& 1323 & d2bbb,d2lo,d2nl,ecut,dtset%ecutsm,dtset%effmass_free,& 1324 & gmet,gprimd,gsqcut,idir,& 1325 & ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,dtset%mband,mband_mem_rbz,mgfftf,& 1326 & mkmem,mk1mem,mpert,mpi_enreg,psps%mpsang,mpw,mpw1,& 1327 & dtset%natom,nband_rbz,nfftf,ngfftf,& 1328 & nkpt_rbz,nkxc,dtset%nloalg,& 1329 & npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,& 1330 & nsym1,dtset%ntypat,occ_rbz,& 1331 & ph1d,dtset%prtbbb,psps,dtset%qptn,rhog,& 1332 & rhor,rhor1,rmet,rprimd,symrc1,dtset%typat,ucvol,& 1333 & wtk_rbz,xred,ylm,ylm1,ylmgr,ylmgr1,& 1334 & rfstrs_ref=dtset%rfstrs_ref) 1335 end if 1336 1337 !Use of NSTPAW3 for NCPP (instead of DFPT_NSELT/DFPT_NSTDY) can be forced with userie=919 1338 !MT oct. 2015: this works perfectly on all automatic tests 1339 if(ipert<=dtset%natom+4)then 1340 if (psps%usepaw==1.or.dtset%userie==919) then 1341 call dfpt_nstpaw(blkflg,cg,cgq,cg1,cplex,cprj,cprjq,docckqde,doccde_rbz,dtfil,dtset,d2lo,d2nl,d2ovl,& 1342 & eigenq,eigen0,eigen1,eovl1,gmet,gprimd,gsqcut,idir,indkpt1,indsy1,ipert,irrzon1,istwfk_rbz,& 1343 & kg,kg1,kpt_rbz,kxc,mgfftf,mkmem,mkqmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,mband_mem_rbz,ncpgr,& 1344 & nfftf,ngfftf,nhat,nhat1,nkpt_rbz,nkxc,npwarr,npwar1,nspden,dtset%nspinor,dtset%nsppol,& 1345 & nsym1,n3xccc,occkq,occ_rbz,paw_an,paw_an1,paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrad,& 1346 & pawrhoij,pawrhoij1,pawtab,phnons1,ph1d,ph1df,psps,rhog,rhor,rhor1,rmet,rprimd,symaf1,symrc1,& 1347 & symrl1,tnons1,ucvol,usecprj,psps%usepaw,usexcnhat,useylmgr1,vectornd,vhartr1,vpsp1,vtrial,vtrial1,vxc,& 1348 & with_vectornd,wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr1) 1349 else 1350 if (dtset%nspden==4) then 1351 call dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,gmet,& 1352 & gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mkmem,mk1mem,mpert,mpi_enreg,& 1353 & mpw,mpw1,nattyp,nband_rbz,nfftf,ngfftf,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,& 1354 & dtset%nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,symrc1,ucvol,& 1355 & wtk_rbz,xred,ylm,ylm1,rhor=rhor,vxc=vxc) 1356 else 1357 call dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,gmet,& 1358 & gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mband_mem_rbz,mkmem,mk1mem,mpert,mpi_enreg,& 1359 & mpw,mpw1,nattyp,nband_rbz,nfftf,ngfftf,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,& 1360 & dtset%nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,symrc1,ucvol,& 1361 & wtk_rbz,xred,ylm,ylm1) 1362 end if 1363 end if 1364 end if 1365 1366 call timab(150,2,tsec) 1367 call timab(160,1,tsec) 1368 1369 !calculate Born effective charge and store it in d2lo 1370 if ((dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1371 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17).and.& 1372 & ipert<=dtset%natom) then 1373 call dfptff_bec(cg,cg1,dtefield,dtset%natom,d2lo,idir,ipert,dtset%mband,mband_mem_rbz,mkmem,& 1374 & mpi_enreg,mpw,mpw1,mpert,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd) 1375 blkflg(:,dtset%natom+2,:,1:dtset%natom)=1 1376 end if 1377 1378 !calculate dielectric tensor and store it in d2lo 1379 if ((dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1380 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17).and.& 1381 & ipert==dtset%natom+2) then 1382 call dfptff_die(cg,cg1,dtefield,d2lo,idir,ipert,dtset%mband,mband_mem_rbz,mkmem,& 1383 & mpi_enreg,mpw,mpw1,mpert,nkpt,npwarr,npwar1,dtset%nsppol,dtset%nspinor,pwindall,qmat,rprimd) 1384 blkflg(:,dtset%natom+2,:,dtset%natom+2)=1 1385 end if 1386 1387 !If SCF convergence was not reached (for nstep>0), 1388 !print a warning to the output file (non-dummy arguments: nstep, 1389 !residm, diffor - infos from tollist have been saved inside ) 1390 !Set also the value of conv_retcode 1391 choice=3 1392 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,& 1393 & etotal,favg,fcart,fermie,fermie,dtfil%fnametmp_eig,dtfil%filnam_ds(1),& 1394 & 1,iscf_mod,istep,istep_fock_outer,istep_mix,kpt_rbz,maxfor,& 1395 & mvdum,mpi_enreg,nband_rbz,nkpt_rbz,& 1396 & nstep,occ_rbz,0,prtfor,0,& 1397 & quit,res2,resid,residm,response,& 1398 & tollist,psps%usepaw,vxcavg,wtk_rbz,xred,conv_retcode) 1399 1400 !Update the content of the header (evolving variables) 1401 bantot_rbz = sum(nband_rbz(1:nkpt_rbz*dtset%nsppol)) 1402 call hdr%update(bantot_rbz,etotal,fermie,fermie,& 1403 & residm,rprimd,occ_rbz,pawrhoij1,xred,dtset%amu_orig(:,1),& 1404 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab ) 1405 1406 !Optionally provide output of charge density and/or potential in real space, 1407 !as well as analysis of geometrical factors (bond lengths and bond angles). 1408 !Warnings : 1409 !- core charge is excluded from the charge density; 1410 !- the potential is the INPUT vtrial. 1411 1412 if(ipert==dtset%natom+5.or.ipert<=dtset%natom)then 1413 prtopt=1 1414 if(ipert==dtset%natom+5) then 1415 prtopt=idir+1; 1416 call calcdenmagsph(mpi_enreg,dtset%natom,nfftf,ngfftf,nspden,& 1417 & dtset%ntypat,dtset%ratsm,dtset%ratsph,rhor1,rprimd,dtset%typat,xred,& 1418 & prtopt,cplex,intgden=intgden,dentot=dentot,rhomag=rhomag) 1419 call prtdenmagsph(cplex,intgden,dtset%natom,nspden,dtset%ntypat,ab_out,prtopt,dtset%ratsm,dtset%ratsph,rhomag,dtset%typat) 1420 end if 1421 end if 1422 1423 if (iwrite_fftdatar(mpi_enreg)) then 1424 if (dtset%prtden>0) then 1425 rdwrpaw=0 1426 call appdig(pertcase,dtfil%fnameabo_den,fi1o) 1427 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1428 call fftdatar_write_from_hdr("first_order_density",fi1o,dtset%iomode,hdr,& 1429 ngfftf,cplex,nfftf,dtset%nspden,rhor1,mpi_enreg) 1430 end if 1431 1432 ! Write first order potentials (needed by EPH) 1433 ! In DFPT, prtpot is automatically set to 1 unless the user set it to 0 explictly in the input 1434 ! See invars2 1435 ! (actually we should avoid writing 1WFK) 1436 if (dtset%prtpot > 0) then 1437 rdwrpaw=0 1438 call appdig(pertcase,dtfil%fnameabo_pot,fi1o) 1439 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1440 call fftdatar_write_from_hdr("first_order_potential",fi1o,dtset%iomode,hdr,& 1441 ngfftf,cplex,nfftf,dtset%nspden,vtrial1,mpi_enreg) 1442 1443 ! Add rhog1(G=0) to file 1444 ! This part is obsolete. I keep it just to maintain compatibility with the fileformat. 1445 if (mpi_enreg%me_g0 == 1) then 1446 if (dtset%iomode == IO_MODE_ETSF) then 1447 #ifdef HAVE_NETCDF 1448 NCF_CHECK(nctk_open_modify(ncid, nctk_ncify(fi1o), xmpi_comm_self)) 1449 ncerr = nctk_def_one_array(ncid, nctkarr_t('rhog1_g0', "dp", "two"), varid=varid) 1450 NCF_CHECK(ncerr) 1451 NCF_CHECK(nctk_set_datamode(ncid)) 1452 NCF_CHECK(nf90_put_var(ncid, varid, rhog1(:,1))) 1453 NCF_CHECK(nf90_close(ncid)) 1454 #endif 1455 else 1456 ! Handle Fortran files. 1457 if (open_file(fi1o, msg, newunit=ncid, form='unformatted', status='old', action="readwrite") /= 0) then 1458 ABI_ERROR(msg) 1459 end if 1460 if (fort_denpot_skip(ncid, msg) /= 0) ABI_ERROR(msg) 1461 write(ncid) rhog1(:,1) 1462 close(ncid) 1463 end if 1464 end if 1465 1466 end if 1467 1468 ! output files for perturbed potential components: vhartr1,vpsp1,vxc 1469 ! NB: only 1 spin for these 1470 if (dtset%prtvha > 0) then 1471 rdwrpaw=0 1472 ABI_MALLOC(vhartr1_tmp, (cplex*nfftf, dtset%nspden)) 1473 vhartr1_tmp = zero 1474 vhartr1_tmp(:,1) = vhartr1(:) 1475 call appdig(pertcase,dtfil%fnameabo_vha,fi1o) 1476 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1477 call fftdatar_write_from_hdr("first_order_vhartree",fi1o,dtset%iomode,hdr,& 1478 ngfftf,cplex,nfftf,dtset%nspden,vhartr1_tmp,mpi_enreg) 1479 ABI_FREE(vhartr1_tmp) 1480 end if 1481 1482 ! vpsp1 needs to be copied to a temp array - intent(inout) in fftdatar_write_from_hdr though I do not know why 1483 ! if (dtset%prtvpsp > 0) then 1484 ! rdwrpaw=0 1485 ! call appdig(pertcase,dtfil%fnameabo_vpsp,fi1o) 1486 ! ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1487 ! call fftdatar_write_from_hdr("first_order_vpsp",fi1o,dtset%iomode,hdr,& 1488 ! ngfftf,cplex,nfftf,1,vpsp1,mpi_enreg) 1489 ! end if 1490 1491 if (dtset%prtvxc > 0) then 1492 rdwrpaw=0 1493 call appdig(pertcase,dtfil%fnameabo_vxc,fi1o) 1494 ! TODO: should we write pawrhoij1 or pawrhoij. Note that ioarr writes hdr%pawrhoij 1495 call fftdatar_write_from_hdr("first_order_vxc",fi1o,dtset%iomode,hdr,& 1496 ngfftf,cplex,nfftf,dtset%nspden,vxc1,mpi_enreg) 1497 end if 1498 1499 end if ! iwrite_fftdatar(mpi_enreg) 1500 1501 !All procs waiting here... 1502 if(mpi_enreg%paral_kgb==1)then 1503 call timab(61,1,tsec) 1504 call xmpi_barrier(spaceComm) 1505 call timab(61,2,tsec) 1506 end if 1507 1508 !Deallocate arrays 1509 ABI_FREE(fcart) 1510 ABI_FREE(vtrial1) 1511 if (.not.kramers_deg) then 1512 ABI_FREE(vhartr1_mq) 1513 ABI_FREE(vxc1_mq) 1514 ABI_FREE(vtrial1_pq) 1515 ABI_FREE(vtrial1_mq) 1516 end if 1517 ABI_FREE(vhartr1) 1518 ABI_FREE(vxc1) 1519 ABI_FREE(pwindall) 1520 ABI_FREE(qmat) 1521 if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.& 1522 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then 1523 call destroy_efield(dtefield) 1524 if(allocated(mpi_enreg%kpt_loc2ibz_sp)) then 1525 ABI_FREE(mpi_enreg%kpt_loc2ibz_sp) 1526 end if 1527 end if 1528 1529 if(ALLOCATED(vectornd)) then 1530 ABI_FREE(vectornd) 1531 end if 1532 1533 if(psps%usepaw==1) then 1534 call paw_an_free(paw_an1) 1535 call paw_ij_free(paw_ij1) 1536 end if 1537 ABI_FREE(paw_an1) 1538 ABI_FREE(paw_ij1) 1539 ABI_FREE(nhat1) 1540 1541 call timab(160,2,tsec) 1542 call timab(120,2,tsec) 1543 1544 DBG_EXIT("COLL") 1545 1546 end subroutine dfpt_scfcv
ABINIT/dfpt_wfkfermi [ Functions ]
NAME
dfpt_wfkfermi
FUNCTION
This routine computes the partial Fermi-level density at a given k-point, and the fixed contribution to the 1st-order Fermi energy (nonlocal and kinetic)
INPUTS
cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)=planewave coefficients of wavefunctions
cgq(2,mcgq)=array for planewave coefficients of wavefunctions.
cplex=1 if rhoaug is real, 2 if rhoaug is complex
cprj(natom,nspinor*mband_mem*mkmem*nsppol*usecprj)= wave functions at k
projected with non-local projectors: cprj=<p_i|Cnk>
cprjq(natom,nspinor*mband_mem*mkqmem*nsppol*usecprj)= wave functions at k+q
projected with non-local projectors: cprjq=<p_i|Cnk+q>
dtfil <type(datafiles_type)>=variables related to files
eig0_k(nband_k)=GS eigenvalues at k (hartree)
gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
ibg=shift to be applied on the location of data in the array cprj
ibgq=shift to be applied on the location of data in the array cprjq
icg=shift to be applied on the location of data in the array cg
icgq=shift to be applied on the location of data in the array cgq
idir=direction of the current perturbation
ikpt=number of the k-point
ipert=type of the perturbation
isppol=1 for unpolarized, 2 for spin-polarized
kptopt=option for the generation of k points
mband_mem=maximum number of bands on this cpu
mcgq=second dimension of the cgq array
mcprjq=second dimension of the cprjq array
mkmem =number of k points treated by this node.
mpi_enreg=information about MPI parallelization
mpw=maximum dimensioned size of npw or wfs at k
nband_k=number of bands at this k point for that spin polarization
ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
npw_k=number of plane waves at this k point
npw1_k=number of plane waves at this k+q point
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
occ_k(nband_k)=occupation number for each band (usually 2) for each k.
prtvol=control print volume and debugging output
rf_hamkq <type(gs_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
rhoaug(cplex*n4,n5,n6)= density in electrons/bohr**3,
on the augmented fft grid. (cumulative, so input as well as output)
rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
if this ratio has been attributed to the band n (second argument), zero otherwise
wtk_k=weight assigned to the k point.
OUTPUT
eig1_k(2*nband_k**2)=first-order eigenvalues (hartree) fe1fixed_k(nband_k)=contribution to 1st-order Fermi energy from changes of occupation from all bands at this k point. fe1norm_k(nband_k)=contribution to normalization for above rhoaug(cplex*n4,n5,n6)= Fermi-level density in electrons/bohr**3, on the augmented fft grid. (cumulative, so input as well as output). ==== if (gs_hamkq%usepaw==1) ==== pawrhoijfermi(natom) <type(pawrhoij_type)>= paw rhoij occupancies at Fermi level (cumulative, so input as well as output)
SOURCE
4233 subroutine dfpt_wfkfermi(cg,cgq,cplex,cprj,cprjq,& 4234 & dtfil,eig0_k,eig1_k,fe1fixed_k,fe1norm_k,gs_hamkq,& 4235 & ibg,ibgq,icg,icgq,idir,ikpt,ipert,isppol,& 4236 & kptopt,mband_mem,mcgq,mcprjq,mkmem,mpi_enreg,mpw,nband_k,ncpgr,& 4237 & npw_k,npw1_k,nspinor,nsppol,occ_k,pawrhoijfermi,prtvol,& 4238 & rf_hamkq,rhoaug,rocceig,wtk_k) 4239 4240 !Arguments ------------------------------------ 4241 !scalars 4242 integer,intent(in) :: cplex,ibg,ibgq,icg,icgq,idir,ikpt 4243 integer,intent(in) :: ipert,isppol,kptopt,mcgq,mcprjq,mkmem,mpw,ncpgr 4244 integer,intent(in) :: mband_mem 4245 integer,intent(in) :: npw1_k,nspinor,nsppol,prtvol 4246 integer,intent(inout) :: nband_k,npw_k 4247 real(dp),intent(in) :: wtk_k 4248 type(MPI_type),intent(in) :: mpi_enreg 4249 type(datafiles_type),intent(in) :: dtfil 4250 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 4251 type(rf_hamiltonian_type),intent(inout) :: rf_hamkq 4252 !arrays 4253 real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol),cgq(2,mcgq) 4254 real(dp),intent(in) :: eig0_k(nband_k),occ_k(nband_k),rocceig(nband_k,nband_k) 4255 real(dp),intent(inout) :: rhoaug(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc) 4256 real(dp),intent(inout) :: eig1_k(2*nband_k**2) 4257 real(dp),intent(out) :: fe1fixed_k(nband_k) 4258 real(dp),intent(out) :: fe1norm_k(nband_k) 4259 !TODO distribute cprj over bands 4260 type(pawcprj_type),intent(in) :: cprj(gs_hamkq%natom,nspinor*mband_mem*mkmem*nsppol*gs_hamkq%usecprj) 4261 type(pawcprj_type),intent(in) :: cprjq(gs_hamkq%natom,mcprjq) 4262 type(pawrhoij_type),intent(inout) :: pawrhoijfermi(gs_hamkq%natom*gs_hamkq%usepaw) 4263 4264 !Local variables------------------------------- 4265 !scalars 4266 integer,parameter :: level=18 4267 integer :: berryopt,iband,ii,indx,iorder_cprj 4268 integer :: iband_me, nband_me 4269 integer :: ipw,me,nkpt_max,optlocal,optnl,opt_accrho,opt_corr 4270 integer :: opt_gvnlx1,sij_opt,tim_fourwf,tim_getgh1c,usevnl 4271 real(dp) :: dotr,lambda,wtband 4272 character(len=500) :: msg 4273 !arrays 4274 real(dp) :: dum_grad_berry(1,1),dum_gvnlx1(1,1),dum_gs1(1,1),tsec(2) 4275 real(dp),allocatable :: cwave0(:,:),cwaveq(:,:),gh1(:,:) 4276 type(pawcprj_type),allocatable :: cwaveprj0(:,:),cwaveprjq(:,:),cwaveprj_tmp(:,:) 4277 4278 ! ********************************************************************* 4279 4280 DBG_ENTER('COLL') 4281 4282 !Check arguments validity 4283 if (ipert>gs_hamkq%natom.and.ipert/=gs_hamkq%natom+3.and.ipert/=gs_hamkq%natom+4.and.ipert/=gs_hamkq%natom+5) then !SPr rfmagn deb 4284 ABI_BUG('wrong ipert argument !') 4285 end if 4286 if (cplex/=1) then 4287 ABI_BUG('wrong cplex/=1 argument !') 4288 end if 4289 4290 !Debugging statements 4291 if(prtvol==-level)then 4292 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,'dfpt_wfkfermi : enter' 4293 call wrtout(std_out,msg,'PERS') 4294 end if 4295 nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1 4296 4297 if(prtvol>2 .or. ikpt<=nkpt_max)then 4298 write(msg, '(a,a,i5,2x,a,3f9.5,2x,a)' ) ch10,& 4299 & ' Non-SCF iterations; k pt #',ikpt,'k=',gs_hamkq%kpt_k(:),' band residuals:' 4300 call wrtout(std_out,msg,'PERS') 4301 end if 4302 4303 !Retrieve parallelism data 4304 me=mpi_enreg%me_kpt 4305 !Initializations and allocations 4306 4307 ABI_MALLOC(gh1,(2,npw1_k*nspinor)) 4308 ABI_MALLOC(cwave0,(2,npw_k*nspinor)) 4309 ABI_MALLOC(cwaveq,(2,npw1_k*nspinor)) 4310 iorder_cprj=0 ; eig1_k(:)=zero 4311 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 4312 ABI_MALLOC(cwaveprj0,(gs_hamkq%natom,nspinor)) 4313 ABI_MALLOC(cwaveprjq,(gs_hamkq%natom,nspinor)) 4314 call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj) 4315 call pawcprj_alloc(cwaveprjq,0,gs_hamkq%dimcprj) 4316 else 4317 ABI_MALLOC(cwaveprj0,(0,0)) 4318 ABI_MALLOC(cwaveprjq,(0,0)) 4319 end if 4320 !Arguments of getgh1c routine (want only (NL+kin) frozen H(1)) 4321 berryopt=0;usevnl=0;sij_opt=-gs_hamkq%usepaw;tim_getgh1c=3 4322 optlocal=0;optnl=1;opt_gvnlx1=0 4323 if(ipert==gs_hamkq%natom+5) optnl=0; ! no 1st order NL in H(1), also no kin, but this will be taken into account later 4324 !if(ipert==gs_hamkq%natom+5) optlocal=0; ! 1st order LOCAL potential present 4325 4326 !Arguments of the dfpt_accrho routine 4327 tim_fourwf=5 ; opt_accrho=1 ; opt_corr=0 4328 !Null potentially unassigned output variables 4329 fe1fixed_k(:)=zero; fe1norm_k(:)=zero 4330 4331 !Read the npw and kg records of wf files 4332 call timab(139,1,tsec) 4333 4334 !Loop over bands 4335 iband_me = 0 4336 nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 4337 do iband=1,nband_k 4338 4339 ! Skip bands not treated by current proc 4340 if(mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me) cycle 4341 iband_me = iband_me + 1 4342 4343 ! Select occupied bands 4344 if(abs(occ_k(iband))>tol8.and.abs(rocceig(iband,iband))>tol8)then 4345 4346 wtband=rocceig(iband,iband)/occ_k(iband) 4347 ! Get ground-state wavefunctions at k 4348 do ipw=1,npw_k*nspinor 4349 cwave0(1,ipw)=cg(1,ipw+(iband_me-1)*npw_k*nspinor+icg) 4350 cwave0(2,ipw)=cg(2,ipw+(iband_me-1)*npw_k*nspinor+icg) 4351 end do 4352 4353 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 4354 ! Read PAW ground state projected WF (cprj) 4355 ! cprj is already distributed in band and k, just get the corresponding cprj in cwaveprj0 4356 call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,gs_hamkq%natom,iband_me,ibg,ikpt,iorder_cprj,& 4357 & isppol,mband_mem,mkmem,gs_hamkq%natom,1,nband_me,nspinor,nsppol,dtfil%unpaw,& 4358 !& mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,& 4359 & icpgr=idir,ncpgr=ncpgr) 4360 end if 4361 4362 ! Read ground-state wavefunctions at k+q 4363 indx=npw1_k*nspinor*(iband_me-1)+icgq 4364 cwaveq(:,1:npw_k*nspinor)=wtband*cgq(:,1+indx:npw_k*nspinor+indx) 4365 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 4366 ! Read PAW ground state projected WF (cprj) 4367 indx=nspinor*(iband_me-1)+ibgq 4368 ! TODO: cprj distributed -> iband_me 4369 !indx=nspinor*(iband-1)+ibgq 4370 call pawcprj_copy(cprjq(:,1+indx:nspinor+indx),cwaveprjq) 4371 call pawcprj_axpby(zero,wtband,cwaveprj_tmp,cwaveprjq) 4372 end if 4373 4374 ! Apply H^(1)-Esp.S^(1) to Psi^(0) (H(^1)=only (NL+kin) frozen part) 4375 lambda=eig0_k(iband) 4376 call getgh1c(berryopt,cwave0,cwaveprj0,gh1,dum_grad_berry,dum_gs1,gs_hamkq,dum_gvnlx1,& 4377 & idir,ipert,lambda,mpi_enreg,optlocal,optnl,opt_gvnlx1,rf_hamkq,sij_opt,& 4378 & tim_getgh1c,usevnl) 4379 ! Compute Eig1=<Psi^(0)|H^(1)-Eps.S^(1)|Psi(0)> 4380 call dotprod_g(dotr,lambda,gs_hamkq%istwf_k,npw_k*nspinor,1,cwave0,gh1,mpi_enreg%me_g0, & 4381 & mpi_enreg%comm_spinorfft) 4382 indx=2*iband-1+(iband-1)*2*nband_k 4383 eig1_k(indx)=dotr 4384 ! Compute the fixed contribution to the 1st-order Fermi energy 4385 fe1fixed_k(iband)=two*wtband*eig1_k(indx) 4386 fe1norm_k(iband) =two*wtband 4387 4388 ! Accumulate contribution to density and PAW occupation matrix 4389 4390 call dfpt_accrho(cplex,cwave0,cwaveq,cwaveq,cwaveprj0,cwaveprjq,dotr,& 4391 gs_hamkq,iband,0,0,isppol,kptopt,mpi_enreg,gs_hamkq%natom,nband_k,ncpgr,& 4392 npw_k,npw1_k,nspinor,occ_k,opt_accrho,pawrhoijfermi,rhoaug,tim_fourwf,& 4393 opt_corr,wtk_k) 4394 end if ! End of non-zero occupation and rocceig 4395 4396 end do ! End loop over bands 4397 4398 call timab(139,2,tsec) 4399 call timab(130,1,tsec) 4400 4401 ABI_FREE(cwave0) 4402 ABI_FREE(cwaveq) 4403 ABI_FREE(gh1) 4404 if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then 4405 call pawcprj_free(cwaveprj0) 4406 call pawcprj_free(cwaveprjq) 4407 end if 4408 ABI_FREE(cwaveprj0) 4409 ABI_FREE(cwaveprjq) 4410 4411 !Structured debugging : if prtvol=-level, stop here. 4412 if(prtvol==-level)then 4413 write(msg,'(a,a1,a,i2,a)')' fermie3 : exit prtvol=-',level,', debugging mode => stop ' 4414 ABI_ERROR(msg) 4415 end if 4416 4417 call timab(130,2,tsec) 4418 4419 DBG_EXIT('COLL') 4420 4421 end subroutine dfpt_wfkfermi
ABINIT/m_dfpt_scfcv [ Modules ]
NAME
m_dfpt_scfcv
FUNCTION
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (XG, DRH, MB, XW, MT, SPr, XW, MV, MM, AR) 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
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_dfpt_scfcv 22 23 use defs_basis 24 use m_ab7_mixing 25 use m_efield 26 use m_errors 27 use m_dtset 28 use m_abicore 29 use m_wfk 30 use m_wffile 31 use m_xmpi 32 use m_nctk 33 use m_hdr 34 use m_dtfil 35 use m_hamiltonian 36 #ifdef HAVE_NETCDF 37 use netcdf 38 #endif 39 40 use defs_datatypes, only : pseudopotential_type 41 use defs_abitypes, only : MPI_type 42 use m_cgtools, only : mean_fftr, overlap_g, dotprod_vn, dotprod_vn, dotprod_g 43 use m_fstrings, only : int2char4, sjoin 44 use m_geometry, only : metric, stresssym 45 use m_time, only : abi_wtime, sec2str, timab 46 use m_io_tools, only : open_file, file_exists, get_unit, iomode_from_fname 47 use m_exit, only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in 48 use m_mpinfo 49 use m_kg, only : getcut, mkkin, kpgstr, mkkpg 50 use m_fft, only : fftpac, fourdp 51 use m_symtk, only : mati3inv 52 use m_dynmat, only : dfpt_sygra 53 use m_occ, only : occeig 54 use m_paw_mkrho, only : pawmkrho 55 use m_mkffnl, only : mkffnl 56 use m_getgh1c, only : getgh1c 57 use m_dfpt_mkrho, only : dfpt_accrho 58 use m_nonlop, only : nonlop 59 use m_ioarr, only : ioarr, fftdatar_write_from_hdr, fort_denpot_skip 60 use m_pawang, only : pawang_type 61 use m_pawrad, only : pawrad_type 62 use m_pawtab, only : pawtab_type 63 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_reset_flags 64 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_reset_flags 65 use m_pawfgrtab,only : pawfgrtab_type 66 use m_pawrhoij, only : pawrhoij_type, pawrhoij_init_unpacked, pawrhoij_gather, pawrhoij_filter, & 67 pawrhoij_alloc, pawrhoij_free, pawrhoij_nullify, & 68 pawrhoij_free_unpacked, pawrhoij_mpisum_unpacked, pawrhoij_inquire_dim 69 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_copy, pawcprj_axpby, pawcprj_free, pawcprj_getdim 70 use m_pawdij, only : pawdij, pawdijfr, symdij 71 use m_pawfgr, only : pawfgr_type 72 use m_paw_denpot, only : pawdenpot 73 use m_paw_dfpt, only : pawdfptenergy 74 use m_paw_nhat, only : pawmknhat,pawnhatfr 75 use m_rf2, only : rf2_getidirs 76 use m_dens, only : calcdenmagsph, prtdenmagsph 77 use m_dfpt_fef, only : dfptff_initberry, qmatrix, dfptff_edie, dfptff_ebp, dfptff_die, dfptff_bec 78 use m_dfpt_vtorho, only : dfpt_vtorho 79 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 80 use m_common, only : scprqt 81 use m_prcref, only : moddiel 82 use m_dfpt_rhotov, only : dfpt_rhotov 83 use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll 84 use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr 85 use m_mklocl, only : dfpt_vlocal, vlocalstr 86 use m_dfpt_nstwf, only : dfpt_nstpaw, dfpt_nstwf 87 use m_mkcore, only : dfpt_mkcore 88 use m_spacepar, only : hartrestr, make_vectornd, symrhg 89 90 implicit none 91 92 private
ABINIT/newfermie1 [ Functions ]
NAME
newfermie1
FUNCTION
This routine computes the derivative of the fermi energy wrt the active perturbation for use in evaluating the edocc term and active subspace contribution to the first-order wavefunctions in the case of metals. This is presently used only for the strain and magnetic field perturbations, and only for Q = 0.
INPUTS
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX fe1fixed=fixed contribution to the first-order Fermi energy ipert=index of perturbation istep=index of the number of steps in the routine scfcv ixc= choice of exchange-correlation scheme mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor natom=number of atoms nfft=(effective) number of FFT grid points (for this processor) nfftot= total number of FFT grid points nhatfermi(nfft,nspden)=fermi-level compensation charge density (PAW only) nspden=number of spin-density components ntypat=number of atom types occopt=option for occupancies paw_an(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh paw_ij1(natom) <type(paw_ij_type)>=(1st-order) paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawnzlm=-- PAW only -- option for the computation of non-zero lm moments of the on-sites densities pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawrhoij1(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies pawrhoijfermi(natom) <type(pawrhoij_type)>=paw rhoij occupancies at Fermi level pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) xclevel= XC functional level prtvol=control print volume and debugging output rhorfermi(nfft,nspden)=fermi-level electronic density ucvol=unit cell volume in bohr**3 usepaw=1 if PAW is activated usexcnhat= -PAW only- flag controling use of compensation density in Vxc vtrial1(cplex*nfft,nspden)=1-st order potential vxc1(cplex*nfft,nspden)=1-st order XC potential
OUTPUT
(see side effects)
SIDE EFFECTS
fermie1=derivative of fermi energy wrt perturbation at input : old value at output : updated value
SOURCE
1743 subroutine newfermie1(cplex,fermie1,fe1fixed,ipert,istep,ixc,my_natom,natom,nfft,nfftot,& 1744 & nhatfermi,nspden,ntypat,occopt,paw_an,paw_an1,paw_ij1,pawang,pawnzlm,pawrad,& 1745 & pawrhoij1,pawrhoijfermi,pawtab,pawxcdev,prtvol,rhorfermi,& 1746 & ucvol,usepaw,usexcnhat,vtrial1,vxc1,xclevel,& 1747 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 1748 1749 !Arguments ------------------------------- 1750 !scalars 1751 integer,intent(in) :: cplex,ipert,istep,ixc,my_natom,natom,nfft,nfftot,nspden,ntypat 1752 integer,intent(in) :: occopt,pawnzlm,pawxcdev,prtvol,usepaw,usexcnhat,xclevel 1753 integer,optional,intent(in) :: comm_atom 1754 real(dp),intent(in) :: fe1fixed,ucvol 1755 real(dp),intent(inout) :: fermie1 1756 type(pawang_type),intent(in) :: pawang 1757 !arrays 1758 integer,optional,target,intent(in) :: mpi_atmtab(:) 1759 real(dp),intent(in) :: rhorfermi(nfft,nspden),vtrial1(cplex*nfft,nspden) 1760 real(dp),intent(in) :: nhatfermi(:,:),vxc1(:,:) 1761 type(paw_an_type),intent(in) :: paw_an(my_natom*usepaw) 1762 type(paw_an_type),intent(inout) :: paw_an1(my_natom*usepaw) 1763 type(paw_ij_type),intent(inout) :: paw_ij1(my_natom*usepaw) 1764 type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw) 1765 type(pawrhoij_type),intent(in) :: pawrhoij1(my_natom*usepaw),pawrhoijfermi(my_natom*usepaw) 1766 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 1767 1768 !Local variables------------------------------- 1769 !scalars 1770 integer :: ipert0,my_comm_atom,nzlmopt,nzlmopt_fermi,option,pawprtvol 1771 logical :: my_atmtab_allocated,paral_atom 1772 real(dp) :: doti,fe1_scf,fe1_tmp,fermie1_new,fermie1rs 1773 character(len=500) :: msg 1774 !arrays 1775 integer, pointer :: my_atmtab(:) 1776 real(dp) :: fe1_paw(2) 1777 real(dp), allocatable :: rhor_nonhat(:,:),vtrial1_novxc(:,:) 1778 1779 ! ********************************************************************* 1780 1781 !Tests 1782 if (cplex==2) then 1783 msg='Not compatible with cplex=2!' 1784 ABI_BUG(msg) 1785 end if 1786 if (usepaw==1.and.usexcnhat==0.and.(size(nhatfermi)<=0.or.size(vxc1)<=0)) then 1787 msg='Should have nhatfermi and vxc1 allocated with usexcnhat=0!' 1788 ABI_BUG(msg) 1789 end if 1790 1791 !Set up parallelism over atoms 1792 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 1793 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 1794 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 1795 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 1796 1797 if(occopt>=3 .and. occopt <=8) then 1798 1799 ! The product of the current trial potential and the so-called Fermi level 1800 ! density is integrated to give the local potential contributions to the 1801 ! first-order Fermi level. 1802 option=1 1803 if (usepaw==1.and.usexcnhat==0) then 1804 ABI_MALLOC(rhor_nonhat,(nfft,nspden)) 1805 ABI_MALLOC(vtrial1_novxc,(nfft,nspden)) 1806 rhor_nonhat(1:nfft,1:nspden)=rhorfermi(1:nfft,1:nspden)-nhatfermi(1:nfft,1:nspden) 1807 vtrial1_novxc(1:nfft,1:nspden)=vtrial1(1:nfft,1:nspden)-vxc1(1:nfft,1:nspden) 1808 call dotprod_vn(cplex,rhor_nonhat,fe1_scf,doti,nfft,nfftot,& 1809 & nspden,option,vtrial1,ucvol) 1810 call dotprod_vn(cplex,nhatfermi,fe1_tmp,doti,nfft,nfftot,& 1811 & nspden,option,vtrial1_novxc,ucvol) 1812 fe1_scf=fe1_scf+fe1_tmp 1813 ABI_FREE(rhor_nonhat) 1814 ABI_FREE(vtrial1_novxc) 1815 else 1816 call dotprod_vn(cplex,rhorfermi,fe1_scf,doti,nfft,nfftot,& 1817 & nspden,option,vtrial1,ucvol) 1818 end if 1819 1820 fe1_paw(:)=zero 1821 ! PAW on-site contribution (use Fermi level occupation matrix) 1822 if (usepaw==1) then 1823 ipert0=0;pawprtvol=0 1824 nzlmopt=0;if (istep>1) nzlmopt=pawnzlm 1825 if (istep==1.and.pawnzlm>0) nzlmopt=-1 1826 nzlmopt_fermi=0;if (pawnzlm>0) nzlmopt_fermi=-1 1827 call pawdfptenergy(fe1_paw,ipert,ipert0,ixc,my_natom,natom,ntypat,nzlmopt,& 1828 & nzlmopt_fermi,paw_an,paw_an1,paw_ij1,pawang,pawprtvol,pawrad,& 1829 & pawrhoij1,pawrhoijfermi,pawtab,pawxcdev,xclevel,& 1830 & mpi_atmtab=my_atmtab, comm_atom=my_comm_atom) 1831 end if 1832 1833 ! The fixed contributions consisting of non-local potential and kinetic terms 1834 ! are added 1835 fermie1_new=fe1fixed+fe1_scf+fe1_paw(1) 1836 fermie1rs=(fermie1-fermie1_new)**2 1837 fermie1=fermie1_new 1838 1839 if(prtvol>=10)then 1840 write(msg, '(a,i5,2es18.8)' ) ' fermie1, residual squared',istep,fermie1,fermie1rs 1841 call wrtout(std_out,msg,'COLL') 1842 end if 1843 1844 else 1845 fermie1=zero 1846 end if 1847 1848 !Destroy atom table used for parallelism 1849 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 1850 1851 end subroutine newfermie1