LCOV - code coverage report
Current view: top level - src/72_response - m_dfptlw_pert.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 84.3 % 248 209
Test Date: 2026-11-10 14:34:30 Functions: 100.0 % 3 3
Branches: 66.4 % 584 388

             Branch data     Line data    Source code
       1                 :             : !{\src2tex{textfont=tt}}
       2                 :             : !!****m* ABINIT/m_dfptlw_pert
       3                 :             : !! NAME
       4                 :             : !!  m_dfptlw_pert
       5                 :             : !!
       6                 :             : !! FUNCTION
       7                 :             : !!
       8                 :             : !! COPYRIGHT
       9                 :             : !!  Copyright (C) 2022-2026 ABINIT group (MR)
      10                 :             : !!  This file is distributed under the terms of the
      11                 :             : !!  GNU General Public License, see ~abinit/COPYING
      12                 :             : !!  or http://www.gnu.org/copyleft/gpl.txt .
      13                 :             : !!
      14                 :             : !! SOURCE
      15                 :             : 
      16                 :             : #if defined HAVE_CONFIG_H
      17                 :             : #include "config.h"
      18                 :             : #endif
      19                 :             : 
      20                 :             : #include "abi_common.h"
      21                 :             : 
      22                 :             : module m_dfptlw_pert
      23                 :             : 
      24                 :             :  use defs_basis
      25                 :             :  use defs_abitypes
      26                 :             :  use defs_datatypes
      27                 :             :  use m_dtset
      28                 :             :  use m_dtfil
      29                 :             :  use m_errors
      30                 :             :  use m_profiling_abi
      31                 :             :  use m_hamiltonian
      32                 :             :  use m_pawcprj
      33                 :             :  use m_pawfgr
      34                 :             :  use m_wfk
      35                 :             :  use m_xmpi
      36                 :             :  use m_getgh1c
      37                 :             :  use m_mklocl
      38                 :             :  use m_initylmg,   only : initylmg
      39                 :             :  use m_fstrings,   only : itoa, sjoin
      40                 :             :  use m_io_tools,   only : file_exists
      41                 :             :  use m_time,       only : cwtime
      42                 :             :  use m_kg,         only : mkkpg
      43                 :             :  use m_mpinfo,     only : proc_distrb_cycle
      44                 :             :  use m_dfptlw_wf
      45                 :             :  use m_dfpt_mkvxc, only : dfpt_mkvxcggadq
      46                 :             :  use m_dfptlw_nv,  only : dfptlw_geom
      47                 :             :  use m_spacepar,   only : hartredq
      48                 :             :  use m_cgtools,    only : dotprod_vn
      49                 :             :  use m_mkffnl,     only : mkffnl
      50                 :             : 
      51                 :             :  implicit none
      52                 :             : 
      53                 :             :  public :: dfptlw_pert
      54                 :             :  public :: preca_ffnl
      55                 :             : 
      56                 :             :  private
      57                 :             : 
      58                 :             : ! *************************************************************************
      59                 :             : 
      60                 :             : contains
      61                 :             : !!***
      62                 :             : 
      63                 :             : !!****f* ABINIT/m_dfptlw_pert/dfptlw_pert
      64                 :             : !! NAME
      65                 :             : !!  dfptlw_pert
      66                 :             : !!
      67                 :             : !! FUNCTION
      68                 :             : !! Compute first-order response function contributions to the spatial-dispersion
      69                 :             : !! 3rd order energy derivatives of the longwave driver.
      70                 :             : !! The main inputs are :
      71                 :             : !!   - GS WFs and Hamiltonian (cg,gs_hamkq)
      72                 :             : !!   - 1st-order WFs for two perturbations i1pert/i1dir,i2pert/i2dir (cg1,cg2)
      73                 :             : !!   - 1st-order Local+SCF potentials for i1pert and i2pert
      74                 :             : !!   - 1st-order WFs DDK and 2nd-order WF D2_DKDK (d2_dkdk_f)
      75                 :             : !!
      76                 :             : !! INPUTS
      77                 :             : !!  cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = array for planewave
      78                 :             : !!                                          coefficients of wavefunctions
      79                 :             : !!  cg1 = first derivative of cg with respect the perturbation i1pert
      80                 :             : !!  cg2 = first derivative of cg with respect the perturbation i2pert
      81                 :             : !!  cplex= if 1, real space 1-order functions on FFT grid are REAL,
      82                 :             : !!          if 2, COMPLEX
      83                 :             : !!  dimffnl= third dimension of ffnl
      84                 :             : !!  dtset <type(dataset_type)>=all input variables for this dataset
      85                 :             : !!  eigen1(2*mband*mband*nkpt*nsppol)=1st-order eigenvalues for i1pert,i1dir (hartree)
      86                 :             : !!  eigen2(2*mband*mband*nkpt*nsppol)=1st-order eigenvalues for i2pert,i2dir (hartree)
      87                 :             : !!  ffnl(dtset%mkmem,dtset%mpw,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives
      88                 :             : !!  gmet(3,3)=reciprocal space metric tensor in bohr**-2
      89                 :             : !!  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
      90                 :             : !!  gsqcut=large sphere cut-off
      91                 :             : !!  i1dir,i2dir,i3dir=directions of the corresponding perturbations
      92                 :             : !!  i1pert,i2pert,i3pert = type of perturbation that has to be computed
      93                 :             : !!  kg(3,mpw*mkmem_rbz)=reduced planewave coordinates
      94                 :             : !!  kxc(nfft,nkxc)=exchange and correlation kernel
      95                 :             : !!  mband = maximum number of bands
      96                 :             : !!  mkmem_rbz = maximum number of k points which can fit in core memory
      97                 :             : !!  mk1mem = maximum number of k points for first-order WF
      98                 :             : !!           which can fit in core memory
      99                 :             : !!  mpert =maximum number of ipert
     100                 :             : !!  mpi_enreg=MPI-parallelisation information
     101                 :             : !!  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
     102                 :             : !!  mpw   = maximum number of planewaves in basis sphere (large number)
     103                 :             : !!  natom = number of atoms in unit cell
     104                 :             : !!  n1dq= third dimension of vlocal1_i1pertdq
     105                 :             : !!  n2dq= third dimension of vlocal1_i2pertdq
     106                 :             : !!  nfft= number of FFT grid points (for this proc)
     107                 :             : !!  ngfft(1:18)=integer array with FFT box dimensions and other
     108                 :             : !!  nkpt = number of k points
     109                 :             : !!  nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
     110                 :             : !!  nspden = number of spin-density components
     111                 :             : !!  nspinor = number of spinorial components of the wavefunctions
     112                 :             : !!  nsppol = number of channels for spin-polarization (1 or 2)
     113                 :             : !!  npwarr(nkpt) = array holding npw for each k point
     114                 :             : !!  nylmgr=second dimension of ylmgr_k
     115                 :             : !!  occ(mband*nkpt*nsppol) = occupation number for each band and k
     116                 :             : !!  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
     117                 :             : !!  psps <type(pseudopotential_type)> = variables related to pseudopotentials
     118                 :             : !!  rho1g1(2,nfft)=G-space RF electron density in electrons/bohr**3 (i1pert)
     119                 :             : !!  rho1r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i1pert)
     120                 :             : !!  rho2r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i2pert)
     121                 :             : !!  rmet(3,3)=real space metric tensor in bohr**2
     122                 :             : !!  rprimd(3,3) = dimensional primitive translations (bohr)
     123                 :             : !!  samepert= .true. if i1pert=i2pert and i1dir=i2dir
     124                 :             : !!  ucvol=volume of the unit cell
     125                 :             : !!  useylmgr= if 1 use the derivative of spherical harmonics
     126                 :             : !!  vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of first-order
     127                 :             : !!          gradient Hamiltonian for i1pert
     128                 :             : !!  vpsp1_i1pertdq(cplex*nfft,nspden,n1dq)= local potential of second-order
     129                 :             : !!          gradient Hamiltonian for i1pert
     130                 :             : !!  vpsp1_i1pertdq_geom(cplex*nfft,nspden,3)= local potential of first-order
     131                 :             : !!          gradient Hamiltonian for i1pert wrp to i3dir and i2dir
     132                 :             : !!  vpsp1_i2pertdq(cplex*nfft,nspden,n2dq)= local potential of first-order
     133                 :             : !!          gradient Hamiltonian for i2pert
     134                 :             : !!  ddk_f = wf files
     135                 :             : !!  d2_dkdk_f = wf files
     136                 :             : !!  ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics
     137                 :             : !!  ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics
     138                 :             : !!
     139                 :             : !! OUTPUT
     140                 :             : !!  d3etot(2,3,mpert,3,mpert,3,mpert) = third derivatives of the energy tensor
     141                 :             : !!  d3etot_t4(2,n2dq)= t4 term which might need to be converted to type-II
     142                 :             : !!  d3etot_t5(2,n1dq)= t5 term which might need to be converted to type-II
     143                 :             : !!  d3etot_tgeom(2,2)= Geometric term which needs to be converted to type-II
     144                 :             : !!
     145                 :             : !! SIDE EFFECTS
     146                 :             : !!  TO DO!
     147                 :             : !!
     148                 :             : !! SOURCE
     149                 :             : 
     150                 :         591 : subroutine dfptlw_pert(cg,cg1,cg2,cplex,d3etot,d3etot_t4,d3etot_t5,d3etot_tgeom,&
     151                 :         591 : & dimffnl,dtset,eigen1,eigen2,ffnl,gmet,gs_hamkq,gsqcut,i1dir,i2dir,i3dir,&
     152                 :         591 : & i1pert,i2pert,i3pert,kg,kxc,mband,mkmem_rbz,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,&
     153                 :             : & n1dq,n2dq,nfft,ngfft,nkpt,nkxc,&
     154                 :         591 : & nspden,nspinor,nsppol,npwarr,nylmgr,occ,pawfgr,psps,rho1g1,rho1r1,rho2r1,rmet,rprimd,samepert,&
     155                 :         591 : & ucvol,useylmgr,vpsp1_i1pertdq,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,vpsp1_i2pertdq,ddk_f,d2_dkdk_f,d2_dkdk_f2,ylm,ylmgr)
     156                 :             : 
     157                 :             : !Arguments ------------------------------------
     158                 :             : !scalars
     159                 :             :  integer,intent(in) :: cplex,dimffnl,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband
     160                 :             :  integer,intent(in) :: mk1mem,mkmem_rbz,mpert,mpsang,mpw,natom,n1dq,n2dq,nfft,nkpt,nkxc,nspden
     161                 :             :  integer,intent(in) :: nspinor,nsppol,nylmgr,useylmgr
     162                 :             :  real(dp),intent(in) :: gsqcut,ucvol
     163                 :             :  logical,intent(in) :: samepert
     164                 :             :  type(MPI_type),intent(inout) :: mpi_enreg
     165                 :             :  type(dataset_type),intent(in) :: dtset
     166                 :             :  type(pseudopotential_type),intent(in) :: psps
     167                 :             :  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
     168                 :             :  type(pawfgr_type),intent(in) :: pawfgr
     169                 :             :  type(wfk_t),intent(inout) :: ddk_f,d2_dkdk_f, d2_dkdk_f2
     170                 :             : 
     171                 :             : !arrays
     172                 :             :  integer,intent(in) :: kg(3,mpw*mkmem_rbz),ngfft(18),npwarr(nkpt)
     173                 :             :  real(dp),intent(in) :: eigen1(2*mband*mband*nkpt*nsppol)
     174                 :             :  real(dp),intent(in) :: eigen2(2*mband*mband*nkpt*nsppol)
     175                 :             :  real(dp),intent(in) :: ffnl(mkmem_rbz,mpw,dimffnl,psps%lmnmax,psps%ntypat)
     176                 :             :  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol)
     177                 :             :  real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
     178                 :             :  real(dp),intent(in) :: cg2(2,mpw*nspinor*mband*mk1mem*nsppol)
     179                 :             :  real(dp),intent(in) :: gmet(3,3),kxc(nfft,nkxc)
     180                 :             :  real(dp),intent(in) :: occ(mband*nkpt*nsppol)
     181                 :             :  real(dp),intent(in) :: rho1g1(2,nfft),rho1r1(cplex*nfft,dtset%nspden)
     182                 :             :  real(dp),intent(in) :: rho2r1(cplex*nfft,dtset%nspden)
     183                 :             :  real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
     184                 :             :  real(dp),intent(in) :: vpsp1_i1pertdq(2*nfft,nspden,n1dq)
     185                 :             :  real(dp),intent(in) :: vpsp1_i1pertdqdq(2*nfft,nspden,n2dq)
     186                 :             :  real(dp),intent(in) :: vpsp1_i1pertdq_geom(2*nfft,nspden,3)
     187                 :             :  real(dp),intent(in) :: vpsp1_i2pertdq(2*nfft,nspden,n2dq)
     188                 :             :  real(dp),intent(inout) :: d3etot(2,3,mpert,3,mpert,3,mpert)
     189                 :             :  real(dp),intent(out) :: d3etot_t4(2,n2dq),d3etot_t5(2,n1dq)
     190                 :             :  real(dp),intent(out) :: d3etot_tgeom(2,n2dq)
     191                 :             :  real(dp),intent(in) :: ylm(mpw*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
     192                 :             :  real(dp),intent(in) :: ylmgr(mpw*mk1mem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
     193                 :             : 
     194                 :             : !Variables ------------------------------------
     195                 :             : !scalars
     196                 :             :  integer :: bandtot,bd2tot,icg,idq,ierr,ii,ikc,ikg,ikpt,ilm,isppol,istwf_k,me,n1,n2,n3,n4,n5,n6
     197                 :             :  integer :: nband_k,npw_k,spaceworld,tim_getgh1c
     198                 :             :  integer :: usepaw
     199                 :             :  real(dp) :: tmpim,tmpre,wtk_k
     200                 :             :  real(dp) :: cpu, wall, gflops
     201                 :             :  character(len=1000) :: msg
     202                 :             :  logical :: with_nonlocal_i1pert, with_nonlocal_i2pert
     203                 :             : !arrays
     204                 :         591 :  integer,allocatable :: kg_k(:,:)
     205                 :             :  real(dp) :: d3etot_t1(2),d3etot_t1_k(2)
     206                 :             :  real(dp) :: d3etot_t2(2),d3etot_t2_k(2)
     207                 :             :  real(dp) :: d3etot_t3(2),d3etot_t3_k(2)
     208                 :        1182 :  real(dp) :: d3etot_t4_k(2,n2dq)
     209                 :        1182 :  real(dp) :: d3etot_t5_k(2,n1dq)
     210                 :        1182 :  real(dp) :: d3etot_tgeom_k(2,n2dq)
     211                 :             :  real(dp) :: d3etot_telec(2)
     212                 :             :  real(dp) :: e3tot(2),kpt(3)
     213                 :         591 :  real(dp),allocatable :: eig1_k(:),eig2_k(:),occ_k(:)
     214                 :         591 :  real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:)
     215                 :         591 :  real(dp),allocatable :: ffnl_k(:,:,:,:)
     216                 :             : 
     217                 :             : ! *************************************************************************
     218                 :             : 
     219                 :             :  DBG_ENTER("COLL")
     220                 :             : 
     221                 :         591 :  write(msg,'(2a,3(a,i2,a,i1))') ch10,'LONGWAVE : ',&
     222                 :        1182 :  ' perts : ',i1pert,'.',i1dir,' / ',i2pert,'.',i2dir,' / ',i3pert,'.',i3dir
     223                 :         591 :  call wrtout(std_out,msg,'COLL')
     224                 :         591 :  call wrtout(ab_out,msg,'COLL')
     225                 :             : 
     226                 :             : !Init parallelism
     227                 :         591 :  spaceworld=mpi_enreg%comm_cell
     228                 :         591 :  me=mpi_enreg%me_kpt
     229                 :             : 
     230                 :             : !Additional definitions
     231                 :         591 :  tim_getgh1c=0
     232                 :         591 :  usepaw=dtset%usepaw
     233                 :         591 :  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
     234                 :         591 :  n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
     235                 :         591 :  with_nonlocal_i1pert=.true. ; if (i1pert==natom+2) with_nonlocal_i1pert=.false.
     236                 :         591 :  with_nonlocal_i2pert=.true. ; if (i2pert==natom+2) with_nonlocal_i2pert=.false.
     237                 :             : 
     238                 :             : !Initialize d3etot parts
     239                 :         591 :  d3etot_t1=zero
     240                 :         591 :  d3etot_t2=zero
     241                 :         591 :  d3etot_t3=zero
     242   [ +  +  +  + ]:        3093 :  d3etot_t4=zero
     243   [ +  +  +  + ]:        2364 :  d3etot_t5=zero
     244                 :             :  d3etot_telec=zero
     245   [ +  +  +  + ]:        3093 :  d3etot_tgeom=zero
     246                 :             : 
     247                 :             : !Calculate the electrostatic contribution
     248                 :             :  call lw_elecstic(cplex,d3etot_telec,gmet,gs_hamkq%gprimd,gsqcut,&
     249                 :         591 : & i3dir,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,rho1g1,rho1r1,rho2r1,ucvol)
     250                 :             : 
     251                 :             : !Loop over spins
     252                 :         591 :  bandtot = 0
     253                 :         591 :  bd2tot = 0
     254                 :         591 :  icg=0
     255         [ +  + ]:        1182 :  do isppol = 1, nsppol
     256                 :             : 
     257                 :             : !  Loop over k-points
     258                 :         591 :    ikg = 0
     259                 :         591 :    ikc = 0
     260         [ +  + ]:       72042 :    do ikpt = 1, nkpt
     261                 :             : 
     262                 :       70860 :      nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
     263                 :       70860 :      npw_k = npwarr(ikpt)
     264                 :       70860 :      istwf_k = dtset%istwfk(ikpt)
     265                 :             : 
     266         [ -  + ]:       70860 :      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,mpi_enreg%me)) then
     267                 :           0 :        bandtot = bandtot + nband_k
     268                 :           0 :        bd2tot = bd2tot + 2*nband_k**2
     269                 :           0 :        cycle ! Skip the rest of the k-point loop
     270                 :             :      end if
     271                 :       70860 :      ikc= ikc + 1
     272                 :             : 
     273   [ -  +  +  -  :      212580 :      ABI_MALLOC(occ_k,(nband_k))
          +  -  +  -  -  
                +  -  + ]
     274         [ +  + ]:      356236 :      occ_k(:) = occ(1+bandtot:nband_k+bandtot)
     275                 :       70860 :      wtk_k    = dtset%wtk(ikpt)
     276         [ +  + ]:      283440 :      kpt(:) = dtset%kptns(:,ikpt)
     277                 :             : 
     278   [ +  -  -  +  :      212580 :      ABI_MALLOC(eig1_k,(2*nband_k**2))
          +  -  +  -  +  
             -  -  +  -  
                      + ]
     279   [ -  +  +  -  :      141720 :      ABI_MALLOC(eig2_k,(2*nband_k**2))
             -  +  -  + ]
     280   [ -  +  +  -  :      212580 :      ABI_MALLOC(kg_k,(3,npw_k))
          +  -  +  -  -  
                +  -  + ]
     281   [ -  +  +  -  :      283440 :      ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
          -  +  +  -  +  
          -  +  -  -  +  
                   -  + ]
     282   [ -  +  +  -  :      354300 :      ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
          -  +  +  -  +  
          +  +  -  +  -  
          +  +  -  +  -  
                      + ]
     283   [ -  +  +  -  :      425160 :      ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
          +  +  +  -  -  
          +  +  -  -  +  
          +  -  +  -  +  
             +  -  +  -  
                      + ]
     284                 :             : 
     285                 :             :      !Get plane-wave vectors and related data at k
     286   [ +  +  +  + ]:    29521220 :      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
     287         [ +  + ]:       70860 :      if (dtset%ffnl_lw==1) then
     288         [ +  - ]:       22400 :        if (psps%useylm==1) then
     289         [ +  + ]:      224000 :          do ilm=1,psps%mpsang*psps%mpsang
     290         [ +  + ]:    20138300 :            ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
     291                 :             :          end do
     292         [ +  - ]:       22400 :          if (useylmgr==1) then
     293         [ +  + ]:      224000 :            do ilm=1,psps%mpsang*psps%mpsang
     294         [ +  + ]:     2038400 :              do ii=1,nylmgr
     295         [ +  + ]:   181244700 :                ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
     296                 :             :              end do
     297                 :             :            end do
     298                 :             :          end if
     299                 :             :        end if
     300         [ +  - ]:       48460 :      else if (dtset%ffnl_lw==0) then
     301   [ +  +  +  +  :   665791688 :        ffnl_k(1:npw_k,:,:,:)=ffnl(ikc,1:npw_k,:,:,:)
             +  +  +  + ]
     302                 :             :      end if
     303                 :             : 
     304                 :             :      !Get matrix elements for uniform perturbations
     305         [ +  + ]:     2434508 :      eig1_k(:)=eigen1(1+bd2tot:2*nband_k**2+bd2tot)
     306         [ +  + ]:     2434508 :      eig2_k(:)=eigen2(1+bd2tot:2*nband_k**2+bd2tot)
     307                 :             : 
     308                 :             :      !Compute the stationary terms of d3etot depending on response functions
     309                 :             :      call dfpt_1wf(cg,cg1,cg2,cplex,ddk_f,d2_dkdk_f,d2_dkdk_f2,d3etot_t1_k,d3etot_t2_k,d3etot_t3_k,&
     310                 :             :      & d3etot_t4_k,d3etot_t5_k,dimffnl,dtset,eig1_k,eig2_k,ffnl_k,gs_hamkq,icg,&
     311                 :             :      & i1dir,i2dir,i3dir,i1pert,i2pert,ikpt,isppol,istwf_k,&
     312                 :             :      & kg_k,kpt,mkmem_rbz,mpi_enreg,mpw,natom,nband_k,&
     313                 :             :      & n1dq,n2dq,nfft,ngfft,npw_k,nspden,nsppol,nylmgr,occ_k,&
     314                 :             :      & pawfgr,psps,rmet,rprimd,samepert,useylmgr,&
     315                 :             :      & vpsp1_i1pertdq,vpsp1_i2pertdq,&
     316                 :       70860 :      & wtk_k,ylm_k,ylmgr_k)
     317                 :             : 
     318                 :             : !    Add the contribution from each k-point.
     319         [ +  + ]:      212580 :      d3etot_t1=d3etot_t1 + d3etot_t1_k
     320         [ +  + ]:      212580 :      d3etot_t2=d3etot_t2 + d3etot_t2_k
     321         [ +  + ]:      212580 :      d3etot_t3=d3etot_t3 + d3etot_t3_k
     322   [ +  +  +  + ]:      376752 :      d3etot_t4=d3etot_t4 + d3etot_t4_k
     323   [ +  +  +  + ]:      283440 :      d3etot_t5=d3etot_t5 + d3etot_t5_k
     324                 :             : 
     325                 :             :      !Compute the nonvariational geometric term
     326                 :       70860 :      call cwtime(cpu, wall, gflops, "start")
     327   [ +  +  +  +  :       70860 :      if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then
                   +  + ]
     328                 :             :        call dfptlw_geom(cg,d3etot_tgeom_k,dimffnl,dtset, &
     329                 :             :        &  ffnl_k,gs_hamkq,icg, &
     330                 :             :        &  i1dir,i2dir,i3dir,i1pert,i2pert,ikpt, &
     331                 :             :        &  isppol,istwf_k,kg_k,kpt,mkmem_rbz,mpi_enreg,natom,mpw,nband_k,n2dq,nfft, &
     332                 :             :        &  ngfft,npw_k,nspden,nsppol,nylmgr,occ_k, &
     333                 :       41472 :        &  psps,rmet,rprimd,useylmgr,vpsp1_i1pertdqdq,vpsp1_i1pertdq_geom,wtk_k,ylm_k,ylmgr_k)
     334                 :             : 
     335                 :             :        !Add the contribution from each k-point
     336   [ +  +  +  + ]:      257484 :        d3etot_tgeom=d3etot_tgeom + d3etot_tgeom_k
     337                 :             :      end if
     338                 :       70860 :      call cwtime(cpu, wall, gflops, "stop")
     339                 :             : 
     340                 :             : !    Keep track of total number of bands
     341                 :       70860 :      bandtot = bandtot + nband_k
     342                 :       70860 :      bd2tot = bd2tot + 2*nband_k**2
     343                 :             : 
     344                 :             : !    Shift arrays memory
     345                 :       70860 :      icg=icg+npw_k*dtset%nspinor*nband_k
     346                 :       70860 :      ikg=ikg+npw_k
     347                 :             : 
     348         [ -  + ]:       70860 :      ABI_FREE(eig1_k)
     349         [ -  + ]:       70860 :      ABI_FREE(eig2_k)
     350         [ -  + ]:       70860 :      ABI_FREE(occ_k)
     351         [ -  + ]:       70860 :      ABI_FREE(kg_k)
     352         [ -  + ]:       70860 :      ABI_FREE(ylm_k)
     353         [ -  + ]:       70860 :      ABI_FREE(ylmgr_k)
     354         [ -  + ]:      213171 :      ABI_FREE(ffnl_k)
     355                 :             : 
     356                 :             :    end do !ikpt
     357                 :             : 
     358                 :             :  end do !isppol
     359                 :             : 
     360                 :             : 
     361                 :             : !=== MPI communications ==================
     362                 :             :  if (xmpi_paral==1) then
     363                 :         591 :    call xmpi_sum(d3etot_t1,spaceworld,ierr)
     364                 :         591 :    call xmpi_sum(d3etot_t2,spaceworld,ierr)
     365                 :         591 :    call xmpi_sum(d3etot_t3,spaceworld,ierr)
     366                 :         591 :    call xmpi_sum(d3etot_t4,spaceworld,ierr)
     367                 :         591 :    call xmpi_sum(d3etot_t5,spaceworld,ierr)
     368                 :         591 :    call xmpi_sum(d3etot_tgeom,spaceworld,ierr)
     369                 :             :  end if
     370                 :             : 
     371                 :             : !Apply +i or -i in case of strain perturbation.
     372   [ +  -  -  + ]:         591 :  if (i1pert==natom+3.or.i1pert==natom+4) then
     373                 :           0 :    tmpre=d3etot_telec(1);tmpim=d3etot_telec(2) ; d3etot_telec(1)=tmpim;d3etot_telec(2)=-tmpre
     374                 :           0 :    tmpre=d3etot_t1(1);tmpim=d3etot_t1(2) ; d3etot_t1(1)=tmpim;d3etot_t1(2)=-tmpre
     375                 :           0 :    tmpre=d3etot_t2(1);tmpim=d3etot_t2(2) ; d3etot_t2(1)=tmpim;d3etot_t2(2)=-tmpre
     376                 :           0 :    tmpre=d3etot_t3(1);tmpim=d3etot_t3(2) ; d3etot_t3(1)=tmpim;d3etot_t3(2)=-tmpre
     377         [ #  # ]:           0 :    do idq=1,n2dq
     378                 :           0 :      tmpre=d3etot_t4(1,idq);tmpim=d3etot_t4(2,idq) ; d3etot_t4(1,idq)=tmpim;d3etot_t4(2,idq)=-tmpre
     379                 :             :    end do
     380         [ #  # ]:           0 :    do idq=1,n1dq
     381                 :         591 :      tmpre=d3etot_t5(1,idq);tmpim=d3etot_t5(2,idq) ; d3etot_t5(1,idq)=tmpim;d3etot_t5(2,idq)=-tmpre
     382                 :             :    end do
     383                 :             :  end if
     384   [ +  +  +  + ]:         591 :  if (i2pert==natom+3.or.i2pert==natom+4) then
     385                 :         486 :    tmpre=d3etot_telec(1);tmpim=d3etot_telec(2) ; d3etot_telec(1)=-tmpim;d3etot_telec(2)=tmpre
     386                 :         486 :    tmpre=d3etot_t1(1);tmpim=d3etot_t1(2) ; d3etot_t1(1)=-tmpim;d3etot_t1(2)=tmpre
     387                 :         486 :    tmpre=d3etot_t2(1);tmpim=d3etot_t2(2) ; d3etot_t2(1)=-tmpim;d3etot_t2(2)=tmpre
     388                 :         486 :    tmpre=d3etot_t3(1);tmpim=d3etot_t3(2) ; d3etot_t3(1)=-tmpim;d3etot_t3(2)=tmpre
     389         [ +  + ]:        1215 :    do idq=1,n2dq
     390                 :         729 :      tmpre=d3etot_t4(1,idq);tmpim=d3etot_t4(2,idq) ; d3etot_t4(1,idq)=-tmpim;d3etot_t4(2,idq)=tmpre
     391         [ +  + ]:        1215 :      if (i1pert<=natom) then
     392                 :         486 :        tmpre=d3etot_tgeom(1,idq);tmpim=d3etot_tgeom(2,idq) ; d3etot_tgeom(1,idq)=-tmpim;d3etot_tgeom(2,idq)=tmpre
     393                 :             :      end if
     394                 :             :    end do
     395         [ +  + ]:         972 :    do idq=1,n1dq
     396                 :        1077 :      tmpre=d3etot_t5(1,idq);tmpim=d3etot_t5(2,idq) ; d3etot_t5(1,idq)=-tmpim;d3etot_t5(2,idq)=tmpre
     397                 :             :    end do
     398                 :             :  end if
     399                 :             : 
     400                 :             : !Join all the contributions in e3tot except t4 and t5 which may need to be
     401                 :             : !converted to type-II in case of strain perturbation.
     402                 :             : !Apply here the two factor to the stationary wf1 contributions
     403                 :             : !(see PRB 105, 064101 (2022))
     404         [ +  + ]:        1773 :  d3etot_t1(:)=two*d3etot_t1(:)
     405         [ +  + ]:        1773 :  d3etot_t2(:)=two*d3etot_t2(:)
     406         [ +  + ]:        1773 :  d3etot_t3(:)=two*d3etot_t3(:)
     407   [ +  +  +  + ]:        3093 :  d3etot_t4(:,:)=two*d3etot_t4(:,:)
     408   [ +  +  +  + ]:        2364 :  d3etot_t5(:,:)=two*d3etot_t5(:,:)
     409         [ +  + ]:        1773 :  e3tot(:)=d3etot_t1(:)+d3etot_t2(:)+d3etot_t3(:)+d3etot_telec(:)
     410                 :             : 
     411                 :             : 
     412                 :             : !Before printing, set small contributions to zero
     413                 :             :  !Real parts
     414         [ +  + ]:         591 :  if (abs(d3etot_t1(1))<tol8) d3etot_t1(1)= zero
     415         [ +  + ]:         591 :  if (abs(d3etot_t2(1))<tol8) d3etot_t2(1)= zero
     416         [ +  + ]:         591 :  if (abs(d3etot_t3(1))<tol8) d3etot_t3(1)= zero
     417         [ +  + ]:        1425 :  do idq=1,n2dq
     418         [ +  + ]:         834 :    if (abs(d3etot_t4(1,idq))<tol8) d3etot_t4(1,idq)= zero
     419         [ +  + ]:        1425 :    if (abs(d3etot_tgeom(1,idq))<tol8) d3etot_tgeom(1,idq)= zero
     420                 :             :  end do
     421         [ +  + ]:        1182 :  do idq=1,n1dq
     422         [ +  + ]:        1182 :    if (abs(d3etot_t5(1,idq))<tol8) d3etot_t5(1,idq)= zero
     423                 :             :  end do
     424         [ +  + ]:         591 :  if (abs(d3etot_telec(1))<tol8) d3etot_telec(1)= zero
     425         [ +  + ]:         591 :  if (abs(e3tot(1))    <tol8)     e3tot(1)= zero
     426                 :             : 
     427                 :             :  !Imaginary parts
     428         [ +  + ]:         591 :  if (abs(d3etot_t1(2))<tol8) d3etot_t1(2)= zero
     429         [ -  + ]:         591 :  if (abs(d3etot_t2(2))<tol8) d3etot_t2(2)= zero
     430         [ -  + ]:         591 :  if (abs(d3etot_t3(2))<tol8) d3etot_t3(2)= zero
     431         [ +  + ]:        1425 :  do idq=1,n2dq
     432         [ -  + ]:         834 :    if (abs(d3etot_t4(2,idq))<tol8) d3etot_t4(2,idq)= zero
     433         [ +  + ]:        1425 :    if (abs(d3etot_tgeom(2,idq))<tol8) d3etot_tgeom(2,idq)= zero
     434                 :             :  end do
     435         [ +  + ]:        1182 :  do idq=1,n1dq
     436         [ +  + ]:        1182 :    if (abs(d3etot_t5(2,idq))<tol8) d3etot_t5(2,idq)= zero
     437                 :             :  end do
     438         [ +  + ]:         591 :  if (abs(d3etot_telec(2))<tol8) d3etot_telec(2)= zero
     439         [ +  + ]:         591 :  if (abs(e3tot(2))    <tol8)     e3tot(2)= zero
     440                 :             : 
     441         [ -  + ]:         591 :  if (dtset%prtvol>=10) then
     442                 :             :    write(msg,'(4(a,2(a,f18.8)),a)') &
     443                 :           0 :    ch10,'          d3etot_telec = ',d3etot_telec(1),  ',',d3etot_telec(2),&
     444                 :           0 :    ch10,'             d3etot_t1 = ',d3etot_t1(1),  ',',d3etot_t1(2),&
     445                 :           0 :    ch10,'             d3etot_t2 = ',d3etot_t2(1),  ',',d3etot_t2(2),&
     446                 :           0 :    ch10,'             d3etot_t3 = ',d3etot_t3(1),  ',',d3etot_t3(2)
     447                 :           0 :    call wrtout(std_out,msg,'COLL')
     448                 :           0 :    call wrtout(ab_out,msg,'COLL')
     449         [ #  # ]:           0 :    if (n2dq==1) then
     450                 :             :      write(msg,'(2(a,f18.8))') &
     451                 :           0 :      '             d3etot_t4 = ',d3etot_t4(1,1),  ',',d3etot_t4(2,1)
     452         [ #  # ]:           0 :    else if (n2dq==2) then
     453                 :             :      write(msg,'(2(2(a,f18.8),a))') &
     454                 :           0 :      '   d3etot_t4(dw shear) = ',d3etot_t4(1,1),  ',',d3etot_t4(2,1),ch10,&
     455                 :           0 :      '   d3etot_t4(up shear) = ',d3etot_t4(1,2),  ',',d3etot_t4(2,2)
     456                 :             :    end if
     457                 :           0 :    call wrtout(std_out,msg,'COLL')
     458                 :           0 :    call wrtout(ab_out,msg,'COLL')
     459         [ #  # ]:           0 :    if (n1dq==1) then
     460                 :             :      write(msg,'(2(a,f18.8))') &
     461                 :           0 :      '             d3etot_t5 = ',d3etot_t5(1,1),  ',',d3etot_t5(2,1)
     462         [ #  # ]:           0 :    else if (n1dq==2) then
     463                 :             :      write(msg,'(2(2(a,f18.8),a))') &
     464                 :           0 :      '   d3etot_t5(dw shear) = ',d3etot_t5(1,1),  ',',d3etot_t5(2,1),ch10,&
     465                 :           0 :      '   d3etot_t5(up shear) = ',d3etot_t5(1,2),  ',',d3etot_t5(2,2)
     466                 :             :    end if
     467                 :           0 :    call wrtout(std_out,msg,'COLL')
     468                 :           0 :    call wrtout(ab_out,msg,'COLL')
     469   [ #  #  #  #  :           0 :    if (i1pert<=natom.and.(i2pert==natom+3.or.i2pert==natom+4)) then
                   #  # ]
     470         [ #  # ]:           0 :      if (n2dq==1) then
     471                 :             :        write(msg,'(2(a,f18.8))') &
     472                 :           0 :        '          d3etot_tgeom = ',d3etot_tgeom(1,1),  ',',d3etot_tgeom(2,1)
     473         [ #  # ]:           0 :      else if (n2dq==2) then
     474                 :             :        write(msg,'(2(2(a,f18.8),a))') &
     475                 :           0 :        'd3etot_tgeom(dw shear) = ',d3etot_tgeom(1,1),  ',',d3etot_tgeom(2,1),ch10,&
     476                 :           0 :        'd3etot_tgeom(up shear) = ',d3etot_tgeom(1,2),  ',',d3etot_tgeom(2,2)
     477                 :             :      end if
     478                 :           0 :      call wrtout(std_out,msg,'COLL')
     479                 :           0 :      call wrtout(ab_out,msg,'COLL')
     480                 :             :    end if
     481                 :             :  end if
     482                 :             : 
     483         [ +  + ]:        1773 :  d3etot(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=e3tot(:)
     484                 :             : 
     485                 :             : !Deallocations
     486                 :             : 
     487                 :             :  DBG_EXIT("COLL")
     488                 :             : 
     489                 :         591 : end subroutine dfptlw_pert
     490                 :             : !!***
     491                 :             : 
     492                 :             : !!****f* ABINIT/m_dfptlw_pert/lw_elecstic
     493                 :             : !! NAME
     494                 :             : !!  lw_elecstic
     495                 :             : !!
     496                 :             : !! FUNCTION
     497                 :             : !!  This routine calculates the electrostatic term of the spatial-dispersion
     498                 :             : !!  third-order energy derivative for a couple of perturbations and a gradient
     499                 :             : !!  direction.
     500                 :             : !!
     501                 :             : !! INPUTS
     502                 :             : !!  cplex= if 1, real space 1-order functions on FFT grid are REAL,
     503                 :             : !!          if 2, COMPLEX
     504                 :             : !!  gmet(3,3)=reciprocal space metric tensor in bohr**-2
     505                 :             : !!  gprimd(3,3)=reciprocal space dimensional primitive translations
     506                 :             : !!  gsqcut=large sphere cut-off
     507                 :             : !!  i3dir= directions of the 3th perturbations
     508                 :             : !!  kxc(nfft,nkxc)=exchange and correlation kernel
     509                 :             : !!  mpi_enreg=information about MPI parallelization
     510                 :             : !!  nfft= number of FFT grid points (for this proc)
     511                 :             : !!  ngfft(1:18)=integer array with FFT box dimensions and other
     512                 :             : !!  nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed.
     513                 :             : !!  nspden = number of spin-density components
     514                 :             : !!  rho1g1(2,nfft)=G-space RF electron density in electrons/bohr**3 (i1pert)
     515                 :             : !!  rho1r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i1pert)
     516                 :             : !!  rho2r1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3 (i2pert)
     517                 :             : !!  ucvol=volume of the unit cell
     518                 :             : !!
     519                 :             : !! OUTPUT
     520                 :             : !!  d3etot_telec(2)= Electrostatic term of the third-order energy derivative
     521                 :             : !!
     522                 :             : !! SOURCE
     523                 :             : 
     524                 :         591 : subroutine lw_elecstic(cplex,d3etot_telec,gmet,gprimd,gsqcut,&
     525                 :         591 : & i3dir,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,rho1g1,rho1r1,rho2r1,ucvol)
     526                 :             : 
     527                 :             : !Arguments ------------------------------------
     528                 :             :  integer,intent(in) :: cplex,i3dir
     529                 :             :  integer,intent(in) :: nfft,nkxc,nspden
     530                 :             :  real(dp),intent(in) :: gsqcut,ucvol
     531                 :             :  type(MPI_type),intent(inout) :: mpi_enreg
     532                 :             : !arrays
     533                 :             :  integer,intent(in) :: ngfft(18)
     534                 :             :  real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
     535                 :             :  real(dp),intent(in) :: rho1g1(2,nfft),rho1r1(cplex*nfft,nspden)
     536                 :             :  real(dp),intent(in) :: rho2r1(cplex*nfft,nspden),kxc(nfft,nkxc)
     537                 :             :  real(dp),intent(out) :: d3etot_telec(2)
     538                 :             : 
     539                 :             : !Local variables-------------------------------
     540                 :             : !scalars
     541                 :             :  integer :: ii,jj,nfftot,qcar
     542                 :             :  real(dp) :: doti,dotr
     543                 :             : 
     544                 :             : !arrays
     545                 :         591 :  real(dp),allocatable :: rhor1_cplx(:,:)
     546                 :         591 :  real(dp),allocatable :: vxc1dq(:,:),vxc1dq_car(:,:,:),vqgradhart(:)
     547                 :             : 
     548                 :             : ! *************************************************************************
     549                 :             : 
     550                 :             :  DBG_ENTER("COLL")
     551                 :             : 
     552                 :             : !If GGA xc first calculate the Cartesian q gradient of the xc potential
     553         [ +  + ]:         591 :  if (nkxc == 7) then
     554   [ +  -  -  +  :         995 :    ABI_MALLOC(vxc1dq,(2*nfft,nspden))
          +  -  -  +  +  
          -  +  -  +  -  
             -  +  -  + ]
     555   [ -  +  +  -  :         995 :    ABI_MALLOC(vxc1dq_car,(2*nfft,nspden,3))
          -  +  +  -  +  
          -  +  -  +  -  
             -  +  -  + ]
     556         [ +  + ]:         796 :    do qcar=1,3
     557                 :         597 :      call dfpt_mkvxcggadq(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nkxc,nspden,qcar,rho1r1,vxc1dq)
     558   [ +  +  +  + ]:     4134967 :      vxc1dq_car(:,:,qcar)=vxc1dq(:,:)
     559                 :             :    end do
     560                 :             :  end if
     561                 :             : 
     562                 :             : !Calculate the q gradient of the Hartree potential
     563   [ +  -  -  +  :        2364 :  ABI_MALLOC(vqgradhart,(2*nfft))
          +  -  +  -  +  
             -  -  +  -  
                      + ]
     564                 :         591 :  call hartredq(2,gmet,gsqcut,mpi_enreg,nfft,ngfft,i3dir,rho1g1,vqgradhart)
     565                 :             : 
     566                 :             : !If GGA convert the gradient of xc kernel to reduced coordinates and incorporate it to the Hartree part
     567         [ +  + ]:         591 :  if (nkxc == 7) then
     568   [ +  +  +  + ]:     1378256 :    vxc1dq=zero
     569         [ +  + ]:         796 :    do qcar=1,3
     570   [ +  +  +  + ]:     4134967 :      vxc1dq(:,:)=vxc1dq(:,:) + gprimd(qcar,i3dir) * vxc1dq_car(:,:,qcar)
     571                 :             :    end do
     572         [ +  + ]:     1378057 :    vqgradhart(:)=vqgradhart(:)+vxc1dq(:,1)
     573         [ -  + ]:         199 :    ABI_FREE(vxc1dq_car)
     574                 :             :  end if
     575                 :             : 
     576                 :             : !Calculate the electrostatic energy term with the i2pert density response
     577                 :             : !I need a complex density for the dotprod_vn
     578   [ -  +  +  -  :        2364 :  ABI_MALLOC(rhor1_cplx,(2*nfft,nspden))
          -  +  +  -  +  
          -  +  -  -  +  
                   +  - ]
     579   [ +  +  +  + ]:     5568228 :  rhor1_cplx=zero
     580         [ +  + ]:     2784114 :  do ii=1,nfft
     581                 :     2783523 :    jj=ii*2
     582         [ +  + ]:     5567637 :    rhor1_cplx(jj-1,:)=rho2r1(ii,:)
     583                 :             :  end do
     584                 :             : 
     585                 :         591 :  nfftot=ngfft(1)*ngfft(2)*ngfft(3)
     586                 :         591 :  call dotprod_vn(2,rhor1_cplx,dotr,doti,nfft,nfftot,nspden,2,vqgradhart,ucvol)
     587                 :             : 
     588                 :         591 :  d3etot_telec(1)=dotr
     589                 :         591 :  d3etot_telec(2)=doti
     590                 :             : 
     591                 :             : !Deallocations
     592         [ +  + ]:         591 :  ABI_SFREE(vxc1dq)
     593         [ -  + ]:         591 :  ABI_FREE(vqgradhart)
     594         [ -  + ]:         591 :  ABI_FREE(rhor1_cplx)
     595                 :             : 
     596                 :             :  DBG_EXIT("COLL")
     597                 :             : 
     598         [ -  + ]:         591 : end subroutine lw_elecstic
     599                 :             : !!***
     600                 :             : 
     601                 :             : !!****f* ABINIT/m_dfptlw_pert/preca_ffnl
     602                 :             : !! NAME
     603                 :             : !!  preca_ffnl
     604                 :             : !!
     605                 :             : !! FUNCTION
     606                 :             : !!  Calculates the nonlocal form factors and derivatives for all the atoms
     607                 :             : !!  and k points.
     608                 :             : !!
     609                 :             : !! INPUTS
     610                 :             : !!  dimffnl= second dimension of ffnl
     611                 :             : !!  gmet(3,3)= reciprocal-space metric tensor
     612                 :             : !!  gprimd(3,3)= dimensional reciprocal space primitive translations (b^-1)
     613                 :             : !!  ider= if 1 first order derivatives of ffnl are calculated
     614                 :             : !!        if 2 first and second order derivatives of ffnl are calculated
     615                 :             : !!  idir0= variable that controls the way in which the derivatives of ffnl are
     616                 :             : !!         calculated and saved
     617                 :             : !!  kg(3,mpw)=integer coordinates of G vectors in basis sphere
     618                 :             : !!  kptns(3,nkpt)=k points in terms of reciprocal translations
     619                 :             : !!  mband= masimum number of bands
     620                 :             : !!  mkmem= maximum number of k points which can fit in core memory
     621                 :             : !!  mpi_enreg=information about MPI parallelization
     622                 :             : !!  mpw   = maximum number of planewaves in basis sphere (large number)
     623                 :             : !!  nkpt = number of k point
     624                 :             : !!  npwarr(nkpt)=array holding npw for each k point
     625                 :             : !!  nylmgr=second dimension of ylmgr
     626                 :             : !!  psps <type(pseudopotential_type)> = variables related to pseudopotentials
     627                 :             : !!  rmet(3,3)= real-space metric tensor
     628                 :             : !!  useylmgr= if 1 use the derivative of spherical harmonics
     629                 :             : !!  ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics
     630                 :             : !!  ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics
     631                 :             : !!
     632                 :             : !! OUTPUT
     633                 :             : !!  ffnl(mkmem,npw_k,dimffnl,psps%lmnmax,psps%ntypat)= Nonlocal projectors and their derivatives
     634                 :             : !!
     635                 :             : !! SIDE EFFECTS
     636                 :             : !!
     637                 :             : !! NOTES
     638                 :             : !!
     639                 :             : !! SOURCE
     640                 :             : 
     641                 :          12 : subroutine preca_ffnl(dimffnl,ffnl,gmet,gprimd,ider,idir0,kg,kptns,mband,mkmem,mpi_enreg,mpw,nkpt, &
     642                 :          12 : & npwarr,nylmgr,psps,rmet,useylmgr,ylm,ylmgr)
     643                 :             : 
     644                 :             : !Arguments ------------------------------------
     645                 :             : !scalars
     646                 :             :  integer , intent(in)  :: dimffnl,ider,idir0,mband,mkmem,mpw,nkpt,nylmgr,useylmgr
     647                 :             :  type(pseudopotential_type),intent(in) :: psps
     648                 :             :  type(MPI_type),intent(in) :: mpi_enreg
     649                 :             : !arrays
     650                 :             :  integer,intent(in) :: kg(3,mpw*mkmem)
     651                 :             :  integer,intent(in) :: npwarr(nkpt)
     652                 :             :  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),kptns(3,nkpt),rmet(3,3)
     653                 :             :  real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
     654                 :             :  real(dp),intent(in) :: ylmgr(mpw*mkmem,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)
     655                 :             :  real(dp),intent(out) :: ffnl(mkmem,mpw,dimffnl,psps%lmnmax,psps%ntypat)
     656                 :             : 
     657                 :             : !Local variables-------------------------------
     658                 :             : !scalars
     659                 :             :  integer :: ii,ikc,ikg,ikpt,ilm,nkpg,npw_k
     660                 :             :  !character(len=500) :: msg
     661                 :             : !arrays
     662                 :          12 :  integer,allocatable :: kg_k(:,:)
     663                 :             :  real(dp) :: kpt(3)
     664                 :          12 :  real(dp),allocatable :: ffnl_k(:,:,:,:),kpg_k(:,:)
     665                 :          12 :  real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:),ylmgr_k_part(:,:,:)
     666                 :             : 
     667                 :             : ! *************************************************************************
     668                 :             : 
     669                 :             :  DBG_ENTER("COLL")
     670                 :             : 
     671                 :             :  !Loop over k-points
     672                 :          12 :  ikg=0
     673                 :          12 :  ikc=0
     674         [ +  + ]:         876 :  do ikpt = 1, nkpt
     675                 :             : 
     676                 :         864 :    npw_k = npwarr(ikpt)
     677         [ -  + ]:         864 :    if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,1,mpi_enreg%me)) then
     678                 :             :      cycle ! Skip the rest of the k-point loop
     679                 :             :    end if
     680                 :         864 :    ikc= ikc + 1
     681                 :             : 
     682   [ -  +  +  -  :        2592 :    ABI_MALLOC(kg_k,(3,npw_k))
          +  -  +  -  -  
                +  -  + ]
     683   [ -  +  +  -  :        3456 :    ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     684   [ -  +  +  -  :        4320 :    ABI_MALLOC(ylmgr_k,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
          -  +  +  -  -  
          +  +  -  +  -  
          +  -  -  +  -  
                      + ]
     685                 :             : 
     686         [ +  + ]:        3456 :    kpt(:)= kptns(:,ikpt)
     687                 :             : 
     688                 :             :    !Get plane-wave vectors and related data at k
     689   [ +  +  +  + ]:      572632 :    kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
     690         [ +  - ]:         864 :    if (psps%useylm==1) then
     691         [ +  + ]:        9496 :      do ilm=1,psps%mpsang*psps%mpsang
     692         [ +  + ]:     1404265 :        ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
     693                 :             :      end do
     694         [ +  - ]:         864 :      if (useylmgr==1) then
     695         [ +  + ]:        9496 :        do ilm=1,psps%mpsang*psps%mpsang
     696         [ +  + ]:       87184 :          do ii=1,nylmgr
     697         [ +  + ]:    12639241 :            ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
     698                 :             :          end do
     699                 :             :        end do
     700                 :             :      end if
     701                 :             :    end if
     702                 :             : 
     703         [ +  + ]:         864 :    if (dimffnl==2.or.dimffnl==4) then
     704   [ -  +  +  -  :        2432 :      ABI_MALLOC(ylmgr_k_part,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
          -  +  +  -  +  
          -  +  -  -  +  
             -  +  -  + ]
     705   [ +  +  +  +  :     3527451 :      ylmgr_k_part(:,:,:)=ylmgr_k(:,1:3,:)
                   +  + ]
     706         [ +  - ]:         256 :    else if (dimffnl==10) then
     707   [ -  +  +  -  :        1024 :      ABI_MALLOC(ylmgr_k_part,(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr))
          -  +  +  -  -  
          +  +  -  -  +  
             -  +  -  + ]
     708   [ +  +  +  +  :     2071624 :      ylmgr_k_part(:,:,:)=ylmgr_k(:,:,:)
                   +  + ]
     709                 :             :    else
     710                 :           0 :      ABI_BUG('wrong size for ffnl via dimffnl!')
     711                 :             :    end if
     712                 :             : 
     713                 :             : 
     714                 :         864 :    nkpg=0
     715   [ -  +  +  -  :        1728 :    ABI_MALLOC(kpg_k,(npw_k,nkpg))
             -  +  -  + ]
     716   [ -  +  +  -  :        5184 :    ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
          -  +  +  -  -  
          +  +  -  -  +  
          +  -  +  -  +  
             -  -  +  -  
                      + ]
     717                 :             :    call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_k,psps%ffspl,gmet,gprimd,ider,idir0,&
     718                 :             :  & psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
     719                 :         864 :  & npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k_part)
     720                 :             : 
     721   [ +  +  +  +  :    12905392 :    ffnl(ikc,1:npw_k,:,:,:)=ffnl_k(1:npw_k,:,:,:)
             +  +  +  + ]
     722                 :             : 
     723         [ -  + ]:         864 :    ABI_FREE(kg_k)
     724         [ -  + ]:         864 :    ABI_FREE(ylm_k)
     725         [ -  + ]:         864 :    ABI_FREE(ylmgr_k)
     726         [ -  + ]:         864 :    ABI_FREE(ylmgr_k_part)
     727         [ -  + ]:         864 :    ABI_FREE(ffnl_k)
     728         [ -  + ]:         864 :    ABI_FREE(kpg_k)
     729                 :             : 
     730                 :             :    !Shift arrays memory
     731                 :         876 :    ikg=ikg+npw_k
     732                 :             : 
     733                 :             :  end do
     734                 :             : 
     735                 :             :  DBG_EXIT("COLL")
     736                 :             : 
     737   [ -  +  -  + ]:          12 : end subroutine preca_ffnl
     738                 :             : !!***
     739                 :             : 
     740                 :             : end module m_dfptlw_pert
     741                 :             : !!***
        

Generated by: LCOV version 2.3-1