TABLE OF CONTENTS
ABINIT/longwave [ Functions ]
NAME
longwave
FUNCTION
Primary routine for conducting DFPT calculations of spatial dispersion properties
INPUTS
codvsn = code version dtfil <type(datafiles_type)> = variables related to files dtset <type(dataset_type)> = all input variables for this dataset etotal = new total energy (no meaning at output) mpi_enreg=information about MPI pnarallelization occ(mband*nkpt*nsppol) = occupation number for each band and k xred(3,natom) = reduced atomic coordinates
OUTPUT
npwtot(nkpt) = total number of plane waves at each k point
SIDE EFFECTS
pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)> = variables related to pseudopotentials
NOTES
PARENTS
m_driver
CHILDREN
check_kxc,ddb_hdr%free,ddb_hdr%open_write,ddb_hdr_init, dfpt_lw_doutput,ebands_free fourdp,hdr%free,hdr%update,hdr_init,inwffil,kpgio,matr3inv,mkcore,mkrho pawfgr_init,pspini,read_rhor,rhotoxc,setsym,setup1,symmetrize_xred wffclose,xcdata_init
SOURCE
116 subroutine longwave(codvsn,dtfil,dtset,etotal,mpi_enreg,npwtot,occ,& 117 & pawrad,pawtab,psps,xred) 118 119 #ifdef FC_INTEL 120 !DEC$ NOOPTIMIZE 121 #endif 122 123 124 implicit none 125 126 !Arguments ------------------------------------ 127 !scalars 128 real(dp),intent(inout) :: etotal 129 character(len=8),intent(in) :: codvsn 130 type(MPI_type),intent(inout) :: mpi_enreg 131 type(datafiles_type),intent(in) :: dtfil 132 type(dataset_type),intent(inout) :: dtset 133 type(pseudopotential_type),intent(inout) :: psps 134 !arrays 135 integer,intent(out) :: npwtot(dtset%nkpt) 136 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),xred(3,dtset%natom) 137 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 138 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 139 140 !Local variables------------------------------- 141 !scalars 142 integer,parameter :: cplex1=1,formeig=0,response=1 143 integer :: ask_accurate,bantot,coredens_method,dimffnl,dimffnl_i 144 integer :: gscase,iatom,ierr,indx,ireadwf0,iscf_eff,itypat 145 integer :: ider,idir0,idir 146 integer :: i1dir,i1pert,i2dir,ii,i2pert,i3dir,i3pert 147 integer :: mcg,mgfftf,natom,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim 148 ! integer :: isym 149 integer :: mpert,my_natom,n1,nkxc,nk3xc,ntypat,n3xccc,nylmgr 150 integer :: option,optorth,psp_gencond,rdwrpaw,spaceworld,timrev,tim_mkrho 151 integer :: usexcnhat,useylmgr 152 real(dp) :: ecore,ecutdg_eff,ecut_eff,enxc,etot,fermie,fermih,gsqcut_eff,gsqcutc_eff,residm 153 real(dp) :: ucvol,vxcavg 154 logical :: non_magnetic_xc 155 character(len=500) :: msg 156 type(ebands_t) :: bstruct 157 type(ddb_hdr_type) :: ddb_hdr 158 type(ddb_type) :: ddb 159 type(paw_dmft_type) :: paw_dmft 160 type(pawfgr_type) :: pawfgr 161 type(hdr_type) :: hdr,hdr_den 162 type(xcdata_type) :: xcdata 163 type(wvl_data) :: wvl 164 type(wffile_type) :: wffgs,wfftgs 165 !arrays 166 integer :: ngfft(18),ngfftf(18),perm(6) 167 real(dp) :: dummy6(6),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3) 168 real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3) 169 real(dp) :: strsxc(6) 170 integer,allocatable :: atindx(:),atindx1(:) 171 integer,allocatable :: blkflg(:,:,:,:,:,:),blkflg_car(:,:,:,:,:,:) 172 integer,allocatable :: d3e_pert1(:),d3e_pert2(:),d3e_pert3(:) 173 integer,allocatable :: indsym(:,:,:),irrzon(:,:,:),kg(:,:) 174 integer,allocatable :: nattyp(:),npwarr(:),pertsy(:,:),symrec(:,:,:) 175 integer,allocatable :: rfpert(:,:,:,:,:,:) 176 real(dp),allocatable :: cg(:,:) 177 real(dp),allocatable :: d3etot(:,:,:,:,:,:,:),d3etot_car(:,:,:,:,:,:,:) 178 real(dp),allocatable :: d3etot_nv(:,:,:,:,:,:,:),doccde(:) 179 real(dp),allocatable :: eigen0(:),ffnl(:,:,:,:,:),ffnl_i(:,:,:,:,:) 180 real(dp),allocatable :: grxc(:,:),kxc(:,:),vxc(:,:),nhat(:,:),nhatgr(:,:,:) 181 real(dp),allocatable :: phnons(:,:,:),rhog(:,:),rhor(:,:),dummy_dyfrx2(:,:,:) 182 ! real(dp),allocatable :: symrel_cart(:,:,:) 183 real(dp),allocatable :: work(:),xccc3d(:) 184 real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:) 185 type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:) 186 ! ************************************************************************* 187 188 DBG_ENTER("COLL") 189 190 !Not valid for PAW 191 if (psps%usepaw==1) then 192 msg='This routine cannot be used for PAW!' 193 ABI_BUG(msg) 194 end if 195 196 !Not valid for finite wave-vector perturbations 197 if (sqrt(sum(dtset%qptn**2))/=0_dp) then 198 msg='This routine cannot be used for q /= 0 ' 199 ABI_BUG(msg) 200 end if 201 202 !Only usable with spherical harmonics 203 if (dtset%useylm/=1.and.(dtset%lw_qdrpl/=0.or.dtset%lw_flexo/=0)) then 204 msg='This routine can only be used with useylm/=1 for lw_natopt=1' 205 ABI_BUG(msg) 206 end if 207 208 !Not valid for spin-dependent calculations 209 if (dtset%nspinor/=1.or.dtset%nsppol/=1.or.dtset%nspden/=1) then 210 msg='This routine cannot be used for spin-dependent calculations' 211 ABI_BUG(msg) 212 end if 213 214 !Not usable with core electron density corrections 215 if (psps%n1xccc/=0) then 216 msg='This routine cannot be used for n1xccc/=0' 217 ABI_BUG(msg) 218 end if 219 220 221 !Define some data 222 ntypat=psps%ntypat 223 natom=dtset%natom 224 timrev=1 225 226 !Init spaceworld 227 spaceworld=mpi_enreg%comm_cell 228 my_natom=mpi_enreg%my_natom 229 230 !Define FFT grid(s) sizes (be careful !) 231 !See NOTES in the comments at the beginning of this file. 232 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 233 nfftot=product(ngfft(1:3)) 234 nfftotf=product(ngfftf(1:3)) 235 236 !Set up for iterations 237 call setup1(dtset%acell_orig(1:3,1),bantot,dtset,& 238 & ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,& 239 & ngfftf,ngfft,dtset%nkpt,dtset%nsppol,& 240 & response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw) 241 242 !Define the set of admitted perturbations taking into account 243 !the possible permutations 244 ! -> natom+8 refers to ddq perturbation 245 mpert=natom+8 246 ABI_MALLOC(blkflg,(3,mpert,3,mpert,3,mpert)) 247 ABI_MALLOC(d3etot,(2,3,mpert,3,mpert,3,mpert)) 248 ABI_MALLOC(d3etot_nv,(2,3,mpert,3,mpert,3,mpert)) 249 ABI_MALLOC(rfpert,(3,mpert,3,mpert,3,mpert)) 250 ABI_MALLOC(d3e_pert1,(mpert)) 251 ABI_MALLOC(d3e_pert2,(mpert)) 252 ABI_MALLOC(d3e_pert3,(mpert)) 253 blkflg(:,:,:,:,:,:) = 0 254 d3etot(:,:,:,:,:,:,:) = zero 255 d3etot_nv(:,:,:,:,:,:,:) = zero 256 rfpert(:,:,:,:,:,:) = 0 257 d3e_pert1(:) = 0 ; d3e_pert2(:) = 0 ; d3e_pert3(:) = 0 258 259 d3e_pert3(natom+8)=1 260 if (dtset%lw_qdrpl==1) then 261 d3e_pert1(natom+2)=1 262 d3e_pert2(1:natom)=1 263 end if 264 265 if (dtset%lw_flexo==2.or.dtset%lw_flexo==1) then 266 d3e_pert1(natom+2)=1 267 d3e_pert2(natom+3:natom+4)=1 268 end if 269 270 if (dtset%lw_flexo==3.or.dtset%lw_flexo==1) then 271 d3e_pert1(natom+2)=1 ; d3e_pert1(1:natom)=1 272 d3e_pert2(1:natom)=1 273 end if 274 275 if (dtset%lw_flexo==4.or.dtset%lw_flexo==1) then 276 d3e_pert1(1:natom)=1 277 d3e_pert2(natom+3:natom+4)=1 278 end if 279 280 if (dtset%lw_natopt==1) then 281 d3e_pert1(natom+2)=1 282 d3e_pert2(natom+2)=1 283 end if 284 285 perm(:)=0 286 do i1pert = 1, mpert 287 do i2pert = 1, mpert 288 do i3pert = 1, mpert 289 perm(1)=d3e_pert1(i1pert)*d3e_pert2(i2pert)*d3e_pert3(i3pert) 290 ! perm(2)=d3e_pert1(i1pert)*d3e_pert2(i3pert)*d3e_pert3(i2pert) 291 ! perm(3)=d3e_pert1(i2pert)*d3e_pert2(i1pert)*d3e_pert3(i3pert) 292 ! perm(4)=d3e_pert1(i2pert)*d3e_pert2(i3pert)*d3e_pert3(i1pert) 293 ! perm(5)=d3e_pert1(i3pert)*d3e_pert2(i2pert)*d3e_pert3(i1pert) 294 ! perm(6)=d3e_pert1(i3pert)*d3e_pert2(i1pert)*d3e_pert3(i2pert) 295 if ( sum(perm(:)) > 0 ) rfpert(:,i1pert,:,i2pert,:,i3pert)=1 296 end do 297 end do 298 end do 299 300 !Do symmetry stuff 301 ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 302 ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 303 ABI_MALLOC(indsym,(4,dtset%nsym,natom)) 304 ABI_MALLOC(symrec,(3,3,dtset%nsym)) 305 irrzon=0;indsym=0;symrec=0;phnons=zero 306 !If the density is to be computed by mkrho, need irrzon and phnons 307 iscf_eff=0;if(dtset%getden==0)iscf_eff=1 308 call setsym(indsym,irrzon,iscf_eff,natom,& 309 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 310 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred) 311 312 !Symmetrize atomic coordinates over space group elements: 313 call symmetrize_xred(natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym) 314 315 ! Get symmetries in cartesian coordinates 316 ! ABI_MALLOC(symrel_cart, (3, 3, dtset%nsym)) 317 ! do isym =1,dtset%nsym 318 ! call symredcart(rprimd, gprimd, symrel_cart(:,:,isym), dtset%symrel(:,:,isym)) 319 ! ! purify operations in cartesian coordinates. 320 ! where (abs(symrel_cart(:,:,isym)) < tol14) 321 ! symrel_cart(:,:,isym) = zero 322 ! end where 323 ! end do 324 325 ! call sylwtens(indsym,mpert,natom,dtset%nsym,rfpert,symrec,dtset%symrel,symrel_cart) 326 call sylwtens(indsym,mpert,natom,dtset%nsym,rfpert,symrec,dtset%symrel) 327 328 write(msg,'(a,a,a,a,a)') ch10, & 329 & ' The list of irreducible elements of the spatial-dispersion third-order energy derivatives is: ', ch10,& 330 & ' (in reduced coordinates except for strain pert.) ', ch10 331 call wrtout(ab_out,msg,'COLL') 332 call wrtout(std_out,msg,'COLL') 333 334 write(msg,'(12x,a)') 'i1dir i1pert i2dir i2pert i3dir i3pert' 335 call wrtout(ab_out,msg,'COLL') 336 call wrtout(std_out,msg,'COLL') 337 n1 = 0 338 do i3pert = 1, mpert 339 do i3dir = 1, 3 340 do i2pert = 1, mpert 341 do i2dir = 1,3 342 do i1pert = 1, mpert 343 do i1dir = 1, 3 344 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 345 n1 = n1 + 1 346 write(msg,'(2x,i4,a,6(5x,i3))') n1,')', & 347 & i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 348 call wrtout(ab_out,msg,'COLL') 349 call wrtout(std_out,msg,'COLL') 350 else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) then 351 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 352 if (dtset%prtvol>=10) then 353 n1 = n1 + 1 354 write(msg,'(2x,i4,a,6(5x,i3),a)') n1,')', & 355 & i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,' => must be zero, not computed' 356 call wrtout(ab_out,msg,'COLL') 357 call wrtout(std_out,msg,'COLL') 358 end if 359 else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-1) then 360 if (dtset%prtvol>=10) then 361 n1 = n1 + 1 362 write(msg,'(2x,i4,a,6(5x,i3),a)') n1,')', & 363 & i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,' => symmetric of another element, not computed' 364 call wrtout(ab_out,msg,'COLL') 365 call wrtout(std_out,msg,'COLL') 366 end if 367 end if 368 end do 369 end do 370 end do 371 end do 372 end do 373 end do 374 write(msg,'(a,a)') ch10,ch10 375 call wrtout(ab_out,msg,'COLL') 376 call wrtout(std_out,msg,'COLL') 377 378 !In some cases (e.g. getcell/=0), the plane wave vectors have 379 !to be generated from the original simulation cell 380 rprimd_for_kg=rprimd 381 if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=dtset%rprimd_orig(:,:,1) 382 call matr3inv(rprimd_for_kg,gprimd_for_kg) 383 gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg) 384 385 !Set up the basis sphere of planewaves 386 ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem)) 387 ABI_MALLOC(npwarr,(dtset%nkpt)) 388 call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg,& 389 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,& 390 & dtset%nsppol) 391 392 !Open and read pseudopotential files 393 ecore=zero 394 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,& 395 & psps,rprimd,comm_mpi=mpi_enreg%comm_cell) 396 397 !Initialize band structure datatype 398 bstruct = ebands_from_dtset(dtset, npwarr) 399 400 !Initialize PAW atomic occupancies to zero 401 ABI_MALLOC(pawrhoij,(0)) 402 403 !Initialize header 404 gscase=0 405 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr, & 406 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 407 408 !Update header, with evolving variables, when available 409 !Here, rprimd, xred and occ are available 410 etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm 411 412 !If parallelism over atom, hdr is distributed 413 call hdr%update(bantot,etot,fermie,fermih,& 414 & residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), & 415 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 416 417 !Clean band structure datatype (should use it more in the future !) 418 call ebands_free(bstruct) 419 420 !Initialize wavefunction files and wavefunctions. 421 ireadwf0=1 422 423 mcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 424 ABI_STAT_MALLOC(cg,(2,mcg), ierr) 425 ABI_CHECK(ierr==0, "out-of-memory in cg") 426 427 ABI_MALLOC(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol)) 428 eigen0(:)=zero ; ask_accurate=1 429 optorth=0 430 431 hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors 432 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,& 433 & formeig,hdr,ireadwf0,dtset%istwfk,kg,dtset%kptns,& 434 & dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,& 435 & dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,& 436 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 437 & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,dtfil%fnamewffk,wvl) 438 hdr%rprimd=rprimd 439 440 !Close wffgs, if it was ever opened (in inwffil) 441 if (ireadwf0==1) then 442 call WffClose(wffgs,ierr) 443 end if 444 445 446 !Generate an index table of atoms, in order for them to be used 447 !type after type. 448 ABI_MALLOC(atindx,(natom)) 449 ABI_MALLOC(atindx1,(natom)) 450 ABI_MALLOC(nattyp,(ntypat)) 451 indx=1 452 do itypat=1,ntypat 453 nattyp(itypat)=0 454 do iatom=1,natom 455 if(dtset%typat(iatom)==itypat)then 456 atindx(iatom)=indx 457 atindx1(indx)=iatom 458 indx=indx+1 459 nattyp(itypat)=nattyp(itypat)+1 460 end if 461 end do 462 end do 463 464 !Derivative of occupations is always zero for non metallic systems 465 ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol)) 466 doccde(:)=zero 467 468 !Read ground-state charge density from diskfile in case getden /= 0 469 !or compute it from wfs that were read previously : rhor 470 471 ABI_MALLOC(rhog,(2,nfftf)) 472 ABI_MALLOC(rhor,(nfftf,dtset%nspden)) 473 474 if (dtset%getden /= 0 .or. dtset%irdden /= 0) then 475 ! Read rho1(r) from a disk file and broadcast data. 476 ! This part is not compatible with MPI-FFT (note single_proc=.True. below) 477 478 rdwrpaw=psps%usepaw 479 ABI_MALLOC(pawrhoij_read,(0)) 480 481 ! MT july 2013: Should we read rhoij from the density file ? 482 call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor, & 483 hdr_den, pawrhoij_read, spaceworld, check_hdr=hdr) 484 etotal = hdr_den%etot; call hdr_den%free() 485 486 ABI_FREE(pawrhoij_read) 487 488 ! Compute up+down rho(G) by fft 489 ABI_MALLOC(work,(nfftf)) 490 work(:)=rhor(:,1) 491 call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,1,ngfftf,0) 492 ABI_FREE(work) 493 else 494 ! Obtain the charge density from read wfs 495 ! Be careful: in PAW, compensation density has to be added ! 496 tim_mkrho=4 497 paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented 498 paw_dmft%use_dmft=0 ! respfn with dmft not implemented 499 500 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 501 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 502 end if ! getden 503 ! ABI_FREE(cg) 504 505 !Pseudo core electron density by method 2 506 !TODO: The code is not still adapted to consider n3xccc in the long-wave 507 !driver. 508 n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf 509 ABI_MALLOC(xccc3d,(n3xccc)) 510 coredens_method=2 511 if (coredens_method==2.and.psps%n1xccc/=0) then 512 option=1 513 ABI_MALLOC(dummy_dyfrx2,(3,3,natom)) ! dummy 514 ABI_MALLOC(vxc,(0,0)) ! dummy 515 ABI_MALLOC(grxc,(3,natom)) 516 call mkcore(dummy6,dummy_dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,& 517 & ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,& 518 & psps%xcccrc,psps%xccc1d,xccc3d,xred) 519 ABI_FREE(dummy_dyfrx2) ! dummy 520 ABI_FREE(vxc) ! dummy 521 ABI_FREE(grxc) ! dummy 522 end if 523 524 !Set up xc potential. Compute kxc here. 525 !TODO: Iclude nonlinear core corrections (see m_respfn_driver.F90) 526 option=2 ; nk3xc=1 527 nkxc=2*min(dtset%nspden,2)-1;if(dtset%xclevel==2)nkxc=12*min(dtset%nspden,2)-5 528 call check_kxc(dtset%ixc,dtset%optdriver) 529 ABI_MALLOC(kxc,(nfftf,nkxc)) 530 ABI_MALLOC(vxc,(nfftf,dtset%nspden)) 531 532 nhatgrdim=0;nhatdim=0 533 ABI_MALLOC(nhat,(0,0)) 534 ABI_MALLOC(nhatgr,(0,0,0)) 535 ! n3xccc=0 536 ! ABI_MALLOC(xccc3d,(n3xccc)) 537 non_magnetic_xc=.false. 538 539 enxc=zero; usexcnhat=0 540 541 call xcdata_init(xcdata,dtset=dtset) 542 call rhotoxc(enxc,kxc,mpi_enreg,nfftf,ngfftf,& 543 & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor,& 544 & rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata) 545 546 !Set up the spherical harmonics (Ylm) and gradients at each k point 547 if (psps%useylm==1) then 548 useylmgr=1; option=2 ; nylmgr=9 549 ABI_MALLOC(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)) 550 ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 551 call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,& 552 & psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,option,& 553 & rprimd,ylm,ylmgr) 554 end if 555 556 !Compute nonlocal form factors ffnl1, for all atoms and all k-points. 557 if (dtset%ffnl_lw == 0) then 558 if (dtset%lw_natopt==1) then 559 ider=1;dimffnl=4;dimffnl_i=2 560 ABI_MALLOC(ffnl,(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)) 561 ABI_MALLOC(ffnl_i,(dtset%mkmem,dtset%mpw,dimffnl_i,psps%lmnmax,psps%ntypat)) 562 do idir=1, 3 563 idir0=idir 564 call preca_ffnl(dimffnl_i,ffnl_i,gmet,gprimd,ider,idir0,kg, & 565 & dtset%kptns,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw, & 566 & dtset%nkpt,npwarr,nylmgr,psps,rmet,useylmgr,ylm,ylmgr) 567 ffnl(:,:,1,:,:)=ffnl_i(:,:,1,:,:) 568 ffnl(:,:,1+idir,:,:)=ffnl_i(:,:,2,:,:) 569 end do 570 ABI_FREE(ffnl_i) 571 if (psps%useylm==1) then 572 useylmgr=0 573 ABI_FREE(ylmgr) 574 ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 575 end if 576 else 577 if (dtset%lw_qdrpl==1.or.dtset%lw_flexo==3) ider=1; idir0=4; dimffnl=4 578 if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) then 579 ider=2; idir0=4; dimffnl=10 580 end if 581 ABI_MALLOC(ffnl,(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)) 582 call preca_ffnl(dimffnl,ffnl,gmet,gprimd,ider,idir0,kg, & 583 & dtset%kptns,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw, & 584 & dtset%nkpt,npwarr,nylmgr,psps,rmet,useylmgr,ylm,ylmgr) 585 useylmgr=0 586 ABI_FREE(ylmgr) 587 ABI_MALLOC(ylmgr,(dtset%mpw*dtset%mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 588 end if 589 else if (dtset%ffnl_lw == 1) then 590 dimffnl=0 591 ABI_MALLOC(ffnl,(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)) 592 end if 593 594 !TODO: This part of the implementation does not work properly to select specific directions 595 ! for each perturbation. This development is temporarily frozen. 596 !Initialize the list of perturbations rfpert 597 ! mpert=natom+11 598 ! ABI_MALLOC(rfpert,(mpert)) 599 ! rfpert(:)=0 600 ! rfpert(natom+1)=1 601 ! if (dtset%lw_qdrpl==1.or.dtset%lw_flexo==1.or.dtset%lw_flexo==3.or.dtset%lw_flexo==4 & 602 !&.or.dtset%d3e_pert1_phon==1.or.dtset%d3e_pert2_phon==1) then 603 ! if (dtset%d3e_pert1_phon==1) rfpert(dtset%d3e_pert1_atpol(1):dtset%d3e_pert1_atpol(2))=1 604 ! if (dtset%d3e_pert2_phon==1) rfpert(dtset%d3e_pert2_atpol(1):dtset%d3e_pert2_atpol(2))=1 605 ! end if 606 ! if (dtset%lw_qdrpl==1.or.dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==3.or.& 607 !& dtset%d3e_pert1_elfd==1) then 608 ! rfpert(natom+2)=1 609 ! rfpert(natom+10)=1 610 ! rfpert(natom+11)=1 611 ! end if 612 ! if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4.or.dtset%d3e_pert2_strs/=0) then 613 ! if (dtset%d3e_pert2_strs==1.or.dtset%d3e_pert2_strs==3) rfpert(natom+3)=1 614 ! if (dtset%d3e_pert2_strs==2.or.dtset%d3e_pert2_strs==3) rfpert(natom+4)=1 615 ! endif 616 ! 617 !!Determine which directions treat for each type of perturbation 618 ! ABI_MALLOC(pertsy,(3,natom+6)) 619 ! pertsy(:,:)=0 620 ! !atomic displacement 621 ! do ipert=1,natom 622 ! if (rfpert(ipert)==1.and.dtset%d3e_pert1_phon==1) then 623 ! do idir=1,3 624 ! if (dtset%d3e_pert1_dir(idir)==1) pertsy(idir,ipert)=1 625 ! end do 626 ! endif 627 ! if (rfpert(ipert)==1.and.dtset%d3e_pert2_phon==1) then 628 ! do idir=1,3 629 ! if (dtset%d3e_pert2_dir(idir)==1) pertsy(idir,ipert)=1 630 ! end do 631 ! end if 632 ! end do 633 ! !ddk 634 ! do idir=1,3 635 ! if (dtset%d3e_pert3_dir(idir)==1) pertsy(idir,natom+1)=1 636 ! end do 637 ! !electric field 638 ! if (rfpert(natom+2)==1) then 639 ! do idir=1,3 640 ! if (dtset%d3e_pert1_dir(idir)==1) pertsy(idir,natom+2)=1 641 ! end do 642 ! end if 643 ! !strain 644 ! if (rfpert(natom+3)==1) pertsy(:,natom+3)=1 645 ! if (rfpert(natom+4)==1) pertsy(:,natom+4)=1 646 647 !TODO:Add perturbation symmetries. See m_respfn_driver.F90. 648 !........ 649 650 ! All perturbations and directions are temporarily activated 651 ABI_MALLOC(pertsy,(3,natom+6)) 652 pertsy(:,:)=1 653 654 655 !############# SPATIAL-DISPERSION PROPERTIES CALCULATION ########################### 656 657 !Anounce start of spatial-dispersion calculation 658 write(msg, '(a,80a,a,a,a)' ) ch10,('=',ii=1,80),ch10,& 659 & ' ==> Compute spatial-dispersion 3rd-order energy derivatives <== ',ch10 660 call wrtout(std_out,msg,'COLL') 661 call wrtout(ab_out,msg,'COLL') 662 663 if (dtset%prtvol>=10) then 664 write(msg,'(5a)') ' CAUTION: Individual contributions to the 3rd-order energy derivatives ',ch10, & 665 & ' are not written in a unified form. Mixed cartesian/reduced coordinates ',ch10, & 666 & ' and/or type-I/type-II forms are used.' 667 call wrtout(std_out,msg,'COLL') 668 call wrtout(ab_out,msg,'COLL') 669 end if 670 671 !Calculate the nonvariational Ewald terms 672 if (dtset%lw_flexo==1.or.dtset%lw_flexo==3.or.dtset%lw_flexo==4) then 673 call dfptlw_nv(d3etot_nv,dtset,gmet,gprimd,mpert,my_natom,rfpert,rmet,rprimd,ucvol,xred,psps%ziontypat, & 674 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 675 end if 676 677 !Main loop over the perturbations to calculate the stationary part 678 call dfptlw_loop(atindx,blkflg,cg,d3e_pert1,d3e_pert2,d3etot,dimffnl,dtfil,dtset,& 679 & ffnl,gmet,gprimd,& 680 & hdr,kg,kxc,dtset%mband,dtset%mgfft,& 681 & dtset%mkmem,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,natom,nattyp,ngfftf,nfftf,& 682 & dtset%nkpt,nkxc,dtset%nspinor,dtset%nsppol,npwarr,nylmgr,occ,& 683 & pawfgr,pawtab,& 684 & psps,rfpert,rhog,rhor,rmet,rprimd,ucvol,useylmgr,xred,ylm,ylmgr) 685 686 !Merge stationay and nonvariational contributions 687 d3etot(:,:,:,:,:,:,:)=d3etot(:,:,:,:,:,:,:) + d3etot_nv(:,:,:,:,:,:,:) 688 689 !Real (imaginary) part of d3etot is zero for first (second) momentum derivatives 690 if (dtset%kptopt /= 3) then 691 do i3pert = 1, mpert 692 do i3dir = 1, 3 693 do i2pert = 1, mpert 694 do i2dir = 1,3 695 do i1pert = 1, mpert 696 do i1dir = 1, 3 697 if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1) then 698 if (i2pert /= natom+3 .and. i2pert /= natom+4) then 699 d3etot(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero 700 else 701 d3etot(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero 702 end if 703 end if 704 end do 705 end do 706 end do 707 end do 708 end do 709 end do 710 end if 711 712 713 !Complete missing elements using symmetry operations 714 ! has_strain=.false. 715 ! if (dtset%lw_flexo==1.or.dtset%lw_flexo==2.or.dtset%lw_flexo==4) has_strain=.true. 716 ! call d3lwsym(blkflg,d3etot,has_strain,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel,symrel_cart) 717 call d3lwsym(blkflg,d3etot,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 718 719 !Deallocate global proc_distrib 720 if(xmpi_paral==1) then 721 ABI_FREE(mpi_enreg%proc_distrb) 722 end if 723 724 ! Write the DDB file 725 call ddb_hdr%init(dtset,psps,pawtab,& 726 dscrpt=' Note : temporary (transfer) database ',& 727 nblok=1,xred=xred,occ=occ) 728 729 call ddb%init(dtset, 1, mpert, with_d3E=.true.) 730 731 call ddb%set_d3matr(1, d3etot, blkflg, lw=.true.) 732 733 call ddb%write(ddb_hdr, dtfil%fnameabo_ddb) 734 735 736 !Calculate spatial-dispersion quantities in Cartesian coordinates and write 737 !them in abi_out 738 ABI_MALLOC(blkflg_car,(3,mpert,3,mpert,3,mpert)) 739 ABI_MALLOC(d3etot_car,(2,3,mpert,3,mpert,3,mpert)) 740 call lwcart(blkflg,blkflg_car,d3etot,d3etot_car,gprimd,mpert,natom,rprimd) 741 call dfptlw_out(blkflg_car,d3etot_car,dtset%lw_flexo,dtset%lw_qdrpl,dtset%lw_natopt,mpert,natom,ucvol) 742 !end if 743 744 call ddb_hdr%free() 745 call ddb%free() 746 747 !Deallocate arrays 748 ABI_FREE(atindx) 749 ABI_FREE(atindx1) 750 ABI_FREE(blkflg) 751 ABI_FREE(doccde) 752 ABI_FREE(eigen0) 753 ABI_FREE(cg) 754 ABI_SFREE(ffnl) 755 ABI_FREE(indsym) 756 ABI_FREE(irrzon) 757 ABI_FREE(nattyp) 758 ABI_FREE(kg) 759 ABI_FREE(kxc) 760 ABI_FREE(npwarr) 761 ABI_FREE(phnons) 762 ABI_FREE(rhog) 763 ABI_FREE(rhor) 764 ABI_FREE(symrec) 765 ! ABI_FREE(symrel_cart) 766 ABI_FREE(vxc) 767 ABI_FREE(d3etot) 768 ABI_FREE(d3etot_nv) 769 ABI_FREE(pertsy) 770 ABI_FREE(rfpert) 771 ABI_FREE(d3e_pert1) 772 ABI_FREE(d3e_pert2) 773 ABI_FREE(d3e_pert3) 774 ABI_SFREE(pawrhoij) 775 ABI_FREE(xccc3d) 776 ABI_SFREE(nhat) 777 ABI_SFREE(nhatgr) 778 ABI_SFREE(ylm) 779 ABI_SFREE(ylmgr) 780 ABI_SFREE(blkflg_car) 781 ABI_SFREE(d3etot_car) 782 783 ! Clean the header 784 call hdr%free() 785 786 DBG_EXIT("COLL") 787 788 end subroutine longwave
ABINIT/m_dfptlw_loop/dfptlw_out [ Functions ]
NAME
dfptlw_out
FUNCTION
Write the relevant spatial-dispersion quantities in Cartesian coordinates
COPYRIGHT
Copyright (C) 2022-2024 ABINIT group (MR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
blkflg_car(3,mpert,3,mpert,3,mpert) =flags for each element of the 3DTE d3etot_car(2,3,mpert,3,mpert,3,mpert) =array with the cartesian thir-order derivatives lw_qdrpl= flag that activates quadrupoles calculation lw_flexo= flag that activates flexoelectric tensor calculation mpert =maximum number of ipert natom = number of atoms in unit cell
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
824 #if defined HAVE_CONFIG_H 825 #include "config.h" 826 #endif 827 828 #include "abi_common.h" 829 830 831 subroutine dfptlw_out(blkflg_car,d3etot_car,lw_flexo,lw_qdrpl,lw_natopt,mpert,natom,ucvol) 832 833 use defs_basis 834 use m_errors 835 use m_profiling_abi 836 837 implicit none 838 839 !Arguments ------------------------------------ 840 !scalars 841 integer,intent(in) :: lw_flexo,lw_qdrpl,lw_natopt,mpert,natom 842 real(dp),intent(in) :: ucvol 843 !arrays 844 integer,intent(in) :: blkflg_car(3,mpert,3,mpert,3,mpert) 845 real(dp),intent(out) :: d3etot_car(2,3,mpert,3,mpert,3,mpert) 846 847 !Local variables------------------------------- 848 !scalar 849 integer :: beta,delta,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,istr 850 !arrays 851 integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 852 real(dp),allocatable :: qdrp(:,:,:,:,:,:,:) 853 real(dp) :: piezoci(2),piezofr(2),celastci(2) 854 855 ! ************************************************************************* 856 857 DBG_ENTER("COLL") 858 859 i3pert=natom+8 860 if (lw_qdrpl==1.or.lw_flexo==3.or.lw_flexo==1) then 861 write(ab_out,'(a)')' First real-space moment of the polarization response ' 862 write(ab_out,'(a)')' to an atomic displacementatom, in cartesian coordinates,' 863 write(ab_out,'(a)')' (1/ucvol factor not included),' 864 write(ab_out,'(a)')' efidir atom atdir qgrdir real part imaginary part' 865 i1pert=natom+2 866 do i3dir=1,3 867 do i1dir=1,3 868 do i2pert=1,natom 869 do i2dir=1,3 870 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 871 write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i2pert,i2dir,i3dir, & 872 & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert), & 873 & -d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 874 end if 875 end do 876 end do 877 end do 878 write(ab_out,*)' ' 879 end do 880 881 !Calculate cuadrupoles (symmetrize i1dir/i3dir) 882 ABI_MALLOC(qdrp,(2,3,mpert,3,mpert,3,mpert)) 883 i1pert=natom+2 884 do i2pert=1,natom 885 do i2dir=1,3 886 do i1dir=1,3 887 do i3dir=1,i1dir-1 888 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 889 !real part 890 qdrp(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=& 891 & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + & 892 & d3etot_car(2,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert) 893 894 qdrp(1,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert)=& 895 & qdrp(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 896 897 !imaginary part 898 qdrp(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=& 899 & -(d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + & 900 & d3etot_car(1,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert) ) 901 902 qdrp(2,i3dir,i1pert,i2dir,i2pert,i1dir,i3pert)=& 903 & qdrp(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 904 end if 905 end do 906 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)==1) then 907 !real part 908 qdrp(1,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)=& 909 & two*d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert) 910 911 !imaginary part 912 qdrp(2,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert)=& 913 &-two*d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i1dir,i3pert) 914 end if 915 end do 916 end do 917 end do 918 919 write(ab_out,'(a)')' Quadrupole tensor, in cartesian coordinates,' 920 write(ab_out,'(a)')' efidir atom atdir qgrdir real part imaginary part' 921 do i3dir=1,3 922 do i1dir=1,3 923 do i2pert=1,natom 924 do i2dir=1,3 925 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 926 write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i2pert,i2dir,i3dir, & 927 & qdrp(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert), & 928 & qdrp(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 929 end if 930 end do 931 end do 932 end do 933 write(ab_out,*)' ' 934 end do 935 ABI_FREE(qdrp) 936 937 write(ab_out,'(a)')' Electronic (clamped-ion) contribution to the piezoelectric tensor,' 938 write(ab_out,'(a)')' in cartesian coordinates, (from sum rule of dynamic quadrupoles or P^1 tensor)' 939 write(ab_out,'(a)')' efidir atdir qgrdir real part imaginary part' 940 do i3dir=1,3 941 do i1dir=1,3 942 do i2dir=1,3 943 piezoci=zero 944 do i2pert=1,natom 945 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 946 piezoci(1)=piezoci(1)+d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 947 piezoci(2)=piezoci(2)-d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 948 end if 949 end do 950 piezoci(1)=-piezoci(1)/ucvol 951 piezoci(2)=-piezoci(2)/ucvol 952 write(ab_out,'(3(i5,3x),2(1x,f20.10))') i1dir,i2dir,i3dir,piezoci(1),piezoci(2) 953 end do 954 end do 955 write(ab_out,*)' ' 956 end do 957 end if 958 959 if (lw_flexo==2.or.lw_flexo==1) then 960 write(ab_out,'(a)')' Clamped-ion flexoelectric tensor (type-II), in cartesian coordinates,' 961 write(ab_out,'(a)')' efidir qgrdir strdir1 strdir2 real part imaginary part' 962 i1pert=natom+2 963 do i3dir=1,3 964 do i2pert=natom+3,natom+4 965 do i2dir=1,3 966 istr=(i2pert-natom-3)*3+i2dir 967 beta=idx(2*istr-1); delta=idx(2*istr) 968 do i1dir=1,3 969 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 970 write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i3dir,beta,delta, & 971 & d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/ucvol, & 972 & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/ucvol 973 end if 974 end do 975 end do 976 write(ab_out,*)' ' 977 end do 978 end do 979 end if 980 981 if (lw_flexo==3.or.lw_flexo==1) then 982 write(ab_out,'(a)')' 1st real-space moment of IFCs, in cartesian coordinates,' 983 write(ab_out,'(a)')' iatdir iatom jatdir jatom qgrdir real part imaginary part' 984 do i3dir=1,3 985 do i1pert=1,natom 986 do i1dir=1,3 987 do i2pert=1,natom 988 do i2dir=1,3 989 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 990 write(ab_out,'(5(i5,4x),2(1x,f20.10))') i1dir,i1pert,i2dir,i2pert,i3dir, & 991 & -d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),& 992 & d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 993 end if 994 end do 995 end do 996 end do 997 end do 998 write(ab_out,*)' ' 999 end do 1000 1001 write(ab_out,'(a)')' Piezoelectric force-response tensor, in cartesian coordinates ' 1002 write(ab_out,'(a)')' (from sum rule of 1st moment of IFCs),' 1003 write(ab_out,'(a)')' (for non-vanishing forces in the cell it lacks an improper contribution),' 1004 write(ab_out,'(a)')' iatom iatddir jatddir qgrdir real part imaginary part' 1005 do i3dir=1,3 1006 do i1pert=1,natom 1007 do i1dir=1,3 1008 do i2dir=1,3 1009 piezofr=zero 1010 do i2pert=1,natom 1011 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 1012 piezofr(1)=piezofr(1)-d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 1013 piezofr(2)=piezofr(2)+d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 1014 end if 1015 end do 1016 write(ab_out,'(4(i5,4x),2(1x,f20.10))') i1pert,i1dir,i2dir,i3dir, & 1017 & piezofr(1), piezofr(2) 1018 end do 1019 end do 1020 end do 1021 write(ab_out,*)' ' 1022 end do 1023 end if 1024 1025 if (lw_flexo==4.or.lw_flexo==1) then 1026 write(ab_out,'(a)')' Clamped-ion flexoelectric force-response tensor (type-II), in cartesian coordinates,' 1027 write(ab_out,'(a)')' atom atdir qgrdir strdir1 strdir2 real part imaginary part' 1028 do i3dir=1,3 1029 do i1pert=1,natom 1030 do i1dir=1,3 1031 do i2pert=natom+3, natom+4 1032 do i2dir=1,3 1033 istr=(i2pert-natom-3)*3+i2dir 1034 beta=idx(2*istr-1); delta=idx(2*istr) 1035 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 1036 write(ab_out,'(5(i5,3x),2(1x,f20.10))') i1pert,i1dir,i3dir,beta,delta, & 1037 & d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),& 1038 & d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 1039 end if 1040 end do 1041 end do 1042 write(ab_out,*)' ' 1043 end do 1044 end do 1045 end do 1046 1047 1048 write(ab_out,'(a)')' Clamped-ion elastic tensor, in cartesian coordinates ' 1049 write(ab_out,'(a)')' (from sum rule of clamped-ion flexoelectric force-response tensor),' 1050 write(ab_out,'(a)')' (for stressed cells it lacks an improper contribution),' 1051 write(ab_out,'(a)')' atdir qgrdir strdir1 strdir2 real part imaginary part' 1052 do i1dir=1,3 1053 do i3dir=1,i1dir 1054 do i2pert=natom+3, natom+4 1055 do i2dir=1,3 1056 istr=(i2pert-natom-3)*3+i2dir 1057 beta=idx(2*istr-1); delta=idx(2*istr) 1058 celastci=zero 1059 do i1pert=1,natom 1060 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 1061 celastci(1)=celastci(1)+d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 1062 celastci(2)=celastci(2)+d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 1063 end if 1064 end do 1065 write(ab_out,'(4(i5,3x),2(1x,f20.10))') i1dir,i3dir,beta,delta, & 1066 & celastci(1)/ucvol,celastci(2)/ucvol 1067 end do 1068 end do 1069 write(ab_out,*)' ' 1070 end do 1071 end do 1072 end if 1073 1074 if (lw_natopt==1) then 1075 write(ab_out,'(a)')' Natural optical activity tensor, in cartesian coordinates,' 1076 write(ab_out,'(a)')' (1/ucvol factor not included),' 1077 write(ab_out,'(a)')' efidir1 efidir2 qgrdir real part imaginary part' 1078 i1pert=natom+2 1079 i2pert=natom+2 1080 do i3dir=1,3 1081 do i1dir=1,3 1082 do i2dir=1,3 1083 if (blkflg_car(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 1084 write(ab_out,'(3(i5,3x),2(1x,f20.10))') i1dir,i2dir,i3dir, & 1085 & -four*pi*d3etot_car(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert), & 1086 & four*pi*d3etot_car(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 1087 end if 1088 end do 1089 end do 1090 write(ab_out,*)' ' 1091 end do 1092 end if 1093 1094 DBG_EXIT("COLL") 1095 1096 end subroutine dfptlw_out
ABINIT/m_longwave [ Modules ]
NAME
m_longwave
FUNCTION
DFPT long-wave calculation of spatial dispersion properties
COPYRIGHT
Copyright (C) 2019-2024 ABINIT group (MR, MS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
SOURCE
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 module m_longwave 25 26 use defs_basis 27 use m_errors 28 use m_xmpi 29 use defs_datatypes 30 use defs_abitypes, only : MPI_type 31 use defs_wvltypes 32 use m_dtset 33 use m_dtfil 34 use m_xcdata 35 use m_hdr 36 use m_ebands 37 use m_wffile 38 39 use m_pspini, only : pspini 40 use m_common, only : setup1 41 use m_pawfgr, only : pawfgr_type, pawfgr_init 42 use m_pawrhoij, only : pawrhoij_type 43 use m_paw_dmft, only : paw_dmft_type 44 use m_pawrad, only : pawrad_type 45 use m_pawtab, only : pawtab_type 46 use m_drivexc, only : check_kxc 47 use m_rhotoxc, only : rhotoxc 48 use m_ioarr, only : read_rhor 49 use m_symtk, only : matr3inv,symmetrize_xred 50 use m_kg, only : kpgio 51 use m_inwffil, only : inwffil 52 use m_spacepar, only : setsym 53 use m_mkrho, only : mkrho 54 use m_fft, only : fourdp 55 use m_ddb, only : ddb_type,lwcart 56 use m_ddb_hdr, only : ddb_hdr_type 57 use m_mkcore, only : mkcore 58 use m_dfptlw_loop, only : dfptlw_loop 59 use m_dfptlw_nv, only : dfptlw_nv 60 use m_dfptlw_pert, only : preca_ffnl 61 use m_initylmg, only : initylmg 62 use m_dynmat, only : d3lwsym, sylwtens 63 use m_geometry, only : symredcart 64 65 implicit none 66 67 private