LCOV - code coverage report
Current view: top level - src/72_response - m_dfptlw_wf.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 98.2 % 276 271
Test Date: 2026-11-09 20:57:54 Functions: 100.0 % 1 1
Branches: 61.7 % 926 571

             Branch data     Line data    Source code
       1                 :             : !{\src2tex{textfont=tt}}
       2                 :             : !!****m* ABINIT/m_dfptlw_wf
       3                 :             : !! NAME
       4                 :             : !!  m_dfptlw_wf
       5                 :             : !!
       6                 :             : !! FUNCTION
       7                 :             : !!  FIXME: add description.
       8                 :             : !!
       9                 :             : !! COPYRIGHT
      10                 :             : !!  Copyright (C) 2022-2026 ABINIT group (MR)
      11                 :             : !!  This file is distributed under the terms of the
      12                 :             : !!  GNU General Public License, see ~abinit/COPYING
      13                 :             : !!  or http://www.gnu.org/copyleft/gpl.txt .
      14                 :             : !!
      15                 :             : !! SOURCE
      16                 :             : 
      17                 :             : #if defined HAVE_CONFIG_H
      18                 :             : #include "config.h"
      19                 :             : #endif
      20                 :             : 
      21                 :             : #include "abi_common.h"
      22                 :             : 
      23                 :             : module m_dfptlw_wf
      24                 :             : 
      25                 :             :  use defs_basis
      26                 :             :  use defs_abitypes
      27                 :             :  use defs_datatypes
      28                 :             :  use m_dtset
      29                 :             :  use m_errors
      30                 :             :  use m_profiling_abi
      31                 :             :  use m_hamiltonian
      32                 :             :  use m_cgtools
      33                 :             :  use m_pawcprj
      34                 :             :  use m_pawfgr
      35                 :             :  use m_wfk
      36                 :             :  use m_xmpi
      37                 :             :  use m_getgh1c
      38                 :             :  use m_mklocl
      39                 :             : 
      40                 :             :  use m_fstrings, only : itoa, sjoin
      41                 :             :  use m_io_tools, only : file_exists
      42                 :             :  use m_time, only : cwtime
      43                 :             :  use m_kg, only : mkkpg
      44                 :             : 
      45                 :             :  implicit none
      46                 :             : 
      47                 :             :  public :: dfpt_1wf
      48                 :             : 
      49                 :             :  private
      50                 :             : 
      51                 :             : ! *************************************************************************
      52                 :             : 
      53                 :             : contains
      54                 :             : !!***
      55                 :             : 
      56                 :             : !!****f* ABINIT/m_dfptlw_wf/dfpt_1wf
      57                 :             : !! NAME
      58                 :             : !!  dfpt_1wf
      59                 :             : !!
      60                 :             : !! FUNCTION
      61                 :             : !!  Compute the spin, band and kpt resolved contributions
      62                 :             : !!  to the spatial-dispersion third-order energy derivatives
      63                 :             : !!  that depend on first-order response functions.
      64                 :             : !!
      65                 :             : !! INPUTS
      66                 :             : !!  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k
      67                 :             : !!  cplex: if 1, several magnitudes are REAL, if 2, COMPLEX
      68                 :             : !!  ddk_f = wf files
      69                 :             : !!  d2_dkdk_f = wf files
      70                 :             : !!  dimffnl= third dimension of ffnl_k
      71                 :             : !!  dtset <type(dataset_type)>=all input variables for this dataset
      72                 :             : !!  eig1_k(2*nband_k**2)=1st-order eigenvalues at k for i1pert,i1dir
      73                 :             : !!  eig2_k(2*nband_k**2)=1st-order eigenvalues at k for i2pert,i2dir
      74                 :             : !!  ffnl_k(dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives for this k point
      75                 :             : !!  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
      76                 :             : !!  cg1 = first derivative of cg with respect the perturbation i1pert
      77                 :             : !!  cg2 = first derivative of cg with respect the perturbation i2pert
      78                 :             : !!  icg=shift to be applied on the location of data in the array cg
      79                 :             : !!  i1dir,i2dir,i3dir=directions of the corresponding perturbations
      80                 :             : !!  i1pert,i2pert,i3pert = type of perturbation that has to be computed
      81                 :             : !!  ikpt=number of the k-point
      82                 :             : !!  isppol=1 for unpolarized, 2 for spin-polarized
      83                 :             : !!  istwf_k=parameter that describes the storage of wfs
      84                 :             : !!  kg_k(3,npw_k)=reduced planewave coordinates.
      85                 :             : !!  kpt(3)=reduced coordinates of k point
      86                 :             : !!  mkmem =number of k points treated by this node
      87                 :             : !!  mpi_enreg=information about MPI parallelization
      88                 :             : !!  mpw=maximum dimensioned size of npw or wfs at k
      89                 :             : !!  natom= number of atoms in the unit cell
      90                 :             : !!  natpert=number of atomic displacement perturbations
      91                 :             : !!  nband_k=number of bands at this k point for that spin polarization
      92                 :             : !!  nfft=(effective) number of FFT grid points (for this proc)
      93                 :             : !!  ngfft(1:18)=integer array with FFT box dimensions and other
      94                 :             : !!  npw_k=number of plane waves at this k point
      95                 :             : !!  nspden=number of spin-density components
      96                 :             : !!  nsppol=1 for unpolarized, 2 for spin-polarized
      97                 :             : !!  nylmgr=second dimension of ylmgr_k
      98                 :             : !!  occ_k(nband_k)=occupation number for each band (usually 2) for each k.
      99                 :             : !!  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
     100                 :             : !!  psps <type(pseudopotential_type)>=variables related to pseudopotentials
     101                 :             : !!  rmet(3,3)=real space metric (bohr**2)
     102                 :             : !!  rprimd(3,3) = dimensional primitive translations (bohr)
     103                 :             : !!  samepert= .true. if i1pert=i2pert and i1dir=i2dir
     104                 :             : !!  useylmgr= if 1 use the derivative of spherical harmonics
     105                 :             : !!  vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of first-order
     106                 :             : !!          gradient Hamiltonian for i1pert
     107                 :             : !!  vpsp1_i2pertdq(cplex*nfft,nspden,n2dq)= local potential of first-order
     108                 :             : !!          gradient Hamiltonian for i2pert
     109                 :             : !!  wtk_k=weight assigned to the k point.
     110                 :             : !!  ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point
     111                 :             : !!  ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical
     112                 :             : !!                                                                      harmonics for the k point
     113                 :             : !!
     114                 :             : !! OUTPUT
     115                 :             : !! d3etot_t(1-5)_k= stationary 1wf contributions to the third-order energy
     116                 :             : !!                  derivatives for kpt
     117                 :             : !!
     118                 :             : !! SOURCE
     119                 :             : 
     120                 :       70860 : subroutine dfpt_1wf(cg,cg1,cg2,cplex,ddk_f,d2_dkdk_f,d2_dkdk_f2,&
     121                 :             :      & d3etot_t1_k,d3etot_t2_k,d3etot_t3_k,&
     122                 :       70860 :      & d3etot_t4_k,d3etot_t5_k,dimffnl,dtset,eig1_k,eig2_k,ffnl_k,gs_hamkq,icg,&
     123                 :             :      & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt,isppol,istwf_k,&
     124                 :       70860 :      & kg_k,kpt,mkmem,mpi_enreg,mpw,natom,nband_k,&
     125                 :       70860 :      & n1dq,n2dq,nfft,ngfft,npw_k,nspden,nsppol,nylmgr,occ_k,&
     126                 :             :      & pawfgr,psps,rmet,rprimd,samepert,useylmgr,&
     127                 :       70860 :      & vpsp1_i1pertdq,vpsp1_i2pertdq,&
     128                 :       70860 :      & wtk_k,ylm_k,ylmgr_k)
     129                 :             : 
     130                 :             : !Arguments ------------------------------------
     131                 :             : !scalars
     132                 :             :  integer,intent(in) :: cplex,dimffnl,i1dir,i1pert,i2dir,i2pert,i3dir
     133                 :             :  integer,intent(in) :: icg,ikpt,isppol,istwf_k
     134                 :             :  integer,intent(in) :: mkmem,mpw,natom,nband_k,n1dq,n2dq,nfft
     135                 :             :  integer,intent(in) :: npw_k,nspden,nsppol,nylmgr
     136                 :             :  integer,intent(in) :: useylmgr
     137                 :             :  real(dp),intent(in) :: wtk_k
     138                 :             :  logical,intent(in) :: samepert
     139                 :             :  type(dataset_type),intent(in) :: dtset
     140                 :             :  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
     141                 :             :  type(MPI_type),intent(in) :: mpi_enreg
     142                 :             :  type(pseudopotential_type),intent(in) :: psps
     143                 :             :  type(wfk_t),intent(inout) :: ddk_f,d2_dkdk_f, d2_dkdk_f2
     144                 :             :  type(pawfgr_type),intent(in) :: pawfgr
     145                 :             : 
     146                 :             : !arrays
     147                 :             :  integer,intent(in) :: kg_k(3,npw_k),ngfft(18)
     148                 :             :  real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
     149                 :             :  real(dp),intent(in) :: cg1(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
     150                 :             :  real(dp),intent(in) :: cg2(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol)
     151                 :             :  real(dp),intent(out) :: d3etot_t1_k(2)
     152                 :             :  real(dp),intent(out) :: d3etot_t2_k(2)
     153                 :             :  real(dp),intent(out) :: d3etot_t3_k(2)
     154                 :             :  real(dp),intent(out) :: d3etot_t4_k(2,n2dq)
     155                 :             :  real(dp),intent(out) :: d3etot_t5_k(2,n1dq)
     156                 :             :  real(dp),intent(in) :: eig1_k(2*nband_k**2),eig2_k(2*nband_k**2)
     157                 :             :  real(dp),intent(in) :: ffnl_k(npw_k,dimffnl,psps%lmnmax,psps%ntypat)
     158                 :             :  real(dp),intent(in) :: kpt(3),occ_k(nband_k)
     159                 :             :  real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
     160                 :             :  real(dp),intent(in) :: vpsp1_i1pertdq(2*nfft,nspden,n1dq)
     161                 :             :  real(dp),intent(in) :: vpsp1_i2pertdq(2*nfft,nspden,n2dq)
     162                 :             :  real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
     163                 :             :  real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
     164                 :             : 
     165                 :             : !Local variables-------------------------------
     166                 :             : !scalars
     167                 :             :  integer :: berryopt,dimffnlk,dimffnl1,iband,idir,idq,ii,jband,nkpg,nkpg1,nylmgrtmp
     168                 :             :  integer :: offset_cgi,opt_gvnl1,optlocal,optnl,reuse_ffnlk,reuse_ffnl1,sij_opt
     169                 :             :  integer :: size_wf,tim_getgh1c,usepaw,usevnl,useylmgr1
     170                 :             :  real(dp) :: cprodi,cprodr,doti,dotr,dum_lambda,fac,tmpim,tmpre
     171                 :             :  real(dp) :: cpu,wall,gflops
     172                 :             :  logical :: with_nonlocal_i1pert,with_nonlocal_i2pert
     173                 :       70860 :  type(rf_hamiltonian_type) :: rf_hamkq
     174                 :             : 
     175                 :             : !arrays
     176                 :       70860 :  real(dp),allocatable :: cg1_aux(:,:),cg1_ddk(:,:,:),cwave0i(:,:),cwave0j(:,:)
     177                 :       70860 :  real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
     178                 :       70860 :  real(dp),allocatable :: dkinpw(:),gv1c(:,:)
     179                 :       70860 :  real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:)
     180                 :       70860 :  real(dp),allocatable :: gvloc1dqc(:,:),gvnl1dqc(:,:)
     181                 :             :  real(dp) :: cj_h1_ci(2),dum_grad_berry(1,1),dum_gs1(1,1),dum_gvnl1(1,1)
     182                 :       70860 :  real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:)
     183                 :       70860 :  real(dp),allocatable :: part_ylmgr_k(:,:,:),ph3d(:,:,:),ph3d1(:,:,:)
     184                 :       70860 :  real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1(:,:,:,:),dum_vpsp(:)
     185                 :       70860 :  real(dp),allocatable :: vpsp1(:)
     186                 :       70860 :  type(pawcprj_type),allocatable :: dum_cwaveprj(:,:)
     187                 :             : 
     188                 :             : ! *************************************************************************
     189                 :             : 
     190                 :             :  DBG_ENTER("COLL")
     191                 :             : 
     192                 :             : !Additional definitions
     193                 :       70860 :  tim_getgh1c=0
     194                 :       70860 :  useylmgr1=useylmgr
     195                 :       70860 :  usepaw=dtset%usepaw
     196                 :       70860 :  size_wf= dtset%nspinor*npw_k
     197         [ +  + ]:       70860 :  with_nonlocal_i1pert=.true. ; if (i1pert==natom+2) with_nonlocal_i1pert=.false.
     198         [ +  + ]:       70860 :  with_nonlocal_i2pert=.true. ; if (i2pert==natom+2) with_nonlocal_i2pert=.false.
     199         [ +  + ]:       70860 :  reuse_ffnlk=1 ; if (dtset%ffnl_lw==1) reuse_ffnlk=0
     200         [ +  + ]:       70860 :  reuse_ffnl1=1 ; if (dtset%ffnl_lw==1) reuse_ffnl1=0
     201                 :             : 
     202                 :             : !Additional allocations
     203   [ -  +  +  -  :      212580 :  ABI_MALLOC(cwave0i,(2,size_wf))
          +  -  +  -  -  
                +  -  + ]
     204   [ -  +  +  -  :      141720 :  ABI_MALLOC(cwave0j,(2,size_wf))
             -  +  -  + ]
     205   [ -  +  +  -  :      141720 :  ABI_MALLOC(cwavef1,(2,size_wf))
             -  +  -  + ]
     206   [ -  +  +  -  :      141720 :  ABI_MALLOC(cwavef2,(2,size_wf))
             -  +  -  + ]
     207   [ -  +  +  -  :      141720 :  ABI_MALLOC(cg1_aux,(2,size_wf))
             -  +  -  + ]
     208   [ -  +  +  -  :      141720 :  ABI_MALLOC(gv1c,(2,size_wf))
             -  +  -  + ]
     209   [ -  +  +  -  :      425160 :  ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
          -  +  +  -  -  
          +  +  -  -  +  
          +  -  +  -  +  
             -  -  +  -  
                      + ]
     210   [ -  +  +  -  :      212580 :  ABI_MALLOC(dum_vpsp,(nfft))
          +  -  +  -  -  
                +  -  + ]
     211   [ -  +  +  -  :      425160 :  ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
          -  +  +  -  -  
          +  +  -  -  +  
          +  -  +  -  +  
             -  -  +  -  
                      + ]
     212   [ -  +  +  -  :      212580 :  ABI_MALLOC(vpsp1,(cplex*nfft))
          +  -  +  -  -  
                +  -  + ]
     213         [ -  + ]:       70860 :  ABI_MALLOC(dum_cwaveprj,(0,0))
     214   [ -  +  +  -  :      283440 :  ABI_MALLOC(part_ylmgr_k,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
          +  +  +  -  +  
          -  +  +  -  +  
                   +  - ]
     215   [ +  +  +  +  :    60620160 :  part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:)
                   +  + ]
     216                 :             : 
     217                 :             : !------------------------------------T1------------------------------------------------
     218                 :             : !q1-gradient of gs Hamiltonian:
     219                 :             : ! < u_{i,k}^{\lambda1}} | \partial_{gamma} H^{(0)} | u_{i,k}^{\lambda2} >
     220                 :             : !--------------------------------------------------------------------------------------
     221                 :             : 
     222                 :       70860 :  call cwtime(cpu, wall, gflops, "start")
     223                 :             : 
     224                 :             : !Specific definitions
     225                 :       70860 :  d3etot_t1_k=zero
     226   [ +  +  +  +  :   264293592 :  vlocal1=zero
             +  +  +  + ]
     227                 :       70860 :  dum_lambda=zero
     228                 :       70860 :  berryopt=0;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0
     229                 :             : 
     230                 :             : !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup)
     231                 :             :  call rf_hamkq%init(cplex,gs_hamkq,natom+1,&
     232                 :             :  & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
     233                 :       70860 :  & mpi_spintab=mpi_enreg%my_isppoltab)
     234                 :             : 
     235                 :       70860 :  call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.)
     236                 :             : 
     237                 :             :  !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
     238         [ +  + ]:       70860 :  if (dtset%ffnl_lw==0) then
     239   [ -  +  +  -  :      193840 :    ABI_MALLOC(ffnlk,(npw_k,0,psps%lmnmax,psps%ntypat))
          -  +  +  -  -  
          +  +  -  -  +  
                   -  + ]
     240   [ -  +  +  -  :      242300 :    ABI_MALLOC(ffnl1,(npw_k,2,psps%lmnmax,psps%ntypat))
          -  +  +  -  -  
          +  +  -  +  -  
          +  -  -  +  -  
                      + ]
     241   [ +  +  +  +  :    75830651 :    ffnl1(:,1,:,:)=ffnl_k(:,1,:,:)
                   +  + ]
     242   [ +  +  +  +  :    75830651 :    ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:)
                   +  + ]
     243                 :             :  end if
     244                 :             :  call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,&                        ! In
     245                 :             :  kpt,kpt,i3dir,natom+1,natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In
     246                 :             :  npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,&               ! In
     247                 :             :  dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,&           ! Out
     248                 :       70860 :  reuse_ffnlk=reuse_ffnlk, reuse_ffnl1=reuse_ffnl1)                        ! Optional
     249                 :             : 
     250                 :             :  !LOOP OVER BANDS
     251         [ +  + ]:      356236 :  do iband=1,nband_k
     252                 :             : 
     253         [ -  + ]:      285376 :    if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
     254                 :             : 
     255                 :             :    !Select bks wf1
     256                 :      285376 :    offset_cgi = (iband-1)*size_wf+icg
     257   [ +  +  +  + ]:    90262492 :    cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi)
     258   [ +  +  +  + ]:    90262492 :    cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi)
     259                 :             : 
     260                 :             :    !Compute < g |\partial_{gamma} H^{(0)} | u_{i,k}^{\lambda2} >
     261                 :             :    call getgh1c(berryopt,cwavef2,dum_cwaveprj,gv1c,dum_grad_berry,&
     262                 :             :  & dum_gs1,gs_hamkq,dum_gvnl1,i3dir,natom+1,(/dum_lambda/),mpi_enreg,1,optlocal,&
     263         [ +  + ]:      570752 :  & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
     264                 :             : 
     265                 :             :    !Apply the dot product with the ket wf (take into account occupation here)
     266                 :             :    ! < u_{i,k}^{\lambda1}} | \partial_{gamma} H^{(0)} | u_{i,k}^{lambda2}} >
     267                 :             :    call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,gv1c, &
     268                 :      285376 :  & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
     269                 :             : 
     270                 :      285376 :    d3etot_t1_k(1)=d3etot_t1_k(1)+occ_k(iband)*dotr
     271                 :      356236 :    d3etot_t1_k(2)=d3etot_t1_k(2)+occ_k(iband)*doti
     272                 :             : 
     273                 :             :  end do !iband
     274                 :             : 
     275                 :             : !Clean rf_hamiltonian
     276                 :       70860 :  call rf_hamkq%free()
     277                 :             : 
     278                 :             :  !Deallocations
     279         [ -  + ]:       70860 :  ABI_FREE(kpg_k)
     280         [ -  + ]:       70860 :  ABI_FREE(kpg1_k)
     281         [ -  + ]:       70860 :  ABI_FREE(dkinpw)
     282         [ -  + ]:       70860 :  ABI_FREE(kinpw1)
     283         [ -  + ]:       70860 :  ABI_FREE(ffnlk)
     284         [ -  + ]:       70860 :  ABI_FREE(ffnl1)
     285         [ -  + ]:       70860 :  ABI_FREE(ph3d)
     286                 :             : 
     287                 :       70860 :  call cwtime(cpu, wall, gflops, "stop")
     288                 :             : 
     289                 :             : !------------------------------------T2------------------------------------------------
     290                 :             : !q-gradient of CB projector x rf Hamiltonian lambda 2:
     291                 :             : ! < u_{i,k}^{\lambda1}} | \partial_{gamma} Q_k H^{\lambda2} | u_{i,k}^{(0)} >
     292                 :             : !--------------------------------------------------------------------------------------
     293                 :             : 
     294                 :             : !Create array for ddk 1wf from file
     295   [ -  +  +  -  :      283440 :  ABI_MALLOC(cg1_ddk,(2,size_wf,nband_k))
          -  +  +  -  +  
          -  +  -  -  +  
                   -  + ]
     296         [ +  + ]:      356236 :  do iband=1,nband_k
     297                 :      285376 :    call ddk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=cg1_aux)
     298   [ +  +  +  + ]:    90333352 :    cg1_ddk(:,:,iband)=cg1_aux(:,:)
     299                 :             :  end do
     300                 :             : 
     301                 :       70860 :  call cwtime(cpu, wall, gflops, "start")
     302                 :             : 
     303                 :             : !For \lambda1=\lambda2 T2 is inferred from the cc of T3
     304         [ +  + ]:       70860 : if (.not.samepert) then
     305                 :             : 
     306                 :             :   !Specific definitions
     307                 :       70088 :   d3etot_t2_k=zero
     308                 :             : 
     309                 :             :   !LOOP OVER BANDS
     310         [ +  + ]:      352296 :   do iband=1,nband_k
     311                 :             : 
     312         [ +  - ]:      282208 :     if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
     313                 :             : 
     314                 :             :     !Select bks wfs
     315                 :      282208 :     offset_cgi = (iband-1)*size_wf+icg
     316   [ +  +  +  + ]:    89198692 :     cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi)
     317                 :             : 
     318                 :             :     !LOOP OVER BANDS
     319         [ +  + ]:     1519528 :     do jband=1,nband_k
     320                 :             : 
     321                 :             :       !Select ddk wf1
     322   [ +  +  +  + ]:   388853008 :       cg1_aux(:,:)=cg1_ddk(:,:,jband)
     323                 :             : 
     324                 :             :       !Load < u_{j,k}^{(0) | H^{\lambda2}+V^{\lambda2}} | u_{i,k}^{(0)} >
     325                 :     1167232 :       ii=2*jband-1+(iband-1)*2*nband_k
     326                 :     1167232 :       cj_h1_ci(1)=eig2_k(ii)
     327                 :     1167232 :       cj_h1_ci(2)=eig2_k(ii+1)
     328                 :             : 
     329                 :             :       !Calculate: < u_{i,k}^{lambda1}} | u_{j,k}^{k_{\gamma}} >
     330                 :             :       call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,cg1_aux, &
     331                 :     1167232 :     & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
     332                 :             : 
     333                 :             :       !Calculate the contribution to T2
     334                 :     1167232 :       cprodr=dotr*cj_h1_ci(1)-doti*cj_h1_ci(2)
     335                 :     1167232 :       cprodi=dotr*cj_h1_ci(2)+doti*cj_h1_ci(1)
     336                 :     1167232 :       d3etot_t2_k(1)=d3etot_t2_k(1)-cprodr*occ_k(iband)
     337                 :     1449440 :       d3etot_t2_k(2)=d3etot_t2_k(2)-cprodi*occ_k(iband)
     338                 :             : 
     339                 :             :     end do !jband
     340                 :             : 
     341                 :             :   end do !iband
     342                 :             : 
     343                 :             : end if !samepert
     344                 :             : 
     345                 :       70860 :  call cwtime(cpu, wall, gflops, "stop")
     346                 :             : 
     347                 :             : !------------------------------------T3------------------------------------------------
     348                 :             : !rf Hamiltonian lambda 1 x q-gradient of CB projector
     349                 :             : ! < u_{i,k}^{(0) | (H^{\lambda1})^{\dagger} \partial_{gamma} Q_k | u_{i,k}^{\lambda2}}  >
     350                 :             : !--------------------------------------------------------------------------------------
     351                 :             : 
     352                 :       70860 :  call cwtime(cpu, wall, gflops, "start")
     353                 :             : !Specific definitions
     354                 :       70860 :  d3etot_t3_k=zero
     355                 :             : 
     356                 :             :  !LOOP OVER BANDS
     357         [ +  + ]:      356236 :  do iband=1,nband_k
     358                 :             : 
     359         [ +  - ]:      285376 :    if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
     360                 :             : 
     361                 :             :    !Select bks wfs
     362                 :      285376 :    offset_cgi = (iband-1)*size_wf+icg
     363   [ +  +  +  + ]:    90262492 :    cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi)
     364                 :             : 
     365                 :             :    !LOOP OVER BANDS
     366         [ +  + ]:     1538060 :    do jband=1,nband_k
     367                 :             : 
     368                 :             :      !Select ddk wf1
     369   [ +  +  +  + ]:   396115408 :      cg1_aux(:,:)=cg1_ddk(:,:,jband)
     370                 :             : 
     371                 :             :      !Load (< u_{j,k}^{(0) | H^{\lambda1}+V^{\lambda1}} | u_{i,k}^{(0)} >)^*
     372                 :     1181824 :      ii=2*jband-1+(iband-1)*2*nband_k
     373                 :     1181824 :      cj_h1_ci(1)=eig1_k(ii)
     374                 :     1181824 :      cj_h1_ci(2)=-eig1_k(ii+1)
     375                 :             : 
     376                 :             :      !Calculate: < u_{j,k}^{k_{\gamma}} | u_{i,k}^{lambda2}} >
     377                 :             :      call dotprod_g(dotr,doti,istwf_k,size_wf,2,cg1_aux,cwavef2, &
     378                 :     1181824 :    & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
     379                 :             : 
     380                 :             :      !Calculate the contribution to T3
     381                 :     1181824 :      cprodr=dotr*cj_h1_ci(1)-doti*cj_h1_ci(2)
     382                 :     1181824 :      cprodi=dotr*cj_h1_ci(2)+doti*cj_h1_ci(1)
     383                 :     1181824 :      d3etot_t3_k(1)=d3etot_t3_k(1)-cprodr*occ_k(iband)
     384                 :     1467200 :      d3etot_t3_k(2)=d3etot_t3_k(2)-cprodi*occ_k(iband)
     385                 :             : 
     386                 :             :    end do !jband
     387                 :             : 
     388                 :             :  end do !iband
     389                 :             : 
     390         [ -  + ]:       70860 :  ABI_FREE(cg1_ddk)
     391         [ -  + ]:       70860 :  ABI_FREE(cg1_aux)
     392         [ -  + ]:       70860 :  ABI_FREE(vpsp1)
     393         [ -  + ]:       70860 :  ABI_FREE(vlocal1)
     394                 :             : 
     395         [ +  + ]:       70860 : if (samepert) then
     396                 :         772 :   d3etot_t2_k(1)=d3etot_t3_k(1)
     397                 :         772 :   d3etot_t2_k(2)=-d3etot_t3_k(2)
     398                 :             : end if
     399                 :             : 
     400                 :       70860 :  call cwtime(cpu, wall, gflops, "stop")
     401                 :             : !------------------------------------T4------------------------------------------------
     402                 :             : !q-gradient of rf Hamiltonian lambda 2
     403                 :             : ! < u_{i,k}^{\lambda1} | H^{\lambda2}_{gamma} | u_{i,k}^{(0)} >
     404                 :             : !--------------------------------------------------------------------------------------
     405                 :             : 
     406                 :             : !For \lambda1=\lambda2 T4 is inferred from the cc of T5
     407         [ +  + ]:       70860 : if (.not.samepert) then
     408                 :             : 
     409                 :             :   !Specific definitions and allocations
     410   [ +  +  +  + ]:      373664 :    d3etot_t4_k=zero
     411                 :       70088 :    optlocal=1;optnl=1
     412                 :       70088 :    dimffnlk=0
     413         [ +  + ]:       70088 :    if (i2pert/=natom+2) then
     414   [ +  -  -  +  :      420432 :      ABI_MALLOC(vlocal1,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
          +  -  -  +  +  
          -  -  +  +  -  
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     415   [ +  -  -  +  :      280288 :      ABI_MALLOC(vpsp1,(2*nfft))
          +  -  +  -  +  
             -  -  +  -  
                      + ]
     416   [ -  +  +  -  :      140144 :      ABI_MALLOC(gvloc1dqc,(2,size_wf))
             -  +  -  + ]
     417   [ -  +  +  -  :      140144 :      ABI_MALLOC(gvnl1dqc,(2,size_wf))
             -  +  -  + ]
     418                 :             :    end if
     419                 :             :    if (i2pert<=natom) fac=-one
     420         [ +  + ]:       70088 :    if (i2pert==natom+2) fac=one
     421   [ +  +  +  + ]:       70088 :    if (i2pert==natom+3.or.i2pert==natom+4) fac=-half
     422         [ +  + ]:       70088 :    if (i2pert<=natom) then
     423                 :        7864 :      nylmgrtmp=3
     424                 :        7864 :      dimffnlk=1
     425                 :        7864 :      dimffnl1=2
     426   [ +  +  +  + ]:       62224 :    else if (i2pert==natom+3.or.i2pert==natom+4) then
     427                 :       62208 :      nylmgrtmp=nylmgr
     428                 :       62208 :      dimffnl1=10
     429         [ -  + ]:       62208 :      ABI_FREE(part_ylmgr_k)
     430   [ -  +  +  -  :      311040 :      ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
          -  +  +  -  +  
          +  +  -  +  -  
          +  +  -  +  +  
                      - ]
     431   [ +  +  +  +  :   167843032 :      part_ylmgr_k(:,:,:)=ylmgr_k(:,:,:)
                   +  + ]
     432                 :             :    end if
     433                 :             : 
     434                 :             :   !Do loop to compute both extradiagonal shear-strain components
     435         [ +  + ]:      171280 :    do idq=1,n2dq
     436                 :             : 
     437                 :      101192 :  call cwtime(cpu, wall, gflops, "start")
     438                 :             : 
     439         [ +  + ]:      101192 :      if (i2pert/=natom+2) then
     440         [ +  + ]:      101176 :        idir=i2dir; if (i2pert==natom+4) idir=idq*3+i2dir
     441                 :             :        !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup)
     442                 :             :        call rf_hamkq%init(2,gs_hamkq,i2pert,&
     443                 :             :        & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
     444                 :      101176 :        & mpi_spintab=mpi_enreg%my_isppoltab)
     445                 :             : 
     446                 :             :        !Set up local potentials with proper dimensioning
     447                 :             :        !and load the spin-dependent part of the Hamiltonians
     448   [ +  -  -  +  :   694024304 :        vpsp1=vpsp1_i2pertdq(:,isppol,idq)
             -  -  +  + ]
     449                 :             :        call rf_transgrid_and_pack(isppol,nspden,usepaw,2,nfft,nfft,ngfft,&
     450                 :      101176 :        & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1)
     451                 :             :        call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,&
     452                 :      101176 :        & with_nonlocal=with_nonlocal_i2pert)
     453                 :             : 
     454                 :             :        !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
     455         [ +  + ]:      101176 :        if (dtset%ffnl_lw==0) then
     456   [ -  +  +  -  :      343320 :          ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat))
          -  +  +  -  -  
          +  +  -  +  -  
          +  +  -  +  -  
                +  -  + ]
     457   [ +  +  +  +  :    18878035 :          if (dimffnlk==1) ffnlk(:,1,:,:)=ffnl_k(:,1,:,:)
             +  +  +  + ]
     458   [ -  +  +  -  :      411984 :          ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat))
          +  -  -  +  +  
          -  -  +  +  -  
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     459         [ +  + ]:       68664 :          if (dimffnl1==2) then
     460   [ +  +  +  +  :    18815827 :            ffnl1(:,1,:,:)=ffnl_k(:,1,:,:)
                   +  + ]
     461   [ +  +  +  +  :    18815827 :            ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:)
                   +  + ]
     462                 :             :          else
     463   [ +  +  +  +  :   838935144 :            ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:)
             +  +  +  + ]
     464                 :             :          end if
     465                 :             :        end if
     466                 :             : 
     467                 :             :        call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,i2pert,i3dir, &
     468                 :             :      & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrtmp,useylmgr1,kg_k, &
     469                 :             :      & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, &
     470   [ -  +  -  +  :      101176 :      & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1)
          -  +  -  +  -  
                +  -  + ]
     471                 :             : 
     472                 :             :      end if
     473                 :             : 
     474                 :             :      !LOOP OVER BANDS
     475         [ +  + ]:      507816 :      do iband=1,nband_k
     476                 :             : 
     477         [ -  + ]:      406624 :        if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
     478                 :             : 
     479                 :             :        !Select bks wfs
     480                 :      406624 :        offset_cgi = (iband-1)*size_wf+icg
     481   [ +  +  +  + ]:   126193012 :        cwavef1(:,:)= cg1(:,1+offset_cgi:size_wf+offset_cgi)
     482                 :             : 
     483                 :             :        !Perturbation-specific part
     484         [ +  + ]:      406624 :        if (i2pert==natom+2) then
     485                 :             :          if (samepert) then
     486                 :             :            ! Read from d2_dkdk_f
     487                 :             :            call d2_dkdk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c)
     488                 :             :          else
     489                 :             :            ! Read from d2_dkdk_f2
     490                 :         384 :            call d2_dkdk_f2%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c)
     491                 :             :          end if
     492                 :             :        else
     493   [ +  +  +  + ]:   125591572 :          cwave0i(:,:)= cg(:,1+offset_cgi:size_wf+offset_cgi)
     494                 :             : 
     495                 :             :          !Compute < g |H^{\lambda2}}_{\gamma} | u_{i,k}^{(0)} >
     496                 :             :          call getgh1dqc(cwave0i,dum_cwaveprj,gv1c,gvloc1dqc,gvnl1dqc,gs_hamkq, &
     497                 :      406240 :          & idir,i2pert,mpi_enreg,optlocal,optnl,i3dir,rf_hamkq)
     498                 :             :        end if
     499                 :             : 
     500                 :             :        !Calculate: < u_{j,k}^{\lambda1} | |H^{\lambda2}}_{\gamma} | u_{i,k}^{(0)} >
     501                 :             :        call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef1,gv1c, &
     502                 :      406624 :      & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
     503                 :             : 
     504                 :             :        !Calculate the contribution to T4
     505                 :      406624 :        d3etot_t4_k(1,idq)=d3etot_t4_k(1,idq)+dotr*occ_k(iband)
     506                 :      507816 :        d3etot_t4_k(2,idq)=d3etot_t4_k(2,idq)+doti*occ_k(iband)
     507                 :             : 
     508                 :             :      end do !iband
     509                 :             : 
     510         [ +  + ]:      101192 :      if (i2pert/=natom+2) then
     511                 :             : 
     512                 :             :        !Clean the rf_hamiltonian
     513                 :      101176 :        call rf_hamkq%free()
     514                 :             : 
     515                 :             :        !Deallocations
     516         [ -  + ]:      101176 :        ABI_FREE(kpg_k)
     517         [ -  + ]:      101176 :        ABI_FREE(kpg1_k)
     518         [ -  + ]:      101176 :        ABI_FREE(dkinpw)
     519         [ -  + ]:      101176 :        ABI_FREE(kinpw1)
     520         [ -  + ]:      101176 :        ABI_FREE(ffnlk)
     521         [ -  + ]:      101176 :        ABI_FREE(ffnl1)
     522         [ -  + ]:      101176 :        ABI_FREE(ph3d)
     523                 :             : 
     524                 :             :      end if
     525                 :             : 
     526                 :      101192 :  call cwtime(cpu, wall, gflops, "stop")
     527                 :             : 
     528                 :             :      !Apply the perturbation-dependent prefactors on T4
     529                 :      101192 :      tmpre=d3etot_t4_k(1,idq); tmpim=d3etot_t4_k(2,idq)
     530   [ +  +  +  + ]:      101192 :      if (i2pert<=natom.or.i2pert==natom+2) then
     531                 :        7880 :        d3etot_t4_k(1,idq)=-tmpim
     532                 :        7880 :        d3etot_t4_k(2,idq)=tmpre
     533                 :             :      end if
     534         [ +  + ]:      576048 :      d3etot_t4_k(:,idq)=d3etot_t4_k(:,idq)*fac
     535                 :             : 
     536                 :             :    end do !idq
     537                 :             : 
     538         [ +  + ]:       70088 :    if (i2pert/=natom+2) then
     539         [ -  + ]:       70072 :      ABI_FREE(gvloc1dqc)
     540         [ -  + ]:       70072 :      ABI_FREE(gvnl1dqc)
     541         [ -  + ]:       70072 :      ABI_FREE(vlocal1)
     542         [ -  + ]:       70072 :      ABI_FREE(vpsp1)
     543                 :             :    end if
     544                 :             : 
     545                 :             :  end if !samepert
     546                 :             : 
     547                 :             : !------------------------------------T5------------------------------------------------
     548                 :             : !q-gradient of rf Hamiltonian lambda 1
     549                 :             : ! < u_{i,k}^{(0)} | (H^{\lambda1}_{gamma})^{\dagger} | u_{i,k}^{\lambda2} >
     550                 :             : !--------------------------------------------------------------------------------------
     551                 :             : 
     552                 :             : !Specific definitions and allocations
     553   [ +  +  +  + ]:      283440 :  d3etot_t5_k=zero
     554                 :       70860 :  optlocal=1;optnl=1
     555                 :       70860 :  dimffnlk=0
     556         [ +  + ]:       70860 :  if (i1pert/=natom+2) then
     557   [ +  -  -  +  :      271872 :    ABI_MALLOC(vlocal1,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc))
          +  -  -  +  +  
          -  -  +  +  -  
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     558   [ +  -  -  +  :      181248 :    ABI_MALLOC(vpsp1,(2*nfft))
          +  -  +  -  +  
          -  -  +  -  +  
                   -  + ]
     559   [ -  +  +  -  :       90624 :    ABI_MALLOC(gvloc1dqc,(2,size_wf))
          -  +  -  +  -  
                      + ]
     560   [ -  +  +  -  :       90624 :    ABI_MALLOC(gvnl1dqc,(2,size_wf))
          -  +  -  +  -  
                      + ]
     561                 :             :  end if
     562         [ +  + ]:       70860 :  if (i1pert<=natom) fac=-one
     563         [ +  + ]:       70860 :  if (i1pert==natom+2) fac=one
     564   [ +  -  -  + ]:       70860 :  if (i1pert==natom+3.or.i1pert==natom+4) fac=-half
     565         [ +  + ]:       70860 :  if (i1pert<=natom) then
     566                 :       45312 :    nylmgrtmp=3
     567                 :       45312 :    dimffnlk=1
     568                 :       45312 :    dimffnl1=2
     569         [ -  + ]:       45312 :    ABI_FREE(part_ylmgr_k)
     570   [ -  +  +  -  :      181248 :    ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
          +  -  +  +  +  
             -  -  +  +  
                      - ]
     571   [ +  +  +  +  :    40872840 :    part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:)
                   +  + ]
     572   [ +  -  -  + ]:       25548 :  else if (i1pert==natom+3.or.i1pert==natom+4) then
     573                 :           0 :    nylmgrtmp=nylmgr
     574                 :           0 :    dimffnl1=10
     575         [ #  # ]:           0 :    ABI_FREE(part_ylmgr_k)
     576   [ #  #  #  #  :           0 :    ABI_MALLOC(part_ylmgr_k,(npw_k,nylmgrtmp,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
          #  #  #  #  #  
          #  #  #  #  #  
          #  #  #  #  #  
                      # ]
     577   [ -  -  -  -  :       25548 :    part_ylmgr_k(:,:,:)=ylmgr_k(:,:,:)
                   -  - ]
     578                 :             :  end if
     579                 :             : 
     580                 :             : !Do loop to compute both extradiagonal shear-strain components
     581         [ +  + ]:      141720 :  do idq=1,n1dq
     582                 :             : 
     583                 :       70860 :  call cwtime(cpu, wall, gflops, "start")
     584         [ +  + ]:       70860 :    if (i1pert/=natom+2) then
     585         [ -  + ]:       45312 :      idir=i1dir; if (i1pert==natom+4) idir=idq*3+i1dir
     586                 :             :      !Initialize rf Hamiltonian (the k-dependent part is prepared in getgh1c_setup)
     587                 :             :      call rf_hamkq%init(2,gs_hamkq,i1pert,&
     588                 :             :      & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
     589                 :       45312 :      & mpi_spintab=mpi_enreg%my_isppoltab)
     590                 :             : 
     591                 :             :      !Set up local potentials with proper dimensioning
     592                 :             :      !and load the spin-dependent part of the Hamiltonians
     593   [ +  -  -  +  :   305946624 :      vpsp1=vpsp1_i1pertdq(:,isppol,idq)
             -  -  +  + ]
     594                 :             :      call rf_transgrid_and_pack(isppol,nspden,usepaw,2,nfft,nfft,ngfft,&
     595                 :       45312 :      & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1)
     596                 :             :      call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,&
     597                 :       45312 :      & with_nonlocal=with_nonlocal_i1pert)
     598                 :             : 
     599                 :             :      !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
     600         [ +  + ]:       45312 :      if (dtset%ffnl_lw==0) then
     601   [ -  +  +  -  :      151040 :        ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,psps%ntypat))
          -  +  +  -  -  
          +  +  -  +  -  
          +  -  -  +  -  
                +  -  + ]
     602   [ +  -  +  +  :    40752008 :        if (dimffnlk==1) ffnlk(:,1,:,:)=ffnl_k(:,1,:,:)
             +  +  +  + ]
     603   [ -  +  +  -  :      181248 :        ABI_MALLOC(ffnl1,(npw_k,dimffnl1,psps%lmnmax,psps%ntypat))
          +  -  -  +  +  
          -  -  +  +  -  
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     604         [ -  + ]:       30208 :        if (dimffnl1==2) then
     605   [ +  +  +  +  :    40752008 :          ffnl1(:,1,:,:)=ffnl_k(:,1,:,:)
                   +  + ]
     606   [ +  +  +  +  :    40752008 :          ffnl1(:,2,:,:)=ffnl_k(:,1+i3dir,:,:)
                   +  + ]
     607                 :             :        else
     608   [ #  #  #  #  :           0 :          ffnl1(:,1:dimffnl1,:,:)=ffnl_k(:,1:dimffnl1,:,:)
             #  #  #  # ]
     609                 :             :        end if
     610                 :             :      end if
     611                 :             :      call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,i1pert,i3dir, &
     612                 :             :    & dtset%natom,rmet,rprimd,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrtmp,useylmgr1,kg_k, &
     613                 :             :    & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, &
     614   [ -  +  -  +  :       45312 :    & reuse_ffnlk=reuse_ffnlk,reuse_ffnl1=reuse_ffnl1)
          -  +  -  +  -  
                +  -  + ]
     615                 :             : 
     616                 :             :    end if
     617                 :             : 
     618                 :             :    !LOOP OVER BANDS
     619         [ +  + ]:      356236 :    do iband=1,nband_k
     620                 :             : 
     621         [ -  + ]:      285376 :      if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle
     622                 :             : 
     623                 :             :      !Select bks wfs
     624                 :      285376 :      offset_cgi = (iband-1)*size_wf+icg
     625   [ +  +  +  + ]:    90262492 :      cwavef2(:,:)= cg2(:,1+offset_cgi:size_wf+offset_cgi)
     626                 :             : 
     627                 :             :      !Perturbation-specific part
     628         [ +  + ]:      285376 :      if (i1pert==natom+2) then
     629                 :      104128 :        call d2_dkdk_f%read_bks(iband,ikpt,isppol,xmpio_single,cg_bks=gv1c)
     630                 :             :      else
     631   [ +  +  +  + ]:    53892960 :        cwave0i(:,:)= cg(:,1+offset_cgi:size_wf+offset_cgi)
     632                 :             : 
     633                 :             :        !Compute < g |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} >
     634                 :             :        call getgh1dqc(cwave0i,dum_cwaveprj,gv1c,gvloc1dqc,gvnl1dqc,gs_hamkq, &
     635                 :      181248 :        & idir,i1pert,mpi_enreg,optlocal,optnl,i3dir,rf_hamkq)
     636                 :             :      end if
     637                 :             : 
     638                 :             :      !Calculate: < u_{j,k}^{\lambda2} | |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} >
     639                 :             :      call dotprod_g(dotr,doti,istwf_k,size_wf,2,cwavef2,gv1c, &
     640                 :      285376 :    & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
     641                 :             : 
     642                 :             :      !Calculate the contribution to T5:
     643                 :      285376 :      d3etot_t5_k(1,idq)=d3etot_t5_k(1,idq)+dotr*occ_k(iband)
     644                 :      356236 :      d3etot_t5_k(2,idq)=d3etot_t5_k(2,idq)+doti*occ_k(iband)
     645                 :             : 
     646                 :             :    end do !iband
     647                 :             : 
     648         [ +  + ]:       70860 :    if (i1pert/=natom+2) then
     649                 :             : 
     650                 :             :      !Clean the rf_hamiltonian
     651                 :       45312 :      call rf_hamkq%free()
     652                 :             : 
     653                 :             :      !Deallocations
     654         [ -  + ]:       45312 :      ABI_FREE(kpg_k)
     655         [ -  + ]:       45312 :      ABI_FREE(kpg1_k)
     656         [ -  + ]:       45312 :      ABI_FREE(dkinpw)
     657         [ -  + ]:       45312 :      ABI_FREE(kinpw1)
     658         [ -  + ]:       45312 :      ABI_FREE(ffnlk)
     659         [ -  + ]:       45312 :      ABI_FREE(ffnl1)
     660         [ -  + ]:       45312 :      ABI_FREE(ph3d)
     661                 :             : 
     662                 :             :    end if
     663                 :             : 
     664                 :       70860 :  call cwtime(cpu, wall, gflops, "stop")
     665                 :             : 
     666                 :             :    !Apply the perturbation-dependent prefactors on T5
     667                 :       70860 :    tmpre=d3etot_t5_k(1,idq); tmpim=d3etot_t5_k(2,idq)
     668   [ +  +  +  - ]:       70860 :    if (i1pert<=natom.or.i1pert==natom+2) then
     669                 :       70860 :      d3etot_t5_k(1,idq)=-tmpim
     670                 :       70860 :      d3etot_t5_k(2,idq)=tmpre
     671                 :             :    end if
     672         [ +  + ]:      212580 :    d3etot_t5_k(:,idq)=d3etot_t5_k(:,idq)*fac
     673                 :             : 
     674                 :             :    !Apply now the conjugate complex:
     675                 :             :    !(< u_{j,k}^{\lambda2} | |H^{\lambda1}}_{\gamma} | u_{i,k}^{(0)} >)*
     676                 :       70860 :    tmpim=d3etot_t5_k(2,idq)
     677                 :      283440 :    d3etot_t5_k(2,idq)=-tmpim
     678                 :             : 
     679                 :             :  end do !idq
     680                 :             : 
     681                 :             : 
     682         [ +  + ]:       70860 :  if (i1pert/=natom+2) then
     683         [ -  + ]:       45312 :    ABI_FREE(gvloc1dqc)
     684         [ -  + ]:       45312 :    ABI_FREE(gvnl1dqc)
     685         [ -  + ]:       45312 :    ABI_FREE(vlocal1)
     686         [ -  + ]:       45312 :    ABI_FREE(vpsp1)
     687                 :             :  end if
     688                 :             : 
     689         [ +  + ]:       70860 : if (samepert) then
     690         [ +  + ]:        1544 :   d3etot_t4_k(1,:)=d3etot_t5_k(1,:)
     691         [ +  + ]:        1544 :   d3etot_t4_k(2,:)=-d3etot_t5_k(2,:)
     692                 :             : end if
     693                 :             : 
     694                 :             : !Scale d3etot_k contributions by the kpt weight
     695         [ +  + ]:      212580 : d3etot_t1_k(:)=d3etot_t1_k(:)*wtk_k
     696         [ +  + ]:      212580 : d3etot_t2_k(:)=d3etot_t2_k(:)*wtk_k
     697         [ +  + ]:      212580 : d3etot_t3_k(:)=d3etot_t3_k(:)*wtk_k
     698   [ +  +  +  + ]:      376752 : d3etot_t4_k(:,:)=d3etot_t4_k(:,:)*wtk_k
     699   [ +  +  +  + ]:      283440 : d3etot_t5_k(:,:)=d3etot_t5_k(:,:)*wtk_k
     700                 :             : 
     701                 :             : !Deallocations
     702         [ -  + ]:       70860 :  ABI_FREE(cwave0i)
     703         [ -  + ]:       70860 :  ABI_FREE(cwave0j)
     704         [ -  + ]:       70860 :  ABI_FREE(cwavef1)
     705         [ -  + ]:       70860 :  ABI_FREE(cwavef2)
     706         [ -  + ]:       70860 :  ABI_FREE(gv1c)
     707         [ -  + ]:       70860 :  ABI_FREE(dum_vpsp)
     708         [ -  + ]:       70860 :  ABI_FREE(dum_vlocal)
     709   [ +  -  -  +  :       70860 :  ABI_FREE(dum_cwaveprj)
          -  -  -  -  -  
                      + ]
     710         [ -  + ]:       70860 :  ABI_FREE(part_ylmgr_k)
     711                 :             : 
     712                 :             : 
     713                 :             :  DBG_EXIT("COLL")
     714                 :             : 
     715   [ -  +  -  +  :      425160 : end subroutine dfpt_1wf
          -  +  -  +  -  
          +  -  +  -  +  
          -  +  -  +  -  
          +  -  +  -  +  
          -  +  -  +  -  
          -  -  -  -  -  
                   -  + ]
     716                 :             : !!***
     717                 :             : 
     718                 :             : end module m_dfptlw_wf
     719                 :             : !!***
        

Generated by: LCOV version 2.3-1