TABLE OF CONTENTS
ABINIT/dfptlw_geom [ Functions ]
NAME
dfptlw_geom
FUNCTION
This routine computes the nonvariational geometric contribution to the third-order energy derivative of the flexoelectric force-response tensor.
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cplex: if 1, several magnitudes are REAL, if 2, COMPLEX dimffnl= third dimension of ffnl_k dtset <type(dataset_type)>=all input variables for this dataset ffnl_k(dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives for this k point gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg i1dir,i2dir,i3dir=directions of the corresponding perturbations i1pert,i2pert = type of perturbation that has to be computed ikpt=number of the k-point isppol=1 for unpolarized, 2 for spin-polarized istwf_k=parameter that describes the storage of wfs kg_k(3,npw_k)=reduced planewave coordinates. kpt(3)=reduced coordinates of k point natom= number of atoms in the cell mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k natpert=number of atomic displacement perturbations nband_k=number of bands at this k point for that spin polarization n2dq= second dimension of d3etot_tgeom_k nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other npw_k=number of plane waves at this k point nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nylmgr=second dimension of ylmgr_k occ_k(nband_k)=occupation number for each band (usually 2) for each k. psps <type(pseudopotential_type)>=variables related to pseudopotentials rmet(3,3)=real space metric (bohr**2) rprimd(3,3) = dimensional primitive translations (bohr) vpsp1_i1pertdqdq(cplex*nfft,nspden,n2dq)= local potential of second-order gradient Hamiltonian for i1pert vpsp1_i1pertdq_geom(cplex*nfft,nspden,3)= local potential of first-order gradient Hamiltonian for i1pert wrt i3dir and i2dir useylmgr= if 1 use the derivative of spherical harmonics wtk_k=weight assigned to the k point. ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics for the k point
OUTPUT
d3etot_tgeom_k(2,n2dq)= nonvariational geometric contribution to d3etot for
SIDE EFFECTS
NOTES
PARENTS
m_dfpt_lw
CHILDREN
dfpt_vlocaldq,dfpt_vlocaldqdq,dotprod_g,getgh1dqc,getgh1dqc_setup init_rf_hamiltonian,mkkpg,rf_hamkq%free,rf_hamkq%load_spin rf_transgrid_and_pack
SOURCE
410 subroutine dfptlw_geom(cg,d3etot_tgeom_k,dimffnl,dtset,ffnl_k, & 411 & gs_hamkq,icg, & 412 & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt, & 413 & isppol,istwf_k,kg_k,kpt,mkmem,mpi_enreg,natom,mpw,nband_k,n2dq,nfft, & 414 & ngfft,npw_k,nspden,nsppol,nylmgr,occ_k, & 415 & psps,rmet,rprimd,useylmgr,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,wtk_k,ylm_k,ylmgr_k) 416 417 !Arguments ------------------------------------ 418 !scalars 419 integer,intent(in) :: dimffnl,icg,ikpt,isppol,istwf_k 420 integer,intent(in) :: i1dir,i1pert,i2dir,i2pert,i3dir 421 integer,intent(in) :: natom,mkmem,mpw,nband_k,nfft 422 integer,intent(in) :: npw_k,n2dq,nspden,nsppol,nylmgr 423 integer,intent(in) :: useylmgr 424 real(dp),intent(in) :: wtk_k 425 type(dataset_type),intent(in) :: dtset 426 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 427 type(MPI_type),intent(in) :: mpi_enreg 428 type(pseudopotential_type),intent(in) :: psps 429 430 !arrays 431 integer,intent(in) :: kg_k(3,npw_k),ngfft(18) 432 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 433 real(dp),intent(in) :: ffnl_k(npw_k,dimffnl,psps%lmnmax,psps%ntypat) 434 real(dp),intent(in) :: kpt(3),occ_k(nband_k) 435 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 436 real(dp),intent(in) :: vpsp1_i1pertdqdq(2*nfft,nspden,n2dq) 437 real(dp),intent(in) :: vpsp1_i1pertdq_geom(2*nfft,nspden,3) 438 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 439 real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 440 real(dp),intent(out) :: d3etot_tgeom_k(2,n2dq) 441 442 !Local variables------------------------------- 443 !scalars 444 integer :: beta,delta,dimffnlk,dimffnl1,gamma,iband,idq,ii,ipw,istr,nkpg,nkpg1,nylmgrpart 445 integer :: optlocal,optnl,q1dir,q2dir,reuse_ffnlk,reuse_ffnl1,tim_getgh1c,useylmgr1 446 real(dp) :: doti,dotr 447 type(pawfgr_type) :: pawfgr 448 type(rf_hamiltonian_type) :: rf_hamkq 449 450 !arrays 451 integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 452 real(dp) :: q1dirs(2),q2dirs(2) 453 real(dp),allocatable :: cwave0i(:,:) 454 real(dp),allocatable :: dkinpw(:) 455 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 456 real(dp),allocatable :: gh1dqc(:,:),gh1dqpkc(:,:),gvloc1dqc(:,:),gvnl1dqc(:,:) 457 real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:),kpg_pk(:,:),ph3d(:,:,:),ph3d1(:,:,:) 458 real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1dq(:,:,:,:), dum_vpsp(:) 459 real(dp),allocatable :: vpsp1dq(:),part_ylmgr_k(:,:,:) 460 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 461 462 ! ************************************************************************* 463 464 DBG_ENTER("COLL") 465 466 !Definitions 467 tim_getgh1c=0 468 useylmgr1=useylmgr;optlocal=1;optnl=1 469 nylmgrpart=3 470 nkpg=3 471 d3etot_tgeom_k(:,:)=zero 472 reuse_ffnlk=1 ; if (dtset%ffnl_lw==1) reuse_ffnlk=0 473 reuse_ffnl1=1 ; if (dtset%ffnl_lw==1) reuse_ffnl1=0 474 475 !Allocations 476 ABI_MALLOC(cwave0i,(2,npw_k*dtset%nspinor)) 477 ABI_MALLOC(dum_vpsp,(nfft)) 478 ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 479 ABI_MALLOC(dum_cwaveprj,(0,0)) 480 ABI_MALLOC(vpsp1dq,(2*nfft)) 481 ABI_MALLOC(vlocal1dq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 482 ABI_MALLOC(gh1dqc,(2,npw_k*dtset%nspinor)) 483 ABI_MALLOC(gvloc1dqc,(2,npw_k*dtset%nspinor)) 484 ABI_MALLOC(gvnl1dqc,(2,npw_k*dtset%nspinor)) 485 ABI_MALLOC(part_ylmgr_k,(npw_k,3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 486 part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:) 487 ABI_MALLOC(gh1dqpkc,(2,npw_k*dtset%nspinor)) 488 ABI_MALLOC(kpg_pk,(npw_k,nkpg)) 489 490 !Generate k+G vectors 491 call mkkpg(kg_k,kpg_pk,kpt,nkpg,npw_k) 492 493 !Since this is a type-I term, it has to be done for both up and down 494 !extradiagonal shear strains 495 gamma=i3dir 496 do idq=1, n2dq 497 if (i2pert==natom+3) then 498 istr=i2dir 499 else 500 istr=idq*3+i2dir 501 endif 502 beta=idx(2*istr-1); delta=idx(2*istr) 503 504 !----------------------------------------------------------------------------------------------- 505 ! q1-gradient of atomic displacement 1st order hamiltonian: 506 ! < u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}_{\q1dir} \delta_{\beta\q2dir}| u_{i,k}^{(0)} > 507 !----------------------------------------------------------------------------------------------- 508 dimffnlk=1 509 dimffnl1=2 510 q1dirs=(/gamma,delta/) 511 q2dirs=(/delta,gamma/) 512 do ii=1,2 513 q1dir=q1dirs(ii) 514 q2dir=q2dirs(ii) 515 516 if (beta==q2dir) then 517 518 !Get q-gradient of first-order local part of the pseudopotential 519 ! call dfpt_vlocaldq(atindx,2,gs_hamkq%gmet,gsqcut,i1dir,i1pert,mpi_enreg, & 520 ! & psps%mqgrid_vl,dtset%natom,& 521 ! & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 522 ! & ph1d,q1dir,psps%qgrid_vl,& 523 ! & dtset%qptn,ucvol,psps%vlspl,vpsp1dq) 524 ! write(300,*) vpsp1dq(:)-vpsp1_i1pertdq_geom(:,isppol,q1dir) 525 526 527 !Set up q-gradient of local potential vlocal1dq with proper dimensioning 528 vpsp1dq(:)=vpsp1_i1pertdq_geom(:,isppol,q1dir) 529 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,nfft,ngfft,& 530 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq) 531 532 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 533 call init_rf_hamiltonian(2,gs_hamkq,i1pert,rf_hamkq,& 534 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 535 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.) 536 537 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 538 if (dtset%ffnl_lw==0) then 539 ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat)) 540 ffnlk(:,1,:,:)=ffnl_k(:,1,:,:) 541 ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat)) 542 ffnl1(:,1,:,:)=ffnl_k(:,1,:,:) 543 ffnl1(:,2,:,:)=ffnl_k(:,1+q1dir,:,:) 544 end if 545 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,i1dir,i1pert,q1dir, & 546 & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrpart,useylmgr1,kg_k, & 547 & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,& 548 & ph3d,ph3d1,reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1) 549 550 !LOOP OVER BANDS 551 do iband=1,nband_k 552 553 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 554 555 !Read ket ground-state wavefunctions 556 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 557 558 !Compute < g |H^{\tau_{\kappa\alpha}}_{\q1dir} | u_{i,k}^{(0)} > 559 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, & 560 & i1dir,i1pert,mpi_enreg,optlocal,optnl,q1dir,rf_hamkq) 561 562 !Calculate: 563 !<u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\q1dir} | u_{i,k}^{(0)} > 564 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0i,gh1dqc, & 565 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 566 567 !Take into account the two pi factor from the term 568 !(\hat{p}_{k\beta + \frac{q_{\beta}}{2}}) appearing before the double q-derivation 569 !Take also into account here the -i factor and the complex conjugate 570 d3etot_tgeom_k(1,idq)=d3etot_tgeom_k(1,idq)-occ_k(iband)*half*doti*two_pi 571 d3etot_tgeom_k(2,idq)=d3etot_tgeom_k(2,idq)-occ_k(iband)*half*dotr*two_pi 572 573 end do !iband 574 575 !Clean the rf_hamiltonian 576 call rf_hamkq%free() 577 578 !Deallocations 579 ABI_FREE(kpg_k) 580 ABI_FREE(kpg1_k) 581 ABI_FREE(dkinpw) 582 ABI_FREE(kinpw1) 583 ABI_FREE(ffnlk) 584 ABI_FREE(ffnl1) 585 ABI_FREE(ph3d) 586 587 end if 588 589 end do !ii 590 591 !----------------------------------------------------------------------------------------------- 592 ! 2nd q-gradient of atomic displacement 1st order hamiltonian * momentum operator : 593 ! <u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\gamma\delta} (k+G)_{\beta} | u_{i,k}^{(0)} > 594 !----------------------------------------------------------------------------------------------- 595 596 !Get q-gradient of first-order local part of the pseudopotential 597 ! call dfpt_vlocaldqdq(atindx,2,gs_hamkq%gmet,gsqcut,i1dir,i1pert,mpi_enreg, & 598 ! & psps%mqgrid_vl,dtset%natom,& 599 ! & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 600 ! & ph1d,gamma,delta,psps%qgrid_vl,& 601 ! & dtset%qptn,ucvol,psps%vlspl,vpsp1dq) 602 603 !Set up q-gradient of local potential vlocal1dq with proper dimensioning 604 vpsp1dq(:)=vpsp1_i1pertdqdq(:,isppol,idq) 605 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 606 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq) 607 608 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 609 call init_rf_hamiltonian(2,gs_hamkq,i1pert,rf_hamkq,& 610 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 611 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.) 612 613 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 614 if (dtset%ffnl_lw==0) then 615 dimffnlk=1 616 ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat)) 617 ffnlk(:,1,:,:)=ffnl_k(:,1,:,:) 618 dimffnl1=10 619 ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat)) 620 ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:) 621 end if 622 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,i1dir,i1pert,gamma, & 623 & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgr,useylmgr1,kg_k, & 624 & ylm_k,kg_k,ylm_k,ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, & 625 & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1,qdir2=delta) 626 627 !LOOP OVER BANDS 628 do iband=1,nband_k 629 630 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 631 632 !Read ket ground-state wavefunctions 633 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 634 635 !Compute < g |H^{\tau_{\kappa\alpha}}_{\gamma\delta} | u_{i,k}^{(0)} > 636 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, & 637 & i1dir,i1pert,mpi_enreg,optlocal,optnl,gamma,rf_hamkq,qdir2=delta) 638 639 640 do ipw=1,npw_k 641 gh1dqpkc(:,ipw)=gh1dqc(:,ipw)*kpg_pk(ipw,beta) 642 end do 643 644 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0i,gh1dqpkc, & 645 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 646 647 !Take into account the two pi factor from the term 648 !(\hat{p}_{k\beta + \frac{q_{\beta}}{2}}) appearing before the double q-derivation 649 d3etot_tgeom_k(1,idq)=d3etot_tgeom_k(1,idq)-occ_k(iband)*doti*two_pi 650 d3etot_tgeom_k(2,idq)=d3etot_tgeom_k(2,idq)-occ_k(iband)*dotr*two_pi 651 652 end do !iband 653 654 !Clean the rf_hamiltonian 655 call rf_hamkq%free() 656 657 !Deallocations 658 ABI_FREE(kpg_k) 659 ABI_FREE(kpg1_k) 660 ABI_FREE(dkinpw) 661 ABI_FREE(kinpw1) 662 ABI_FREE(ffnlk) 663 ABI_FREE(ffnl1) 664 ABI_FREE(ph3d) 665 666 end do !idq 667 668 !scale by the k-point weight 669 d3etot_tgeom_k(:,:)=d3etot_tgeom_k(:,:)*wtk_k 670 671 !Deallocations 672 ABI_FREE(dum_cwaveprj) 673 ABI_FREE(gh1dqc) 674 ABI_FREE(gh1dqpkc) 675 ABI_FREE(gvloc1dqc) 676 ABI_FREE(gvnl1dqc) 677 ABI_FREE(vpsp1dq) 678 ABI_FREE(vlocal1dq) 679 ABI_FREE(dum_vpsp) 680 ABI_FREE(dum_vlocal) 681 ABI_FREE(kpg_pk) 682 ABI_FREE(cwave0i) 683 ABI_FREE(part_ylmgr_k) 684 685 DBG_EXIT("COLL") 686 687 end subroutine dfptlw_geom
ABINIT/m_dfptlw_nv [ Modules ]
NAME
m_dfptlw_nv
FUNCTION
FIXME: add description.
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 .
NOTES
PARENTS
CHILDREN
SOURCE
23 #if defined HAVE_CONFIG_H 24 #include "config.h" 25 #endif 26 27 #include "abi_common.h" 28 29 module m_dfptlw_nv 30 31 use defs_basis 32 use defs_abitypes 33 use defs_datatypes 34 use m_abicore 35 use m_xmpi 36 use m_errors 37 use m_mpinfo 38 use m_dtset 39 use m_hamiltonian 40 use m_cgtools 41 use m_wfk 42 use m_xmpi 43 use m_getgh1c 44 use m_mklocl 45 use m_pawcprj 46 use m_pawfgr 47 48 use m_dfpt_elt, only : dfpt_ewalddq, dfpt_ewalddqdq 49 use m_kg, only : mkkpg 50 use m_dynmat, only : cart39 51 52 implicit none 53 54 private
ABINIT/m_dfptlw_nv/dfptlw_nv [ Functions ]
NAME
dfptlw_nv
FUNCTION
This routine calculates the nonvariational Ewald contributions to the spatial-dispersion third-order energy derivatives.
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset gmet(3,3)=reciprocal space metric tensor in bohr**-2 gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1) mpert=maximum number of ipert my_natom=number of atoms treated by current processor rmet(3,3)=metric tensor in real space (length units squared) rprimd(3,3)=dimensional primitive translations (bohr) rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations that have to be computed 1 -> element has to be computed explicitely -1 -> use symmetry operations to obtain the corresponding element ucvol=unit cell volume in (whatever length scale units)**3 xred(3,natom)=relative coords of atoms in unit cell (dimensionless) zion(ntypat)=charge on each type of atom (real number) mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms
OUTPUT
d3etot_nv(2,3,mpert,3,mpert,3,mpert)= array with the nonvariational
contributions of d3etot
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
106 subroutine dfptlw_nv(d3etot_nv,dtset,gmet,gprimd,mpert,my_natom,rfpert,rmet,rprimd,ucvol,xred,zion, & 107 & mpi_atmtab,comm_atom ) ! optional arguments (parallelism)) 108 109 implicit none 110 111 !Arguments ------------------------------------ 112 !scalars 113 integer , intent(in) :: mpert,my_natom 114 real(dp) :: ucvol 115 type(dataset_type),intent(in) :: dtset 116 integer,optional,intent(in) :: comm_atom 117 118 !arrays 119 integer,optional,target,intent(in) :: mpi_atmtab(:) 120 integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert) 121 real(dp), intent(inout) :: d3etot_nv(2,3,mpert,3,mpert,3,mpert) 122 real(dp), intent(in) :: gmet(3,3),rmet(3,3),xred(3,dtset%natom),zion(*) 123 real(dp), intent(in) :: gprimd(3,3),rprimd(3,3) 124 125 !Local variables------------------------------- 126 !scalars 127 integer :: alpha,beta,delta,gamma,i1dir,i2dir,i3dir,ii,i1pert,i2pert,i3pert,istr,natom,sumg0 128 real(dp) :: fac,tmpim,tmpre 129 character(len=500) :: msg 130 131 !arrays 132 integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 133 integer :: flg1(3),flg2(3) 134 real(dp),allocatable :: dyewdq(:,:,:,:,:,:),dyewdqdq(:,:,:,:,:,:) 135 real(dp),allocatable :: dyewdqdq_tII(:,:,:,:,:,:) 136 real(dp) :: qphon(3),vec1(3),vec2(3) 137 138 ! ************************************************************************* 139 140 DBG_ENTER("COLL") 141 142 !Initialiations 143 natom=dtset%natom 144 145 if (dtset%lw_flexo==1.or.dtset%lw_flexo==3) then 146 147 !1st q-gradient of Ewald contribution to the IFCs 148 ABI_MALLOC(dyewdq,(2,3,natom,3,natom,3)) 149 sumg0=0;qphon(:)=zero 150 call dfpt_ewalddq(dyewdq,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,zion,& 151 & mpi_atmtab=mpi_atmtab,comm_atom=comm_atom) 152 153 i3pert=natom+8 154 do i1pert=1,natom 155 do i1dir=1,3 156 do i2pert=1,natom 157 do i2dir=1,3 158 do i3dir=1,3 159 tmpre=dyewdq(1,i1dir,i1pert,i2dir,i2pert,i3dir) 160 tmpim=dyewdq(2,i1dir,i1pert,i2dir,i2pert,i3dir) 161 if (abs(tmpre)>=tol8) d3etot_nv(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= tmpre 162 if (abs(tmpim)>=tol8) d3etot_nv(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= tmpim 163 end do 164 end do 165 end do 166 end do 167 end do 168 ABI_FREE(dyewdq) 169 170 end if 171 172 if (dtset%lw_flexo==1.or.dtset%lw_flexo==4) then 173 174 !2nd q-gradient of Ewald contribution to the IFCs 175 ABI_MALLOC(dyewdqdq,(2,3,natom,3,3,3)) 176 ABI_MALLOC(dyewdqdq_tII,(2,3,natom,3,3,3)) 177 sumg0=1;qphon(:)=zero 178 call dfpt_ewalddqdq(dyewdqdq,gmet,my_natom,natom,qphon,rmet,sumg0,dtset%typat,ucvol,xred,zion,& 179 & mpi_atmtab=mpi_atmtab,comm_atom=comm_atom) 180 181 182 !Convert the indexes labelling the strain perturbation into cartesian coordinates 183 !Transform the metric perturbation direction 184 !(treat it as an atomic displacement) 185 flg1(:)=1 186 do i1pert=1,natom 187 do i1dir=1,3 188 do gamma=1,3 189 do ii=1,2 190 do delta=1,3 191 do beta=1,3 192 vec1(beta)=dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma) 193 end do 194 call cart39(flg1,flg2,gprimd,i1pert,natom,rprimd,vec1,vec2) 195 do beta=1,3 196 dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)=vec2(beta) 197 end do 198 end do 199 end do 200 end do 201 end do 202 end do 203 204 !Transform the second q-gradient direction 205 !(treat it as an electric field) 206 do i1pert=1,natom 207 do i1dir=1,3 208 do gamma=1,3 209 do ii=1,2 210 do beta=1,3 211 do delta=1,3 212 vec1(delta)=dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma) 213 end do 214 call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2) 215 do delta=1,3 216 dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)=vec2(delta) 217 end do 218 end do 219 end do 220 end do 221 end do 222 end do 223 224 !Transform the first q-gradient direction 225 !(treat it as an electric field) 226 do i1pert=1,natom 227 do i1dir=1,3 228 do ii=1,2 229 do beta=1,3 230 do delta=1,3 231 do gamma=1,3 232 vec1(gamma)=dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma) 233 end do 234 call cart39(flg1,flg2,gprimd,natom+2,natom,rprimd,vec1,vec2) 235 do gamma=1,3 236 dyewdqdq(ii,i1dir,i1pert,beta,delta,gamma)=vec2(gamma) 237 end do 238 end do 239 end do 240 end do 241 end do 242 end do 243 244 !Convert to a type-II quantity 245 dyewdqdq_tII(:,:,:,:,:,:)=zero 246 do i1pert=1,natom 247 do alpha=1,3 248 do gamma=1,3 249 do beta=1,3 250 do delta=1,3 251 dyewdqdq_tII(:,alpha,i1pert,gamma,beta,delta)= & 252 & dyewdqdq(:,alpha,i1pert,beta,delta,gamma) + & 253 & dyewdqdq(:,alpha,i1pert,delta,gamma,beta) - & 254 & dyewdqdq(:,alpha,i1pert,gamma,beta,delta) 255 end do 256 end do 257 end do 258 end do 259 end do 260 ABI_FREE(dyewdqdq) 261 262 !Transform back the first q-gradient direction to reduced coordinates 263 !(treat it as an electric field) 264 fac=two_pi**2 265 do i1pert=1,natom 266 do alpha=1,3 267 do ii=1,2 268 do delta=1,3 269 do beta=1,3 270 do gamma=1,3 271 vec1(gamma)=dyewdqdq_tII(ii,alpha,i1pert,gamma,beta,delta) 272 end do 273 call cart39(flg1,flg2,transpose(rprimd),natom+2,natom,transpose(gprimd),vec1,vec2) 274 do gamma=1,3 275 dyewdqdq_tII(ii,alpha,i1pert,gamma,beta,delta)=vec2(gamma)*fac 276 end do 277 end do 278 end do 279 end do 280 end do 281 end do 282 283 i3pert=natom+8 284 do i1pert=1,natom 285 do i1dir=1,3 286 do i2pert=natom+3,natom+4 287 do i2dir=1,3 288 istr=(i2pert-natom-3)*3+i2dir 289 beta=idx(2*istr-1); delta=idx(2*istr) 290 do i3dir=1,3 291 tmpre=dyewdqdq_tII(1,i1dir,i1pert,i3dir,beta,delta) 292 tmpim=dyewdqdq_tII(2,i1dir,i1pert,i3dir,beta,delta) 293 if (abs(tmpre)>=tol8) d3etot_nv(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= half*tmpre 294 if (abs(tmpim)>=tol8) d3etot_nv(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)= half*tmpim 295 end do 296 end do 297 end do 298 end do 299 end do 300 ABI_FREE(dyewdqdq_tII) 301 302 end if 303 304 !Print results 305 if (dtset%prtvol>=10) then 306 write(msg,'(3a)') ch10,' LONGWAVE NONVARIATIONAL EWALD D3ETOT: ',ch10 307 call wrtout(std_out,msg,'COLL') 308 call wrtout(ab_out,msg,'COLL') 309 do i1pert=1,mpert 310 do i1dir=1,3 311 do i2pert=1,mpert 312 do i2dir=1,3 313 do i3pert=1,mpert 314 do i3dir=1,3 315 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 316 tmpre=d3etot_nv(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 317 tmpim=d3etot_nv(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 318 if (abs(tmpre)>zero.or.abs(tmpim)>zero) then 319 write(msg,'(3(a,i2,a,i1),2f18.8)') & 320 ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir,& 321 & tmpre, tmpim 322 call wrtout(std_out,msg,'COLL') 323 call wrtout(ab_out,msg,'COLL') 324 end if 325 end if 326 end do 327 end do 328 end do 329 end do 330 end do 331 end do 332 write(msg,'(a)') ch10 333 call wrtout(std_out,msg,'COLL') 334 call wrtout(ab_out,msg,'COLL') 335 end if 336 337 DBG_EXIT("COLL") 338 339 end subroutine dfptlw_nv