TABLE OF CONTENTS
ABINIT/dfpt_nstdy [ Functions ]
NAME
dfpt_nstdy
FUNCTION
This routine compute the non-stationary expression for the second derivative of the total energy, for a whole row of mixed derivatives. Only for norm-conserving pseudopotentials (no PAW)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG, DCA, GMR, MM, AR, MV, MB, MT) 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
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q. cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree) eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues gmet(3,3)=reciprocal space metric tensor in bohr**-2. gsqcut=cutoff on (k+G)^2 (bohr^-2) idir=direction of the perturbation indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points indsy1(4,nsym1,natom)=indirect indexing array for atom labels ipert=type of the perturbation istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates. kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ kxc(nfft,nkxc)=exchange and correlation kernel mkmem =number of k points treated by this node (GS data) mk1mem =number of k points treated by this node (RF data) mpert =maximum number of ipert mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs). nattyp(ntypat)= # atoms of each type. nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt=number of k points in the full BZ nkpt_rbz=number of k points in the reduced BZ for this perturbation nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym1=number of symmetry elements in space group consistent with i perturbation occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k in the reduced Brillouin zone (usually =2) ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information psps <type(pseudopotential_type)>=variables related to pseudopotentials rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3. rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) symrc1(3,3,nsym1)=symmetry operations in reciprocal space ucvol=unit cell volume in bohr**3. wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point
OUTPUT
blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed) d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some second order derivatives d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs
NOTES
Note that the ddk perturbation should not be treated here.
PARENTS
dfpt_scfcv
CHILDREN
appdig,destroy_hamiltonian,dfpt_mkcore,dfpt_mkvxc,dfpt_mkvxc_noncoll dfpt_nstwf,dfpt_sygra,dfpt_vlocal,dotprod_vn,init_hamiltonian load_spin_hamiltonian,mati3inv,timab,wfk_close,wfk_open_read,wrtout xmpi_sum
SOURCE
92 #if defined HAVE_CONFIG_H 93 #include "config.h" 94 #endif 95 96 #include "abi_common.h" 97 98 subroutine dfpt_nstdy(atindx,blkflg,cg,cg1,cplex,dtfil,dtset,d2bbb,d2lo,d2nl,eigen0,eigen1,& 99 & gmet,gsqcut,idir,indkpt1,indsy1,ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mkmem,mk1mem,& 100 & mpert,mpi_enreg,mpw,mpw1,nattyp,nband_rbz,nfft,ngfft,nkpt,nkpt_rbz,nkxc,& 101 & npwarr,npwar1,nspden,nsppol,nsym1,occ_rbz,ph1d,psps,rhor1,rmet,rprimd,& 102 & symrc1,ucvol,wtk_rbz,xred,ylm,ylm1,rhor,vxc) 103 104 use defs_basis 105 use defs_datatypes 106 use defs_abitypes 107 use m_profiling_abi 108 use m_xmpi 109 use m_errors 110 use m_wfk 111 use m_nctk 112 use m_hamiltonian 113 114 use m_io_tools, only : file_exists 115 use m_cgtools, only : dotprod_vn 116 use m_hdr, only : hdr_skip 117 use m_pawtab, only : pawtab_type 118 119 !This section has been created automatically by the script Abilint (TD). 120 !Do not modify the following lines by hand. 121 #undef ABI_FUNC 122 #define ABI_FUNC 'dfpt_nstdy' 123 use interfaces_14_hidewrite 124 use interfaces_18_timing 125 use interfaces_32_util 126 use interfaces_56_xc 127 use interfaces_72_response, except_this_one => dfpt_nstdy 128 !End of the abilint section 129 130 implicit none 131 132 !Arguments ------------------------------- 133 !scalars 134 integer,intent(in) :: cplex,idir,ipert,mk1mem,mkmem,mpert,mpw,mpw1,nfft,nkpt,nkpt_rbz,nkxc,nspden,nsppol,nsym1 135 real(dp),intent(in) :: gsqcut,ucvol 136 type(MPI_type),intent(in) :: mpi_enreg 137 type(datafiles_type),intent(in) :: dtfil 138 type(dataset_type),intent(in) :: dtset 139 type(pseudopotential_type),intent(in) :: psps 140 !arrays 141 integer,intent(in) :: atindx(dtset%natom),indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom) 142 integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem) 143 integer,intent(in) :: nattyp(dtset%ntypat),nband_rbz(nkpt_rbz*nsppol),ngfft(18) 144 integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz),symrc1(3,3,nsym1) 145 integer,intent(inout) :: blkflg(3,mpert,3,mpert) !vz_i 146 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 147 real(dp),intent(in) :: cg1(2,mpw1*dtset%nspinor*dtset%mband*mk1mem*nsppol) 148 real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*nsppol) 149 real(dp),intent(in) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*nsppol) 150 real(dp),intent(in) :: gmet(3,3),kpt_rbz(3,nkpt_rbz) 151 real(dp),intent(in) :: kxc(nfft,nkxc),occ_rbz(dtset%mband*nkpt_rbz*nsppol) 152 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 153 real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3) 154 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,dtset%natom) 155 real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm) 156 real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm) 157 real(dp),intent(inout) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*dtset%prtbbb)!vz_i 158 real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert),d2nl(2,3,mpert,3,mpert) !vz_i 159 ! optional 160 real(dp),optional,intent(in) :: rhor(nfft,nspden) 161 real(dp),optional,intent(in) :: vxc(cplex*nfft,nspden) 162 163 !Local variables------------------------------- 164 !scalars 165 integer,parameter :: formeig1=1 166 integer :: ban2tot,bantot,bdtot_index,ddkcase,iband,icg,icg1,idir1 167 integer :: ierr,ifft,ii,ikg,ikg1,ikpt,ilm,ipert1,ispden,isppol 168 integer :: istwf_k,isym,jj,master,me,n1,n2,n3,n3xccc,n4,n5,n6 169 integer :: nband_k,nfftot,npw1_k,npw_k,nspinor_,option,spaceworld,optnc 170 real(dp) :: doti,dotr,wtk_k 171 logical :: t_exist 172 character(len=500) :: msg 173 character(len=fnlen) :: fiwfddk 174 type(gs_hamiltonian_type) :: gs_hamkq 175 !arrays 176 integer :: ddkfil(3) 177 integer,allocatable :: kg1_k(:,:),kg_k(:,:),symrl1(:,:,:) 178 real(dp) :: d2nl_elfd(2,3),d2nl_mgfd(2,3),kpoint(3),kpq(3),sumelfd(2),summgfd(2),tsec(2) 179 real(dp),allocatable :: buffer1(:),buffer2(:),d2bbb_k(:,:,:,:),d2nl_k(:,:,:) 180 real(dp),allocatable :: eig1_k(:),eig_k(:),occ_k(:) 181 real(dp) :: rhodummy(0,0) 182 real(dp),allocatable :: vpsp1(:),vxc1(:,:),work1(:,:,:),xccc3d1(:),ylm1_k(:,:),ylm_k(:,:) 183 type(pawtab_type) :: pawtab(dtset%ntypat*psps%usepaw) 184 type(wfk_t) :: ddks(3) 185 186 ! ********************************************************************* 187 188 ABI_UNUSED(nkpt) 189 190 DBG_ENTER("COLL") 191 192 !Not valid for PAW 193 if (psps%usepaw==1) then 194 msg='This routine cannot be used for PAW (use pawnst3 instead) !' 195 MSG_BUG(msg) 196 end if 197 198 !Keep track of total time spent in dfpt_nstdy 199 call timab(101,1,tsec) 200 201 !Init parallelism 202 spaceworld=mpi_enreg%comm_cell 203 me=mpi_enreg%me_kpt 204 205 master =0 206 207 !Zero only portion of nonlocal matrix to be computed here 208 d2nl(:,:,1:dtset%natom+2,idir,ipert)=zero 209 210 ABI_ALLOCATE(d2bbb_k,(2,3,dtset%mband,dtset%mband*dtset%prtbbb)) 211 ABI_ALLOCATE(d2nl_k,(2,3,mpert)) 212 ABI_ALLOCATE(eig_k,(nsppol*dtset%mband)) 213 ABI_ALLOCATE(eig1_k,(2*nsppol*dtset%mband**2)) 214 ABI_ALLOCATE(kg_k,(3,mpw)) 215 ABI_ALLOCATE(kg1_k,(3,mpw1)) 216 217 !Do not try to open electric field file 218 ddkfil(:)=0 219 !The treatment of homogeneous electric field potential need the existence of d/dk files. 220 do idir1=1,3 221 ddkcase=idir1+dtset%natom*3 222 call appdig(ddkcase,dtfil%fnamewffddk,fiwfddk) 223 224 ! Check that ddk file exists 225 t_exist = file_exists(fiwfddk) 226 if (.not. t_exist) then 227 ! Try netcdf file. 228 t_exist = file_exists(nctk_ncify(fiwfddk)) 229 if (t_exist) then 230 fiwfddk = nctk_ncify(fiwfddk) 231 write(msg,"(3a)")"- File: ",trim(fiwfddk)," does not exist but found netcdf file with similar name." 232 call wrtout(std_out,msg,'COLL') 233 end if 234 end if 235 236 if (t_exist) then 237 ! Note the use of unit numbers 21, 22 and 23 238 ddkfil(idir1)=20+idir1 239 write(msg, '(a,a)') '-open ddk wf file :',trim(fiwfddk) 240 call wrtout(std_out,msg,'COLL') 241 call wrtout(ab_out,msg,'COLL') 242 call wfk_open_read(ddks(idir1),fiwfddk,formeig1,dtset%iomode,ddkfil(idir1),spaceworld) 243 end if 244 end do 245 246 !Update list of computed matrix elements 247 if (ipert /= dtset%natom + 1) then 248 do ipert1=1,mpert 249 do idir1=1,3 250 if(ipert1 <= dtset%natom .or. ipert1==dtset%natom+2 & 251 & .and. ddkfil(idir1)/=0) then 252 blkflg(idir1,ipert1,idir,ipert)=1 253 end if 254 end do 255 end do 256 else 257 ipert1 = dtset%natom + 1 258 do idir1=1,3 259 ! If was already computed in another run or dataset, or if is to be computed in the present one 260 if ((ddkfil(idir1) /= 0).or. (dtset%rfdir(idir1)/=0.and. idir1<=idir) ) then 261 ! if ((ddkfil(idir1) /= 0).or. (idir1==idir) ) then 262 blkflg(idir1,ipert1,idir,ipert)=1 263 end if 264 end do 265 end if 266 267 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 268 n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6) 269 nspinor_=dtset%nspinor 270 271 bantot = 0 272 ban2tot = 0 273 274 !==== Initialize most of the Hamiltonian ==== 275 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 276 !2) Perform the setup needed for the non-local factors: 277 !3) Constant kleimann-Bylander energies are copied from psps to gs_hamk. 278 call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,dtset%natom,& 279 & dtset%typat,xred,nfft,dtset%mgfft,ngfft,rprimd,dtset%nloalg,ph1d=ph1d,& 280 & use_gpu_cuda=dtset%use_gpu_cuda) 281 282 !LOOP OVER SPINS 283 bdtot_index=0 284 icg=0;icg1=0 285 do isppol=1,nsppol 286 287 ikg=0;ikg1=0 288 289 ! Continue to initialize the Hamiltonian 290 call load_spin_hamiltonian(gs_hamkq,isppol,with_nonlocal=.true.) 291 292 ! BIG FAT k POINT LOOP 293 do ikpt=1,nkpt_rbz 294 295 nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz) 296 istwf_k=istwfk_rbz(ikpt) 297 npw_k=npwarr(ikpt) 298 npw1_k=npwar1(ikpt) 299 300 eig_k(1:nband_k) = eigen0(1+bantot:nband_k+bantot) 301 eig1_k(1:2*nband_k**2) = eigen1(1+ban2tot:2*nband_k**2+ban2tot) 302 bantot = bantot + nband_k 303 ban2tot = ban2tot + 2*nband_k**2 304 305 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 306 bdtot_index=bdtot_index+nband_k 307 ! The wavefunction blocks for ddk file is skipped elsewhere in the loop 308 ! Skip the rest of the k-point loop 309 cycle 310 end if 311 312 ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 313 ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)) 314 315 ! In case of electric field pert1, read ddk wfs file 316 ! Note that the symmetries are not used for ddk, so read each k point 317 ! Also take into account implicitely the parallelism over k points 318 319 do idir1=1,3 320 if (ddkfil(idir1)/=0) then 321 ii = wfk_findk(ddks(idir1), kpt_rbz(:, ikpt)) 322 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1") 323 end if 324 end do 325 326 ABI_ALLOCATE(occ_k,(nband_k)) 327 occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index) 328 kpoint(:)=kpt_rbz(:,ikpt) 329 kpq(:)=kpoint(:)+dtset%qptn(:) 330 wtk_k=wtk_rbz(ikpt) 331 d2nl_k(:,:,:)=zero 332 if(dtset%prtbbb==1)d2bbb_k(:,:,:,:)=zero 333 334 ! Get plane-wave vectors and related data at k 335 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 336 if (psps%useylm==1) then 337 do ilm=1,psps%mpsang*psps%mpsang 338 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 339 end do 340 end if 341 342 ! Get plane-wave vectors and related data at k+q 343 kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1) 344 if (psps%useylm==1) then 345 do ilm=1,psps%mpsang*psps%mpsang 346 ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm) 347 end do 348 end if 349 350 ! Compute the eigenvalues, wavefunction, 351 ! contributions to kinetic energy, nonlocal energy, forces, 352 ! and update of rhor1 to this k-point and this spin polarization. 353 ! Note that dfpt_nstwf is called with kpoint, while kpt is used inside dfpt_vtowfk 354 call dfpt_nstwf(cg,cg1,ddkfil,dtset,d2bbb_k,d2nl_k,eig_k,eig1_k,gs_hamkq,& 355 & icg,icg1,idir,ikpt,ipert,isppol,istwf_k,kg_k,kg1_k,kpoint,kpq,mkmem,mk1mem,mpert,& 356 & mpi_enreg,mpw,mpw1,nband_k,npw_k,npw1_k,nsppol,& 357 & occ_k,psps,rmet,ddks,wtk_k,ylm_k,ylm1_k) 358 359 d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:) 360 if(dtset%prtbbb==1)d2bbb(:,:,idir,ipert,:,:)=d2bbb(:,:,idir,ipert,:,:)+d2bbb_k(:,:,:,:) 361 362 ! Keep track of total number of bands 363 bdtot_index=bdtot_index+nband_k 364 365 ! Shift arrays memory 366 if (mkmem/=0) then 367 icg=icg+npw_k*dtset%nspinor*nband_k 368 ikg=ikg+npw_k 369 end if 370 if (mk1mem/=0) then 371 icg1=icg1+npw1_k*dtset%nspinor*nband_k 372 ikg1=ikg1+npw1_k 373 end if 374 375 ABI_DEALLOCATE(occ_k) 376 ABI_DEALLOCATE(ylm_k) 377 ABI_DEALLOCATE(ylm1_k) 378 379 ! End big k point loop 380 end do 381 382 ! End loop over spins 383 end do 384 385 !if(xmpi_paral==1)then 386 !call timab(161,1,tsec) 387 !call wrtout(std_out,' dfpt_nstdy: loop on k-points and spins done in parallel','COLL') 388 !call xmpi_barrier(spaceworld) 389 !call timab(161,2,tsec) 390 !end if 391 392 call destroy_hamiltonian(gs_hamkq) 393 394 !Treat fixed occupation numbers (as in vtorho) 395 if(xmpi_paral==1)then 396 ABI_ALLOCATE(buffer1,(2*3*mpert)) 397 ABI_ALLOCATE(buffer2,(2*3*mpert)) 398 ! Pack d2nl 399 buffer1(1:2*3*mpert)=reshape(d2nl(:,:,:,idir,ipert),(/2*3*mpert/)) 400 ! Build sum of everything 401 call timab(48,1,tsec) 402 call xmpi_sum(buffer1,buffer2,2*3*mpert,spaceworld,ierr) 403 call timab(48,2,tsec) 404 ! Unpack the final result 405 d2nl(:,:,:,idir,ipert)=reshape(buffer2(:),(/2,3,mpert/)) 406 ABI_DEALLOCATE(buffer1) 407 ABI_DEALLOCATE(buffer2) 408 if(dtset%prtbbb==1)then 409 ABI_ALLOCATE(buffer1,(2*3*dtset%mband*dtset%mband)) 410 ABI_ALLOCATE(buffer2,(2*3*dtset%mband*dtset%mband)) 411 ! Pack d2bbb 412 buffer1(1:2*3*dtset%mband*dtset%mband)=reshape(d2bbb(:,:,idir,ipert,:,:),(/2*3*dtset%mband*dtset%mband/)) 413 ! Build sum of everything 414 call timab(48,1,tsec) 415 call xmpi_sum(buffer1,buffer2,2*3*dtset%mband*dtset%mband,spaceworld,ierr) 416 call timab(48,2,tsec) 417 ! Unpack the final result 418 d2bbb(:,:,idir,ipert,:,:)=reshape(buffer2(:),(/2,3,dtset%mband,dtset%mband/)) 419 ABI_DEALLOCATE(buffer1) 420 ABI_DEALLOCATE(buffer2) 421 end if 422 end if ! xmpi_paral==1 423 424 !In the case of the strain perturbation time-reversal symmetry will always 425 !be true so imaginary part of d2nl will be must be set to zero here since 426 !the symmetry-reduced kpt set will leave a non-zero imaginary part. 427 if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) d2nl(2,:,:,idir,ipert)=zero 428 429 !In case of electric field ipert1, close the ddk wf files 430 do idir1=1,3 431 if (ddkfil(idir1)/=0)then 432 call wfk_close(ddks(idir1)) 433 end if 434 end do 435 436 !Symmetrize the non-local contributions, 437 !as was needed for the forces in a ground-state calculation 438 !However, here the quantity is complex, and there are phases ! 439 440 !Do the transform 441 ABI_ALLOCATE(work1,(2,3,dtset%natom)) 442 do ipert1=1,dtset%natom 443 do idir1=1,3 444 work1(1,idir1,ipert1)=d2nl(1,idir1,ipert1,idir,ipert) 445 work1(2,idir1,ipert1)=d2nl(2,idir1,ipert1,idir,ipert) 446 end do 447 end do 448 call dfpt_sygra(dtset%natom,d2nl(:,:,:,idir,ipert),work1,indsy1,ipert,nsym1,dtset%qptn,symrc1) 449 ABI_DEALLOCATE(work1) 450 451 !Must also symmetrize the electric/magnetic field perturbation response ! 452 !(XG 000803 This was not implemented until now) 453 if(sum(ddkfil(:))/=0)then 454 ! Get the symmetry matrices in terms of real space basis 455 ABI_ALLOCATE(symrl1,(3,3,nsym1)) 456 do isym=1,nsym1 457 call mati3inv(symrc1(:,:,isym),symrl1(:,:,isym)) 458 end do 459 ! There should not be any imaginary part, but stay general (for debugging) 460 d2nl_elfd(:,:)=d2nl(:,:,dtset%natom+2,idir,ipert) 461 do ii=1,3 462 sumelfd(:)=zero 463 summgfd(:)=zero 464 do isym=1,nsym1 465 do jj=1,3 466 if(symrl1(ii,jj,isym)/=0)then 467 if(ddkfil(jj)==0)then 468 blkflg(ii,dtset%natom+2,idir,ipert)=0 469 end if 470 end if 471 end do 472 sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+& 473 & dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+& 474 & dble(symrl1(ii,3,isym))*d2nl_elfd(:,3) 475 summgfd(:)=summgfd(:)+dble(symrl1(ii,1,isym))*d2nl_mgfd(:,1)+& 476 & dble(symrl1(ii,2,isym))*d2nl_mgfd(:,2)+& 477 & dble(symrl1(ii,3,isym))*d2nl_mgfd(:,3) 478 end do 479 d2nl(:,ii,dtset%natom+2,idir,ipert)=sumelfd(:)/dble(nsym1) 480 end do 481 482 if ((dtset%prtbbb==1).and.(ipert<=dtset%natom)) then 483 do iband = 1,dtset%mband 484 d2nl_elfd(:,:)=d2bbb(:,:,idir,ipert,iband,iband) 485 do ii=1,3 486 sumelfd(:)=zero 487 do isym=1,nsym1 488 sumelfd(:)=sumelfd(:)+dble(symrl1(ii,1,isym))*d2nl_elfd(:,1)+& 489 & dble(symrl1(ii,2,isym))*d2nl_elfd(:,2)+& 490 & dble(symrl1(ii,3,isym))*d2nl_elfd(:,3) 491 end do 492 d2bbb(:,ii,idir,ipert,iband,iband)=sumelfd(:)/dble(nsym1) 493 end do 494 end do !iband 495 end if 496 497 ABI_DEALLOCATE(symrl1) 498 end if 499 500 !---------------------------------------------------------------------------- 501 !Now, treat the local contribution 502 503 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 504 ABI_ALLOCATE(vpsp1,(cplex*nfft)) 505 if (ipert /= dtset%natom + 1) then 506 n3xccc=0;if(psps%n1xccc/=0) n3xccc=nfft 507 ABI_ALLOCATE(xccc3d1,(cplex*n3xccc)) 508 ABI_ALLOCATE(vxc1,(cplex*nfft,nspden)) 509 510 do ipert1=1,mpert 511 do idir1=1,3 512 if(ipert1 <= dtset%natom)then 513 514 ! Get first-order local potential and first-order pseudo core density 515 call dfpt_vlocal(atindx,cplex,gmet,gsqcut,idir1,ipert1,mpi_enreg,psps%mqgrid_ff,dtset%natom,& 516 & nattyp,nfft,ngfft,dtset%ntypat,n1,n2,n3,dtset%paral_kgb,ph1d,psps%qgrid_ff,& 517 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 518 if(psps%n1xccc/=0)then 519 call dfpt_mkcore(cplex,idir1,ipert1,dtset%natom,dtset%ntypat,n1,psps%n1xccc,& 520 & n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred) 521 end if 522 523 ! Get first-order exchange-correlation potential (core-correction contribution only !) 524 if(psps%n1xccc/=0)then 525 option=0 526 !FR SPr EB non-collinear magnetism 527 if (nspden==4.and.present(rhor).and.present(vxc)) then 528 optnc=1 529 call dfpt_mkvxc_noncoll(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,rhodummy,0,& 530 & nkxc,nspden,n3xccc,optnc,option,dtset%paral_kgb,dtset%qptn,rhor,rhor1,& 531 & rprimd,0,vxc,vxc1,xccc3d1) 532 else 533 call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,ngfft,rhodummy,0,rhodummy,0,& 534 & nkxc,nspden,n3xccc,option,dtset%paral_kgb,dtset%qptn,rhodummy,& 535 & rprimd,0,vxc1,xccc3d1) 536 end if 537 else 538 vxc1(:,:)=zero 539 end if 540 541 ! Norm-conserving pseudpopotential case: 542 ! Combines density j2 with local potential j1 (vpsp1 and vxc1) 543 ! XG030514 : this is a first possible coding, however, each dotprod contains 544 ! a parallel section (reduction), so it is better to use only one dotprod ... 545 ! call dotprod_vn(cplex,rhor1,dr_psp1,di_psp1,mpi_enreg,nfft,nfftot,1,2,vpsp1,ucvol) 546 ! call dotprod_vn(cplex,rhor1,dr_xc1,di_xc1,mpi_enreg,nfft,nfftot,nspden,2,vxc1,ucvol) 547 ! dotr=dr_psp1+dr_xc1;doti=di_psp1+di_xc1... but then, one needs to overload vxc1 548 do ispden=1,min(nspden,2) 549 do ifft=1,cplex*nfft 550 vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft) 551 end do 552 end do 553 call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol) 554 555 ! MVeithen 021212 : in case ipert = 2, these lines compute the local part 556 ! of the Born effective charges from phonon and electric 557 ! field type perturbations, see eq. 43 of 558 ! X. Gonze and C. Lee, PRB 55, 10355 (1997) 559 ! The minus sign is due to the fact that the effective charges 560 ! are minus the second derivatives of the energy 561 if (ipert == dtset%natom+2) then 562 d2lo(1,idir1,ipert1,idir,ipert)=-dotr 563 d2lo(2,idir1,ipert1,idir,ipert)=-doti 564 else 565 d2lo(1,idir1,ipert1,idir,ipert)=dotr 566 d2lo(2,idir1,ipert1,idir,ipert)=doti 567 end if 568 ! Endif ipert1<=natom 569 end if 570 end do 571 end do 572 573 ABI_DEALLOCATE(vxc1) 574 ABI_DEALLOCATE(xccc3d1) 575 576 end if ! ipert /= natom +1 577 578 ABI_DEALLOCATE(d2bbb_k) 579 ABI_DEALLOCATE(d2nl_k) 580 ABI_DEALLOCATE(kg_k) 581 ABI_DEALLOCATE(kg1_k) 582 ABI_DEALLOCATE(vpsp1) 583 ABI_DEALLOCATE(eig_k) 584 ABI_DEALLOCATE(eig1_k) 585 586 call timab(101,2,tsec) 587 588 DBG_EXIT("COLL") 589 590 end subroutine dfpt_nstdy