TABLE OF CONTENTS
ABINIT/dfptnl_resp [ Functions ]
NAME
dfptnl_resp
FUNCTION
Compute the linear response part to the 3dte
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (MVeithen) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave
coefficients of wavefunctions
cg1 = first-order wavefunction relative to the perturbations i1pert
cg3 = first-order wavefunction relative to the perturbations i3pert
cplex= if 1, real space 1-order functions on FFT grid are REAL,
if 2, COMPLEX
dtfil <type(datafiles_type)>=variables related to files
dtset <type(dataset_type)>=all input variables for this dataset
i1dir,i2dir,i3dir=directions of the corresponding perturbations
i1pert,i2pert,i3pert = type of perturbation that has to be computed
kg(3,mpw*mkmem)=reduced planewave coordinates
mband = maximum number of bands
mgfft=maximum size of 1D FFTs
mkmem = maximum number of k points which can fit in core memory
mk1mem = maximum number of k points for first-order WF
which can fit in core memory
mpert =maximum number of ipert
mpi_enreg=MPI-parallelisation information
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
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
nspinor = number of spinorial components of the wavefunctions
nsppol = number of channels for spin-polarization (1 or 2)
npwarr(nkpt) = array holding npw for each k point
occ(mband*nkpt*nsppol) = occupation number for each band and k
ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
psps <type(pseudopotential_type)> = variables related to pseudopotentials
rprimd(3,3) = dimensional primitive translations (bohr)
vtrial1(cplex*nfft,nspden)=firs-order local potential
xred(3,natom) = reduced atomic coordinates
ylm(mpw*mkmem,mpsang*mpsang*useylm)= spherical harmonics for
each G and k point
OUTPUT
d3lo(2,3,mpert,3,mpert,3,mpert) = matrix of the 3DTEs
PARENTS
dfptnl_loop
CHILDREN
destroy_hamiltonian,dotprod_g,fftpac,fourwf,init_hamiltonian load_k_hamiltonian,mkffnl,mkkpg,nonlop,status,xmpi_sum
SOURCE
65 #if defined HAVE_CONFIG_H 66 #include "config.h" 67 #endif 68 69 #include "abi_common.h" 70 71 72 subroutine dfptnl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,& 73 & i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,& 74 & kg,mband,mgfft,mkmem,mk1mem,& 75 & mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,nspinor,nsppol,& 76 & npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm) 77 78 use defs_basis 79 use defs_datatypes 80 use defs_abitypes 81 use m_profiling_abi 82 use m_xmpi 83 84 use m_cgtools, only : dotprod_g 85 use m_pawtab, only : pawtab_type 86 use m_pawcprj, only : pawcprj_type 87 use m_hamiltonian,only : init_hamiltonian,destroy_hamiltonian,& 88 & load_k_hamiltonian,gs_hamiltonian_type 89 90 !This section has been created automatically by the script Abilint (TD). 91 !Do not modify the following lines by hand. 92 #undef ABI_FUNC 93 #define ABI_FUNC 'dfptnl_resp' 94 use interfaces_32_util 95 use interfaces_53_ffts 96 use interfaces_66_nonlocal 97 !End of the abilint section 98 99 implicit none 100 101 !Arguments ------------------------------------ 102 !scalars 103 integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft 104 integer,intent(in) :: mk1mem,mkmem,mpert,mpsang,mpw,natom,nfft,nkpt,nspden 105 integer,intent(in) :: nspinor,nsppol 106 type(MPI_type),intent(in) :: mpi_enreg 107 type(datafiles_type),intent(in) :: dtfil 108 type(dataset_type),intent(in) :: dtset 109 type(pseudopotential_type),intent(in) :: psps 110 !arrays 111 integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt) 112 real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol) 113 real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol) 114 real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol) 115 real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom),rprimd(3,3) 116 real(dp),intent(in) :: xred(3,natom),ylm(mpw*mkmem,mpsang*mpsang*psps%useylm) 117 real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden) 118 real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert) 119 120 !Local variables------------------------------- 121 !scalars 122 integer,parameter :: level=52 123 integer :: bantot,choice,counter,cpopt,dimffnl,iband,icg0,ider,ierr,iexit 124 integer :: ii,ikg,ikpt,ilm,ipw,isppol,istwf_k,jband,jj 125 integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nnlout,npw_k 126 integer :: option,paw_opt,signs,spaceComm,tim_fourwf,tim_nonlop 127 real(dp) :: dot1i,dot1r,dot2i,dot2r,doti,dotr,lagi,lagr,sumi,sumr,weight 128 type(gs_hamiltonian_type) :: gs_hamk 129 !arrays 130 integer,allocatable :: kg_k(:,:) 131 real(dp) :: buffer(2),enlout(3),kpq(3),kpt(3) 132 real(dp) :: dum_svectout(1,1),dum(1),rmet(3,3),ylmgr_dum(1,1,1) 133 real(dp),allocatable :: cwave0(:,:),cwavef3(:,:),ffnlk(:,:,:,:) 134 real(dp),allocatable :: gh0(:,:),gh1(:,:),gvnl(:,:),kpg_k(:,:) 135 real(dp),allocatable :: vlocal1(:,:,:),wfraug(:,:,:,:),ylm_k(:,:) 136 type(pawcprj_type) :: cprj_dum(1,1) 137 type(pawtab_type) :: pawtab_dum(0) 138 139 !*********************************************************************** 140 141 me = mpi_enreg%me 142 spaceComm=mpi_enreg%comm_cell 143 144 call status(0,dtfil%filstat,iexit,level,'enter ') 145 146 bantot = 0 147 icg0 = 0 148 option = 2 149 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 150 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 151 152 ABI_ALLOCATE(vlocal1,(cplex*n4,n5,n6)) 153 ABI_ALLOCATE(wfraug,(2,n4,n5,n6)) 154 155 !Initialize Hamiltonian (k-independent terms) - NCPP only 156 call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,& 157 & dtset%typat,xred,nfft,mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 158 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 159 !& paw_ij=paw_ij) 160 rmet = MATMUL(TRANSPOSE(rprimd),rprimd) 161 162 sumr = zero ; sumi = zero 163 164 !Loop over spins 165 166 do isppol = 1, nsppol 167 168 call status(0,dtfil%filstat,iexit,level,'call fftpac ') 169 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,vtrial1,vlocal1,option) 170 171 ! Loop over k-points 172 173 ikg = 0 174 do ikpt = 1, nkpt 175 176 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me))cycle 177 178 counter = 100*ikpt 179 180 nband_k = dtset%nband(ikpt+(isppol-1)*nkpt) 181 npw_k = npwarr(ikpt) 182 istwf_k = dtset%istwfk(ikpt) 183 184 kpt(:) = dtset%kptns(:,ikpt) 185 kpq(:) = dtset%kptns(:,ikpt) ! In case of non zero q, kpt = kpt + q 186 187 ABI_ALLOCATE(cwave0,(2,npw_k*dtset%nspinor)) 188 ABI_ALLOCATE(cwavef3,(2,npw_k*dtset%nspinor)) 189 ABI_ALLOCATE(gh0,(2,npw_k*dtset%nspinor)) 190 ABI_ALLOCATE(gvnl,(2,npw_k*dtset%nspinor)) 191 ABI_ALLOCATE(gh1,(2,npw_k*dtset%nspinor)) 192 193 ABI_ALLOCATE(kg_k,(3,npw_k)) 194 ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang*psps%useylm)) 195 kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg) 196 if (psps%useylm==1) then 197 do ilm=1,mpsang*mpsang 198 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 199 end do 200 end if 201 202 ! Compute (k+G) and (k+q+G) vectors (only if useylm=1) 203 nkpg=0;if (i2pert<natom+1) nkpg=3*dtset%nloalg(3) 204 ABI_ALLOCATE(kpg_k,(npw_k,nkpg)) 205 if (nkpg>0) then 206 call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k) 207 end if 208 209 ! Compute nonlocal form factors ffnl at (k+G), for all atoms 210 dimffnl=1 211 ABI_ALLOCATE(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 212 if (i2pert<natom+1) then 213 ider=0 214 call status(counter,dtfil%filstat,iexit,level,'call mkffnl ') 215 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,& 216 & ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,& 217 & psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,& 218 & psps%usepaw,psps%useylm,ylm_k,ylmgr_dum) 219 end if 220 221 ! Load k-dependent part in the Hamiltonian datastructure 222 call load_k_hamiltonian(gs_hamk,kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,& 223 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk,compute_gbound=.true.) 224 ! Load k+q-dependent part in the Hamiltonian datastructure 225 ! call load_kprime_hamiltonian... !! To be activated when q/=0 226 227 ! Loop over bands 228 229 do iband = 1,nband_k 230 231 cwave0(:,:)=cg(:,1+(iband - 1)*npw_k*dtset%nspinor+icg0:& 232 & iband*npw_k*dtset%nspinor+icg0) 233 cwavef3(:,:)=cg3(:,1+(iband-1)*npw_k*dtset%nspinor+icg0:& 234 & iband*npw_k*dtset%nspinor+icg0) 235 236 ! Compute vtrial1 | cwafef3 > 237 tim_fourwf = 0 ; weight = one 238 call status(counter,dtfil%filstat,iexit,level,'call fourwf ') 239 call fourwf(cplex,vlocal1,cwavef3,gh1,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,& 240 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,& 241 & dtset%paral_kgb,tim_fourwf,weight,weight,& 242 & use_gpu_cuda=dtset%use_gpu_cuda) 243 244 ! In case i2pert = phonon-type perturbation 245 ! add first-order change in the nonlocal potential 246 if (i2pert<natom+1) then 247 signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1 248 call status(counter,dtfil%filstat,iexit,level,'call nonlop ') 249 call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,& 250 & signs,dum_svectout,tim_nonlop,cwavef3,gvnl,iatom_only=i2pert) 251 gh1(:,:) = gh1(:,:) + gvnl(:,:) 252 end if 253 254 ii = (iband-1)*npw_k*dtset%nspinor + icg0 255 call dotprod_g(dotr,doti,istwf_k,npw_k,2,cg1(:,ii+1:ii+npw_k),gh1,mpi_enreg%me_g0,xmpi_comm_self) 256 257 ! Compute vtrial1 | cwave0 > 258 tim_fourwf = 0 ; weight = one 259 call status(counter,dtfil%filstat,iexit,level,'call fourwf ') 260 call fourwf(cplex,vlocal1,cwave0,gh0,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,& 261 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,& 262 & dtset%paral_kgb,tim_fourwf,weight,weight,use_gpu_cuda=dtset%use_gpu_cuda) 263 264 ! In case i2pert = phonon-type perturbation 265 ! add first-order change in the nonlocal potential 266 if (i2pert<natom+1) then 267 signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1 268 call status(counter,dtfil%filstat,iexit,level,'call nonlop ') 269 call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,& 270 & signs,dum_svectout,tim_nonlop,cwave0,gvnl,iatom_only=i2pert) 271 gh0(:,:) = gh0(:,:) + gvnl(:,:) 272 end if 273 274 ! Compute the dft contribution to the Lagrange multiplier 275 ! cwavef3 and cwave0 have been transferred to gh1 and gh0 276 ! these vectors will be used to store the wavefunctions of band iband 277 ! cg1 and gh0 contain the wavefunctions of band jband 278 279 lagr = zero ; lagi = zero 280 do jband = 1, nband_k 281 282 ii = (jband - 1)*npw_k*dtset%nspinor + icg0 283 jj = (iband - 1)*npw_k*dtset%nspinor + icg0 284 285 ! dot1r and dot1i contain < u_mk | v^(1) | u_nk > 286 ! dot2r and dot2i contain < u_nk^(1) | u_mk^(1) > 287 ! m -> jband and n -> iband 288 289 dot1r = zero ; dot1i = zero 290 dot2r = zero ; dot2i = zero 291 do ipw = 1, npw_k 292 ii = ii + 1 ; jj = jj + 1 293 dot1r = dot1r + cg(1,ii)*gh0(1,ipw) + cg(2,ii)*gh0(2,ipw) 294 dot1i = dot1i + cg(1,ii)*gh0(2,ipw) - cg(2,ii)*gh0(1,ipw) 295 dot2r = dot2r + cg1(1,jj)*cg3(1,ii) + & 296 & cg1(2,jj)*cg3(2,ii) 297 dot2i = dot2i + cg1(1,jj)*cg3(2,ii) - & 298 & cg1(2,jj)*cg3(1,ii) 299 end do ! ipw 300 301 lagr = lagr + dot1r*dot2r - dot1i*dot2i 302 lagi = lagi + dot1r*dot2i + dot1i*dot2r 303 304 end do ! jband 305 306 sumr = sumr + & 307 & dtset%wtk(ikpt)*occ(bantot+iband)*(dotr-lagr) 308 sumi = sumi + & 309 & dtset%wtk(ikpt)*occ(bantot+iband)*(doti-lagi) 310 311 end do ! end loop over bands 312 313 bantot = bantot + nband_k 314 icg0 = icg0 + npw_k*dtset%nspinor*nband_k 315 ikg = ikg + npw_k 316 317 ABI_DEALLOCATE(cwave0) 318 ABI_DEALLOCATE(cwavef3) 319 ABI_DEALLOCATE(gh0) 320 ABI_DEALLOCATE(gh1) 321 ABI_DEALLOCATE(gvnl) 322 ABI_DEALLOCATE(kg_k) 323 ABI_DEALLOCATE(ylm_k) 324 ABI_DEALLOCATE(ffnlk) 325 ABI_DEALLOCATE(kpg_k) 326 327 end do ! end loop over k-points 328 329 end do ! end loop over spins 330 331 if (xmpi_paral == 1) then 332 buffer(1) = sumr ; buffer(2) = sumi 333 call xmpi_sum(buffer,spaceComm,ierr) 334 sumr = buffer(1) ; sumi = buffer(2) 335 end if 336 337 338 d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr 339 !d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi 340 341 !In some cases, the imaginary part is /= 0 because of the 342 !use of time reversal symmetry 343 d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero 344 345 call destroy_hamiltonian(gs_hamk) 346 347 ABI_DEALLOCATE(vlocal1) 348 ABI_DEALLOCATE(wfraug) 349 350 call status(0,dtfil%filstat,iexit,level,'exit ') 351 352 end subroutine dfptnl_resp