TABLE OF CONTENTS
- ABINIT/dfpt_dyfro
- ABINIT/dfpt_dyout
- ABINIT/dfpt_dyxc1
- ABINIT/dfpt_gatherdy
- ABINIT/m_respfn_driver
- m_respfn_driver/respfn
- m_respfn_driver/wrtloctens
ABINIT/dfpt_dyfro [ Functions ]
NAME
dfpt_dyfro
FUNCTION
Compute the different parts of the frozen-wavefunction part of the dynamical matrix, except the non-local one, computed previously. Also (when installed) symmetrize the different part and their sum.
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx dyfr_cplex=1 if dyfrnl is real, 2 if it is complex dyfr_nondiag=1 if dyfrnl and dyfrwf are non diagonal with respect to atoms; 0 otherwise gmet(3,3)=reciprocal space metric (bohr^-2) gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1) gsqcut=cutoff on G^2 based on ecut indsym(4,nsym,natom)=index showing transformation of atom labels under symmetry operations (computed in symatm) mgfft=maximum size of 1D FFTs mpi_enreg=information about MPI parallelization mqgrid=dimensioned number of q grid points for local psp spline natom=number of atoms in unit cell nattyp(ntypat)=number of atoms of each type 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 nspden=number of spin-density components nsym=number of symmetries in space group ntypat=number of types of atoms n1xccc=dimension of xccc1d ; 0 if no XC core correction is used n3xccc=dimension of the xccc3d array (0 or nfft). psps <type(pseudopotential_type)>=variables related to pseudopotentials pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information qgrid(mqgrid)=q point array for local psp spline fits qphon(3)=wavevector of the phonon rhog(2,nfft)=electron density in G space rprimd(3,3)=dimensional primitive translation vectors (bohr) symq(4,2,nsym)=1 if symmetry preserves present qpoint. From littlegroup_q symrec(3,3,nsym)=symmetries in reciprocal space typat(natom)=integer type for each atom in cell ucvol=unit cell volume (bohr**3). usepaw= 0 for non paw calculation; =1 for paw calculation vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom. vxc(nfft,nspden)=exchange-correlation potential (hartree) in real space--only used when n1xccc/=0 xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives, for each type of atom, from psp xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xred(3,natom)=reduced coordinates for atoms in unit cell
OUTPUT
dyfrlo(3,3,natom)=frozen wavefunctions part of the dynamical matrix
(local only)
dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
frozen wavefunctions part of the dynamical matrix
(local + non-local)
If NCPP, it depends on one atom
If PAW, it depends on two atoms
dyfrxc(3,3,natom)=frozen wavefunctions part of the dynamical matrix
(non-linear xc core correction)
SIDE EFFECTS
Input/Output
dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=
frozen wavefunctions part of the dynamical matrix
(non-local only)
If NCPP, it depends on one atom
If PAW, it depends on two atoms
SOURCE
3939 subroutine dfpt_dyfro(atindx1,dyfrnl,dyfrlo,dyfrwf,dyfrxc,dyfr_cplex,dyfr_nondiag,& 3940 & dtset,gmet,gprimd,gsqcut,indsym,mgfft,mpi_enreg,mqgrid,natom,nattyp,& 3941 & nfft,ngfft,nspden,nsym,ntypat,n1xccc,n3xccc,psps,pawtab,ph1d,qgrid,& 3942 & qphon,rhog,rprimd,symq,symrec,typat,ucvol,usepaw,vlspl,vxc,& 3943 & xcccrc,xccc1d,xccc3d,xred) 3944 3945 !Arguments ------------------------------------ 3946 !scalars 3947 integer,intent(in) :: dyfr_cplex,dyfr_nondiag,mgfft,mqgrid,n1xccc,n3xccc,natom,nfft,nspden 3948 integer,intent(in) :: nsym,ntypat,usepaw 3949 real(dp),intent(in) :: gsqcut,ucvol 3950 type(dataset_type),intent(in) :: dtset 3951 type(pseudopotential_type),intent(in) :: psps 3952 type(MPI_type),intent(in) :: mpi_enreg 3953 !arrays 3954 integer,intent(in) :: atindx1(natom),indsym(4,nsym,natom),nattyp(ntypat) 3955 integer,intent(in) :: ngfft(18),symq(4,2,nsym),symrec(3,3,nsym),typat(natom) 3956 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom) 3957 real(dp),intent(in) :: qgrid(mqgrid),qphon(3),rhog(2,nfft) 3958 real(dp),intent(in) :: vlspl(mqgrid,2,ntypat),vxc(nfft,nspden) 3959 real(dp),intent(in) :: xccc1d(n1xccc,6,ntypat),xcccrc(ntypat),xred(3,natom) 3960 real(dp),intent(inout) :: rprimd(3,3) 3961 real(dp),intent(inout) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag),xccc3d(n3xccc) 3962 real(dp),intent(out) :: dyfrlo(3,3,natom),dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag) 3963 real(dp),intent(inout) :: dyfrxc(3,3,natom) !vz_i 3964 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 3965 3966 !Local variables------------------------------- 3967 !scalars 3968 logical, parameter :: do_final_sym=.true. 3969 integer :: iatom,jatom,n1,n2,n3,optatm,optdyfr,opteltfr,optgr,option 3970 integer :: optn,optn2,optstr,optv 3971 real(dp) :: eei 3972 !arrays 3973 integer :: qprtrb(3) 3974 real(dp) :: dummy6(6),dum_strn(6),dum_strv(6) 3975 real(dp) :: tsec(2),vprtrb(2) 3976 real(dp) :: dum_atmrho(0),dum_atmvloc(0),dum_gauss(0),dum_grn(0),dum_grv(0),dum_eltfrxc(0) 3977 real(dp),allocatable :: dyfrlo_tmp1(:,:,:),dyfrlo_tmp2(:,:,:,:,:),dyfrsym_tmp(:,:,:,:,:) 3978 real(dp),allocatable :: gr_dum(:,:),v_dum(:),vxctotg(:,:) 3979 3980 ! ************************************************************************* 3981 3982 if(nspden==4 .and. usepaw==1)then 3983 ABI_WARNING('dfpt_dyfro : DFPT with nspden=4 works at the moment just for norm-conserving psp! (no paw support yet with nspden=4)') 3984 end if 3985 3986 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 3987 3988 if (usepaw==1 .or. psps%nc_xccc_gspace==1) then 3989 3990 ! PAW or NC with nc_xccc_gspace: compute local psp and core charge contribs together 3991 ! in reciprocal space 3992 ! ----------------------------------------------------------------------- 3993 call timab(563,1,tsec) 3994 if (n3xccc>0) then 3995 ABI_MALLOC(v_dum,(nfft)) 3996 ABI_MALLOC(vxctotg,(2,nfft)) 3997 v_dum(:)=vxc(:,1);if (nspden>=2) v_dum(:)=0.5_dp*(v_dum(:)+vxc(:,2)) 3998 call fourdp(1,vxctotg,v_dum,-1,mpi_enreg,nfft,1,ngfft,0) 3999 call zerosym(vxctotg,2,ngfft(1),ngfft(2),ngfft(3),& 4000 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 4001 ABI_FREE(v_dum) 4002 else 4003 ABI_MALLOC(vxctotg,(0,0)) 4004 end if 4005 optatm=0;optdyfr=1;optgr=0;optstr=0;optv=1;optn=n3xccc/nfft;optn2=1;opteltfr=0 4006 call atm2fft(atindx1,dum_atmrho,dum_atmvloc,dyfrxc,dyfrlo,dum_eltfrxc,& 4007 & dum_gauss,gmet,gprimd,dum_grn,dum_grv,gsqcut,mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,& 4008 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,qgrid,qprtrb,dtset%rcut,rhog,& 4009 & rprimd,dum_strn,dum_strv,ucvol,usepaw,vxctotg,vxctotg,vxctotg,vprtrb,vlspl) 4010 ABI_FREE(vxctotg) 4011 if (n3xccc==0) dyfrxc=zero 4012 else 4013 4014 ! Norm-conserving: compute local psp contribution in reciprocal space 4015 ! and core charge contribution in real space 4016 ! ----------------------------------------------------------------------- 4017 option=4 4018 ABI_MALLOC(dyfrlo_tmp1,(3,3,natom)) 4019 ABI_MALLOC(gr_dum,(3,natom)) 4020 ABI_MALLOC(v_dum,(nfft)) 4021 call mklocl_recipspace(dyfrlo_tmp1,eei,gmet,gprimd,& 4022 & gr_dum,gsqcut,dtset%icutcoul,dummy6,mgfft,mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,dtset%nkpt,& 4023 & ntypat,option,ph1d,qgrid,qprtrb,dtset%rcut,rhog,rprimd,ucvol,dtset%vcutgeo,vlspl,vprtrb,v_dum) 4024 do iatom=1,natom 4025 ! Reestablish correct order of atoms 4026 dyfrlo(1:3,1:3,atindx1(iatom))=dyfrlo_tmp1(1:3,1:3,iatom) 4027 end do 4028 ABI_FREE(dyfrlo_tmp1) 4029 ABI_FREE(v_dum) 4030 if(n1xccc/=0)then 4031 call mkcore(dummy6,dyfrxc,gr_dum,mpi_enreg,natom,nfft,nspden,ntypat,& 4032 & n1,n1xccc,n2,n3,option,rprimd,typat,ucvol,vxc,xcccrc,xccc1d,xccc3d,xred) 4033 end if 4034 ABI_FREE(gr_dum) 4035 end if 4036 4037 !Symmetrize dynamical matrix explicitly for given space group: 4038 4039 !Symmetrize local part of the dynamical matrix dyfrlo: 4040 ABI_MALLOC(dyfrsym_tmp,(1,3,3,natom,1)) 4041 ABI_MALLOC(dyfrlo_tmp2,(1,3,3,natom,1)) 4042 dyfrsym_tmp(1,:,:,:,1)=dyfrlo(:,:,:) 4043 call dfpt_sydy(1,dyfrsym_tmp,indsym,natom,0,nsym,qphon,dyfrlo_tmp2,symq,symrec) 4044 dyfrlo(:,:,:)=dyfrlo_tmp2(1,:,:,:,1) 4045 if (do_final_sym) then 4046 dyfrsym_tmp(1,:,:,:,1)=dyfrxc(:,:,:) 4047 call dfpt_sydy(1,dyfrsym_tmp,indsym,natom,0,nsym,qphon,dyfrlo_tmp2,symq,symrec) 4048 dyfrxc(:,:,:)=dyfrlo_tmp2(1,:,:,:,1) 4049 end if 4050 ABI_FREE(dyfrsym_tmp) 4051 ABI_FREE(dyfrlo_tmp2) 4052 4053 !Symmetrize nonlocal part of the dynamical matrix dyfrnl: 4054 !atindx1 is used to reestablish the correct order of atoms 4055 if (dyfr_nondiag==0) then 4056 ABI_MALLOC(dyfrsym_tmp,(dyfr_cplex,3,3,natom,1)) 4057 do iatom=1,natom 4058 dyfrsym_tmp(:,:,:,atindx1(iatom),1)=dyfrnl(:,:,:,iatom,1) 4059 end do 4060 else 4061 ABI_MALLOC(dyfrsym_tmp,(dyfr_cplex,3,3,natom,natom)) 4062 do jatom=1,natom 4063 do iatom=1,natom 4064 dyfrsym_tmp(:,:,:,atindx1(iatom),atindx1(jatom))=dyfrnl(:,:,:,iatom,jatom) 4065 end do 4066 end do 4067 end if 4068 call dfpt_sydy(dyfr_cplex,dyfrsym_tmp,indsym,natom,dyfr_nondiag,nsym,qphon,dyfrnl,symq,symrec) 4069 ABI_FREE(dyfrsym_tmp) 4070 4071 !Collect local, nl xc core, and non-local part 4072 !of the frozen wf dynamical matrix. 4073 dyfrwf(:,:,:,:,:)=dyfrnl(:,:,:,:,:) 4074 if (dyfr_nondiag==0) then 4075 dyfrwf(1,:,:,:,1)=dyfrwf(1,:,:,:,1)+dyfrlo(:,:,:)+dyfrxc(:,:,:) 4076 else 4077 do iatom=1,natom 4078 dyfrwf(1,:,:,iatom,iatom)=dyfrwf(1,:,:,iatom,iatom)+dyfrlo(:,:,iatom)+dyfrxc(:,:,iatom) 4079 end do 4080 end if 4081 4082 end subroutine dfpt_dyfro
ABINIT/dfpt_dyout [ Functions ]
NAME
dfpt_dyout
FUNCTION
Output of all quantities related to the 2nd-order matrix: Ewald part, local and non-local frozen wf part, core contributions, local and non-local variational part, 2nd-order matrix itself, and, for the phonon part, eigenfrequencies, in Hartree, meV and cm-1. Also output unformatted 2nd-order matrix for later use in the Brillouin-zone interpolation
INPUTS
becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only) blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical matrix has been calculated ; 0 otherwise ) carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian 2DTE matrix has been calculated correctly ; 0 otherwise ) ddkfil(3)=components are 1 if corresponding d/dk file exists, otherwise 0 (in what follows, DYMX means dynamical matrix, and D2MX means 2nd-order matrix) dyew(2,3,natom,3,natom)=Ewald part of the DYMX dyfrlo(3,3,natom)=frozen wf local part of the DYMX dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wf nonloc part of the DYMX dyfrx1(2,3,natom,3,natom)=frozen wf nonlin. xc core corr.(1) part of the DYMX dyfrx2(3,3,natom)=frozen wf nonlin. xc core corr.(2) part of the DYMX dyfr_cplex=1 if dyfrnl is real, 2 if it is complex dyfr_nondiag=1 if dyfrnl is non diagonal with respect to atoms; 0 otherwise dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix d2cart(2,3,mpert,3,mpert)=D2MX in cartesian coordinates d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)= band by band decomposition of Born effective charges (calculated from phonon-type perturbation) in cartesian coordinates d2eig0(2,3,mpert,3,mpert)=0-order eigen. station. part of the D2MX d2k0(2,3,mpert,3,mpert)=0-order kinet. station. part of the D2MX d2lo(2,3,mpert,3,mpert)=nonstation. local part of the D2MX d2loc0(2,3,mpert,3,mpert)=0-order loc station. part of the D2MX d2matr(2,3,mpert,3,mpert)=D2MX in non-cartesian coordinates d2nl(2,3,mpert,3,mpert)=nonstation. nonloc part of the D2MX d2nl0(2,3,mpert,3,mpert)=0-order nonloc station. part of the D2MX d2nl1(2,3,mpert,3,mpert)=1-order nonloc station. part of the D2MX d2ovl(2,mpert,3,mpert*usepaw)=1st-order change of WF overlap contributions to the 2DTEs (PAW) d2vn(2,3,mpert,3,mpert)=potential*dens station. part of the D2MX and without masses included) eltcore(6,6)=core contribution to the elastic tensor elteew(6+3*natom,6)=Ewald contribution to the elastic tsenor eltfrhar(6,6)=hartree contribution to the elastic tensor eltfrkin(6,6)=kinetic contribution to the elastic tensor eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor has_full_piezo=the full calculation of the piezoelectric tensor from electric field perturbation is only available if nsym=1 (strain perturbation is not symmetrized) has_allddk= True if all ddk file are present on disk, so the effective charge or piezzo electric tensor are correctly computed (PAW ONLY) iout=unit number for long write-up mband=maximum number of bands mpert =maximum number of ipert natom=number of atoms ntypat=number of atom types outd2=option for the output of the 2nd-order matrix : if outd2=1, non-stationary part if outd2=2, stationary part. pawbec= flag for the computation of frozen part of Born Effective Charges (PAW only) pawpiezo= flag for the computation of frozen part of Piezoelectric tensor (PAW only) prtbbb=if 1, print the band-by-band decomposition prtvol=print volume qzero=1 if zero phonon wavevector rfdir(3)=defines the directions for the perturbations rfpert(mpert)=defines the perturbations rfphon=if 1, there are phonon perturbations rfstrs=if 1,2,3 there are strain perturbations typat(natom)=integer label of each type of atom (1,2,...) usepaw=1 if PAW, 0 otherwise usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use zion(ntypat)=charge corresponding to the atom type
SIDE EFFECTS
d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)
NOTES
This routine is called only by the processor me==0 . In consequence, no use of message and wrtout routine.
SOURCE
2165 subroutine dfpt_dyout(becfrnl,berryopt,blkflg,carflg,ddkfil,dyew,dyfrlo,dyfrnl,& 2166 & dyfrx1,dyfrx2,dyfr_cplex,dyfr_nondiag,dyvdw,d2cart,d2cart_bbb,& 2167 & d2eig0,d2k0,d2lo,d2loc0,d2matr,d2nl,d2nl0,d2nl1,d2ovl,d2vn,& 2168 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 2169 & has_full_piezo,has_allddk,iout,mband,mpert,natom,ntypat,& 2170 & outd2,pawbec,pawpiezo,piezofrnl,prtbbb,prtvol,qzero,typat,rfdir,& 2171 & rfpert,rfphon,rfstrs,usepaw,usevdw,zion) 2172 2173 !Arguments ------------------------------- 2174 !scalars 2175 integer,intent(in) :: berryopt,dyfr_cplex,dyfr_nondiag,iout,mband,mpert 2176 integer,intent(in) :: natom,ntypat,outd2,pawbec,pawpiezo,prtbbb,prtvol,qzero 2177 integer, intent(in) :: rfphon,rfstrs,usepaw,usevdw 2178 !arrays 2179 integer,intent(in) :: blkflg(3,mpert,3,mpert),carflg(3,mpert,3,mpert) 2180 integer,intent(in) :: ddkfil(3),rfdir(3),rfpert(mpert),typat(natom) 2181 real(dp),intent(in) :: becfrnl(3,natom,3*pawbec) 2182 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert),d2eig0(2,3,mpert,3,mpert) 2183 real(dp),intent(in) :: d2k0(2,3,mpert,3,mpert),d2lo(2,3,mpert,3,mpert) 2184 real(dp),intent(in) :: d2loc0(2,3,mpert,3,mpert),d2matr(2,3,mpert,3,mpert) 2185 real(dp),intent(in) :: d2nl(2,3,mpert,3,mpert),d2nl0(2,3,mpert,3,mpert) 2186 real(dp),intent(in) :: d2nl1(2,3,mpert,3,mpert),d2ovl(2,3,mpert,3,mpert*usepaw) 2187 real(dp),intent(in) :: d2vn(2,3,mpert,3,mpert) 2188 real(dp),intent(in) :: dyew(2,3,natom,3,natom),dyfrlo(3,3,natom) 2189 real(dp),intent(in) :: dyfrnl(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag) 2190 real(dp),intent(in) :: dyfrx1(2,3,natom,3,natom),dyfrx2(3,3,natom) 2191 real(dp),intent(in) :: dyvdw(2,3,natom,3,natom*usevdw) 2192 real(dp),intent(in) :: eltcore(6,6),elteew(6+3*natom,6) 2193 real(dp),intent(in) :: eltfrhar(6,6),eltfrkin(6,6),eltfrloc(6+3*natom,6) 2194 real(dp),intent(in) :: eltfrnl(6+3*natom,6),eltfrxc(6+3*natom,6) 2195 real(dp),intent(in) :: eltvdw(6+3*natom,6*usevdw),piezofrnl(6,3*pawpiezo) 2196 real(dp),intent(in) :: zion(ntypat) 2197 real(dp),intent(inout) :: d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb) 2198 logical,intent(in) :: has_allddk,has_full_piezo 2199 2200 !Local variables ------------------------- 2201 !scalars 2202 integer :: iband,idir1,idir2,ii,ipert1,ipert2,jj,nelmts,nline 2203 real(dp) :: zi,zr 2204 !arrays 2205 real(dp) :: delta(3,3) 2206 2207 ! ********************************************************************* 2208 2209 ! GA: As much as I can tell, the option outd2 is always set to 1. 2210 ! This variable should be removed 2211 2212 2213 !Long print : includes detail of every part of the 2nd-order energy 2214 if(prtvol>=10)then 2215 2216 ! In case of phonon 2217 if (rfphon==1)then 2218 2219 ! write the Ewald part of the dynamical matrix 2220 write(iout,*)' ' 2221 write(iout,*)' Ewald part of the dynamical matrix' 2222 write(iout,*)' j1 j2 matrix element' 2223 write(iout,*)' dir pert dir pert real part imaginary part' 2224 do ipert1=1,natom 2225 do idir1=1,3 2226 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2227 & .or. outd2==1 )then 2228 write(iout,*)' ' 2229 do ipert2=1,natom 2230 do idir2=1,3 2231 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2232 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2233 & dyew(1,idir1,ipert1,idir2,ipert2),& 2234 & dyew(2,idir1,ipert1,idir2,ipert2) 2235 end if 2236 end do 2237 end do 2238 end if 2239 end do 2240 end do 2241 2242 ! Now the local frozen wf part 2243 write(iout,*)' ' 2244 write(iout,*)' Frozen wf local part of the dynamical matrix' 2245 write(iout,*)' j1 j2 matrix element' 2246 write(iout,*)' dir pert dir pert real part imaginary part' 2247 do ipert1=1,natom 2248 do idir1=1,3 2249 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2250 & .or. outd2==1 )then 2251 write(iout,*)' ' 2252 do ipert2=1,natom 2253 do idir2=1,3 2254 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2255 if(ipert1==ipert2)then 2256 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2257 & dyfrlo(idir1,idir2,ipert2),zero 2258 else 2259 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2260 & zero,zero 2261 end if 2262 end if 2263 end do 2264 end do 2265 end if 2266 end do 2267 end do 2268 2269 ! Now the nonlo frozen wf part 2270 write(iout,*)' ' 2271 write(iout,*)' Frozen wf non-local part of the dynamical matrix' 2272 write(iout,*)' j1 j2 matrix element' 2273 write(iout,*)' dir pert dir pert real part imaginary part' 2274 do ipert1=1,natom 2275 do idir1=1,3 2276 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2277 & .or. outd2==1 )then 2278 write(iout,*)' ' 2279 do ipert2=1,natom 2280 do idir2=1,3 2281 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2282 if(ipert1==ipert2.or.dyfr_nondiag==1)then 2283 if (dyfr_cplex==1) then 2284 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2285 & dyfrnl(1,idir1,idir2,ipert1,1+(ipert2-1)*dyfr_nondiag),zero 2286 else 2287 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2288 & dyfrnl(:,idir1,idir2,ipert1,1+(ipert2-1)*dyfr_nondiag) 2289 end if 2290 else 2291 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2292 & zero,zero 2293 end if 2294 end if 2295 end do 2296 end do 2297 end if 2298 end do 2299 end do 2300 2301 ! Now the nonlinear xc core correction(1) part 2302 write(iout,*)' ' 2303 write(iout,*)' Frozen wf xc core (1) part',& 2304 & ' of the dynamical matrix' 2305 write(iout,*)' j1 j2 matrix element' 2306 write(iout,*)' dir pert dir pert real part imaginary part' 2307 do ipert1=1,natom 2308 do idir1=1,3 2309 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2310 & .or. outd2==1 )then 2311 write(iout,*)' ' 2312 do ipert2=1,natom 2313 do idir2=1,3 2314 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2315 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2316 & dyfrx1(1,idir1,ipert1,idir2,ipert2),& 2317 & dyfrx1(2,idir1,ipert1,idir2,ipert2) 2318 end if 2319 end do 2320 end do 2321 end if 2322 end do 2323 end do 2324 2325 ! Now the nonlinear xc core correction(2) part 2326 write(iout,*)' ' 2327 write(iout,*)' Frozen wf xc core (2) part',& 2328 & ' of the dynamical matrix' 2329 write(iout,*)' j1 j2 matrix element' 2330 write(iout,*)' dir pert dir pert real part imaginary part' 2331 do ipert1=1,natom 2332 do idir1=1,3 2333 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2334 & .or. outd2==1 )then 2335 write(iout,*)' ' 2336 do ipert2=1,natom 2337 do idir2=1,3 2338 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2339 if(ipert1==ipert2)then 2340 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2341 & dyfrx2(idir1,idir2,ipert2),zero 2342 else 2343 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2344 & zero,zero 2345 end if 2346 end if 2347 end do 2348 end do 2349 end if 2350 end do 2351 end do 2352 2353 ! Now the DFT-D vdw part of the dynamical matrix 2354 if (usevdw==1) then 2355 write(iout,*)' ' 2356 write(iout,*)' DFT-D van der Waals part of the dynamical matrix' 2357 write(iout,*)' j1 j2 matrix element' 2358 write(iout,*)' dir pert dir pert real part imaginary part' 2359 do ipert1=1,natom 2360 do idir1=1,3 2361 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2362 & .or. outd2==1 )then 2363 write(iout,*)' ' 2364 do ipert2=1,natom 2365 do idir2=1,3 2366 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2367 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2368 & dyvdw(1,idir1,ipert1,idir2,ipert2),& 2369 & dyvdw(2,idir1,ipert1,idir2,ipert2) 2370 end if 2371 end do 2372 end do 2373 end if 2374 end do 2375 end do 2376 end if 2377 2378 ! End of the phonon condition 2379 end if 2380 2381 ! In case of atom. strain/electric field perturbation (piezoelectric tensor) 2382 if (pawpiezo==1.and.(rfpert(natom+2)==1.or.rfstrs/=0).and.outd2==1)then 2383 write(iout,*)' ' 2384 write(iout,*)' Frozen wf part of the piezoelectric tensor' 2385 write(iout,*)' j1 j2 matrix element' 2386 write(iout,*)' dir pert dir pert real part imaginary part' 2387 ipert1=natom+2 2388 do idir1=1,3 2389 write(iout,*)' ' 2390 ii=1 2391 do ipert2=natom+3,natom+4 2392 do idir2=1,3 2393 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2394 & piezofrnl(ii,idir1),zero 2395 ii=ii+1 2396 end do 2397 end do 2398 end do 2399 end if 2400 2401 ! In case of atom. displ/electric field perturbation (Born Effective Charges) 2402 if (pawbec==1.and.(rfpert(natom+2)==1.or.rfphon==1).and.outd2==1)then 2403 write(iout,*)' ' 2404 write(iout,*)' Frozen wf part of the Born Effective Charges' 2405 write(iout,*)' j1 j2 matrix element' 2406 write(iout,*)' dir pert dir pert real part imaginary part' 2407 ipert1 = natom+2 2408 do ipert2=1,natom 2409 do idir1=1,3 2410 write(iout,*)' ' 2411 do idir2=1,3 2412 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2413 & becfrnl(idir2,ipert2,idir1),zero 2414 end do 2415 end do 2416 end do 2417 end if 2418 2419 ! In case of strain 2420 if (rfstrs/=0)then 2421 2422 ! Write the Ewald part of the elastic tensor 2423 write(iout,*)' ' 2424 write(iout,*)' Ewald part of the elastic tensor in cartesian coordinates' 2425 write(iout,*)' j1 j2 matrix element' 2426 write(iout,*)' dir pert dir pert real part imaginary part' 2427 do ipert1=natom+3,natom+4 2428 do idir1=1,3 2429 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2430 & .or. outd2==1 )then 2431 ii=idir1+3*(ipert1-natom-3) 2432 write(iout,*)' ' 2433 do ipert2=natom+3,natom+4 2434 do idir2=1,3 2435 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2436 jj=idir2+3*(ipert2-natom-3) 2437 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2438 & elteew(ii,jj),zero 2439 end if 2440 end do 2441 end do 2442 end if 2443 end do 2444 end do 2445 2446 ! Write the Ewald part of the internal strain coupling parameters 2447 write(iout,*)' ' 2448 write(iout,*)' Ewald part of the internal strain coupling parameters' 2449 write(iout,*)' (cartesian strain, reduced atomic coordinates)' 2450 write(iout,*)' j1 j2 matrix element' 2451 write(iout,*)' dir pert dir pert real part imaginary part' 2452 do ipert1=1,natom 2453 do idir1=1,3 2454 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2455 & .or. outd2==1 )then 2456 ii=idir1+6+3*(ipert1-1) 2457 write(iout,*)' ' 2458 do ipert2=natom+3,natom+4 2459 do idir2=1,3 2460 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2461 jj=idir2+3*(ipert2-natom-3) 2462 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2463 & elteew(ii,jj),zero 2464 end if 2465 end do 2466 end do 2467 end if 2468 end do 2469 end do 2470 2471 ! Now the local frozen wf part 2472 write(iout,*)' ' 2473 write(iout,*)' Frozen wf local part of the elastic tensor in cartesian coordinates' 2474 write(iout,*)' j1 j2 matrix element' 2475 write(iout,*)' dir pert dir pert real part imaginary part' 2476 do ipert1=natom+3,natom+4 2477 do idir1=1,3 2478 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2479 & .or. outd2==1 )then 2480 ii=idir1+3*(ipert1-natom-3) 2481 write(iout,*)' ' 2482 do ipert2=natom+3,natom+4 2483 do idir2=1,3 2484 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2485 jj=idir2+3*(ipert2-natom-3) 2486 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2487 & eltfrloc(ii,jj),zero 2488 end if 2489 end do 2490 end do 2491 end if 2492 end do 2493 end do 2494 2495 write(iout,*)' ' 2496 write(iout,*)' Frozen wf local part of the internal strain coupling parameters ' 2497 write(iout,*)' (cartesian strain, reduced atomic coordinates)' 2498 write(iout,*)' j1 j2 matrix element' 2499 write(iout,*)' dir pert dir pert real part imaginary part' 2500 do ipert1=1,natom 2501 do idir1=1,3 2502 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2503 & .or. outd2==1 )then 2504 ii=idir1+6+3*(ipert1-1) 2505 write(iout,*)' ' 2506 do ipert2=natom+3,natom+4 2507 do idir2=1,3 2508 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2509 jj=idir2+3*(ipert2-natom-3) 2510 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2511 & eltfrloc(ii,jj),zero 2512 end if 2513 end do 2514 end do 2515 end if 2516 end do 2517 end do 2518 2519 ! Now the nonlo frozen wf part 2520 write(iout,*)' ' 2521 write(iout,*)' Frozen wf nonlocal part of the elastic tensor in cartesian coordinates' 2522 write(iout,*)' j1 j2 matrix element' 2523 write(iout,*)' dir pert dir pert real part imaginary part' 2524 do ipert1=natom+3,natom+4 2525 do idir1=1,3 2526 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2527 & .or. outd2==1 )then 2528 ii=idir1+3*(ipert1-natom-3) 2529 write(iout,*)' ' 2530 do ipert2=natom+3,natom+4 2531 do idir2=1,3 2532 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2533 jj=idir2+3*(ipert2-natom-3) 2534 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2535 & eltfrnl(ii,jj),zero 2536 end if 2537 end do 2538 end do 2539 end if 2540 end do 2541 end do 2542 2543 write(iout,*)' ' 2544 write(iout,*)' Frozen wf nonlocal part of the internal strain coupling parameters ' 2545 write(iout,*)' (cartesian strain, reduced atomic coordinates)' 2546 write(iout,*)' j1 j2 matrix element' 2547 write(iout,*)' dir pert dir pert real part imaginary part' 2548 do ipert1=1,natom 2549 do idir1=1,3 2550 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2551 & .or. outd2==1 )then 2552 ii=idir1+6+3*(ipert1-1) 2553 write(iout,*)' ' 2554 do ipert2=natom+3,natom+4 2555 do idir2=1,3 2556 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2557 jj=idir2+3*(ipert2-natom-3) 2558 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2559 & eltfrnl(ii,jj),zero 2560 end if 2561 end do 2562 end do 2563 end if 2564 end do 2565 end do 2566 2567 ! Now the xc part 2568 write(iout,*)' ' 2569 write(iout,*)' Frozen wf xc part of the elastic tensor in cartesian coordinates' 2570 write(iout,*)' j1 j2 matrix element' 2571 write(iout,*)' dir pert dir pert real part imaginary part' 2572 do ipert1=natom+3,natom+4 2573 do idir1=1,3 2574 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2575 & .or. outd2==1 )then 2576 ii=idir1+3*(ipert1-natom-3) 2577 write(iout,*)' ' 2578 do ipert2=natom+3,natom+4 2579 do idir2=1,3 2580 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2581 jj=idir2+3*(ipert2-natom-3) 2582 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2583 & eltfrxc(ii,jj),zero 2584 end if 2585 end do 2586 end do 2587 end if 2588 end do 2589 end do 2590 2591 write(iout,*)' ' 2592 write(iout,*)' Frozen wf xc part of the internal strain coupling parameters ' 2593 write(iout,*)' (cartesian strain, reduced atomic coordinates)' 2594 write(iout,*)' j1 j2 matrix element' 2595 write(iout,*)' dir pert dir pert real part imaginary part' 2596 do ipert1=1,natom 2597 do idir1=1,3 2598 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2599 & .or. outd2==1 )then 2600 ii=idir1+6+3*(ipert1-1) 2601 write(iout,*)' ' 2602 do ipert2=natom+3,natom+4 2603 do idir2=1,3 2604 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2605 jj=idir2+3*(ipert2-natom-3) 2606 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2607 & eltfrxc(ii,jj),zero 2608 end if 2609 end do 2610 end do 2611 end if 2612 end do 2613 end do 2614 2615 ! Now the kinetic frozen wf part 2616 write(iout,*)' ' 2617 write(iout,*)' Frozen wf kinetic part of the elastic tensor in cartesian coordinates' 2618 write(iout,*)' j1 j2 matrix element' 2619 write(iout,*)' dir pert dir pert real part imaginary part' 2620 do ipert1=natom+3,natom+4 2621 do idir1=1,3 2622 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2623 & .or. outd2==1 )then 2624 ii=idir1+3*(ipert1-natom-3) 2625 write(iout,*)' ' 2626 do ipert2=natom+3,natom+4 2627 do idir2=1,3 2628 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2629 jj=idir2+3*(ipert2-natom-3) 2630 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2631 & eltfrkin(ii,jj),zero 2632 end if 2633 end do 2634 end do 2635 end if 2636 end do 2637 end do 2638 2639 ! Now the hartree frozen wf part 2640 write(iout,*)' ' 2641 write(iout,*)' Frozen wf hartree part of the elastic tensor in cartesian coordinates' 2642 write(iout,*)' j1 j2 matrix element' 2643 write(iout,*)' dir pert dir pert real part imaginary part' 2644 do ipert1=natom+3,natom+4 2645 do idir1=1,3 2646 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2647 & .or. outd2==1 )then 2648 ii=idir1+3*(ipert1-natom-3) 2649 write(iout,*)' ' 2650 do ipert2=natom+3,natom+4 2651 do idir2=1,3 2652 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2653 jj=idir2+3*(ipert2-natom-3) 2654 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2655 & eltfrhar(ii,jj),zero 2656 end if 2657 end do 2658 end do 2659 end if 2660 end do 2661 end do 2662 2663 ! Now the psp core part 2664 write(iout,*)' ' 2665 write(iout,*)' Psp core part of the elastic tensor in cartesian coordinates' 2666 write(iout,*)' j1 j2 matrix element' 2667 write(iout,*)' dir pert dir pert real part imaginary part' 2668 do ipert1=natom+3,natom+4 2669 do idir1=1,3 2670 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1)& 2671 & .or. outd2==1 )then 2672 ii=idir1+3*(ipert1-natom-3) 2673 write(iout,*)' ' 2674 do ipert2=natom+3,natom+4 2675 do idir2=1,3 2676 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2677 jj=idir2+3*(ipert2-natom-3) 2678 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2679 & eltcore(ii,jj),zero 2680 end if 2681 end do 2682 end do 2683 end if 2684 end do 2685 end do 2686 2687 ! Now the DFT-D vdw part 2688 if (usevdw==1) then 2689 write(iout,*)' ' 2690 write(iout,*)' DFT-D van der Waals part of the elastic tensor in cartesian coordinates' 2691 write(iout,*)' j1 j2 matrix element' 2692 write(iout,*)' dir pert dir pert real part imaginary part' 2693 do ipert1=natom+3,natom+4 2694 do idir1=1,3 2695 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1)then 2696 ii=idir1+3*(ipert1-natom-3) 2697 write(iout,*)' ' 2698 do ipert2=natom+3,natom+4 2699 do idir2=1,3 2700 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2701 jj=idir2+3*(ipert2-natom-3) 2702 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltvdw(ii,jj),zero 2703 end if 2704 end do 2705 end do 2706 end if 2707 end do 2708 end do 2709 2710 write(iout,*)' ' 2711 write(iout,*)' DFT-D2 van der Waals part of the internal strain coupling parameters' 2712 write(iout,*)' (cartesian strain, reduced atomic coordinates)' 2713 write(iout,*)' j1 j2 matrix element' 2714 write(iout,*)' dir pert dir pert real part imaginary part' 2715 do ipert1=1,natom 2716 do idir1=1,3 2717 if ( (rfpert(ipert1)==1.and.rfdir(idir1)==1) .or. outd2==1 )then 2718 ii=idir1+6+3*(ipert1-1) 2719 write(iout,*)' ' 2720 do ipert2=natom+3,natom+4 2721 do idir2=1,3 2722 if (rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2723 jj=idir2+3*(ipert2-natom-3) 2724 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,eltvdw(ii,jj),zero 2725 end if 2726 end do 2727 end do 2728 end if 2729 end do 2730 end do 2731 end if ! usevdw 2732 2733 end if ! strain condition 2734 2735 ! Now the local nonstationary nonfrozenwf part 2736 if (outd2==1)then 2737 write(iout,*)' ' 2738 write(iout,*)' Non-stationary local part of the 2-order matrix' 2739 write(iout,*)' j1 j2 matrix element' 2740 write(iout,*)' dir pert dir pert real part imaginary part' 2741 do ipert1=1,mpert 2742 do idir1=1,3 2743 if ((ipert1<=natom .or.& 2744 & (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.& 2745 & ((ipert1==natom+3.or.ipert1==natom+4).and.& 2746 & (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then 2747 write(iout,*)' ' 2748 do ipert2=1,mpert 2749 do idir2=1,3 2750 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2751 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2752 & d2lo(1,idir1,ipert1,idir2,ipert2),& 2753 & d2lo(2,idir1,ipert1,idir2,ipert2) 2754 end if 2755 end do 2756 end do 2757 end if 2758 end do 2759 end do 2760 end if 2761 2762 ! Now the nonlocal nonstationary nonfrozenwf part 2763 if (outd2==1)then 2764 write(iout,*)' ' 2765 write(iout,*)' Non-stationary non-local part of the 2nd-order matrix' 2766 write(iout,*)' j1 j2 matrix element' 2767 write(iout,*)' dir pert dir pert real part imaginary part' 2768 do ipert1=1,mpert 2769 do idir1=1,3 2770 if ((ipert1<=natom .or.& 2771 & (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.& 2772 & ((ipert1==natom+3.or.ipert1==natom+4).and.& 2773 & (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then 2774 write(iout,*)' ' 2775 do ipert2=1,mpert 2776 do idir2=1,3 2777 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2778 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2779 & d2nl(1,idir1,ipert1,idir2,ipert2),& 2780 & d2nl(2,idir1,ipert1,idir2,ipert2) 2781 end if 2782 end do 2783 end do 2784 end if 2785 end do 2786 end do 2787 end if 2788 2789 ! Now the overlap change nonstationnary nonfrozenwf part (PAW only) 2790 if (outd2==1.and.usepaw==1)then 2791 write(iout,*)' ' 2792 write(iout,*)' PAW: Non-stationary WF-overlap part of the 2nd-order matrix' 2793 write(iout,*)' j1 j2 matrix element' 2794 write(iout,*)' dir pert dir pert real part imaginary part' 2795 do ipert1=1,mpert 2796 do idir1=1,3 2797 if ((ipert1<=natom .or.& 2798 & (ipert1==natom+2.and.qzero==1.and.ddkfil(idir1)/=0)).or.& 2799 & ((ipert1==natom+3.or.ipert1==natom+4).and.& 2800 & (rfpert(natom+3)==1.or.rfpert(natom+4)==1)))then 2801 write(iout,*)' ' 2802 do ipert2=1,mpert 2803 do idir2=1,3 2804 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2805 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2806 & d2ovl(1,idir1,ipert1,idir2,ipert2),& 2807 & d2ovl(2,idir1,ipert1,idir2,ipert2) 2808 end if 2809 end do 2810 end do 2811 end if 2812 end do 2813 end do 2814 end if 2815 2816 ! Now the 0-order local stationary nonfrozenwf part 2817 if (outd2==2)then 2818 write(iout,*)' ' 2819 write(iout,*)' Stationary 0-order local part of the 2nd-order matrix' 2820 write(iout,*)' j1 j2 matrix element' 2821 write(iout,*)' dir pert dir pert real part imaginary part' 2822 do ipert1=1,mpert 2823 do idir1=1,3 2824 if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then 2825 write(iout,*)' ' 2826 do ipert2=1,mpert 2827 do idir2=1,3 2828 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2829 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2830 & d2loc0(1,idir1,ipert1,idir2,ipert2),& 2831 & d2loc0(2,idir1,ipert1,idir2,ipert2) 2832 end if 2833 end do 2834 end do 2835 end if 2836 end do 2837 end do 2838 end if 2839 2840 ! Now the stationary 0-order kinetic nonfrozenwf part 2841 if (outd2==2)then 2842 write(iout,*)' ' 2843 write(iout,*)' Stationary 0-order kinetic part of the 2nd-order matrix' 2844 write(iout,*)' j1 j2 matrix element' 2845 write(iout,*)' dir pert dir pert real part imaginary part' 2846 do ipert1=1,mpert 2847 do idir1=1,3 2848 if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then 2849 write(iout,*)' ' 2850 do ipert2=1,mpert 2851 do idir2=1,3 2852 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2853 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2854 & d2k0(1,idir1,ipert1,idir2,ipert2),& 2855 & d2k0(2,idir1,ipert1,idir2,ipert2) 2856 end if 2857 end do 2858 end do 2859 end if 2860 end do 2861 end do 2862 end if 2863 2864 ! Now the stationary 0-order eigenvalue nonfrozenwf part 2865 if (outd2==2)then 2866 write(iout,*)' ' 2867 write(iout,*)' Stationary 0-order eigenvalue part of the' ,' 2nd-order matrix' 2868 write(iout,*)' j1 j2 matrix element' 2869 write(iout,*)' dir pert dir pert real part imaginary part' 2870 do ipert1=1,mpert 2871 do idir1=1,3 2872 if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then 2873 write(iout,*)' ' 2874 do ipert2=1,mpert 2875 do idir2=1,3 2876 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2877 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2878 & d2eig0(1,idir1,ipert1,idir2,ipert2),& 2879 & d2eig0(2,idir1,ipert1,idir2,ipert2) 2880 end if 2881 end do 2882 end do 2883 end if 2884 end do 2885 end do 2886 end if 2887 2888 ! Now the stationary potential-density nonfrozenwf part 2889 if (outd2==2)then 2890 write(iout,*)' ' 2891 write(iout,*)' Station. potential-density part of the ',& 2892 & ' 2nd-order matrix' 2893 write(iout,*)' (Note : include some xc core-correction) ' 2894 write(iout,*)' j1 j2 matrix element' 2895 write(iout,*)' dir pert dir pert real part imaginary part' 2896 do ipert1=1,mpert 2897 do idir1=1,3 2898 if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then 2899 write(iout,*)' ' 2900 do ipert2=1,mpert 2901 do idir2=1,3 2902 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2903 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2904 & d2vn(1,idir1,ipert1,idir2,ipert2),& 2905 & d2vn(2,idir1,ipert1,idir2,ipert2) 2906 end if 2907 end do 2908 end do 2909 end if 2910 end do 2911 end do 2912 end if 2913 2914 ! Now the stationary 0-order nonloc nonfrozenwf part 2915 if (outd2==2)then 2916 write(iout,*)' ' 2917 write(iout,*)' Stationary 0-order nonlocal part of the 2-order'& 2918 & ,' matrix' 2919 write(iout,*)' j1 j2 matrix element' 2920 write(iout,*)' dir pert dir pert real part imaginary part' 2921 do ipert1=1,mpert 2922 do idir1=1,3 2923 if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then 2924 write(iout,*)' ' 2925 do ipert2=1,mpert 2926 do idir2=1,3 2927 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2928 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2929 & d2nl0(1,idir1,ipert1,idir2,ipert2),& 2930 & d2nl0(2,idir1,ipert1,idir2,ipert2) 2931 end if 2932 end do 2933 end do 2934 end if 2935 end do 2936 end do 2937 end if 2938 2939 ! Now the stationary 1-order nonloc nonfrozenwf part 2940 if (outd2==2)then 2941 write(iout,*)' ' 2942 write(iout,*)' Stationary 1-order nonlocal part of the'& 2943 & ,' 2nd-order matrix' 2944 write(iout,*)' (or the ddk wf part of it, in case of',& 2945 & ' an electric field perturbation )' 2946 write(iout,*)' j1 j2 matrix element' 2947 write(iout,*)' dir pert dir pert real part imaginary part' 2948 do ipert1=1,mpert 2949 do idir1=1,3 2950 if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then 2951 write(iout,*)' ' 2952 do ipert2=1,mpert 2953 do idir2=1,3 2954 if(rfpert(ipert2)==1.and.rfdir(idir2)==1)then 2955 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 2956 & d2nl1(1,idir1,ipert1,idir2,ipert2),& 2957 & d2nl1(2,idir1,ipert1,idir2,ipert2) 2958 end if 2959 end do 2960 end do 2961 end if 2962 end do 2963 end do 2964 end if 2965 2966 ! End of the long print out condition 2967 end if 2968 2969 2970 !Derivative database initialisation 2971 2972 !Calculation of the number of elements 2973 nelmts=0 2974 do ipert1=1,mpert 2975 do idir1=1,3 2976 do ipert2=1,mpert 2977 do idir2=1,3 2978 nelmts=nelmts+blkflg(idir1,ipert1,idir2,ipert2) 2979 end do 2980 end do 2981 end do 2982 end do 2983 2984 !Now the whole 2nd-order matrix, but not in cartesian coordinates, 2985 !and masses not included 2986 write(iout,*)' ' 2987 write(iout,*)' 2nd-order matrix (non-cartesian coordinates,',' masses not included,' 2988 write(iout,*)' asr not included )' 2989 if(rfstrs/=0) then 2990 write(iout,*)' cartesian coordinates for strain terms (1/ucvol factor ' 2991 write(iout,*)' for elastic tensor components not included) ' 2992 end if 2993 write(iout,*)' j1 j2 matrix element' 2994 write(iout,*)' dir pert dir pert real part imaginary part' 2995 nline=1 2996 do ipert1=1,mpert 2997 do idir1=1,3 2998 if(nline/=0)write(iout,*)' ' 2999 nline=0 3000 do ipert2=1,mpert 3001 do idir2=1,3 3002 if(blkflg(idir1,ipert1,idir2,ipert2)==1)then 3003 nline=nline+1 3004 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3005 & d2matr(1,idir1,ipert1,idir2,ipert2),& 3006 & d2matr(2,idir1,ipert1,idir2,ipert2) 3007 end if 3008 end do 3009 end do 3010 end do 3011 end do 3012 3013 !Now the dynamical matrix 3014 if(rfphon==1)then 3015 write(iout,*)' ' 3016 write(iout,*)' Dynamical matrix, in cartesian coordinates,' 3017 write(iout,*)' if specified in the inputs, asr has been imposed' 3018 write(iout,*)' j1 j2 matrix element' 3019 write(iout,*)' dir pert dir pert real part imaginary part' 3020 nline=1 3021 do ipert1=1,natom 3022 do idir1=1,3 3023 if(nline/=0)write(iout,*)' ' 3024 nline=0 3025 do ipert2=1,natom 3026 do idir2=1,3 3027 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3028 nline=nline+1 3029 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3030 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3031 & d2cart(2,idir1,ipert1,idir2,ipert2) 3032 end if 3033 end do 3034 end do 3035 end do 3036 end do 3037 end if 3038 3039 !Now the dielectric tensor ! normal case 3040 if(rfpert(natom+2)==1)then 3041 3042 write(iout,*)' ' 3043 write(iout,*)' Dielectric tensor, in cartesian coordinates,' 3044 write(iout,*)' j1 j2 matrix element' 3045 write(iout,*)' dir pert dir pert real part imaginary part' 3046 ipert1=natom+2 3047 ipert2=natom+2 3048 nline=1 3049 do idir1=1,3 3050 if(nline/=0)write(iout,*)' ' 3051 nline=0 3052 do idir2=1,3 3053 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3054 nline=nline+1 3055 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3056 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3057 & d2cart(2,idir1,ipert1,idir2,ipert2) 3058 end if 3059 end do 3060 end do 3061 3062 if (prtbbb == 1) then 3063 3064 delta(:,:) = zero 3065 delta(1,1) = one ; delta(2,2) = one ; delta(3,3) = one 3066 3067 write(iout,*) 3068 write(iout,*)'Band by band decomposition of the dielectric tensor' 3069 write(iout,*)' ' 3070 3071 write(iout,*)' Vacuum polarization' 3072 write(iout,*)' j1 j2 matrix element' 3073 write(iout,*)' dir pert dir pert real part imaginary part' 3074 nline=1 3075 do idir1=1,3 3076 if(nline/=0)write(iout,*)' ' 3077 nline=0 3078 do idir2=1,3 3079 nline=nline+1 3080 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3081 & delta(idir2,idir1),zero 3082 end do 3083 end do 3084 3085 do iband = 1,mband 3086 write(iout,*)' ' 3087 write(iout,*)' Dielectric tensor, in cartesian coordinates, for band',iband 3088 write(iout,*)' j1 j2 matrix element' 3089 write(iout,*)' dir pert dir pert real part imaginary part' 3090 ipert1 = natom + 2 3091 ipert2 = natom + 2 3092 nline=1 3093 do idir1=1,3 3094 if(nline/=0)write(iout,*)' ' 3095 nline=0 3096 do idir2=1,3 3097 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3098 ! substract vacuum polarization 3099 if (idir1 == idir2) then 3100 d2cart_bbb(1,idir1,idir2,ipert2,iband,iband) = & 3101 & d2cart_bbb(1,idir1,idir2,ipert2,iband,iband) - 1 3102 end if 3103 nline=nline+1 3104 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3105 & d2cart_bbb(1,idir1,idir2,ipert2,iband,iband),& 3106 & d2cart_bbb(2,idir1,idir2,ipert2,iband,iband) 3107 end if 3108 end do 3109 end do 3110 end do !iband 3111 3112 end if !prtbbb 3113 3114 end if ! end natom+2 dielectric output 3115 3116 !Now the effective charges 3117 !In case of the stationary calculation 3118 if(outd2==2 .and. rfpert(natom+2)==1 .and.rfphon==1)then 3119 write(iout,*)' ' 3120 write(iout,*)' Effective charges, in cartesian coordinates,' 3121 write(iout,*)' if specified in the inputs, charge neutrality has been imposed' 3122 write(iout,*)' j1 j2 matrix element' 3123 write(iout,*)' dir pert dir pert real part imaginary part' 3124 ipert1=natom+2 3125 nline=1 3126 do idir1=1,3 3127 if(nline/=0)write(iout,*)' ' 3128 nline=0 3129 do ipert2=1,natom 3130 do idir2=1,3 3131 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3132 nline=nline+1 3133 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3134 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3135 & d2cart(2,idir1,ipert1,idir2,ipert2) 3136 end if 3137 end do 3138 end do 3139 end do 3140 end if 3141 3142 !Now in case of the non-stationary calculation 3143 if(outd2==1 .and. rfpert(natom+2)==1)then 3144 write(iout,*)' ' 3145 if(usepaw==1.and..not.(has_allddk))then 3146 write(iout,*)' Warning: Born effectives charges are not correctly computed' 3147 write(iout,*)' you need all ddk perturbations!' 3148 end if 3149 write(iout,*)' Effective charges, in cartesian coordinates,' 3150 write(iout,*)' (from electric field response) ' 3151 write(iout,*)' if specified in the inputs, charge neutrality has been imposed' 3152 write(iout,*)' j1 j2 matrix element' 3153 write(iout,*)' dir pert dir pert real part imaginary part' 3154 ipert2=natom+2 3155 nline=1 3156 do idir2=1,3 3157 if(nline/=0)write(iout,*)' ' 3158 nline=0 3159 do ipert1=1,natom 3160 do idir1=1,3 3161 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3162 nline=nline+1 3163 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3164 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3165 & d2cart(2,idir1,ipert1,idir2,ipert2) 3166 end if 3167 end do 3168 end do 3169 end do 3170 end if 3171 3172 if(outd2==1 .and. rfphon==1 .and. qzero==1& 3173 & .and. ( (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0) .or. & 3174 & berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17 ) )then !!HONG need to test for fixed E and D 3175 write(iout,*)' ' 3176 if(usepaw==1.and..not.(has_allddk))then 3177 write(iout,*)' Warning: Born effectives charges are not correctly computed' 3178 write(iout,*)' you need all ddk perturbations!' 3179 end if 3180 write(iout,*)' Effective charges, in cartesian coordinates,' 3181 write(iout,*)' (from phonon response) ' 3182 write(iout,*)' if specified in the inputs, charge neutrality has been imposed' 3183 write(iout,*)' j1 j2 matrix element' 3184 write(iout,*)' dir pert dir pert real part imaginary part' 3185 nline=1 3186 do ipert2=1,natom 3187 do idir2=1,3 3188 if(nline/=0)write(iout,*)' ' 3189 nline=0 3190 ipert1=natom+2 3191 do idir1=1,3 3192 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3193 nline=nline+1 3194 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3195 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3196 & d2cart(2,idir1,ipert1,idir2,ipert2) 3197 end if 3198 end do 3199 end do 3200 end do 3201 write(iout,*)' ' 3202 write(iout,*)' ' 3203 write(iout,*)' ' 3204 3205 if (prtbbb == 1) then 3206 3207 write(iout,*)'Band by band decomposition of the Born effective charges' 3208 write(iout,*)' ' 3209 write(iout,*)'Ionic charges in cartesian coordinates' 3210 write(iout,*)' j1 j2 matrix element' 3211 write(iout,*)' dir pert dir pert real part imaginary part' 3212 zr = zero 3213 zi = zero 3214 do ipert2=1,natom 3215 do idir2=1,3 3216 if(nline/=0)write(iout,*)' ' 3217 nline=0 3218 ipert1=natom+2 3219 do idir1=1,3 3220 zr = zero 3221 if (idir1 == idir2) then 3222 zr = zion(typat(ipert2)) 3223 end if 3224 nline=nline+1 3225 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3226 & zr,zi 3227 end do 3228 end do 3229 end do 3230 3231 do iband = 1,mband 3232 write(iout,*)' ' 3233 write(iout,*)' Effective charges, in cartesian coordinates, for band',iband 3234 write(iout,*)' (from phonon response) ' 3235 write(iout,*)' if specified in the inputs, charge neutrality has been imposed' 3236 write(iout,*)' j1 j2 matrix element' 3237 write(iout,*)' dir pert dir pert real part imaginary part' 3238 nline=1 3239 do ipert2=1,natom 3240 do idir2=1,3 3241 if(nline/=0)write(iout,*)' ' 3242 nline=0 3243 ipert1=natom+2 3244 do idir1=1,3 3245 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3246 nline=nline+1 3247 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3248 & d2cart_bbb(1,idir1,idir2,ipert2,iband,iband),& 3249 & d2cart_bbb(2,idir1,idir2,ipert2,iband,iband) 3250 end if 3251 end do 3252 end do 3253 end do 3254 end do !iband 3255 end if !prtbbb 3256 end if ! end of print effective charges 3257 3258 !Now the elastic tensor 3259 if(rfstrs/=0) then 3260 write(iout,*)' ' 3261 write(iout,*)' Rigid-atom elastic tensor , in cartesian coordinates,' 3262 write(iout,*)' j1 j2 matrix element' 3263 write(iout,*)' dir pert dir pert real part imaginary part' 3264 nline=1 3265 do ipert1=natom+3,natom+4 3266 do idir1=1,3 3267 if(nline/=0)write(iout,*)' ' 3268 nline=0 3269 do ipert2=natom+3,natom+4 3270 do idir2=1,3 3271 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3272 nline=nline+1 3273 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3274 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3275 & d2cart(2,idir1,ipert1,idir2,ipert2) 3276 end if 3277 end do 3278 end do 3279 end do 3280 end do 3281 end if 3282 3283 !Now the internal strain coupling parameters 3284 if(rfstrs/=0) then 3285 write(iout,*)' ' 3286 write(iout,*)' Internal strain coupling parameters, in cartesian coordinates,' 3287 write(iout,*)' zero average net force deriv. has been imposed ' 3288 write(iout,*)' j1 j2 matrix element' 3289 write(iout,*)' dir pert dir pert real part imaginary part' 3290 nline=1 3291 do ipert1=1,natom 3292 do idir1=1,3 3293 if(nline/=0)write(iout,*)' ' 3294 nline=0 3295 do ipert2=natom+3,natom+4 3296 do idir2=1,3 3297 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3298 nline=nline+1 3299 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3300 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3301 & d2cart(2,idir1,ipert1,idir2,ipert2) 3302 end if 3303 end do 3304 end do 3305 end do 3306 end do 3307 end if 3308 3309 !Now the piezoelectric tensor 3310 if(rfstrs/=0 .and. (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0))then 3311 write(iout,*)' ' 3312 if(usepaw==1.and..not.(has_allddk))then 3313 write(iout,*)' Warning: Rigid-atom proper piezoelectric tensor is not correctly computed' 3314 write(iout,*)' you need all ddk perturbations!' 3315 end if 3316 write(iout,*)' Rigid-atom proper piezoelectric tensor, in cartesian coordinates,' 3317 write(iout,*)' (from strain response)' 3318 write(iout,*)' j1 j2 matrix element' 3319 write(iout,*)' dir pert dir pert real part imaginary part' 3320 nline=1 3321 ipert1=natom+2 3322 do idir1=1,3 3323 if(nline/=0)write(iout,*)' ' 3324 nline=0 3325 do ipert2=natom+3,natom+4 3326 do idir2=1,3 3327 if(carflg(idir1,ipert1,idir2,ipert2)==1)then 3328 nline=nline+1 3329 write(iout,'(2(i4,i5),2(1x,f20.10))')idir1,ipert1,idir2,ipert2,& 3330 & d2cart(1,idir1,ipert1,idir2,ipert2),& 3331 & d2cart(2,idir1,ipert1,idir2,ipert2) 3332 end if 3333 end do 3334 end do 3335 end do 3336 end if 3337 3338 !Now the piezoelectric tensor 3339 if(outd2==1 .and. (pawpiezo==1.and.rfpert(natom+2)==1)& 3340 & .and. (ddkfil(1)/=0.or.ddkfil(2)/=0.or.ddkfil(3)/=0)) then 3341 write(iout,*)' ' 3342 if(usepaw==1.and..not.(has_allddk))then 3343 write(iout,*)' Warning: Rigid-atom proper piezoelectric tensor is not correctly computed' 3344 write(iout,*)' you need all ddk perturbations!' 3345 end if 3346 if(usepaw==1.and..not.has_full_piezo)then 3347 write(iout,*)' Warning: The rigid-atom proper piezoelectric tensor' 3348 write(iout,*)' from electric field response requires nsym=1' 3349 end if 3350 if (has_full_piezo) then 3351 write(iout,*)' Rigid-atom proper piezoelectric tensor, in cartesian coordinates,' 3352 write(iout,*)' (from electric field response)' 3353 write(iout,*)' j1 j2 matrix element' 3354 write(iout,*)' dir pert dir pert real part imaginary part' 3355 nline=1 3356 ipert1=natom+2 3357 do idir1=1,3 3358 if(nline/=0)write(iout,*)' ' 3359 nline=0 3360 do ipert2=natom+3,natom+4 3361 do idir2=1,3 3362 if(carflg(idir2,ipert2,idir1,ipert1)==1)then 3363 nline=nline+1 3364 write(iout,'(2(i4,i5),2(1x,f20.10))')idir2,ipert2,idir1,ipert1,& 3365 & d2cart(1,idir2,ipert2,idir1,ipert1),& 3366 & d2cart(2,idir2,ipert2,idir1,ipert1) 3367 end if 3368 end do 3369 end do 3370 end do 3371 end if 3372 end if 3373 3374 end subroutine dfpt_dyout
ABINIT/dfpt_dyxc1 [ Functions ]
NAME
dfpt_dyxc1
FUNCTION
Compute 2nd-order non-linear xc core-correction (part1) to the dynamical matrix. In case of derivative with respect to k or electric field perturbation, the 1st-order local potential vanishes.
INPUTS
atindx(natom)=index table for atoms gmet(3,3)=metrix tensor in G space in Bohr**-2. gsqcut=cutoff value on G**2 for sphere inside fft box. ixc= choice of exchange-correlation scheme kxc(nfft,nkxc)=first-order derivative of the xc potential if (nkxc=1) LDA kxc(:,1)= d2Exc/drho2 if (nkxc=2) LDA kxc(:,1)= d2Exc/drho_up drho_up kxc(:,2)= d2Exc/drho_up drho_dn kxc(:,3)= d2Exc/drho_dn drho_dn if (nkxc=7) GGA kxc(:,1)= d2Exc/drho2 kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)| kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| ) kxc(:,5)= gradx(rho) kxc(:,6)= grady(rho) kxc(:,7)= gradz(rho) if (nkxc=19) spin-polarized GGA case (same as nkxc=7 with up and down components) mgfft=maximum size of 1D FFTs mpert=maximum number of ipert mpi_enreg=information about MPI parallelization mqgrid=number of grid pts in q array for f(q) spline. natom=number of atoms in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(3)=fft grid dimensions. nkxc=second dimension of the kxc array (=1 for non-spin-polarized case, =3 for spin-polarized case) nmxc= if true, handle density/potential as non-magnetic (even if it is) nspden=number of spin-density components ntypat=number of types of atoms in cell. n1xccc=dimension of xccc1d ; 0 if no XC core correction is used psps <type(pseudopotential_type)>=variables related to pseudopotentials pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information qgrid(mqgrid)=q grid for spline from 0 to qmax. qphon(3)=wavevector of the phonon rfdir(3)=array that define the directions of perturbations rfpert(mpert)=array defining the type of perturbations that have to be computed rprimd(3,3)=dimensional primitive translation vectors (bohr) timrev=1 if time-reversal preserves the q wavevector; 0 otherwise. typat(natom)=integer type for each atom in cell ucvol=unit cell volume (bohr**3). usepaw= 0 for non paw calculation; =1 for paw calculation xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives, for each type of atom, from psp xred(3,natom)=fractional coordinates for atoms in unit cell OUTPUT! non_magnetic_xc= if true, handle density/potential as non-magnetic (even if it is) blkflgfrx1(3,natom,3,natom)=flag to indicate whether an element has been computed or not dyfrx1(2,3,natom,3,natom)=2nd-order non-linear xc core-correction (part1) part of the dynamical matrix
SOURCE
4149 subroutine dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,ixc,kxc,mgfft,mpert,mpi_enreg,mqgrid,& 4150 & natom,nfft,ngfft,nkxc,nmxc,nspden,ntypat,n1xccc,psps,pawtab,& 4151 & ph1d,qgrid,qphon,rfdir,rfpert,rprimd,timrev,typat,ucvol,usepaw,xcccrc,xccc1d,xred,rhor,vxc) 4152 4153 use m_cgtools, only : dotprod_vn 4154 use m_atm2fft, only : dfpt_atm2fft 4155 use m_dfpt_mkvxc, only : dfpt_mkvxc, dfpt_mkvxc_noncoll 4156 4157 !Arguments ------------------------------------ 4158 !scalars 4159 integer,intent(in) :: ixc,mgfft,mpert,mqgrid,n1xccc,natom,nfft,nkxc,nspden,ntypat 4160 integer,intent(in) :: timrev,usepaw 4161 logical,intent(in) :: nmxc 4162 real(dp),intent(in) :: gsqcut,ucvol 4163 type(pseudopotential_type),intent(in) :: psps 4164 type(MPI_type),intent(in) :: mpi_enreg 4165 !arrays 4166 integer,intent(in) :: atindx(natom),ngfft(18),rfdir(3),rfpert(mpert),typat(natom) 4167 real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc) 4168 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3) 4169 real(dp),intent(in) :: rprimd(3,3),xccc1d(n1xccc,6,ntypat),xcccrc(ntypat) 4170 real(dp),intent(in) :: xred(3,natom) 4171 integer,intent(out) :: blkflgfrx1(3,natom,3,natom) 4172 real(dp),intent(out) :: dyfrx1(2,3,natom,3,natom) 4173 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 4174 !optional 4175 real(dp),optional,intent(in) :: rhor(nfft,nspden) 4176 real(dp),optional,intent(in) :: vxc(nfft,nspden) 4177 4178 !Local variables------------------------------- 4179 !scalars 4180 integer :: cplex,iat1,iatom1,iatom2,idir1,idir2,ierr,ifft,my_natom,comm_atom 4181 integer :: n1,n2,n3,n3xccc,nfftot,option,upperdir,optnc 4182 logical :: paral_atom 4183 real(dp) :: valuei,valuer 4184 !arrays 4185 integer,pointer :: my_atmtab(:) 4186 real(dp) :: tsec(2),gprimd_dummy(3,3) 4187 real(dp) :: dum_nhat(0) 4188 real(dp),allocatable :: rhor1(:,:),vxc10(:,:),xcccwk1(:),xcccwk2(:) 4189 ! ********************************************************************* 4190 4191 call timab(182,1,tsec) 4192 4193 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 4194 nfftot=n1*n2*n3 4195 4196 !Set up parallelism over atoms 4197 my_natom=mpi_enreg%my_natom 4198 my_atmtab=>mpi_enreg%my_atmtab 4199 comm_atom=mpi_enreg%comm_atom 4200 paral_atom=(my_natom/=natom) 4201 4202 !Zero out the output arrays : 4203 blkflgfrx1(:,:,:,:)=0 4204 dyfrx1(:,:,:,:,:)=zero 4205 4206 cplex=2-timrev ; n3xccc=nfft 4207 ABI_MALLOC(vxc10,(cplex*nfft,nspden)) 4208 4209 4210 !Loop on the perturbation j1 4211 do iat1=1,my_natom 4212 iatom1=iat1; if(paral_atom)iatom1=my_atmtab(iat1) 4213 do idir1=1,3 4214 4215 ! Compute the derivative of the core charge with respect to j1 4216 ABI_MALLOC(xcccwk1,(cplex*n3xccc)) 4217 4218 ! PAW or NC with nc_xccc_gspace: 1st-order core charge in reciprocal space 4219 if (usepaw==1 .or. psps%nc_xccc_gspace==1) then 4220 call dfpt_atm2fft(atindx,cplex,gmet,gprimd_dummy,gsqcut,idir1,iatom1,& 4221 & mgfft,mqgrid,natom,1,nfft,ngfft,ntypat,ph1d,qgrid,& 4222 & qphon,typat,ucvol,usepaw,xred,psps,pawtab,atmrhor1=xcccwk1,optn2_in=1) 4223 4224 ! Norm-conserving psp: 1st-order core charge in real space 4225 else 4226 call dfpt_mkcore(cplex,idir1,iatom1,natom,ntypat,n1,n1xccc,& 4227 & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xcccwk1,xred) 4228 end if 4229 4230 ! Compute the corresponding potential 4231 option=0 4232 ABI_MALLOC(rhor1,(cplex*nfft,nspden)) 4233 rhor1=zero 4234 !FR SPr EB Non-collinear magnetism 4235 if (nspden==4.and.present(rhor).and.present(vxc)) then 4236 optnc=1 4237 call dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,dum_nhat,0,dum_nhat,0,dum_nhat,0,nkxc,& 4238 & nmxc,nspden,n3xccc,optnc,option,qphon,rhor,rhor1,rprimd,0,vxc,vxc10,xcccwk1) 4239 else 4240 call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,dum_nhat,0,dum_nhat,0,nkxc,& 4241 & nmxc,nspden,n3xccc,option,qphon,rhor1,rprimd,0,vxc10,xcccwk1) 4242 end if 4243 ABI_FREE(rhor1) 4244 ABI_FREE(xcccwk1) 4245 4246 ! vxc10 will couple with xcccwk2, that behaves like 4247 ! a total density (ispden=1). Only the spin-up + spin-down 4248 ! average of vxc10 is needed. 4249 if (nspden/=1)then 4250 do ifft=1,cplex*nfft 4251 vxc10(ifft,1)=(vxc10(ifft,1)+vxc10(ifft,2))*half 4252 end do 4253 end if 4254 4255 ! Loop on the perturbation j2 4256 do iatom2=1,iatom1 4257 upperdir=3 4258 if(iatom1==iatom2)upperdir=idir1 4259 do idir2=1,upperdir 4260 if( (rfpert(iatom1)==1 .and. rfdir(idir1) == 1) .or. & 4261 & (rfpert(iatom2)==1 .and. rfdir(idir2) == 1) )then 4262 4263 ! Compute the derivative of the core charge with respect to j2 4264 ABI_MALLOC(xcccwk2,(cplex*n3xccc)) 4265 4266 ! PAW or NC with nc_xccc_gspace: 1st-order core charge in reciprocal space 4267 if (usepaw==1 .or. psps%nc_xccc_gspace==1) then 4268 call dfpt_atm2fft(atindx,cplex,gmet,gprimd_dummy,gsqcut,idir2,iatom2,& 4269 & mgfft,mqgrid,natom,1,nfft,ngfft,ntypat,ph1d,qgrid,& 4270 & qphon,typat,ucvol,usepaw,xred,psps,pawtab,atmrhor1=xcccwk2,optn2_in=1) 4271 4272 ! Norm-conserving psp: 1st-order core charge in real space 4273 else 4274 call dfpt_mkcore(cplex,idir2,iatom2,natom,ntypat,n1,n1xccc,& 4275 & n2,n3,qphon,rprimd,typat,ucvol,xcccrc,xccc1d,xcccwk2,xred) 4276 end if 4277 4278 ! Get the matrix element j1,j2 4279 4280 call dotprod_vn(cplex,xcccwk2,valuer,valuei,nfft,nfftot,1,2,vxc10,ucvol) 4281 4282 ABI_FREE(xcccwk2) 4283 4284 dyfrx1(1,idir1,iatom1,idir2,iatom2)= valuer 4285 dyfrx1(2,idir1,iatom1,idir2,iatom2)= valuei 4286 dyfrx1(1,idir2,iatom2,idir1,iatom1)= valuer 4287 dyfrx1(2,idir2,iatom2,idir1,iatom1)=-valuei 4288 blkflgfrx1(idir1,iatom1,idir2,iatom2)=1 4289 blkflgfrx1(idir2,iatom2,idir1,iatom1)=1 4290 end if 4291 end do 4292 end do 4293 end do 4294 end do 4295 4296 if (paral_atom) then 4297 call timab(48,1,tsec) 4298 call xmpi_sum(dyfrx1,comm_atom,ierr) 4299 call xmpi_sum(blkflgfrx1,comm_atom,ierr) 4300 call timab(48,2,tsec) 4301 end if 4302 4303 ABI_FREE(vxc10) 4304 4305 call timab(182,2,tsec) 4306 4307 end subroutine dfpt_dyxc1
ABINIT/dfpt_gatherdy [ Functions ]
NAME
dfpt_gatherdy
FUNCTION
Sum (gather) the different parts of the 2nd-order matrix, to get the matrix of second-order derivatives (d2matr) Then, generates the dynamical matrix, not including the masses, but the correct non-cartesian coordinates ( => d2cart)
INPUTS
asr= (0=> no acoustic sum rule [asr] imposed), (1 or 2=> asr is imposed) only for dynamical matrix at Gamma becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only) berryopt=option for berry phase treatment blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical matrix has been calculated ; 0 otherwise ) chneut= (0=> no charge neutrality sum rule imposed), (1=> charge neutrality is imposed, with equal repartition of the charge neutrality correction for the effective charges), (2=> charge neutrality is imposed, with weighted repartition of the charge neutrality correction for the effective charges), dyew(2,3,natom,3,natom)=Ewald part of the dyn.matrix dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wf part of the dyn.matrix (except xc1) dyfrx1(2,3,natom,3,natom)=xc core correction (1) part of the frozen-wf part of the dynamical matrix. dyfr_cplex=1 if dyfrnl is real, 2 if it is complex dyfr_nondiag=1 if dyfrwf is non diagonal with respect to atoms; 0 otherwise dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some second order derivatives d2nfr(2,3,mpert,3,mpert)=non-frozen wf part of the 2nd-order matr eltcore(6,6)=core contribution to the elastic tensor elteew(6+3*natom,6)=Ewald contribution to the elastic tsenor eltfrhar(6,6)=hartree contribution to the elastic tensor eltfrkin(6,6)=kinetic contribution to the elastic tensor eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor gprimd(3,3)=basis vector in the reciprocal space mband=maximum number of bands mpert =maximum number of ipert natom=number of atoms in unit cell ntypat=number of atom types outd2=option for the output of the 2nd-order matrix : if outd2=1, non-stationary part if outd2=2, stationary part. pawbec= flag for the computation of frozen part of Born Effective Charges (PAW only) prtbbb=if 1, print the band-by-band decomposition, otherwise, prtbbb=0 rfpert(mpert)=define the perturbations rprimd(3,3)=dimensional primitive translations (bohr) typat(natom)=integer label of each type of atom (1,2,...) ucvol=unit cell volume usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use zion(ntypat)=charge corresponding to the atom type
OUTPUT
carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian 2DTE matrix has been calculated correctly ; 0 otherwise ) d2cart(2,3,mpert,3,mpert)= dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)= band by band decomposition of Born effective charges (calculated from phonon-type perturbation) in cartesian coordinates d2matr(2,3,mpert,3,mpert)=2nd-order matrix (masses non included, no cartesian coordinates : simply second derivatives)
SOURCE
3446 subroutine dfpt_gatherdy(asr,becfrnl,berryopt,blkflg,carflg,chneut,dyew,dyfrwf,dyfrx1,& 3447 & dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,& 3448 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 3449 & gprimd,mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,prtbbb,& 3450 & rfpert,rprimd,typat,ucvol,usevdw,zion) 3451 3452 !Arguments ------------------------------- 3453 !scalars 3454 integer,intent(in) :: asr,berryopt,chneut,dyfr_cplex,dyfr_nondiag,mband,mpert,natom,ntypat,outd2 3455 integer,intent(in) :: pawbec,pawpiezo,prtbbb,usevdw 3456 real(dp),intent(in) :: ucvol 3457 !arrays 3458 integer,intent(in) :: rfpert(mpert),typat(natom) 3459 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 3460 integer,intent(out) :: carflg(3,mpert,3,mpert) 3461 real(dp),intent(in) :: becfrnl(3,natom,3*pawbec) 3462 real(dp),intent(in) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb) 3463 real(dp),intent(in) :: d2nfr(2,3,mpert,3,mpert),dyew(2,3,natom,3,natom) 3464 real(dp),intent(in) :: dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag) 3465 real(dp),intent(in) :: dyfrx1(2,3,natom,3,natom),dyvdw(2,3,natom,3,natom*usevdw) 3466 real(dp),intent(in) :: eltcore(6,6),elteew(6+3*natom,6) 3467 real(dp),intent(in) :: eltfrhar(6,6),eltfrkin(6,6),eltfrloc(6+3*natom,6) 3468 real(dp),intent(in) :: eltfrnl(6+3*natom,6),eltfrxc(6+3*natom,6) 3469 real(dp),intent(in) :: eltvdw(6+3*natom,6*usevdw),gprimd(3,3) 3470 real(dp),intent(in) :: piezofrnl(6,3*pawpiezo),rprimd(3,3),zion(ntypat) 3471 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert) 3472 real(dp),intent(out) :: d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb) 3473 real(dp),intent(out) :: d2matr(2,3,mpert,3,mpert) 3474 3475 !Local variables ------------------------- 3476 !scalars 3477 integer :: iband,iblok,idir,idir1,idir2,ii,ipert,ipert1,ipert2 3478 integer :: jj,nblok,selectz 3479 character(len=500) :: message 3480 !arrays 3481 integer :: flg1(3),flg2(3) 3482 real(dp) :: vec1(3),vec2(3) 3483 ! real(dp) :: ter(3,3) ! this variable appears commented out below 3484 real(dp),allocatable :: d2tmp(:,:,:,:,:),d2work(:,:,:,:,:),elfrtot(:,:) 3485 3486 ! ********************************************************************* 3487 3488 3489 if(outd2/=3)then 3490 3491 ! Initialise the 2nd-derivative matrix 3492 d2matr(:,:,:,:,:)=0.0_dp 3493 3494 ! Add the non-frozen-part, the 3495 ! Ewald part and the xc1 part of the frozen-wf part 3496 ! Add the vdw part (if any) 3497 do ipert2=1,mpert 3498 do idir2=1,3 3499 do ipert1=1,mpert 3500 do idir1=1,3 3501 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 3502 do ii=1,2 3503 d2matr(ii,idir1,ipert1,idir2,ipert2)=& 3504 & d2nfr(ii,idir1,ipert1,idir2,ipert2) 3505 if(ipert1<=natom .and. ipert2<=natom) then 3506 d2matr(ii,idir1,ipert1,idir2,ipert2)=& 3507 & d2matr(ii,idir1,ipert1,idir2,ipert2)+& 3508 & dyew(ii,idir1,ipert1,idir2,ipert2) +& 3509 & dyfrx1(ii,idir1,ipert1,idir2,ipert2) 3510 if (usevdw==1) then 3511 d2matr(ii,idir1,ipert1,idir2,ipert2)=& 3512 & d2matr(ii,idir1,ipert1,idir2,ipert2)+& 3513 & dyvdw(ii,idir1,ipert1,idir2,ipert2) 3514 end if 3515 end if 3516 end do 3517 end if 3518 end do 3519 end do 3520 end do 3521 end do 3522 3523 ! Add the frozen-wavefunction part 3524 if (dyfr_nondiag==0) then 3525 do ipert2=1,natom 3526 do idir2=1,3 3527 do idir1=1,3 3528 if( blkflg(idir1,ipert2,idir2,ipert2)==1 ) then 3529 d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)=& 3530 & d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)& 3531 & +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert2,1) 3532 end if 3533 end do 3534 end do 3535 end do 3536 else 3537 do ipert2=1,natom 3538 do ipert1=1,natom 3539 do idir2=1,3 3540 do idir1=1,3 3541 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 3542 d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)=& 3543 & d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)& 3544 & +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert1,ipert2) 3545 end if 3546 end do 3547 end do 3548 end do 3549 end do 3550 end if 3551 3552 ! Add the frozen-wavefunction part of Born Effective Charges 3553 if (pawbec==1) then 3554 ipert2=natom+2 3555 do idir2=1,3 ! Direction of electric field 3556 do ipert1=1,natom ! Atom 3557 do idir1=1,3 ! Direction of atom 3558 if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 3559 d2matr(1,idir1,ipert1,idir2,ipert2)=& 3560 & d2matr(1,idir1,ipert1,idir2,ipert2)+becfrnl(idir1,ipert1,idir2) 3561 end if 3562 if(blkflg(idir2,ipert2,idir1,ipert1)==1 ) then 3563 d2matr(1,idir2,ipert2,idir1,ipert1)=& 3564 & d2matr(1,idir2,ipert2,idir1,ipert1)+becfrnl(idir1,ipert1,idir2) 3565 end if 3566 end do 3567 end do 3568 end do 3569 end if 3570 3571 ! Section for piezoelectric tensor (from electric field response only for PAW) 3572 if(rfpert(natom+2)==1.and.pawpiezo==1) then 3573 ipert2=natom+2 3574 do idir2=1,3 ! Direction of electric field 3575 do ipert1=natom+3,natom+4 ! Strain 3576 do idir1=1,3 3577 ii=idir1+3*(ipert1-natom-3) 3578 if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 3579 d2matr(1,idir1,ipert1,idir2,ipert2)=& 3580 & d2nfr(1,idir1,ipert1,idir2,ipert2)+piezofrnl(ii,idir2) 3581 end if 3582 end do 3583 end do 3584 end do 3585 end if 3586 3587 ! Section for strain perturbation 3588 if(rfpert(natom+3)==1 .or. rfpert(natom+4)==1) then 3589 ! Make sure relevant columns of output are nulled 3590 d2matr(:,:,:,:,natom+3:natom+4)=0.0_dp 3591 ! Accumulate all frozen parts of the elastic tensor 3592 ABI_MALLOC(elfrtot,(6+3*natom,6)) 3593 elfrtot(:,:)=elteew(:,:)+eltfrloc(:,:)+eltfrnl(:,:)+eltfrxc(:,:) 3594 elfrtot(1:6,1:6)=elfrtot(1:6,1:6)+eltcore(:,:)+eltfrhar(:,:)+eltfrkin(:,:) 3595 if (usevdw==1) elfrtot(:,:)=elfrtot(:,:)+eltvdw(:,:) 3596 3597 do ipert2=natom+3,natom+4 3598 do idir2=1,3 3599 ! Internal strain components first 3600 do ipert1=1,natom 3601 do idir1=1,3 3602 if( blkflg(1,ipert1,idir2,ipert2)==1 ) then 3603 ii=idir1+6+3*(ipert1-1) 3604 jj=idir2+3*(ipert2-natom-3) 3605 d2matr(1,idir1,ipert1,idir2,ipert2)=& 3606 & d2nfr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj) 3607 end if 3608 end do 3609 end do 3610 ! Now, electric field - strain mixed derivative (piezoelectric tensor) 3611 ipert1=natom+2 3612 do idir1=1,3 3613 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 3614 d2matr(1,idir1,ipert1,idir2,ipert2)=& 3615 & d2nfr(1,idir1,ipert1,idir2,ipert2) 3616 if (pawpiezo==1) then 3617 ii=idir2+3*(ipert2-natom-3) 3618 d2matr(1,idir1,ipert1,idir2,ipert2)=& 3619 & d2matr(1,idir1,ipert1,idir2,ipert2)+piezofrnl(ii,idir1) 3620 end if 3621 end if 3622 end do 3623 ! Now, strain-strain 2nd derivatives 3624 do ipert1=natom+3,natom+4 3625 do idir1=1,3 3626 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 3627 ii=idir1+3*(ipert1-natom-3) 3628 jj=idir2+3*(ipert2-natom-3) 3629 d2matr(1,idir1,ipert1,idir2,ipert2)=& 3630 & d2nfr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj) 3631 end if 3632 end do 3633 end do 3634 end do 3635 end do 3636 ABI_FREE(elfrtot) 3637 end if 3638 ! End section for strain perturbation 3639 3640 ! The second-order matrix has been computed. 3641 3642 ! Filter now components smaller in absolute value than 1.0d-20, 3643 ! for automatic testing reasons 3644 do ipert2=1,mpert 3645 do idir2=1,3 3646 do ipert1=1,mpert 3647 do idir1=1,3 3648 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 3649 do ii=1,2 3650 if( d2matr(ii,idir1,ipert1,idir2,ipert2)**2 < 1.0d-40)then 3651 d2matr(ii,idir1,ipert1,idir2,ipert2)=zero 3652 end if 3653 end do 3654 end if 3655 end do 3656 end do 3657 end do 3658 end do 3659 3660 ! Cartesian coordinates transformation 3661 iblok=1 ; nblok=1 3662 3663 ! In the case of finite electric field, the convention for the 3664 ! direction of the electric field perturbation was NOT the usual convention ... 3665 ! So, there is a transformation to the usual conventions 3666 ! to be done first ... 3667 if((berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17 ) & 3668 & .and. minval(abs(blkflg(:,natom+2,:,natom+2)))/=0)then !!HONG need to check for fixed D and E calculation 3669 if(minval(abs(blkflg(:,natom+2,:,natom+2)-1))/=0)then 3670 write(message,'(5a)')& 3671 & ' In case of finite electric field, and electric field perturbation,',ch10,& 3672 & ' the three directions for the perturbations must be treated.',ch10,& 3673 & ' Action : set idir to 1 1 1, or forget about finite electric field.' 3674 ABI_ERROR(message) 3675 end if 3676 do ipert=1,mpert 3677 do idir=1,3 3678 do ii=1,2 3679 vec1(:)=d2matr(ii,idir,ipert,:,natom+2) 3680 flg1(:)=blkflg(idir,ipert,:,natom+2) 3681 call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2) 3682 d2matr(ii,idir,ipert,:,natom+2)=vec2(:)*two_pi 3683 blkflg(idir,ipert,:,natom+2)=flg2(:) 3684 end do 3685 end do 3686 end do 3687 do ipert=1,mpert 3688 do idir=1,3 3689 do ii=1,2 3690 vec1(:)=d2matr(ii,:,natom+2,idir,ipert) 3691 flg1(:)=blkflg(:,natom+2,idir,ipert) 3692 call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2) 3693 d2matr(ii,:,natom+2,idir,ipert)=vec2(:)*two_pi 3694 blkflg(:,natom+2,idir,ipert)=flg2(:) 3695 end do 3696 end do 3697 end do 3698 ! Also to be done, a change of sign, that I do not understand (XG071110) 3699 ! Perhaps due to d/dk replacing id/dk ? ! 3700 d2matr(:,:,natom+2,:,natom+2)=-d2matr(:,:,natom+2,:,natom+2) 3701 end if 3702 3703 call cart29(blkflg,d2matr,carflg,d2cart,& 3704 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion) 3705 3706 ! Band by band decomposition of the Born effective charges 3707 if(prtbbb==1)then 3708 ABI_MALLOC(d2work,(2,3,mpert,3,mpert)) 3709 ABI_MALLOC(d2tmp,(2,3,mpert,3,mpert)) 3710 do iband=1,mband 3711 d2work(:,:,:,:,:)=0.0_dp 3712 d2tmp(:,:,:,:,:)=0.0_dp 3713 d2work(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband) 3714 call cart29(blkflg,d2work,carflg,d2tmp,& 3715 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion) 3716 3717 ! Remove the ionic part 3718 do ipert1=1,natom 3719 do idir1=1,3 3720 d2tmp(1,idir1,natom+2,idir1,ipert1) = & 3721 & d2tmp(1,idir1,natom+2,idir1,ipert1) - zion(typat(ipert1)) 3722 end do 3723 end do 3724 3725 d2cart_bbb(:,:,:,:,iband,iband)=d2tmp(:,:,natom+2,:,:) 3726 3727 end do 3728 ABI_FREE(d2tmp) 3729 ABI_FREE(d2work) 3730 end if ! prtbbb==1 3731 3732 ! 3733 ! Now, the cartesian elements are ready for output 3734 ! carflg give the information on their correctness 3735 end if 3736 3737 ! Imposition of the ASR on the analytical part of the DynMat 3738 ! Assume that if asr/=0, the whole cartesian matrix is correct 3739 if(asr/=0)then 3740 3741 ABI_MALLOC(d2work,(2,3,mpert,3,mpert)) 3742 call asria_calc(asr,d2work,d2cart,mpert,natom) 3743 ! The following line imposes ASR: 3744 call asria_corr(asr,d2work,d2cart,mpert,natom) 3745 3746 ABI_FREE(d2work) 3747 3748 ! Imposition of the charge neutrality on the effective charges. 3749 if(rfpert(natom+2)==1)then 3750 selectz=0 3751 call chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion) 3752 end if 3753 3754 end if 3755 3756 !Additional operations on cartesian strain derivatives 3757 if(rfpert(natom+3)==1 .or. rfpert(natom+4)==1) then 3758 ! Impose zero-net-force condition on internal strain tensor 3759 do ipert2=natom+3,natom+4 3760 do idir2=1,3 3761 vec1(:)=0.0_dp 3762 do ipert1=1,natom 3763 do idir1=1,3 3764 if(carflg(idir1,ipert1,idir2,ipert2)==1) then 3765 vec1(idir1)=vec1(idir1)+d2cart(1,idir1,ipert1,idir2,ipert2) 3766 end if 3767 end do 3768 end do 3769 vec1(:)=vec1(:)/dble(natom) 3770 do ipert1=1,natom 3771 do idir1=1,3 3772 if(carflg(idir1,ipert1,idir2,ipert2)==1) then 3773 ! Note minus sign to convert gradients to forces 3774 d2cart(1,idir1,ipert1,idir2,ipert2)=& 3775 & -(d2cart(1,idir1,ipert1,idir2,ipert2)-vec1(idir1)) 3776 end if 3777 end do 3778 end do 3779 end do 3780 end do 3781 ! Divide strain 2nd deriviative by ucvol to give elastic tensor 3782 do ipert2=natom+3,natom+4 3783 do idir2=1,3 3784 do ipert1=natom+3,natom+4 3785 do idir1=1,3 3786 if(carflg(idir1,ipert1,idir2,ipert2)==1) then 3787 d2cart(1,idir1,ipert1,idir2,ipert2)=& 3788 & d2cart(1,idir1,ipert1,idir2,ipert2)/ucvol 3789 end if 3790 end do 3791 end do 3792 end do 3793 end do 3794 end if 3795 3796 !calculate Born effective charges from electric field perturbation 3797 !do ipert1=1,natom 3798 !ter(:,:)=zero 3799 !do idir1=1,3 3800 !do ii=1,3 3801 !do jj=1,3 3802 !if(abs(gprimd(idir1,ii))>1.0d-10)then 3803 !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,idir1,natom+2,jj,ipert1)*gprimd(jj,ii) 3804 !endif 3805 !enddo 3806 !enddo 3807 !add zion to bec 3808 !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1)) 3809 !enddo 3810 !d2cart(1,:,ipert1,:,natom+2)=ter(:,:) 3811 !enddo 3812 !carflg(:,1:natom,:,natom+2)=1 3813 3814 !Born effective charges from phonon perturbation 3815 !do ipert1=1,natom 3816 !ter(:,:)=zero 3817 !do idir1=1,3 3818 !do ii=1,3 3819 !do jj=1,3 3820 !if(abs(gprimd(idir1,ii))>1.0d-10)then 3821 !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,jj,ipert1,idir1,natom+2)*gprimd(jj,ii) 3822 !endif 3823 !enddo 3824 !enddo 3825 !! add zion to bec 3826 !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1)) 3827 !enddo 3828 !d2cart(1,:,natom+2,:,ipert1)=ter(:,:) 3829 !enddo 3830 !carflg(:,natom+2,:,1:natom)=1 3831 3832 3833 !!Dielectric constant 3834 !do ii=1,3 3835 !do jj=1,3 3836 !ter(ii,jj)=d2nfr(1,ii,natom+2,jj,natom+2) 3837 !end do 3838 !end do 3839 !ter(:,:)=pi*four*ter(:,:)/ucvol 3840 ! 3841 !do ii=1,3 3842 !ter(ii,ii)=ter(ii,ii)+one 3843 !end do 3844 !d2cart(1,:,natom+2,:,natom+2)=ter(:,:) 3845 !carflg(:,natom+2,:,1:natom+2)=1 3846 3847 !DEBUG 3848 !Debugging, but needed to make the stuff work on the IBM Dirac ? ! 3849 !write(std_out,*)' d2cart ' 3850 !ipert2=natom+2 3851 !do idir2=1,3 3852 !ipert1=natom+2 3853 !do idir1=1,3 3854 !write(std_out,'(5i4,2d20.10)' )idir1,ipert1,idir2,ipert2,& 3855 !& carflg(idir1,ipert1,idir2,ipert2),& 3856 !& d2cart(1,idir1,ipert1,idir2,ipert2),& 3857 !& d2cart(2,idir1,ipert1,idir2,ipert2) 3858 !end do 3859 !end do 3860 !ENDDEBUG 3861 3862 3863 end subroutine dfpt_gatherdy
ABINIT/m_respfn_driver [ Modules ]
NAME
m_respfn_driver
FUNCTION
Subdriver for DFPT calculations.
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (XG, DRH, MT, MKV, GA) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_respfn_driver 23 24 use defs_basis 25 use defs_wvltypes 26 use m_efmas_defs 27 use m_abicore 28 use m_xmpi 29 use m_exit 30 use m_wffile 31 use m_errors 32 use m_ebands 33 use m_results_respfn 34 use m_hdr 35 use m_crystal 36 use m_xcdata 37 use m_dtset 38 use m_dtfil 39 40 use defs_datatypes, only : pseudopotential_type, ebands_t 41 use defs_abitypes, only : MPI_type 42 use m_time, only : timab 43 use m_fstrings, only : strcat 44 use m_symtk, only : matr3inv, littlegroup_q, symmetrize_xred 45 use m_fft, only : zerosym, fourdp 46 use m_kpts, only : symkchk 47 use m_geometry, only : irreducible_set_pert, symredcart 48 use m_dynmat, only : chkph3, d2sym3, q0dy3_apply, q0dy3_calc, wings3, dfpt_phfrq, sytens, sylwtens, dfpt_prtph, & 49 asria_calc, asria_corr, cart29, cart39, chneu9, dfpt_sydy 50 use m_ddb, only : ddb_type 51 use m_ddb_hdr, only : ddb_hdr_type 52 use m_occ, only : newocc 53 use m_efmas, only : efmasdeg_free_array, efmasval_free_array 54 use m_wfk, only : wfk_read_eigenvalues, wfk_read_my_kptbands 55 use m_ioarr, only : read_rhor 56 use m_pawang, only : pawang_type 57 use m_pawrad, only : pawrad_type 58 use m_pawtab, only : pawtab_type, pawtab_get_lsize 59 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify 60 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify 61 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 62 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, & 63 pawrhoij_bcast, pawrhoij_nullify, pawrhoij_inquire_dim, & 64 pawrhoij_print_rhoij, pawrhoij_io 65 66 use m_pawdij, only : pawdij, symdij, pawdij_print_dij 67 use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy 68 use m_paw_finegrid,only : pawexpiqr 69 use m_pawxc, only : pawxc_get_nkxc 70 use m_paw_dmft, only : paw_dmft_type 71 use m_paw_sphharm, only : setsym_ylm 72 use m_paw_nhat, only : nhatgrid,pawmknhat 73 use m_paw_tools, only : chkpawovlp 74 use m_paw_denpot, only : pawdenpot 75 use m_paw_init, only : pawinit,paw_gencond 76 use m_kg, only : getcut, getph, kpgio 77 use m_eig2d, only : eig2tot, elph2_fanddw 78 use m_inwffil, only : inwffil 79 use m_spacepar, only : hartre, setsym 80 use m_mkrho, only : mkrho 81 use m_vdw_dftd2, only : vdw_dftd2 82 use m_vdw_dftd3, only : vdw_dftd3 83 use m_initylmg, only : initylmg 84 use m_pspini, only : pspini 85 use m_atm2fft, only : atm2fft 86 use m_dfpt_loopert,only : dfpt_looppert, eigen_meandege 87 use m_rhotoxc, only : rhotoxc 88 use m_drivexc, only : check_kxc 89 use m_xc_tb09, only : xc_tb09_update_c 90 use m_mklocl, only : mklocl, mklocl_recipspace 91 use m_common, only : setup1, prteigrs 92 use m_fourier_interpol, only : transgrid 93 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 94 use m_paw_occupancies, only : initrhoij 95 use m_paw_correlations,only : pawpuxinit 96 use m_mkcore, only : mkcore, dfpt_mkcore 97 use m_dfpt_elt, only : dfpt_eltfrxc, dfpt_eltfrloc, dfpt_eltfrkin, dfpt_eltfrhar, elt_ewald, dfpt_ewald 98 use m_d2frnl, only : d2frnl 99 100 #if defined HAVE_GPU_CUDA 101 use m_alloc_hamilt_gpu 102 #endif 103 104 implicit none 105 106 private
m_respfn_driver/respfn [ Functions ]
[ Top ] [ m_respfn_driver ] [ Functions ]
NAME
respfn
FUNCTION
Primary routine for conducting DFT calculations of Response functions.
INPUTS
codvsn=code version cpui=initial cpu time dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset | mband=maximum number of bands | mgfft=maximum single fft dimension | mkmem=Number of k points treated by this node | mpw=maximum number of planewaves in basis sphere (large number) | natom=number of atoms in unit cell | nfft=(effective) number of FFT grid points (for this processor) | nkpt=number of k points | nspden=number of spin-density components | nsppol=number of channels for spin-polarization (1 or 2) | nsym=number of symmetry elements in space group mkmems(3)=array containing the tree values of mkmem (see above) (k-GS, k+q-GS and RF) mpi_enreg=information about MPI parallelization npwtot(nkpt)=number of planewaves in basis and boundary at each k point xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
etotal=total energy (sum of 7 or 8 contributions) (hartree)
SIDE EFFECTS
iexit=index of "exit" on first line of file (0 if not found)
occ(mband*nkpt*nsppol)=occup number for each band (often 2) at each k point
Occupations number may have been read from a previous dataset...
pawang <type(pawang_type)>=paw angular mesh and related data
pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
Some dimensions in pawrad have been set in driver.f
pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
Some dimensions in pawtab have been set in driver.f
psps <type(pseudopotential_type)>=variables related to pseudopotentials
Before entering the first time in respfn, a significant part of psps
has been initialized: the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,
mpsso,mgrid,ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp
All the remaining components of psps are to be initialized in the call
to pspini. The next time the code enters respfn, psps might be identical
to the one of the previous dtset, in which case, no reinitialisation
is scheduled in pspini.f .
results_respfn <type(results_respfn_type)>=stores some results of respfn calls
NOTES
USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut)
for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ...
are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg)
for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ...
Total density, potentials, ...
are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used.
It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf)
are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
189 subroutine respfn(codvsn,cpui,dtfil,dtset,etotal,iexit,& 190 & mkmems,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,results_respfn,xred) 191 192 !Arguments ------------------------------------ 193 integer,intent(inout) :: iexit 194 real(dp),intent(in) :: cpui 195 real(dp),intent(inout) :: etotal !vz_i 196 character(len=8),intent(in) :: codvsn 197 type(MPI_type),intent(inout) :: mpi_enreg 198 type(datafiles_type),intent(in) :: dtfil 199 type(dataset_type),intent(in) :: dtset 200 type(pawang_type),intent(inout) :: pawang 201 type(pseudopotential_type),intent(inout) :: psps 202 integer,intent(in) :: mkmems(3) 203 integer,intent(inout) :: npwtot(dtset%nkpt) 204 real(dp),intent(inout) :: xred(3,dtset%natom) 205 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 206 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 207 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 208 type(results_respfn_type),intent(inout) :: results_respfn 209 210 !Local variables------------------------------- 211 integer,parameter :: formeig=0,level=10 212 integer,parameter :: response=1,syuse=0,master=0,cplex1=1 213 integer :: nk3xc 214 integer :: analyt,ask_accurate,asr,bantot,bdeigrf,chneut,coredens_method,coretau_method,cplex,cplex_rhoij 215 !integer :: nkpt_eff, band_index, ikpt, isppol, nkpt_max, nband_k, 216 integer :: dim_eig2nkq,dim_eigbrd,dyfr_cplex,dyfr_nondiag,gnt_option 217 integer :: gscase,has_dijnd,has_diju,has_vhartree,has_kxc,iatom,iatom_tot,iband,idir,ider,ierr,ifft,ii,indx 218 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 219 integer :: initialized,ipert,ipert2,ireadwf0,iscf,iscf_eff,ispden,isym 220 integer :: itypat,izero,mcg,me,mgfftf,mk1mem,mkqmem,mpert,mu 221 integer :: my_natom,n1,natom,n3xccc,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim 222 integer :: nkpt_rbz,nkxc,nkxc1,nspden_rhoij,ntypat,nzlmopt,openexit 223 integer :: optcut,option,optgr0,optgr1,optgr2,optorth,optrad 224 integer :: optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv 225 integer :: outd2,pawbec,pawpiezo,prtbbb,psp_gencond,qzero,rdwrpaw 226 integer :: rfddk,rfelfd,rfphon,rfstrs,rfuser,rf2_dkdk,rf2_dkde,rfmagn 227 integer :: spaceworld,sumg0,sz1,sz2,tim_mkrho,timrev,usecprj,usevdw 228 integer :: usexcnhat,use_sym,vloc_method,zero_by_symm 229 logical :: has_full_piezo,has_allddk,is_dfpt=.true.,non_magnetic_xc 230 logical :: paral_atom,qeq0,use_nhat_gga,call_pawinit 231 real(dp) :: boxcut,compch_fft,compch_sph,cpus,ecore,ecut_eff,ecutdg_eff,ecutf 232 real(dp) :: eei,eew,ehart,eii,ek,enl,entropy,enxc 233 real(dp) :: epaw,epawdc,etot,evdw,fermie,fermih,gsqcut,gsqcut_eff,gsqcutc_eff,qphnrm,residm 234 real(dp) :: ucvol,vxcavg 235 character(len=500) :: message 236 type(ebands_t) :: bstruct 237 type(hdr_type) :: hdr,hdr_fine,hdr0,hdr_den 238 type(ddb_type) :: ddb 239 type(ddb_hdr_type) :: ddb_hdr 240 type(paw_dmft_type) :: paw_dmft 241 type(pawfgr_type) :: pawfgr 242 type(wvl_data) :: wvl 243 type(xcdata_type) :: xcdata 244 integer :: ddkfil(3),ngfft(18),ngfftf(18),rfdir(3),rf2_dirs_from_rfpert_nl(3,3) 245 integer,allocatable :: atindx(:),atindx1(:),blkflg(:,:,:,:),blkflgfrx1(:,:,:,:),blkflg1(:,:,:,:) 246 integer,allocatable :: blkflg2(:,:,:,:),carflg(:,:,:,:),clflg(:,:),indsym(:,:,:) 247 integer,allocatable :: irrzon(:,:,:),kg(:,:),l_size_atm(:),nattyp(:),npwarr(:) 248 integer,allocatable :: pertsy(:,:),rfpert(:) 249 integer,allocatable :: rfpert_lw(:,:,:,:,:,:),rfpert_nl(:,:,:,:,:,:),symq(:,:,:),symrec(:,:,:) 250 logical,allocatable :: distrb_flags(:,:,:) 251 real(dp) :: dum_gauss(0),dum_dyfrn(0),dum_dyfrv(0),dum_eltfrxc(0) 252 real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0) 253 real(dp) :: dummy6(6),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3),qphon(3) 254 real(dp) :: dummy_in(0) 255 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0) 256 real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),strn_dummy6(6),strv_dummy6(6),strsxc(6),tsec(2) 257 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 258 real(dp),allocatable :: becfrnl(:,:,:),cg(:,:),d2bbb(:,:,:,:,:,:),d2cart(:,:,:,:,:) 259 real(dp),allocatable :: d2cart_bbb(:,:,:,:,:,:),d2eig0(:,:,:,:,:) 260 real(dp),allocatable :: d2k0(:,:,:,:,:),d2lo(:,:,:,:,:),d2loc0(:,:,:,:,:) 261 real(dp),allocatable :: d2matr(:,:,:,:,:),d2nfr(:,:,:,:,:),d2nl(:,:,:,:,:),d2ovl(:,:,:,:,:) 262 real(dp),allocatable :: d2nl0(:,:,:,:,:),d2nl1(:,:,:,:,:),d2tmp(:,:,:,:,:),d2vn(:,:,:,:,:) 263 real(dp),allocatable :: displ(:),doccde(:) 264 real(dp),allocatable :: dyew(:,:,:,:,:),dyewq0(:,:,:),dyfrlo(:,:,:),dyfrlo_indx(:,:,:) 265 real(dp),allocatable :: dyfrnl(:,:,:,:,:),dyfrwf(:,:,:,:,:),dyfrx1(:,:,:,:,:),dyvdw(:,:,:,:,:) 266 real(dp),allocatable :: dyfrx2(:,:,:),eigen0(:),eigval(:),eigvec(:) 267 real(dp),allocatable :: eig2nkq(:,:,:,:,:,:,:),eigbrd(:,:,:,:,:,:,:) 268 real(dp),allocatable :: eigen_fan(:),eigen_ddw(:),eigen_fanddw(:) 269 real(dp),allocatable :: eigen_fan_mean(:),eigen_ddw_mean(:) 270 real(dp),allocatable :: eltcore(:,:),elteew(:,:),eltfrhar(:,:),eltfrkin(:,:) 271 real(dp),allocatable :: eltfrloc(:,:),eltfrnl(:,:),eltfrxc(:,:),eltvdw(:,:),grtn_indx(:,:) 272 real(dp),allocatable :: grxc(:,:),kxc(:,:),nhat(:,:),nhatgr(:,:,:) 273 real(dp),allocatable :: ph1d(:,:),ph1df(:,:),phfrq(:),phnons(:,:,:),piezofrnl(:,:) 274 real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:) 275 real(dp),allocatable :: symrel_cart(:,:,:),taug(:,:),taur(:,:) 276 real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:) 277 real(dp),allocatable :: vxc(:,:),vxctau(:,:,:),xccc3d(:),xcctau3d(:),ylm(:,:),ylmgr(:,:,:) 278 real(dp),pointer :: eigenq_fine(:,:,:),eigen1_pert(:,:,:) 279 real(dp),allocatable :: eigen0_pert(:),eigenq_pert(:),occ_rbz_pert(:) 280 type(efmasdeg_type),allocatable :: efmasdeg(:) 281 type(efmasval_type),allocatable :: efmasval(:,:) 282 type(paw_an_type),allocatable :: paw_an(:) 283 type(paw_ij_type),allocatable :: paw_ij(:) 284 type(pawfgrtab_type),allocatable,save :: pawfgrtab(:) 285 type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:) 286 287 ! *********************************************************************** 288 289 DBG_ENTER("COLL") 290 291 call timab(132,1,tsec) 292 call timab(133,1,tsec) 293 294 !Some data for parallelism 295 296 my_natom=mpi_enreg%my_natom 297 paral_atom=(my_natom/=dtset%natom) 298 !Define FFT grid(s) sizes (be careful !) 299 !See NOTES in the comments at the beginning of this file. 300 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 301 302 !Structured debugging if dtset%prtvol==-level 303 if(dtset%prtvol==-level)then 304 write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' respfn : enter , debug mode ' 305 call wrtout(std_out,message) 306 end if 307 308 !Option input variables 309 iscf=dtset%iscf 310 311 !Respfn input variables 312 asr=dtset%asr ; chneut=dtset%chneut ; rfdir(1:3)=dtset%rfdir(1:3) 313 rfddk=dtset%rfddk ; rfelfd=dtset%rfelfd ; rfmagn=dtset%rfmagn 314 rfphon=dtset%rfphon ; rfstrs=dtset%rfstrs 315 rfuser=dtset%rfuser ; rf2_dkdk=dtset%rf2_dkdk ; rf2_dkde=dtset%rf2_dkde 316 317 pawbec=0 ; if(psps%usepaw==1.and.(rfphon==1.or.(rfelfd==1.or.rfelfd==3))) pawbec=1 318 pawpiezo=0; if(psps%usepaw==1.and.(rfstrs/=0.or.(rfelfd==1.or.rfelfd==3))) pawpiezo=1 319 !AM 10152015 -- WARNING --- the full calculation of the piezoelectric tensor 320 !from electric field perturbation is only available 321 !if nsym/=1 (strain perturbation is not symmetrized): 322 has_full_piezo=.False. ; if(pawpiezo==1.and.dtset%nsym==1) has_full_piezo=.True. 323 usevdw=0;if (dtset%vdw_xc>=5.and.dtset%vdw_xc<=7) usevdw=1 324 !mkmem variables (mkmem is already argument) 325 mkqmem=mkmems(2) ; mk1mem=mkmems(3) 326 327 ntypat=psps%ntypat 328 natom=dtset%natom 329 nfftot=product(ngfft(1:3)) 330 nfftotf=product(ngfftf(1:3)) 331 332 !LIKELY TO BE TAKEN AWAY 333 initialized=0 334 ek=zero ; ehart=zero ; enxc=zero ; eei=zero ; enl=zero 335 eii=zero ; eew=zero ; ecore=zero 336 337 !Set up for iterations 338 call setup1(dtset%acell_orig(1:3,1),bantot,dtset,& 339 ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,& 340 ngfftf,ngfft,dtset%nkpt,dtset%nsppol,& 341 response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw) 342 343 !In some cases (e.g. getcell/=0), the plane wave vectors have 344 ! to be generated from the original simulation cell 345 rprimd_for_kg=rprimd 346 if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=dtset%rprimd_orig(:,:,1) 347 call matr3inv(rprimd_for_kg,gprimd_for_kg) 348 gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg) 349 350 !Define the set of admitted perturbations 351 ! Note that we have a global parameter (mpert=natom+MPERT_MAX) 352 ! with MPERT_MAX=8, but we use a smaller value here. 353 mpert=natom+7 354 if (rf2_dkdk>0.or.rf2_dkde>0) mpert=natom+11 355 356 !Initialize the list of perturbations rfpert 357 ABI_MALLOC(rfpert,(mpert)) 358 rfpert(:)=0 359 if(rfphon==1)rfpert(dtset%rfatpol(1):dtset%rfatpol(2))=1 360 361 if(rfddk==1)rfpert(natom+1)=1 362 if(rfddk==2)rfpert(natom+6)=1 363 364 if(rf2_dkdk/=0)rfpert(natom+10)=1 365 if(rf2_dkde/=0)rfpert(natom+11)=1 366 367 if(rfelfd==1.or.rfelfd==2)rfpert(natom+1)=1 368 if(rfelfd==1.or.rfelfd==3)rfpert(natom+2)=1 369 370 if(rfstrs==1.or.rfstrs==3)rfpert(natom+3)=1 371 if(rfstrs==2.or.rfstrs==3)rfpert(natom+4)=1 372 373 if(rfuser==1.or.rfuser==3)rfpert(natom+6)=1 374 if(rfuser==2.or.rfuser==3)rfpert(natom+7)=1 375 376 if(rfmagn==1) rfpert(natom+5)=1 377 378 qeq0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2<1.d-14) 379 380 !Init spaceworld 381 spaceworld=mpi_enreg%comm_cell 382 me = xmpi_comm_rank(spaceworld) 383 384 !Set up the basis sphere of planewaves 385 ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem)) 386 ABI_MALLOC(npwarr,(dtset%nkpt)) 387 call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg,& 388 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,& 389 & dtset%nsppol) 390 391 !Set up the Ylm for each k point 392 ABI_MALLOC(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 393 if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) then 394 ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,9,psps%mpsang*psps%mpsang*psps%useylm)) 395 else 396 ABI_MALLOC(ylmgr,(0,0,psps%useylm)) 397 end if 398 if (psps%useylm==1) then 399 option=0 400 if (rfstrs/=0.or.pawbec==1.or.pawpiezo==1.or.dtset%efmas>0) option=2 401 call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,& 402 & npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr) 403 end if 404 405 call timab(133,2,tsec) 406 call timab(134,1,tsec) 407 408 !Open and read pseudopotential files 409 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,& 410 & psps,rprimd,comm_mpi=mpi_enreg%comm_cell) 411 412 call timab(134,2,tsec) 413 call timab(135,1,tsec) 414 415 !Initialize band structure datatype 416 bstruct = ebands_from_dtset(dtset, npwarr) 417 418 !Initialize PAW atomic occupancies 419 if (psps%usepaw==1) then 420 ABI_MALLOC(pawrhoij,(my_natom)) 421 call pawrhoij_nullify(pawrhoij) 422 call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu, & 423 & my_natom,natom,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,& 424 & pawrhoij,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,& 425 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 426 else 427 ABI_MALLOC(pawrhoij,(0)) 428 end if 429 430 !Initialize header 431 gscase=0 432 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr, & 433 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 434 435 !Update header, with evolving variables, when available 436 !Here, rprimd, xred and occ are available 437 etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm 438 439 !If parallelism over atom, hdr is distributed 440 call hdr%update(bantot,etot,fermie,fermih,& 441 residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), & 442 comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 443 444 !Clean band structure datatype (should use it more in the future !) 445 call ebands_free(bstruct) 446 447 !Initialize wavefunction files and wavefunctions. 448 ireadwf0=1 449 450 451 mcg=dtset%mpw*dtset%nspinor*dtset%mband_mem*dtset%mkmem*dtset%nsppol 452 ABI_MALLOC(cg,(2,mcg)) 453 !ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr) 454 455 ABI_MALLOC(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol)) 456 eigen0(:)=zero ; ask_accurate=1 457 optorth=0 458 459 ! Initialize the wave function type and read GS WFK 460 ABI_MALLOC(distrb_flags,(dtset%nkpt,dtset%mband,dtset%nsppol)) 461 distrb_flags = (mpi_enreg%proc_distrb == mpi_enreg%me_kpt) 462 call wfk_read_my_kptbands(dtfil%fnamewffk, distrb_flags, spaceworld, dtset%ecut*(dtset%dilatmx)**2, & 463 & formeig, dtset%istwfk, dtset%kptns, mcg, dtset%mband, dtset%mband_mem,dtset%mkmem,dtset%mpw,& 464 & dtset%natom, dtset%nkpt, npwarr, dtset%nspinor, dtset%nsppol, dtset%usepaw,& 465 & cg, eigen=eigen0, pawrhoij=hdr%pawrhoij) 466 ABI_FREE(distrb_flags) 467 468 if (psps%usepaw==1 .and. ireadwf0==1) then 469 ! if parallelism, pawrhoij is distributed, hdr%pawrhoij is not 470 call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,& 471 & mpi_atmtab=mpi_enreg%my_atmtab) 472 end if 473 474 call timab(135,2,tsec) 475 call timab(136,1,tsec) 476 477 ! Report on eigen0 values ! Should use prteigrs.F90 478 !write(message, '(a,a)' ) 479 !call wrtout(std_out,ch10//' respfn : eigen0 array') 480 !nkpt_eff=dtset%nkpt 481 !nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1 482 !if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. dtset%nkpt>nkpt_max ) nkpt_eff=nkpt_max 483 !band_index=0 484 !do isppol=1,dtset%nsppol 485 ! do ikpt=1,dtset%nkpt 486 ! nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 487 ! if(ikpt<=nkpt_eff)then 488 ! write(message, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt 489 ! call wrtout(std_out,message) 490 ! do iband=1,nband_k,4 491 ! write(message, '(a,4es16.6)')' ',eigen0(iband+band_index:min(iband+3,nband_k)+band_index) 492 ! call wrtout(std_out,message) 493 ! end do 494 ! else if(ikpt==nkpt_eff+1)then 495 ! write(message,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigen0.',ch10 496 ! call wrtout(std_out,message) 497 ! end if 498 ! band_index=band_index+nband_k 499 ! end do 500 !end do 501 502 !Allocation for forces and atomic positions (should be taken away, also argument ... ) 503 ABI_MALLOC(grxc,(3,natom)) 504 505 !Do symmetry stuff 506 ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 507 ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 508 ABI_MALLOC(indsym,(4,dtset%nsym,natom)) 509 ABI_MALLOC(symrec,(3,3,dtset%nsym)) 510 irrzon=0;indsym=0;symrec=0;phnons=zero 511 !If the density is to be computed by mkrho, need irrzon and phnons 512 iscf_eff=0;if(dtset%getden==0)iscf_eff=1 513 call setsym(indsym,irrzon,iscf_eff,natom,& 514 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 515 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred) 516 517 !Symmetrize atomic coordinates over space group elements: 518 call symmetrize_xred(natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym) 519 520 !Examine the symmetries of the q wavevector 521 ABI_MALLOC(symq,(4,2,dtset%nsym)) 522 timrev=1 523 524 ! By default use symmetries. 525 use_sym = 1 526 if (dtset%prtgkk == 1)then 527 use_sym = 0 528 call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol,use_sym=use_sym) 529 else 530 call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol) 531 end if 532 533 !Generate an index table of atoms, in order for them to be used 534 !type after type. 535 ABI_MALLOC(atindx,(natom)) 536 ABI_MALLOC(atindx1,(natom)) 537 ABI_MALLOC(nattyp,(ntypat)) 538 indx=1 539 do itypat=1,ntypat 540 nattyp(itypat)=0 541 do iatom=1,natom 542 if(dtset%typat(iatom)==itypat)then 543 atindx(iatom)=indx 544 atindx1(indx)=iatom 545 indx=indx+1 546 nattyp(itypat)=nattyp(itypat)+1 547 end if 548 end do 549 end do 550 551 !Here allocation of GPU for fft calculations 552 #if defined HAVE_GPU_CUDA 553 if (dtset%gpu_option==ABI_GPU_LEGACY .or. dtset%gpu_option==ABI_GPU_KOKKOS) then 554 call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,0,psps,dtset%gpu_option) 555 end if 556 #endif 557 558 !Compute structure factor phases for current atomic pos: 559 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*natom)) 560 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*natom)) 561 call getph(atindx,natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 562 563 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 564 call getph(atindx,natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 565 else 566 ph1df(:,:)=ph1d(:,:) 567 end if 568 569 !Compute occupation numbers and fermi energy, in case occupation scheme is metallic. 570 ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 571 if( dtset%occopt>=3.and.dtset%occopt<=9 ) then 572 573 call newocc(doccde,eigen0,entropy,fermie,fermih,dtset%ivalence,& 574 & dtset%spinmagntarget,dtset%mband,dtset%nband,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,& 575 & dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,dtset%occopt,dtset%prtvol,& 576 & dtset%tphysel,dtset%tsmear,dtset%wtk) 577 578 ! Update fermie and occ 579 etot=hdr%etot ; residm=hdr%residm 580 call hdr%update(bantot,etot,fermie,fermih,residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1)) 581 582 else 583 ! doccde is irrelevant in this case 584 doccde(:)=zero 585 end if 586 587 !Recompute first large sphere cut-off gsqcut, without taking into account dilatmx 588 ecutf=dtset%ecut 589 if (psps%usepaw==1) then 590 ecutf=dtset%pawecutdg 591 call wrtout(std_out,ch10//' FFT (fine) grid used in SCF cycle:') 592 end if 593 594 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf) 595 596 !PAW: 1- Initialize values for several arrays depending only on atomic data 597 !2- Check overlap 598 !3- Identify FFT points in spheres and compute g_l(r).Y_lm(r) (and exp(-i.q.r) if needed) 599 !4- Allocate PAW specific arrays 600 !5- Compute perturbed local potential inside spheres 601 !6- Eventually open temporary storage files 602 if(psps%usepaw==1) then 603 ! 1-Initialize values for several arrays depending only on atomic data 604 gnt_option=1 605 if (dtset%pawxcdev==2.or.dtset%rfphon/=0.or.dtset%rfstrs/=0.or.dtset%rfelfd==1.or.& 606 dtset%rfelfd==3.or.dtset%rf2_dkde==1) gnt_option=2 607 608 ! Test if we have to call pawinit 609 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 610 611 if (psp_gencond==1.or.call_pawinit) then 612 ! Some gen-cond have to be added... 613 call timab(553,1,tsec) 614 call pawinit(dtset%effmass_free,gnt_option,zero,zero,dtset%pawlcutd,dtset%pawlmix,& 615 & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,& 616 & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero) 617 call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,& 618 & rprimd,symrec,pawang%zarot) 619 620 ! Update internal values 621 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 622 623 call timab(553,2,tsec) 624 else 625 if (pawtab(1)%has_kij ==1) pawtab(1:psps%ntypat)%has_kij =2 626 if (pawtab(1)%has_nabla==1) pawtab(1:psps%ntypat)%has_nabla=2 627 end if 628 psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore) 629 call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot) 630 call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,& 631 & is_dfpt,dtset%jpawu,dtset%lexexch,dtset%lpawu,dtset%nspinor,ntypat,dtset%optdcmagpawu,pawang,dtset%pawprtvol,pawrad,& 632 & pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu) 633 compch_fft=-1.d5;compch_sph=-1.d5 634 usexcnhat=maxval(pawtab(:)%usexcnhat) 635 usecprj=dtset%pawusecp 636 ! 2-Check overlap 637 call chkpawovlp(natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred) 638 ! 3-Identify FFT points in spheres and compute g_l(r).Y_lm(r) and exp(-i.q.r) 639 ABI_MALLOC(pawfgrtab,(my_natom)) 640 if (my_natom>0) then 641 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,& 642 & mpi_atmtab=mpi_enreg%my_atmtab) 643 call pawfgrtab_init(pawfgrtab,1,l_size_atm,pawrhoij(1)%nspden,dtset%typat,& 644 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 645 ABI_FREE(l_size_atm) 646 end if 647 use_nhat_gga=(dtset%xclevel==2.and.dtset%pawnhatxc>0.and.usexcnhat>0) 648 optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm 649 if (use_nhat_gga) then 650 optgr1=dtset%pawstgylm 651 if (rfphon==1) optgr2=1 652 end if 653 if (rfphon==1.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1) then ! LB2016-11-28 : Why not rfelfd==1? 654 if (optgr1==0) optgr1=dtset%pawstgylm 655 if (optgr2==0) optgr2=dtset%pawstgylm 656 if (optrad==0.and.(.not.qeq0.or.rfstrs/=0.or.rfelfd==3.or.rf2_dkde==1)) optrad=1 ! LB2016-11-28 : Why not rfelfd==1? 657 end if 658 if (rfelfd==1.or.rfelfd==3.or.rf2_dkde==1) then 659 if (optgr1==0) optgr1=dtset%pawstgylm 660 end if 661 call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfftf,psps%ntypat,& 662 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,& 663 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab ) 664 ! Compute exp(iq.r) factors around the atoms 665 if (.not.qeq0) then 666 do iatom=1,my_natom 667 iatom_tot=iatom; if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom) 668 if (allocated(pawfgrtab(iatom)%expiqr)) then 669 ABI_FREE(pawfgrtab(iatom)%expiqr) 670 end if 671 ABI_MALLOC(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd)) 672 call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,dtset%qptn,& 673 & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot)) 674 pawfgrtab(iatom)%expiqr_allocated=1 675 end do 676 end if 677 ! 4-Allocate PAW specific arrays 678 ABI_MALLOC(paw_an,(my_natom)) 679 ABI_MALLOC(paw_ij,(my_natom)) 680 call paw_an_nullify(paw_an) 681 call paw_ij_nullify(paw_ij) 682 has_kxc=0;nkxc1=0;cplex=1 683 has_dijnd=0; if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1 684 has_diju=merge(0,1,dtset%usepawu==0) 685 if (rfphon/=0.or.rfelfd==1.or.rfelfd==3.or.rfstrs/=0.or.rf2_dkde/=0) then 686 has_kxc=1 687 call pawxc_get_nkxc(nkxc1,dtset%nspden,dtset%xclevel) 688 end if 689 has_vhartree=0 690 if(dtset%orbmag>0 .AND. dtset%pawspnorb > 0) has_vhartree=1 691 call paw_an_init(paw_an,dtset%natom,dtset%ntypat,nkxc1,0,dtset%nspden,& 692 & cplex,dtset%pawxcdev,dtset%typat,pawang,pawtab,has_vxc=1,has_vhartree=has_vhartree,has_vxctau=dtset%usekden,& 693 & has_vxc_ex=1,has_kxc=has_kxc,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 694 call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%pawspnorb,& 695 & natom,dtset%ntypat,dtset%typat,pawtab,has_dij=1,has_dijhartree=1,has_dijnd=has_dijnd,& 696 & has_dijso=1,has_dijU=has_diju,has_pawu_occ=1,has_exexch_pot=1,& 697 & nucdipmom=dtset%nucdipmom,& 698 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 699 700 else ! PAW vs NCPP 701 usexcnhat=0;usecprj=0 702 use_nhat_gga=.false. 703 ABI_MALLOC(paw_an,(0)) 704 ABI_MALLOC(paw_ij,(0)) 705 ABI_MALLOC(pawfgrtab,(0)) 706 end if ! paw 707 708 ABI_MALLOC(rhog,(2,nfftf)) 709 ABI_MALLOC(rhor,(nfftf,dtset%nspden)) 710 ABI_MALLOC(taug,(2,nfftf*dtset%usekden)) 711 ABI_MALLOC(taur,(nfftf,dtset%nspden*dtset%usekden)) 712 713 !>>> Initialize charge density 714 715 if (dtset%getden/=0.or.dtset%irdden/=0) then 716 ! Choice 1: read charge density from a disk file and broadcast data 717 ! This part is not compatible with MPI-FFT (note single_proc=.True. below) 718 rdwrpaw=psps%usepaw ; if(ireadwf0/=0) rdwrpaw=0 719 if (rdwrpaw/=0) then 720 ABI_MALLOC(pawrhoij_read,(natom)) 721 call pawrhoij_nullify(pawrhoij_read) 722 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,& 723 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc) 724 call pawrhoij_alloc(pawrhoij_read,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 725 & dtset%nsppol,dtset%typat,pawtab=pawtab) 726 else 727 ABI_MALLOC(pawrhoij_read,(0)) 728 end if 729 ! Note MT july 2013: should we read rhoij from the density file? 730 call read_rhor(dtfil%fildensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw,& 731 & mpi_enreg,rhor,hdr_den,pawrhoij_read,spaceworld,check_hdr=hdr) 732 etotal = hdr_den%etot 733 call hdr_den%free() 734 if (rdwrpaw/=0) then 735 call pawrhoij_bcast(pawrhoij_read,hdr%pawrhoij,0,spaceworld) 736 call pawrhoij_free(pawrhoij_read) 737 end if 738 ABI_FREE(pawrhoij_read) 739 ! Compute up+down rho(G) by fft 740 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 741 742 else 743 ! Choice 2: obtain the charge density from read wfs 744 ! Warning: in PAW, compensation density has to be added ! 745 tim_mkrho=4 746 paw_dmft%use_dmft=0 ; paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented 747 if (psps%usepaw==1) then 748 ABI_MALLOC(rhowfg,(2,dtset%nfft)) 749 ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden)) 750 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 751 & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 752 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor) 753 ABI_FREE(rhowfg) 754 ABI_FREE(rhowfr) 755 else 756 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 757 & rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 758 end if 759 760 end if ! choice for charge density initialization 761 762 !>>> Initialize kinetic energy density 763 if (dtset%usekden==1) then 764 765 if (dtset%getkden/=0.or.dtset%irdkden/=0) then 766 ! Choice 1: read kinetic energy density from a disk file and broadcast data 767 ! This part is not compatible with MPI-FFT (note single_proc=.True. below) 768 rdwrpaw=0 769 ABI_MALLOC(pawrhoij_read,(0)) 770 call read_rhor(dtfil%filkdensin,cplex1,dtset%nspden,nfftf,ngfftf,rdwrpaw,& 771 & mpi_enreg,taur,hdr_den,pawrhoij_read,spaceworld,check_hdr=hdr) 772 call hdr_den%free() 773 ABI_FREE(pawrhoij_read) 774 ! Compute up+down tau(G) by fft 775 call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 776 777 else 778 ! Choice 2: obtain the kinetic energy density from read wfs 779 tim_mkrho=4 780 paw_dmft%use_dmft=0 ; paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented 781 if (psps%usepaw==1) then 782 ABI_MALLOC(rhowfg,(2,dtset%nfft)) 783 ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden)) 784 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 785 & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 786 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,taug,rhowfr,taur) 787 ABI_FREE(rhowfg) 788 ABI_FREE(rhowfr) 789 else 790 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,& 791 & taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1) 792 end if 793 794 end if ! choice for kinetic energy density initialization 795 end if ! usekden 796 797 !In PAW, compensation density has eventually to be added 798 nhatgrdim=0;nhatdim=0 799 ABI_MALLOC(nhatgr,(0,0,0)) 800 if (psps%usepaw==1.and. ((usexcnhat==0).or.(dtset%getden==0).or.dtset%xclevel==2)) then 801 nhatdim=1 802 803 ABI_MALLOC(nhat,(nfftf,dtset%nspden)) 804 call timab(558,1,tsec) 805 nhatgrdim=0;if (dtset%xclevel==2.and.dtset%pawnhatxc>0) nhatgrdim=usexcnhat 806 ider=2*nhatgrdim 807 if (nhatgrdim>0) then 808 ABI_FREE(nhatgr) 809 ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3)) 810 end if 811 izero=0;cplex=1;ipert=0;idir=0;qphon(:)=zero 812 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,& 813 & nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,& 814 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, & 815 & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom) 816 if (dtset%getden==0) then 817 rhor(:,:)=rhor(:,:)+nhat(:,:) 818 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 819 end if 820 call timab(558,2,tsec) 821 else 822 ABI_MALLOC(nhat,(0,0)) 823 end if 824 825 !The GS irrzon and phnons were only needed to symmetrize the GS density 826 ABI_FREE(irrzon) 827 ABI_FREE(phnons) 828 829 !Will compute now the total potential 830 831 !Compute local ionic pseudopotential vpsp and core electron density xccc3d 832 n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf 833 ABI_MALLOC(xccc3d,(n3xccc)) 834 ABI_MALLOC(vpsp,(nfftf)) 835 if(psps%usepaw==1) then 836 ABI_MALLOC(xcctau3d,(n3xccc*dtset%usekden)) 837 else 838 ABI_MALLOC(xcctau3d,(0)) 839 end if 840 841 842 !Determine by which method the local ionic potential and/or 843 ! the pseudo core charge density have to be computed 844 !Local ionic potential: 845 ! Method 1: PAW ; Method 2: Norm-conserving PP 846 vloc_method=1;if (psps%usepaw==0) vloc_method=2 847 !Pseudo core charge density: 848 ! Method 1: PAW, nc_xccc_gspace ; Method 2: Norm-conserving PP 849 coredens_method=1;if (psps%usepaw==0) coredens_method=2 850 if (psps%nc_xccc_gspace==1) coredens_method=1 851 if (psps%nc_xccc_gspace==0) coredens_method=2 852 853 ! Core kinetic energy density method 854 coretau_method = 0 855 if (dtset%usekden==1.and.psps%usepaw==1) then 856 coretau_method=1 857 if (psps%nc_xccc_gspace==0) then 858 coretau_method=2 859 end if 860 end if 861 862 863 !Local ionic potential and/or pseudo core charge by method 1 864 if (vloc_method==1.or.coredens_method==1) then 865 call timab(562,1,tsec) 866 optv=0;if (vloc_method==1) optv=1 867 optn=0;if (coredens_method==1) optn=n3xccc/nfftf 868 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1 869 call atm2fft(atindx1,xccc3d,vpsp,dum_dyfrn,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,& 870 & dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,& 871 & ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,& 872 & dtset%qprtrb,dtset%rcut,dum_rhog,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl) 873 call timab(562,2,tsec) 874 end if 875 876 if (coretau_method==1) then 877 optv=0;optn=1 878 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=4 879 call atm2fft(atindx1,xcctau3d,dummy_out6,dummy_out1,dummy_out2,dummy_out3,dummy_in,& 880 & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfftf,psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,ntypat,& 881 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,dtset%qprtrb,& 882 & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,& 883 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 884 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 885 end if 886 887 !Local ionic potential by method 2 888 if (vloc_method==2) then 889 option=1 890 ABI_MALLOC(dyfrlo_indx,(3,3,natom)) 891 ABI_MALLOC(grtn_indx,(3,natom)) 892 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,& 893 & grtn_indx,gsqcut,dummy6,mgfftf,mpi_enreg,natom,nattyp,& 894 & nfftf,ngfftf,dtset%nspden,ntypat,option,pawtab,ph1df,psps,& 895 & dtset%qprtrb,rhog,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred) 896 ABI_FREE(dyfrlo_indx) 897 ABI_FREE(grtn_indx) 898 end if 899 900 !Pseudo core electron density by method 2 901 if (coredens_method==2.and.psps%n1xccc/=0) then 902 option=1 903 ABI_MALLOC(dyfrx2,(3,3,natom)) 904 ABI_MALLOC(vxc,(0,0)) ! dummy 905 call mkcore(dummy6,dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,& 906 & ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,& 907 & psps%xcccrc,psps%xccc1d,xccc3d,xred) 908 ABI_FREE(dyfrx2) 909 ABI_FREE(vxc) ! dummy 910 end if 911 912 !Set up hartree and xc potential. Compute kxc here. 913 ABI_MALLOC(vhartr,(nfftf)) 914 915 call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfftf,ngfftf,& 916 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 917 918 option=2 ; nk3xc=1 919 nkxc=2*min(dtset%nspden,2)-1;if(dtset%xclevel==2)nkxc=12*min(dtset%nspden,2)-5 920 call check_kxc(dtset%ixc,dtset%optdriver) 921 ABI_MALLOC(kxc,(nfftf,nkxc)) 922 ABI_MALLOC(vxc,(nfftf,dtset%nspden)) 923 ABI_MALLOC(vxctau,(nfftf,dtset%nspden,4*dtset%usekden)) 924 925 call xcdata_init(xcdata,dtset=dtset) 926 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 927 !If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value 928 if (dtset%xc_tb09_c>99._dp) then 929 call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom,nfftf,ngfftf, & 930 & nhat,psps%usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, & 931 & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,psps%usepaw, & 932 & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 933 end if 934 935 call rhotoxc(enxc,kxc,mpi_enreg,nfftf,ngfftf,& 936 & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor,& 937 & rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 938 & taur=taur,vhartr=vhartr,vxctau=vxctau,xcctau3d=xcctau3d) 939 940 !Compute local + Hxc potential, and subtract mean potential. 941 ABI_MALLOC(vtrial,(nfftf,dtset%nspden)) 942 do ispden=1,min(dtset%nspden,2) 943 do ifft=1,nfftf 944 vtrial(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft) 945 end do 946 end do 947 if (dtset%nspden==4) then 948 do ispden=3,4 949 do ifft=1,nfftf 950 vtrial(ifft,ispden)=vxc(ifft,ispden) 951 end do 952 end do 953 end if 954 955 ABI_FREE(vhartr) 956 957 if(dtset%prtvol==-level) call wrtout(std_out,' respfn: ground-state density and potential set up.') 958 959 !PAW: compute Dij quantities (psp strengths) 960 if (psps%usepaw==1)then 961 cplex=1;ipert=0;option=1 962 nzlmopt=0;if (dtset%pawnzlm>0) nzlmopt=-1 963 964 call pawdenpot(compch_sph,epaw,epawdc,ipert,dtset%ixc,my_natom,natom,dtset%nspden,& 965 & ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,& 966 & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,& 967 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp, & 968 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 969 970 call timab(561,1,tsec) 971 call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,natom,nfftf,nfftotf,& 972 & dtset%nspden,ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,& 973 & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,& 974 & dtset%spnorbscl,ucvol,dtset%cellcharge(1),vtrial,vxc,xred,nucdipmom=dtset%nucdipmom,& 975 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 976 call symdij(gprimd,indsym,ipert,my_natom,natom,dtset%nsym,ntypat,0,& 977 & paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,& 978 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 979 call timab(561,2,tsec) 980 981 end if 982 983 !-----2. Frozen-wavefunctions and Ewald(q=0) parts of 2DTE 984 985 dyfr_nondiag=0;if (psps%usepaw==1.and.rfphon==1) dyfr_nondiag=1 986 dyfr_cplex=1;if (psps%usepaw==1.and.rfphon==1.and.(.not.qeq0)) dyfr_cplex=2 987 ABI_MALLOC(dyew,(2,3,natom,3,natom)) 988 ABI_MALLOC(dyewq0,(3,3,natom)) 989 ABI_MALLOC(dyfrlo,(3,3,natom)) 990 ABI_MALLOC(dyfrx2,(3,3,natom)) 991 ABI_MALLOC(dyfrnl,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)) 992 ABI_MALLOC(dyfrwf,(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)) 993 ABI_MALLOC(becfrnl,(3,natom,3*pawbec)) 994 ABI_MALLOC(piezofrnl,(6,3*pawpiezo)) 995 ABI_MALLOC(dyvdw,(2,3,natom,3,natom*usevdw)) 996 dyew(:,:,:,:,:)=zero 997 dyewq0(:,:,:)=zero 998 dyfrnl(:,:,:,:,:)=zero 999 dyfrwf(:,:,:,:,:)=zero 1000 dyfrlo(:,:,:)=zero 1001 dyfrx2(:,:,:)=zero 1002 if (usevdw==1) dyvdw(:,:,:,:,:)=zero 1003 if (pawbec==1) becfrnl(:,:,:)=zero 1004 if (pawpiezo==1) piezofrnl(:,:)=zero 1005 1006 ABI_MALLOC(eltcore,(6,6)) 1007 ABI_MALLOC(elteew,(6+3*natom,6)) 1008 ABI_MALLOC(eltfrhar,(6,6)) 1009 ABI_MALLOC(eltfrnl,(6+3*natom,6)) 1010 ABI_MALLOC(eltfrloc,(6+3*natom,6)) 1011 ABI_MALLOC(eltfrkin,(6,6)) 1012 ABI_MALLOC(eltfrxc,(6+3*natom,6)) 1013 ABI_MALLOC(eltvdw,(6+3*natom,6*usevdw)) 1014 eltcore(:,:)=zero 1015 elteew(:,:)=zero 1016 eltfrnl(:,:)=zero 1017 eltfrloc(:,:)=zero 1018 eltfrkin(:,:)=zero 1019 eltfrhar(:,:)=zero 1020 eltfrxc(:,:)=zero 1021 if (usevdw==1) eltvdw(:,:)=zero 1022 1023 !Section common to all perturbations 1024 !Compute the nonlocal part of the elastic tensor and/or dynamical matrix 1025 if (rfstrs/=0.or.rfphon==1.or.dtset%efmas>0.or.pawbec==1.or.pawpiezo==1)then 1026 call d2frnl(becfrnl,cg,dtfil,dtset,dyfrnl,dyfr_cplex,& 1027 & dyfr_nondiag,efmasdeg,efmasval,eigen0,eltfrnl,gsqcut,has_allddk,indsym,kg,& 1028 & dtset%mband_mem,dtset%mkmem,mgfftf,& 1029 & mpi_enreg,psps%mpsang,my_natom,natom,nfftf,ngfft,ngfftf,& 1030 & npwarr,occ,paw_ij,pawang,pawbec,pawfgrtab,pawpiezo,pawrad,& 1031 & pawrhoij,pawtab,ph1d,ph1df,piezofrnl,psps,rprimd,rfphon,& 1032 & rfstrs,symrec,vtrial,vxc,xred,ylm,ylmgr) 1033 end if 1034 1035 !No more need of these local derivatives 1036 if (rfphon==1.and.psps%usepaw==1.and.(.not.use_nhat_gga)) then 1037 do iatom=1,my_natom 1038 if (allocated(pawfgrtab(iatom)%gylmgr2)) then 1039 ABI_FREE(pawfgrtab(iatom)%gylmgr2) 1040 end if 1041 pawfgrtab(iatom)%gylmgr2_allocated=0 1042 end do 1043 end if 1044 1045 !Section for the atomic displacement/electric field perturbations 1046 if (rfphon==1) then 1047 1048 ! Compute the local of the dynamical matrix 1049 ! dyfrnl has not yet been symmetrized, but will be in the next routine 1050 call dfpt_dyfro(atindx1,dyfrnl,dyfrlo,dyfrwf,dyfrx2,dyfr_cplex,dyfr_nondiag,& 1051 & dtset,gmet,gprimd,gsqcut,indsym,mgfftf,mpi_enreg,psps%mqgrid_vl,& 1052 & natom,nattyp, nfftf,ngfftf,dtset%nspden,dtset%nsym,ntypat,& 1053 & psps%n1xccc,n3xccc,psps,pawtab,ph1df,psps%qgrid_vl,& 1054 & dtset%qptn,rhog,rprimd,symq,symrec,dtset%typat,ucvol,& 1055 & psps%usepaw,psps%vlspl,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred) 1056 1057 ! Compute Ewald (q=0) contribution 1058 sumg0=0;qphon(:)=zero 1059 call dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,psps%ziontypat,& 1060 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1061 option=1 1062 call q0dy3_calc(natom,dyewq0,dyew,option) 1063 ! Calculate the DFT-D2 vdw part of the dynamical matrix 1064 if (usevdw==1.and.dtset%vdw_xc==5) then 1065 call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,dtset%vdw_tol,& 1066 & xred,psps%znuclpsp,dyn_vdw_dftd2=dyvdw,qphon=dtset%qptn) 1067 end if 1068 ! Calculate the DFT-D3/D3(BJ) vdw part of the dynamical matrix 1069 if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then 1070 call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,1,dtset%typat,rprimd,& 1071 & dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,& 1072 & dyn_vdw_dftd3=dyvdw,qphon=dtset%qptn) 1073 end if 1074 !The frozen-wavefunction part of the dynamical matrix is now: 1075 ! d2frnl: non-local contribution 1076 ! dyfrlo: local contribution 1077 ! dyfrx2: 2nd order xc core correction contribution 1078 ! dyew : Ewald contribution 1079 ! dyvdw : vdw DFT-D contribution 1080 ! dyfrwf: all contributions 1081 ! In case of PAW, it misses a term coming from the perturbed overlap operator 1082 end if 1083 1084 !Section for the strain perturbation 1085 if(rfstrs/=0) then 1086 1087 ! Verify that k-point set has full space-group symmetry; otherwise exit 1088 timrev=1 1089 if (symkchk(dtset%kptns,dtset%nkpt,dtset%nsym,symrec,timrev,message) /= 0) then 1090 ABI_ERROR(message) 1091 end if 1092 1093 ! Calculate the kinetic part of the elastic tensor 1094 call dfpt_eltfrkin(cg,eltfrkin,dtset%ecut,dtset%ecutsm,dtset%effmass_free,& 1095 & dtset%istwfk,kg,dtset%kptns,dtset%mband,dtset%mband_mem,dtset%mgfft,dtset%mkmem,mpi_enreg,& 1096 & dtset%mpw,dtset%nband,dtset%nkpt,ngfft,npwarr,& 1097 & dtset%nspinor,dtset%nsppol,occ,rprimd,dtset%wtk) 1098 1099 ! Calculate the hartree part of the elastic tensor 1100 call dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfftf,ngfftf,rhog) 1101 1102 ! Calculate the xc part of the elastic tensor 1103 call dfpt_eltfrxc(atindx,dtset,eltfrxc,enxc,gsqcut,kxc,mpi_enreg,mgfftf,& 1104 & nattyp,nfftf,ngfftf,ngfftf,nhat,nkxc,n3xccc,pawtab,ph1df,psps,rhor,rprimd,& 1105 & usexcnhat,vxc,xccc3d,xred) 1106 1107 ! Calculate the local potential part of the elastic tensor 1108 call dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfftf,mpi_enreg,psps%mqgrid_vl,& 1109 & natom,nattyp,nfftf,ngfftf,ntypat,ph1df,psps%qgrid_vl,rhog,psps%vlspl) 1110 1111 ! Calculate the Ewald part of the elastic tensor 1112 call elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,& 1113 & dtset%typat,ucvol,xred,psps%ziontypat,& 1114 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1115 1116 ! Calculate the DFT-D2 vdw part of the elastic tensor 1117 if (usevdw==1.and.dtset%vdw_xc==5) then 1118 option=1; if (rfphon==1) option=0 1119 call vdw_dftd2(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,dtset%vdw_tol,& 1120 & xred,psps%znuclpsp,elt_vdw_dftd2=eltvdw) 1121 end if 1122 ! Calculate the DFT-D3/D3(BJ) vdw part of the elastic tensor 1123 if (usevdw==1.and.(dtset%vdw_xc==6.or.dtset%vdw_xc==7)) then 1124 option=1; if (rfphon==1) option=0 1125 call vdw_dftd3(evdw,dtset%ixc,natom,ntypat,option,dtset%typat,rprimd,& 1126 & dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,xred,psps%znuclpsp,& 1127 & elt_vdw_dftd3=eltvdw) 1128 end if 1129 ! Calculate the psp core energy part of elastic tensor (trivial) 1130 eltcore(1:3,1:3)=ecore/ucvol 1131 1132 !The frozen-wavefunction part of the elastic tensor is now: 1133 ! eltfrnl: non-local contribution 1134 ! eltfrloc: local contribution 1135 ! eltfrkin: kinetic contribution 1136 ! eltfrhar: Hartree contribution 1137 ! eltfrx: XC contribution 1138 ! eltcore: psps core contribution 1139 ! elteew: Ewald contribution 1140 ! eltvdw: vdw DFT-D contribution 1141 ! In case of PAW, it misses a term coming from the perturbed overlap operator 1142 end if 1143 1144 ABI_FREE(vpsp) 1145 ABI_FREE(xccc3d) 1146 if(allocated(xcctau3d)) then 1147 ABI_FREE(xcctau3d) 1148 end if 1149 1150 if(dtset%prtvol==-level) call wrtout(std_out,' respfn: frozen wavef. and Ewald(q=0) part of 2DTE done.') 1151 1152 call timab(136,2,tsec) 1153 1154 !-----3. Initialisation of 1st response, taking into account the q vector. 1155 1156 call timab(137,1,tsec) 1157 1158 write(message,'(3a)')ch10,' ==> initialize data related to q vector <== ',ch10 1159 call wrtout([std_out, ab_out] ,message) 1160 1161 qphon(:)=dtset%qptn(:) 1162 sumg0=1 1163 1164 !Treat the case of q=0 or q too close to 0 1165 qzero=0 1166 if(qeq0)then 1167 qphon(:)=zero 1168 write(message,'(3a)')& 1169 & ' respfn : the norm of the phonon wavelength (as input) was small (<1.d-7).',ch10,& 1170 & ' q has been set exactly to (0 0 0)' 1171 call wrtout(std_out,message) 1172 sumg0=0 1173 qzero=1 1174 else 1175 if(rfelfd/=0 .or. rfstrs/=0 .or. rfddk /= 0 .or. rf2_dkdk /= 0 .or. rf2_dkde /= 0) then 1176 ! Temporarily, ... 1177 write(message, '(a,a,a,3es16.6,a,a,5(a,i2),a,a,a)' )ch10,& 1178 & 'The treatment of non-zero wavevector q is restricted to phonons.',& 1179 & 'However, the input normalized qpt is',qphon(:),',',ch10,& 1180 & 'while rfelfd=',rfelfd,', rfddk=',rfddk,', rf2_dkdk=',rf2_dkdk,', rf2_dkde=',rf2_dkde,& 1181 & ' and rfstrs=',rfstrs,'.',ch10,& 1182 & 'Action: change qpt, or rfelfd, or rfstrs in the input file.' 1183 ABI_ERROR(message) 1184 end if 1185 if(chneut/=0)then 1186 write(message,'(2a)')ch10,' chneut/=0 not allowed with q/=0 => chneut reset to 0 locally.' 1187 ABI_WARNING(message) 1188 chneut=0 1189 end if 1190 if(asr/=0)then 1191 write(message,'(2a)')ch10,' asr/=0 not allowed with q/=0 => asr reset to 0 locally.' 1192 ABI_WARNING(message) 1193 asr=0 1194 end if 1195 1196 end if 1197 1198 !Determine the symmetrical perturbations 1199 ABI_MALLOC(pertsy,(3,natom+6)) 1200 call irreducible_set_pert(indsym,natom+6,natom,dtset%nsym,pertsy,rfdir,rfpert,symq,symrec,dtset%symrel) 1201 1202 write(message,'(a)') ' The list of irreducible perturbations for this q vector is:' 1203 call wrtout([std_out, ab_out] ,message) 1204 ii=1 1205 do ipert=1,natom+6 ! GA: Why natom+6 instead of natom+MPERT_MAX ? 1206 do idir=1,3 1207 if(rfpert(ipert)==1.and.rfdir(idir)==1)then 1208 if( pertsy(idir,ipert)==1 )then 1209 write(message, '(i5,a,i2,a,i4)' )ii,') idir=',idir,' ipert=',ipert 1210 call wrtout([std_out, ab_out] ,message) 1211 ii=ii+1 1212 end if 1213 end if 1214 end do 1215 end do 1216 1217 !test if the user left default rfdir 0 0 0 1218 if (ii==1 .and. rf2_dkdk==0 .and. rf2_dkde==0) then 1219 write(message,'(5a)')ch10,& 1220 & ' WARNING: no perturbations to be done at this q-point.',ch10,& 1221 & ' You may have forgotten to set the rfdir or rfatpol variables. Continuing normally.',ch10 1222 call wrtout(ab_out,message) 1223 ABI_WARNING(message) 1224 end if 1225 1226 if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then 1227 ABI_MALLOC(rfpert_nl,(3,natom+2,3,natom+2,3,natom+2)) 1228 rfpert_nl = 0 1229 rfpert_nl(:,natom+2,:,natom+2,:,natom+2) = 1 1230 rfpert_nl(:,1:natom,:,natom+2,:,natom+2) = 1 1231 rfpert_nl(:,natom+2,:,1:natom,:,natom+2) = 1 1232 rfpert_nl(:,natom+2,:,natom+2,:,1:natom) = 1 1233 call sytens(indsym,natom+2,natom,dtset%nsym,rfpert_nl,symrec,dtset%symrel) 1234 write(message, '(a,a,a,a,a)' ) ch10, & 1235 & ' The list of irreducible elements of the Raman and non-linear',& 1236 & ch10,' optical susceptibility tensors is:',ch10 1237 call wrtout(std_out,message) 1238 1239 write(message,'(12x,a)')'i1pert i1dir i2pert i2dir i3pert i3dir' 1240 call wrtout(std_out,message) 1241 n1 = 0 1242 rf2_dirs_from_rfpert_nl(:,:) = 0 1243 do i1pert = 1, natom + 2 1244 do i1dir = 1, 3 1245 do i2pert = 1, natom + 2 1246 do i2dir = 1, 3 1247 do i3pert = 1, natom + 2 1248 do i3dir = 1,3 1249 if (rfpert_nl(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 1250 n1 = n1 + 1 1251 write(message,'(2x,i4,a,6(5x,i3))') n1,')', i1pert,i1dir,i2pert,i2dir,i3pert,i3dir 1252 call wrtout(std_out,message) 1253 if (i2pert==natom+2) then 1254 if (i3pert==natom+2) then 1255 rf2_dirs_from_rfpert_nl(i3dir,i2dir) = 1 1256 else if (i1pert==natom+2) then 1257 rf2_dirs_from_rfpert_nl(i1dir,i2dir) = 1 1258 end if 1259 end if 1260 end if 1261 end do 1262 end do 1263 end do 1264 end do 1265 end do 1266 end do 1267 write(message,'(a,a)') ch10,ch10 1268 call wrtout(std_out,message) 1269 1270 call wrtout(std_out,'rf2_dirs_from_rfpert_nl :') 1271 do i1dir = 1, 3 1272 do i2dir = 1, 3 1273 write(message,'(3(a,i1))') ' ',i1dir,' ',i2dir,' : ',rf2_dirs_from_rfpert_nl(i1dir,i2dir) 1274 call wrtout(std_out,message) 1275 end do 1276 end do 1277 end if 1278 1279 !For longwave calculation: 1280 !Get symmetries in cartesian coordinates 1281 ABI_MALLOC(symrel_cart, (3, 3, dtset%nsym)) 1282 do isym =1,dtset%nsym 1283 call symredcart(rprimd, gprimd, symrel_cart(:,:,isym), dtset%symrel(:,:,isym)) 1284 ! purify operations in cartesian coordinates. 1285 where (abs(symrel_cart(:,:,isym)) < tol14) 1286 symrel_cart(:,:,isym) = zero 1287 end where 1288 end do 1289 1290 if (dtset%prepalw/=0) then 1291 ABI_MALLOC(rfpert_lw,(3,natom+8,3,natom+8,3,natom+8)) 1292 rfpert_lw=0 1293 if (dtset%prepalw==1) then 1294 rfpert_lw(:,1:natom+2,:,1:natom,:,natom+8)=1 1295 rfpert_lw(:,1:natom+2,:,natom+3:natom+4,:,natom+8)=1 1296 else if (dtset%prepalw==2) then 1297 rfpert_lw(:,natom+2,:,1:natom,:,natom+8)=1 1298 else if (dtset%prepalw==3) then 1299 rfpert_lw(:,1:natom+2,:,1:natom,:,natom+8)=1 1300 else if (dtset%prepalw==4) then 1301 rfpert_lw(:,natom+2,:,natom+2,:,natom+8)=1 1302 end if 1303 1304 ! call sylwtens(indsym,natom+8,natom,dtset%nsym,rfpert_lw,symrec,dtset%symrel,symrel_cart) 1305 call sylwtens(indsym,natom+8,natom,dtset%nsym,rfpert_lw,symrec,dtset%symrel) 1306 1307 write(message,'(7a)') ch10, ' The following reducible perturbations will also be ', ch10, & 1308 & ' explicitly calculated for a correct subsequent ', ch10, & 1309 & ' execution of the longwave driver:', ch10 1310 call wrtout(ab_out,message,'COLL') 1311 call wrtout(std_out,message,'COLL') 1312 do i3pert = 1,natom+8 1313 do i3dir = 1, 3 1314 do i2pert = 1, natom+8 1315 do i2dir = 1,3 1316 do i1pert = 1,natom+8 1317 do i1dir = 1, 3 1318 if (rfpert_lw(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 1319 if (pertsy(i1dir,i1pert)==-1) then 1320 pertsy(i1dir,i1pert)=1 1321 write(message,'(a,i2,a,i4)' )' idir=',i1dir,' ipert=',i1pert 1322 call wrtout(ab_out,message,'COLL') 1323 call wrtout(std_out,message,'COLL') 1324 end if 1325 if (pertsy(i2dir,i2pert)==-1) then 1326 pertsy(i2dir,i2pert)=1 1327 write(message,'(a,i2,a,i4)' )' idir=',i2dir,' ipert=',i2pert 1328 call wrtout(ab_out,message,'COLL') 1329 call wrtout(std_out,message,'COLL') 1330 end if 1331 end if 1332 end do 1333 end do 1334 end do 1335 end do 1336 end do 1337 end do 1338 write(message,'(a,a)') ch10,ch10 1339 call wrtout(std_out,message,'COLL') 1340 ABI_FREE(rfpert_lw) 1341 end if 1342 1343 !Contribution to the dynamical matrix from ion-ion energy 1344 if(rfphon==1)then 1345 call dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,psps%ziontypat, & 1346 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1347 call q0dy3_apply(natom,dyewq0,dyew) 1348 end if 1349 1350 !1-order contribution of the xc core correction to the dynamical matrix 1351 ABI_MALLOC(dyfrx1,(2,3,natom,3,natom)) 1352 dyfrx1(:,:,:,:,:)=zero 1353 if(rfphon==1.and.psps%n1xccc/=0)then 1354 ABI_MALLOC(blkflgfrx1,(3,natom,3,natom)) 1355 !FR non-collinear magnetism 1356 if (dtset%nspden==4) then 1357 call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,& 1358 & psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,non_magnetic_xc,dtset%nspden,& 1359 & ntypat,psps%n1xccc,psps,pawtab,ph1df,psps%qgrid_vl,qphon,& 1360 & rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred,rhor=rhor,vxc=vxc) 1361 else 1362 call dfpt_dyxc1(atindx,blkflgfrx1,dyfrx1,gmet,gsqcut,dtset%ixc,kxc,mgfftf,mpert,mpi_enreg,& 1363 & psps%mqgrid_vl,natom,nfftf,ngfftf,nkxc,non_magnetic_xc,dtset%nspden,& 1364 & ntypat,psps%n1xccc,psps,pawtab,ph1df,psps%qgrid_vl,qphon,& 1365 & rfdir,rfpert,rprimd,timrev,dtset%typat,ucvol,psps%usepaw,psps%xcccrc,psps%xccc1d,xred) 1366 end if 1367 end if 1368 1369 !Deallocate the arrays that were needed only for the frozen wavefunction part 1370 ABI_FREE(ph1d) 1371 ABI_FREE(ph1df) 1372 ABI_FREE(cg) 1373 ABI_FREE(kg) 1374 ABI_FREE(npwarr) 1375 if(xmpi_paral==1) then 1376 ABI_FREE(mpi_enreg%proc_distrb) 1377 end if 1378 1379 ABI_MALLOC(blkflg,(3,mpert,3,mpert)) 1380 ABI_MALLOC(d2eig0,(2,3,mpert,3,mpert)) 1381 ABI_MALLOC(d2k0,(2,3,mpert,3,mpert)) 1382 ABI_MALLOC(d2lo,(2,3,mpert,3,mpert)) 1383 ABI_MALLOC(d2loc0,(2,3,mpert,3,mpert)) 1384 ABI_MALLOC(d2nfr,(2,3,mpert,3,mpert)) 1385 ABI_MALLOC(d2nl,(2,3,mpert,3,mpert)) 1386 ABI_MALLOC(d2nl0,(2,3,mpert,3,mpert)) 1387 ABI_MALLOC(d2nl1,(2,3,mpert,3,mpert)) 1388 ABI_MALLOC(d2vn,(2,3,mpert,3,mpert)) 1389 ABI_MALLOC(d2ovl,(2,3,mpert,3,mpert*psps%usepaw)) 1390 blkflg(:,:,:,:)=0 1391 d2eig0(:,:,:,:,:)=zero ; d2k0(:,:,:,:,:)=zero 1392 d2lo(:,:,:,:,:)=zero ; d2loc0(:,:,:,:,:)=zero 1393 d2nfr(:,:,:,:,:)=zero ; d2nl(:,:,:,:,:)=zero 1394 d2nl0(:,:,:,:,:)=zero ; d2nl1(:,:,:,:,:)=zero 1395 d2vn(:,:,:,:,:)=zero 1396 if (psps%usepaw==1) d2ovl(:,:,:,:,:)=zero 1397 1398 prtbbb=dtset%prtbbb 1399 ABI_MALLOC(d2bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)) 1400 ABI_MALLOC(d2cart_bbb,(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)) 1401 if(prtbbb==1)then 1402 d2cart_bbb(:,:,:,:,:,:)=zero 1403 d2bbb(:,:,:,:,:,:)=zero 1404 end if 1405 1406 dim_eig2nkq = 0 1407 if(dtset%ieig2rf /= 0) dim_eig2nkq = 1 1408 ABI_MALLOC(eig2nkq,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eig2nkq)) 1409 dim_eigbrd=0 1410 if(dtset%ieig2rf /= 0 .and. dtset%smdelta>0 ) dim_eigbrd = 1 1411 ABI_MALLOC(eigbrd,(2,dtset%mband*dtset%nsppol,dtset%nkpt,3,natom,3,natom*dim_eigbrd)) 1412 1413 call timab(137,2,tsec) 1414 1415 1416 !Check whether exiting was required by the user. 1417 !If found then do not start minimization steps 1418 !At this first call to exit_check, initialize cpus 1419 cpus=dtset%cpus 1420 if(abs(cpus)>1.0d-5)cpus=cpus+cpui 1421 openexit=1; if(dtset%chkexit==0) openexit=0 1422 call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit) 1423 1424 !TEMPORARY: for testing purpose only 1425 ! if (rfstrs/=0.and.dtset%usepaw==1) iexit=1 1426 1427 if (iexit==0) then 1428 ! ####################################################################### 1429 write(message,'(a,80a)')ch10,('=',mu=1,80) 1430 call wrtout([std_out, ab_out], message) 1431 1432 ddkfil(:)=0 1433 1434 ! MGNAG: WHY THIS? all these variables are declared as optional pointers in dfpt_looppert! 1435 ! but they are allocated here so why using pointers! Moreover OPTIONAL arguments MUST 1436 ! be passed by keyword for better clarity and robustness! 1437 ! People should learn how to program Fort90 before being allowed to change the code! 1438 ! v5[26] crashes in dfpt_looppert 1439 ! The best solution would be using a datatype to gather the results! 1440 ABI_MALLOC(clflg,(3,mpert)) 1441 if(dtset%ieig2rf > 0) then 1442 ABI_MALLOC(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1443 ABI_MALLOC(eigenq_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1444 ABI_MALLOC(occ_rbz_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1445 end if 1446 if(dtset%efmas > 0) then 1447 ABI_MALLOC(eigen0_pert,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1448 end if 1449 ! Note that kg, cg, eigen0, mpw and npwarr are NOT passed to dfpt_looppert : 1450 ! they will be reinitialized for each perturbation, with an eventually 1451 ! reduced set of k point, thanks to the use of symmetry operations. 1452 call dfpt_looppert(atindx,blkflg,codvsn,cpus,dim_eigbrd,dim_eig2nkq,doccde,& 1453 & ddkfil,dtfil,dtset,dyew,dyfrlo,dyfrnl,dyfrx1,dyfrx2,dyvdw,dyfr_cplex,dyfr_nondiag,& 1454 & d2bbb,d2lo,d2nl,d2ovl,efmasdeg,efmasval,eigbrd,eig2nkq,& 1455 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 1456 & etotal,fermie,iexit,indsym,kxc,& 1457 & dtset%mkmem,mkqmem,mk1mem,mpert,mpi_enreg,my_natom,nattyp,& 1458 & nfftf,nhat,dtset%nkpt,nkxc,dtset%nspden,dtset%nsym,occ,& 1459 & paw_an,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,& 1460 & pertsy,prtbbb,psps,rfpert,rf2_dirs_from_rfpert_nl,rhog,rhor,symq,symrec,timrev,& 1461 & usecprj,usevdw,vtrial,vxc,vxcavg,vxctau,xred,clflg,occ_rbz_pert,eigen0_pert,eigenq_pert,& 1462 & eigen1_pert,nkpt_rbz,eigenq_fine,hdr_fine,hdr0) 1463 1464 ! ##################################################################### 1465 end if 1466 1467 call timab(138,1,tsec) 1468 1469 write(message, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,& 1470 ' ---- first-order wavefunction calculations are completed ----',ch10 1471 call wrtout([std_out, ab_out], message) 1472 1473 ABI_FREE(vxc) 1474 ABI_FREE(vxctau) 1475 1476 if (dtset%prepanl==1.and.(rf2_dkdk/=0 .or. rf2_dkde/=0)) then 1477 ABI_FREE(rfpert_nl) 1478 end if 1479 1480 !Output of the localization tensor 1481 if ( rfpert(natom+1) /= 0 .and. (me == 0) .and. dtset%occopt<=2) then 1482 call wrtloctens(blkflg,d2bbb,d2nl,dtset%mband,mpert,natom,dtset%prtbbb,rprimd,psps%usepaw) 1483 end if 1484 1485 !The perturbation natom+1 was only an auxiliary perturbation, 1486 !needed to construct the electric field response, so its flag is now set to 0. 1487 !rfpert(natom+1)=0 1488 1489 !Were 2DTE computed ? 1490 if(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and. rfuser==0 .and. rfmagn==0)then 1491 1492 write(message,'(a,a)' )ch10,' respfn : d/dk was computed, but no 2DTE, so no DDB output.' 1493 call wrtout([std_out, ab_out], message) 1494 1495 ! If 2DTE were computed, only one processor must output them and compute 1496 ! frequencies. 1497 else if(me==0)then 1498 1499 write(message,'(a,a)' )ch10,' ==> Compute Derivative Database <== ' 1500 call wrtout([std_out, ab_out], message) 1501 1502 ! In the RESPFN code, dfpt_nstdy and stady3 were called here 1503 d2nfr(:,:,:,:,:)=d2lo(:,:,:,:,:)+d2nl(:,:,:,:,:) 1504 if (psps%usepaw==1) d2nfr(:,:,:,:,:)=d2nfr(:,:,:,:,:)+d2ovl(:,:,:,:,:) 1505 1506 zero_by_symm=1 1507 if(dtset%rfmeth<0)zero_by_symm=0 1508 1509 ! In case of bbb decomposition 1510 if(prtbbb==1)then 1511 ABI_MALLOC(blkflg1,(3,mpert,3,mpert)) 1512 ABI_MALLOC(blkflg2,(3,mpert,3,mpert)) 1513 blkflg2(:,:,:,:) = blkflg(:,:,:,:) 1514 do ipert = 1, mpert 1515 do ipert2 = 1, mpert 1516 if ((ipert /= natom + 2).and.(ipert>natom).and.(ipert2/=natom+2)) then 1517 blkflg2(:,ipert2,:,ipert) = 0 1518 end if 1519 end do 1520 end do 1521 ABI_MALLOC(d2tmp,(2,3,mpert,3,mpert)) 1522 do iband = 1,dtset%mband 1523 d2tmp(:,:,:,:,:)=zero 1524 blkflg1(:,:,:,:) = blkflg2(:,:,:,:) 1525 d2tmp(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband) 1526 call d2sym3(blkflg1,d2tmp,indsym,mpert,natom,dtset%nsym,qphon,symq,& 1527 & symrec,dtset%symrel,timrev,zero_by_symm) 1528 d2bbb(:,:,:,:,iband,iband) = d2tmp(:,:,natom+2,:,:) 1529 end do 1530 ABI_FREE(blkflg1) 1531 ABI_FREE(blkflg2) 1532 ABI_FREE(d2tmp) 1533 end if 1534 1535 ! Complete the d2nfr matrix by symmetrization of the existing elements 1536 !write(std_out,*)"blkflg before d2sym3: ", blkflg 1537 call d2sym3(blkflg,d2nfr,indsym,mpert,natom,dtset%nsym,qphon,symq,symrec,dtset%symrel,timrev,zero_by_symm) 1538 !write(std_out,*)"blkflg after d2sym3: ", blkflg 1539 1540 if(rfphon==1.and.psps%n1xccc/=0)then 1541 ! Complete the dyfrx1 matrix by symmetrization of the existing elements 1542 call d2sym3(blkflgfrx1,dyfrx1,indsym,natom,natom,dtset%nsym,qphon,symq,symrec,dtset%symrel,timrev,zero_by_symm) 1543 end if 1544 1545 ! Note that d2sym3 usually complete the 2nd-order matrix 1546 ! with elements that are zero by symmetry, automatically, 1547 ! unless it has been explicitly asked not to do so. 1548 ! blkflg is then set to 1 for these matrix elements, even if there has be no calculation. 1549 1550 ! Add the frozen-wf (dyfrwf) part to the ewald part (dyew), 1551 ! the part 1 of the frozen wf part of the xc core correction 1552 ! (dyfrx1) and the non-frozen part (dynfr) to get the second-order 1553 ! derivative matrix (d2matr), then 1554 ! take account of the non-cartesian coordinates (d2cart). 1555 ABI_MALLOC(d2cart,(2,3,mpert,3,mpert)) 1556 ABI_MALLOC(carflg,(3,mpert,3,mpert)) 1557 ABI_MALLOC(d2matr,(2,3,mpert,3,mpert)) 1558 outd2=1 1559 call dfpt_gatherdy(asr,becfrnl,dtset%berryopt,blkflg,carflg,& 1560 & chneut,dyew,dyfrwf,dyfrx1,dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,& 1561 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 1562 & gprimd,dtset%mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,& 1563 & rfpert,rprimd,dtset%typat,ucvol,usevdw,psps%ziontypat) 1564 1565 ! Output of the dynamical matrix (master only) 1566 call dfpt_dyout(becfrnl,dtset%berryopt,blkflg,carflg,ddkfil,dyew,dyfrlo,& 1567 & dyfrnl,dyfrx1,dyfrx2,dyfr_cplex,dyfr_nondiag,dyvdw,d2cart,d2cart_bbb,d2eig0,& 1568 & d2k0,d2lo,d2loc0,d2matr,d2nl,d2nl0,d2nl1,d2ovl,d2vn,& 1569 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 1570 & has_full_piezo,has_allddk,ab_out,dtset%mband,mpert,natom,ntypat,& 1571 & outd2,pawbec,pawpiezo,piezofrnl,dtset%prtbbb,dtset%prtvol,qzero,& 1572 & dtset%typat,rfdir,rfpert,rfphon,rfstrs,psps%usepaw,usevdw,psps%ziontypat) 1573 1574 1575 ! Initialize ddb header object 1576 call ddb_hdr%init(dtset,psps,pawtab,& 1577 & dscrpt=' Note : temporary (transfer) database ',& 1578 & nblok=1,xred=xred,occ=occ,ngfft=ngfft) 1579 1580 ! Initialize ddb object 1581 call ddb%init(dtset, nblok=1, mpert=mpert, with_d2E=.true.) 1582 1583 ! Set the values for the 2nd order derivatives 1584 call ddb%set_qpt(iblok=1, qpt=qphon(1:3)) 1585 call ddb%set_d2matr(1, d2matr, blkflg) 1586 1587 ! Output dynamical matrix 1588 call ddb%write(ddb_hdr, dtfil%fnameabo_ddb) 1589 1590 ! Deallocate ddb object 1591 call ddb_hdr%free() 1592 call ddb%free() 1593 1594 ! In case of phonons, diagonalize the dynamical matrix 1595 if(rfphon==1)then 1596 1597 ! First, suppress the 'wings' elements, 1598 ! for which the diagonal element is not known 1599 call wings3(carflg,d2cart,mpert) 1600 1601 ! Check the analyticity of the dynamical matrix 1602 analyt=0 1603 if (rfpert(natom+2)==0 .or. rfpert(natom+2)==2 .or. sumg0==1 ) analyt=1 1604 1605 ! Diagonalize the analytic part 1606 ABI_MALLOC(displ,(2*3*natom*3*natom)) 1607 ABI_MALLOC(eigval,(3*natom)) 1608 ABI_MALLOC(eigvec,(2*3*natom*3*natom)) 1609 ABI_MALLOC(phfrq,(3*natom)) 1610 qphnrm=one 1611 call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,& 1612 & dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,& 1613 & dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol) 1614 1615 ! Print the phonon frequencies 1616 call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon) 1617 1618 ! Check the completeness of the dynamical matrix and eventually send a warning 1619 call chkph3(carflg,0,mpert,natom) 1620 end if ! end case of phonons 1621 end if !end me == 0 1622 1623 !Compute the other terms for AHC dynamic and AHC full 1624 if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0.or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and. rfuser==0 & 1625 & .and. rfmagn==0)) then 1626 if(rfphon==1) then ! AHC can only be computed in case of phonons 1627 1628 ! Stuff for parallelism 1629 if(master /= me) then 1630 ABI_MALLOC(phfrq,(3*natom)) 1631 ABI_MALLOC(displ,(2*3*natom*3*natom)) 1632 end if 1633 call xmpi_bcast (phfrq,master,mpi_enreg%comm_cell,ierr) !Broadcast phfrq and displ 1634 call xmpi_bcast (displ,master,mpi_enreg%comm_cell,ierr) !to all processes 1635 1636 if(dtset%ieig2rf == 3 .or. dtset%ieig2rf == 4 .or. dtset%ieig2rf == 5 ) then 1637 bdeigrf = dtset%bdeigrf 1638 if(dtset%bdeigrf == -1) bdeigrf = dtset%mband 1639 ! if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1640 ! & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1641 ! write(std_out,*)'Reading the dense grid WF file' 1642 ! call wfk_read_eigenvalues(dtfil%fnameabi_wfkfine,eigenq_fine,hdr_fine,mpi_enreg%comm_world) 1643 ! ABI_CHECK(SIZE(eigenq_fine,DIM=1)==Dtset%mband,"Size eigenq_fine != mband") 1644 ! endif 1645 if(dtset%kptopt==3 .or. dtset%kptopt==0 .or. dtset%nsym==1)then 1646 write(std_out,*) 'Entering: eig2tot' 1647 if(dtset%smdelta>0)then 1648 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1649 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1650 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1651 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1652 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1653 & occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine) 1654 else 1655 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1656 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1657 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1658 & occ_rbz_pert,hdr0,eigbrd) 1659 end if 1660 else 1661 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1662 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1663 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1664 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1665 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1666 & occ_rbz_pert,hdr0,eigbrd,eigenq_fine,hdr_fine) 1667 else 1668 call eig2tot(dtfil,xred,psps,pawtab,natom,bdeigrf,clflg,dim_eig2nkq,eigen0_pert,& 1669 & eigenq_pert,eigen1_pert,eig2nkq,dtset%elph2_imagden,dtset%esmear,dtset%ieig2rf,dtset%mband,mpert,& 1670 & dtset%natom,mpi_enreg,doccde,nkpt_rbz,dtset%nsppol,dtset%smdelta,rprimd,dtset,& 1671 & occ_rbz_pert,hdr0) 1672 end if 1673 end if 1674 write(std_out,*) 'Leaving: eig2tot' 1675 end if 1676 end if 1677 if (dtset%ieig2rf > 0) then 1678 ABI_FREE(eigen0_pert) 1679 ABI_FREE(eigenq_pert) 1680 ABI_FREE(occ_rbz_pert) 1681 ABI_FREE(eigen1_pert) 1682 call hdr0%free() 1683 if ((dtset%getwfkfine /= 0 .and. dtset%irdwfkfine ==0) .or.& 1684 & (dtset%getwfkfine == 0 .and. dtset%irdwfkfine /=0) ) then 1685 call hdr_fine%free() 1686 ABI_FREE(eigenq_fine) 1687 end if 1688 end if ! ieig2rf == 3 or %ieig2rf == 4 or %ieig2rf == 5 1689 end if ! rfphon==1 1690 end if 1691 if(dtset%efmas>0) then 1692 ABI_FREE(eigen0_pert) 1693 ABI_FREE(eigen1_pert) 1694 end if 1695 1696 ABI_FREE(doccde) 1697 1698 if(me==0)then 1699 if (.not.(rfphon==0 .and. (rf2_dkdk/=0 .or. rf2_dkde/=0 .or. rfddk/=0 .or. rfelfd==2) .and. rfstrs==0 .and.rfuser==0 & 1700 & .and. rfmagn==0) )then 1701 if(rfphon==1)then 1702 ! Compute and print the T=0 Fan, and possibly DDW contributions to the eigenenergies. 1703 if(dtset%ieig2rf > 0) then 1704 write(message, '(80a,9a)' ) ('=',mu=1,80),ch10,ch10,& 1705 & ' ---- T=0 shift of eigenenergies due to electron-phonon interation at q ---- ',ch10,& 1706 & ' Warning : the total shift must be computed through anaddb, ',ch10,& 1707 & ' here, only the contribution of one q point is printed. ',ch10,& 1708 & ' Print first the electronic eigenvalues, then the q-dependent Fan shift of eigenvalues.' 1709 call wrtout([std_out, ab_out], message) 1710 1711 if(qeq0)then 1712 write(message, '(a)' )' Phonons at gamma, also compute the Diagonal Debye-Waller shift of eigenvalues.' 1713 call wrtout([std_out, ab_out], message) 1714 end if 1715 1716 write(message, '(a)' ) ' ' 1717 call wrtout([std_out, ab_out], message) 1718 1719 call prteigrs(eigen0,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1720 & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,3,0,dtset%prtvol,& 1721 & eigen0,zero,zero,dtset%wtk) 1722 1723 write(message, '(a)' ) ch10 1724 call wrtout([std_out, ab_out], message) 1725 1726 ! Compute and print Fan contribution 1727 ABI_MALLOC(eigen_fan,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1728 ABI_MALLOC(eigen_fan_mean,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1729 call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_fan,gprimd,& 1730 & dtset%mband,natom,dtset%nkpt,dtset%nsppol,1,phfrq,dtset%prtvol) 1731 call eigen_meandege(eigen0,eigen_fan,eigen_fan_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2) 1732 call prteigrs(eigen_fan_mean,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1733 & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,5,0,dtset%prtvol,& 1734 & eigen0,zero,zero,dtset%wtk) 1735 1736 if(qeq0 .or. dtset%getgam_eig2nkq>0)then 1737 1738 write(message, '(a)' ) ch10 1739 call wrtout([std_out, ab_out], message) 1740 1741 ! Compute and print Diagonal Debye-Waller contribution 1742 ABI_MALLOC(eigen_ddw,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1743 ABI_MALLOC(eigen_ddw_mean,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1744 if(qeq0)then 1745 call elph2_fanddw(dim_eig2nkq,displ,eig2nkq,eigen_ddw,gprimd,& 1746 & dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol) 1747 if(results_respfn%gam_jdtset == -dtset%jdtset)then 1748 sz1=dtset%mband*dtset%nsppol 1749 sz2=natom*dim_eig2nkq 1750 ABI_MALLOC(results_respfn%gam_eig2nkq,(2,sz1,dtset%nkpt,3,natom,3,sz2)) 1751 results_respfn%gam_eig2nkq(:,:,:,:,:,:,:)=eig2nkq(:,:,:,:,:,:,:) 1752 results_respfn%gam_jdtset=dtset%jdtset 1753 end if 1754 else if(dtset%getgam_eig2nkq>0)then 1755 if(results_respfn%gam_jdtset==dtset%getgam_eig2nkq)then 1756 call elph2_fanddw(dim_eig2nkq,displ,results_respfn%gam_eig2nkq,eigen_ddw,& 1757 & gprimd,dtset%mband,natom,dtset%nkpt,dtset%nsppol,2,phfrq,dtset%prtvol) 1758 else 1759 write(message,'(a,i0,2a,i0,2a)')& 1760 & 'results_respfn%gam_jdtset=',results_respfn%gam_jdtset,ch10,& 1761 & 'dtset%getgam_eig2nkq=',dtset%getgam_eig2nkq,ch10,& 1762 & 'So, it seems eig2nkq at gamma has not yet been computed, while it is needed now.' 1763 ABI_BUG(message) 1764 end if 1765 end if 1766 call eigen_meandege(eigen0,eigen_ddw,eigen_ddw_mean,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol,2) 1767 call prteigrs(eigen_ddw_mean,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1768 & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,6,0,dtset%prtvol,& 1769 & eigen0,zero,zero,dtset%wtk) 1770 write(message, '(a)' ) ch10 1771 call wrtout([std_out, ab_out], message) 1772 1773 ! Print sum of mean Fan and DDW 1774 ABI_MALLOC(eigen_fanddw,(dtset%mband*dtset%nkpt*dtset%nsppol)) 1775 eigen_fanddw=eigen_fan_mean+eigen_ddw_mean 1776 call prteigrs(eigen_fanddw,dtset%enunit,fermie,fermih,dtfil%fnameabo_eig,ab_out,-1,dtset%kptns,dtset%kptopt,& 1777 & dtset%mband,dtset%nband,dtset%nbdbuf,dtset%nkpt,1,dtset%nsppol,occ,dtset%occopt,7,0,dtset%prtvol,& 1778 & eigen0,zero,zero,dtset%wtk) 1779 1780 ABI_FREE(eigen_ddw) 1781 ABI_FREE(eigen_ddw_mean) 1782 ABI_FREE(eigen_fanddw) 1783 1784 end if 1785 1786 ABI_FREE(eigen_fan) 1787 ABI_FREE(eigen_fan_mean) 1788 end if 1789 1790 ! In case of a non-analytical part, 1791 ! get the phonon frequencies for three different directions (in cartesian coordinates) 1792 if(analyt==0)then 1793 qphnrm=zero 1794 do idir=1,3 1795 ! Need to know the corresponding dielectric constant 1796 if(carflg(idir,natom+2,idir,natom+2)==1)then 1797 qphon(:)=zero ; qphon(idir)=one 1798 ! Get the phonon frequencies 1799 call dfpt_phfrq(dtset%amu_orig(:,1),displ,d2cart,eigval,eigvec,indsym,mpert,& 1800 & dtset%nsym,natom,dtset%nsym,ntypat,phfrq,qphnrm,qphon,& 1801 & dtset%rprimd_orig(1:3,1:3,1),0,dtset%symrel,dtset%symafm,dtset%typat,ucvol) 1802 ! Print the phonon frequencies 1803 call dfpt_prtph(displ,0,dtset%enunit,ab_out,natom,phfrq,qphnrm,qphon) 1804 ! Check the completeness of the dynamical matrix 1805 ! and eventually send a warning 1806 call chkph3(carflg,idir,mpert,natom) 1807 end if 1808 end do 1809 if (idir < 4) then 1810 qphon(idir)=zero 1811 end if 1812 end if 1813 1814 ABI_FREE(displ) 1815 ABI_FREE(eigval) 1816 ABI_FREE(eigvec) 1817 ABI_FREE(phfrq) 1818 end if ! rfphon == 1 1819 ABI_FREE(carflg) 1820 ABI_FREE(d2cart) 1821 ABI_FREE(d2matr) 1822 end if ! End condition on if.not. 1823 end if ! master node 1824 1825 !Deallocate arrays 1826 ABI_SFREE(displ) 1827 ABI_SFREE(eigval) 1828 ABI_SFREE(eigvec) 1829 ABI_SFREE(phfrq) 1830 1831 ABI_FREE(clflg) 1832 ABI_FREE(atindx) 1833 ABI_FREE(atindx1) 1834 ABI_FREE(blkflg) 1835 ABI_FREE(dyew) 1836 ABI_FREE(dyewq0) 1837 ABI_FREE(dyfrlo) 1838 ABI_FREE(dyfrnl) 1839 ABI_FREE(dyfrwf) 1840 ABI_FREE(dyfrx1) 1841 ABI_FREE(dyfrx2) 1842 ABI_FREE(dyvdw) 1843 ABI_FREE(d2bbb) 1844 ABI_FREE(d2cart_bbb) 1845 ABI_FREE(d2eig0) 1846 ABI_FREE(d2k0) 1847 ABI_FREE(d2lo) 1848 ABI_FREE(d2loc0) 1849 ABI_FREE(d2nfr) 1850 ABI_FREE(d2nl) 1851 ABI_FREE(d2nl0) 1852 ABI_FREE(d2nl1) 1853 ABI_FREE(d2ovl) 1854 ABI_FREE(d2vn) 1855 ABI_FREE(eigen0) 1856 ABI_FREE(eig2nkq) 1857 ABI_FREE(eigbrd) 1858 ABI_FREE(eltcore) 1859 ABI_FREE(elteew) 1860 ABI_FREE(eltfrhar) 1861 ABI_FREE(eltfrnl) 1862 ABI_FREE(eltfrloc) 1863 ABI_FREE(eltfrkin) 1864 ABI_FREE(eltfrxc) 1865 ABI_FREE(eltvdw) 1866 ABI_FREE(becfrnl) 1867 ABI_FREE(piezofrnl) 1868 call efmasdeg_free_array(efmasdeg) 1869 call efmasval_free_array(efmasval) 1870 ABI_FREE(grxc) 1871 ABI_FREE(indsym) 1872 ABI_FREE(kxc) 1873 ABI_FREE(nattyp) 1874 ABI_FREE(pertsy) 1875 ABI_FREE(rfpert) 1876 ABI_FREE(rhog) 1877 ABI_FREE(rhor) 1878 ABI_FREE(taug) 1879 ABI_FREE(taur) 1880 ABI_FREE(symq) 1881 ABI_FREE(symrec) 1882 ABI_FREE(symrel_cart) 1883 ABI_FREE(vtrial) 1884 ABI_FREE(ylm) 1885 ABI_FREE(ylmgr) 1886 call pawfgr_destroy(pawfgr) 1887 if (psps%usepaw==1) then 1888 call pawrhoij_free(pawrhoij) 1889 call paw_an_free(paw_an) 1890 call paw_ij_free(paw_ij) 1891 call pawfgrtab_free(pawfgrtab) 1892 end if 1893 ABI_FREE(nhat) 1894 ABI_FREE(nhatgr) 1895 ABI_FREE(pawrhoij) 1896 ABI_FREE(paw_an) 1897 ABI_FREE(paw_ij) 1898 ABI_FREE(pawfgrtab) 1899 if(rfphon==1.and.psps%n1xccc/=0)then 1900 ABI_FREE(blkflgfrx1) 1901 end if 1902 1903 ! Clean the header 1904 call hdr%free() 1905 1906 !Clean GPU data 1907 #if defined HAVE_GPU_CUDA 1908 if (dtset%gpu_option==ABI_GPU_LEGACY .or. dtset%gpu_option==ABI_GPU_KOKKOS) then 1909 call dealloc_hamilt_gpu(0,dtset%gpu_option) 1910 end if 1911 #endif 1912 1913 call timab(138,2,tsec) 1914 call timab(132,2,tsec) 1915 1916 DBG_EXIT("COLL") 1917 1918 end subroutine respfn
m_respfn_driver/wrtloctens [ Functions ]
[ Top ] [ m_respfn_driver ] [ Functions ]
NAME
wrtloctens
FUNCTION
Output of the localisation tensor
INPUTS
blkflg = flags for each element of the 2DTE (=1 if computed) d2bbb = band by band decomposition of second order derivatives d2nl = non-local contributions to the 2DTEs mband = maximum number of bands mpert = maximum number of ipert natom = number of atoms in unit cell prtbbb = if = 1, write the band by band decomposition of the localization tensor rprimd = dimensional primitive translations for real space (bohr) usepaw = flag for PAW
OUTPUT
(only writing)
TODO
The localization tensor cannot be defined in the metallic case. It should not be computed.
SOURCE
1947 subroutine wrtloctens(blkflg,d2bbb,d2nl,mband,mpert,natom,prtbbb,rprimd,usepaw) 1948 1949 !Arguments ------------------------------------ 1950 !scalars 1951 integer,intent(in) :: mband,mpert,natom,prtbbb,usepaw 1952 !arrays 1953 integer,intent(in) :: blkflg(3,mpert,3,mpert) 1954 real(dp),intent(in) :: rprimd(3,3) 1955 real(dp),intent(inout) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb) 1956 real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert) 1957 1958 !Local variables ------------------------------ 1959 !scalars 1960 integer :: flag,iband,idir,idir2,jband 1961 character(len=500) :: message 1962 !arrays 1963 real(dp) :: loctenscart(2,3,3) 1964 real(dp),allocatable :: loctenscart_bbb(:,:,:,:,:) 1965 1966 ! ********************************************************************* 1967 1968 !This feature is disabled in the PAW case 1969 if (usepaw==1) return 1970 1971 if(prtbbb==1)then 1972 ABI_MALLOC(loctenscart_bbb,(2,3,3,mband,mband*prtbbb)) 1973 loctenscart_bbb(:,:,:,:,:)=zero 1974 end if 1975 1976 !complete missing elements 1977 flag = 0 1978 do idir2 = 1,3 1979 do idir = 1,3 1980 1981 if (blkflg(idir2,natom+1,idir,natom+1)==0) then 1982 if (blkflg(idir,natom+1,idir2,natom+1)==0) then 1983 flag = 1 1984 else 1985 d2nl(1,idir2,natom+1,idir,natom+1) = d2nl(1,idir,natom+1,idir2,natom+1) 1986 d2nl(2,idir2,natom+1,idir,natom+1) =-d2nl(2,idir,natom+1,idir2,natom+1) 1987 if(prtbbb==1)then 1988 d2bbb(1,idir2,idir,natom+1,:,:) = d2bbb(1,idir,idir2,natom+1,:,:) 1989 d2bbb(2,idir2,idir,natom+1,:,:) =-d2bbb(2,idir,idir2,natom+1,:,:) 1990 end if 1991 1992 end if 1993 end if 1994 1995 end do ! idir=1,3 1996 end do ! idir2=1,3 1997 1998 !Transform the tensor to cartesian coordinates 1999 2000 loctenscart(1,:,:) = matmul(rprimd,d2nl(1,:,natom+1,:,natom+1)) 2001 loctenscart(2,:,:) = matmul(rprimd,d2nl(2,:,natom+1,:,natom+1)) 2002 2003 loctenscart(1,:,:) = matmul(loctenscart(1,:,:),transpose(rprimd)) 2004 loctenscart(2,:,:) = matmul(loctenscart(2,:,:),transpose(rprimd)) 2005 2006 loctenscart(:,:,:) = loctenscart(:,:,:)/(two_pi**2) 2007 2008 if (prtbbb == 1) then 2009 2010 write(message,'(a,a)')ch10, ' Band by band decomposition of the localisation tensor (bohr^2)' 2011 call wrtout([std_out, ab_out], message) 2012 2013 do iband = 1,mband 2014 do jband = 1,mband 2015 2016 loctenscart_bbb(1,:,:,iband,jband) = matmul(rprimd,d2bbb(1,:,:,natom+1,iband,jband)) 2017 loctenscart_bbb(2,:,:,iband,jband) = matmul(rprimd,d2bbb(2,:,:,natom+1,iband,jband)) 2018 2019 loctenscart_bbb(1,:,:,iband,jband) = matmul(loctenscart_bbb(1,:,:,iband,jband),transpose(rprimd)) 2020 loctenscart_bbb(2,:,:,iband,jband) = matmul(loctenscart_bbb(2,:,:,iband,jband),transpose(rprimd)) 2021 2022 loctenscart_bbb(:,:,:,iband,jband) = loctenscart_bbb(:,:,:,iband,jband)/(two_pi**2) 2023 2024 2025 write(message,'(a,a,i5,a,i5,a)')ch10, & 2026 ' Localisation tensor (bohr^2) for band ',iband,',',jband, & 2027 ' in cartesian coordinates' 2028 call wrtout([std_out, ab_out], message) 2029 2030 write(ab_out,*)' direction matrix element' 2031 write(ab_out,*)' alpha beta real part imaginary part' 2032 do idir2 = 1,3 2033 do idir = 1,3 2034 write(ab_out,'(5x,i1,8x,i1,3x,2f16.10)')idir2,idir,& 2035 loctenscart_bbb(1,idir2,idir,iband,jband),& 2036 loctenscart_bbb(2,idir2,idir,iband,jband) 2037 end do 2038 end do 2039 2040 end do !jband 2041 end do !iband 2042 2043 end if !prtbbb 2044 2045 if (usepaw==0) then 2046 write(message,'(a,a,a,a)')ch10, & 2047 & ' Total localisation tensor (bohr^2) in cartesian coordinates',ch10,& 2048 & ' WARNING : still subject to testing - especially symmetries.' 2049 else 2050 write(message,'(a,a,a,a,a,a)')ch10, & 2051 & ' Total localisation tensor (bohr^2) in cartesian coordinates',ch10,& 2052 & ' WARNING : probably wrong for PAW (printing for testing purpose)',ch10,& 2053 & ' WARNING : still subject to testing - especially symmetries.' 2054 end if 2055 call wrtout([std_out, ab_out], message) 2056 2057 write(ab_out,*)' direction matrix element' 2058 write(ab_out,*)' alpha beta real part imaginary part' 2059 do idir2 = 1,3 2060 do idir = 1,3 2061 write(ab_out,'(5x,i1,8x,i1,3x,2f16.10)')idir2,idir,& 2062 & loctenscart(1,idir2,idir),& 2063 & loctenscart(2,idir2,idir) 2064 end do 2065 end do 2066 2067 if (flag == 1) then 2068 write(message,'(6a)')ch10,& 2069 & ' WARNING : Localization tensor calculation (this does not apply to other properties).',ch10,& 2070 & ' Not all d/dk perturbations were computed. So the localization tensor in reciprocal space is incomplete,',ch10,& 2071 & ' and transformation to cartesian coordinates may be wrong. Check input variable rfdir.' 2072 call wrtout([std_out, ab_out], message) 2073 end if 2074 2075 ABI_SFREE(loctenscart_bbb) 2076 2077 end subroutine wrtloctens