TABLE OF CONTENTS


m_paw_dfpt/dsdr_k_paw [ Functions ]

[ Top ] [ m_paw_dfpt ] [ Functions ]

NAME

 dsdr_k_paw

FUNCTION

 compute on-site terms for forces and stresses for finite electric fields with PAW

INPUTS

  cprj_k (pawcprj_type) :: cprj for occupied bands at point k
  cprj_kb :: cprj for occupied bands at point k+b
  dtefield :: structure referring to all efield and berry's phase variables
  kdir :: integer giving direction along which overlap is computed for ket
  kfor :: integer indicating whether to compute forward (1) or backward (2)
    along kpt string
  natom :: number of atoms in cell
  typat :: typat(natom) type of each atom

OUTPUT

SIDE EFFECTS

 dsdr :: array of the on-site PAW parts of the derivatives with respect to atm
               positions and/or strains of the overlaps between Bloch states at points
               k and k+b, for the various pairs of bands

NOTES

 This routine assumes that the cprj are not explicitly ordered by
 atom type.

PARENTS

      berryphase_new

CHILDREN

SOURCE

2145  subroutine dsdr_k_paw(cprj_k,cprj_kb,dsdr,dtefield,kdir,kfor,mband,natom,ncpgr,typat)
2146 
2147 
2148 !This section has been created automatically by the script Abilint (TD).
2149 !Do not modify the following lines by hand.
2150 #undef ABI_FUNC
2151 #define ABI_FUNC 'dsdr_k_paw'
2152 !End of the abilint section
2153 
2154  implicit none
2155 
2156 !Arguments---------------------------
2157 !scalars
2158  integer,intent(in) :: kdir,kfor,mband,natom,ncpgr
2159  character(len=500) :: message
2160  type(efield_type),intent(in) :: dtefield
2161  type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%nspinor*mband)
2162  type(pawcprj_type),intent(in) :: cprj_kb(natom,dtefield%nspinor*mband)
2163 
2164 !arrays
2165  integer,intent(in) :: typat(natom)
2166  real(dp),intent(inout) :: dsdr(2,natom,ncpgr,dtefield%mband_occ,dtefield%mband_occ)
2167 
2168 !Local variables---------------------------
2169 !scalars
2170  integer :: iatom,iband,ibs,icpgr,ilmn,ispinor,itypat
2171  integer :: jband,jbs,jlmn,klmn,nspinor
2172  complex(dpc) :: cpk,cpkb,dcpk,dcpkb,cterm,paw_onsite
2173 
2174 ! *************************************************************************
2175 
2176 !initialize dsdr
2177  dsdr(:,:,:,:,:) = zero
2178 
2179 ! if 3 gradients we are in the ctocprj choice 2 case
2180 ! and the 3 gradients are due to the atomic displacements
2181 ! if 6 gradients we are in the ctocprj choice 3 case
2182 ! and the 6 gradients are due to the strains
2183 ! if 9 gradients we are in the ctocprj choice 23 case
2184 ! and the first six are due to strain, last three due to displacements
2185  if (ncpgr /= 3 .and. ncpgr /= 6 .and. ncpgr /= 9) then
2186    message = ' dsdr_k_paw called with ncpgr /= 3, 6, or 9 (no gradients) '
2187    MSG_BUG(message)
2188  end if
2189 
2190  nspinor = dtefield%nspinor
2191 
2192  do iatom = 1, natom
2193    itypat = typat(iatom)
2194 
2195    do ilmn=1,dtefield%lmn_size(itypat)
2196      do jlmn=1,dtefield%lmn_size(itypat)
2197        klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
2198        paw_onsite = cmplx(dtefield%qijb_kk(1,klmn,iatom,kdir),&
2199 &       dtefield%qijb_kk(2,klmn,iatom,kdir))
2200        if (kfor > 1) paw_onsite = conjg(paw_onsite)
2201        do iband = 1, dtefield%mband_occ
2202          do jband = 1, dtefield%mband_occ
2203            do ispinor = 1, nspinor
2204              do icpgr = 1, ncpgr
2205                ibs = nspinor*(iband-1) + ispinor
2206                jbs = nspinor*(jband-1) + ispinor
2207                cpk=cmplx(cprj_k(iatom,ibs)%cp(1,ilmn),cprj_k(iatom,ibs)%cp(2,ilmn))
2208                dcpk=cmplx(cprj_k(iatom,ibs)%dcp(1,icpgr,ilmn),cprj_k(iatom,ibs)%dcp(2,icpgr,ilmn))
2209                cpkb=cmplx(cprj_kb(iatom,jbs)%cp(1,jlmn),cprj_kb(iatom,jbs)%cp(2,jlmn))
2210                dcpkb=cmplx(cprj_kb(iatom,jbs)%dcp(1,icpgr,jlmn),cprj_kb(iatom,jbs)%dcp(2,icpgr,jlmn))
2211                cterm=paw_onsite*(conjg(dcpk)*cpkb+conjg(cpk)*dcpkb)
2212                dsdr(1,iatom,icpgr,iband,jband) = dsdr(1,iatom,icpgr,iband,jband)+real(cterm)
2213                dsdr(2,iatom,icpgr,iband,jband) = dsdr(2,iatom,icpgr,iband,jband)+aimag(cterm)
2214              end do ! end loop over icpgr
2215            end do ! end loop over ispinor
2216          end do ! end loop over jband
2217        end do ! end loop over iband
2218      end do ! end loop over ilmn
2219    end do ! end loop over jlmn
2220 
2221  end do ! end loop over atoms
2222 
2223  end subroutine    dsdr_k_paw

m_paw_dfpt/m_paw_dfpt [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_dfpt

FUNCTION

  This module contains several routines related to the 1st and 2nd order derivatives
    (in the DFPT approach) of PAW on-site quantities.

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (MT,AM,FJ,JWZ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_paw_dfpt
25 
26  use defs_basis
27  use defs_datatypes
28  use defs_abitypes
29  use m_abicore
30  use m_xmpi
31  use m_errors
32  use m_time, only : timab
33 
34  use m_pawang,       only : pawang_type
35  use m_pawrad,       only : pawrad_type
36  use m_pawtab,       only : pawtab_type
37  use m_paw_an,       only : paw_an_type
38  use m_paw_ij,       only : paw_ij_type
39  use m_pawcprj,      only : pawcprj_type
40  use m_pawdij,       only : pawdijhartree,pawdiju_euijkl
41  use m_pawrhoij,     only : pawrhoij_type,pawrhoij_free,pawrhoij_gather,pawrhoij_nullify
42  use m_pawfgrtab,    only : pawfgrtab_type, pawfgrtab_free, pawfgrtab_nullify, pawfgrtab_gather
43  use m_paw_finegrid, only : pawgylm, pawrfgd_fft, pawexpiqr
44  use m_pawxc,        only : pawxc_dfpt, pawxcm_dfpt
45  use m_paw_denpot,   only : pawdensities
46  use m_paral_atom,   only : get_my_atmtab,free_my_atmtab
47 
48  use m_atm2fft,      only : dfpt_atm2fft
49  use m_distribfft,   only : distribfft_type,init_distribfft_seq,destroy_distribfft
50  use m_geometry,     only : metric, stresssym
51 
52  use m_efield,       only : efield_type
53 
54  implicit none
55 
56  private
57 
58 !public procedures.
59  public :: pawdfptenergy ! Compute Hartree+XC PAW on-site contrib. to a 1st or 2nd-order energy
60  public :: pawgrnl       ! Compute derivatives of total energy due to NL terms (PAW Dij derivatives)
61  public :: dsdr_k_paw    ! Compute PAW on-site terms for forces/stresses for finite electric fields
62 
63 CONTAINS  !========================================================================================

m_paw_dfpt/pawdfptenergy [ Functions ]

[ Top ] [ m_paw_dfpt ] [ Functions ]

NAME

 pawdfptenergy

FUNCTION

 This routine compute the Hartree+XC PAW on-site contributions to a 1st-order or 2nd-order energy.
  These contributions are equal to:
    E_onsite=
       Int{ VHxc[n1_a^(1);nc^(1)].n1_b }
      -Int{ VHxc[tild_n1_a^(j1)+hat_n1_a^(j1);tild_n_c^(j1)].(tild_n1_b+n1_b)^(j2) }
 Some typical uses:
  A-Contribution to non-stationary expression of the 2nd-order total energy:
    In that case, n1_a^(1)[r]=n1^(j1)[r] and n1_b[r]=delta_n1^(j2)[r]
    where j1 and j2 are two given perturbations,
    and delta_n1^(j)[r] is the 1s-order density only due to change of WF overlap.
    See PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq.(80).
    E_onsite=
       Int{ VHxc[n1^(j1);nc^(j1)].delta_n1^(j2) }
      -Int{ VHxc[tild_n1^(j1)+hat_n1^(j1);tild_n_c^(j1)].delta_(tild_n1+hat_n1)^(j2) }
  B-Contribution to first-order Fermi energy:
    In that case, n1_a^(1)[r]=n1^(j1)[r] and n1_b[r]=n1[r,EFermi]
    where j1 is the current perturbation, and n1[r,EFermi] is the density at Fermi level.
    E_onsite=
       Int{ VHxc[n1^(j1);nc^(j1)].n1[r,EFermi] }
      -Int{ VHxc[tild_n1^(j1)+hat_n1^(j1);tild_n_c^(j1)].(tild_n1+hat_n1)[r,EFermi] }

INPUTS

  ipert1,ipert2=indexes of perturbations (j1) and (j2)
                if ipert2<=0, we compute a first-order energy
                if ipert2> 0, we compute a second-order energy
  ixc= choice of exchange-correlation scheme
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  ntypat=number of types of atoms in unit cell.
  nzlmopt_a= For the n1_a density:
            if -1, compute all LM-moments of the density and use non-zero LM-moments
            if  0, compute all LM-moments of the density and use all LM-moments
            if +1, compute only non-zero LM-moments of the density (stored before)
  nzlmopt_b= For the n1_b density:
            if -1, compute all LM-moments of the density and use non-zero LM-moments
            if  0, compute all LM-moments of the density and use all LM-moments
            if +1, compute only non-zero LM-moments of the density (stored before)
  paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
  paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh
                                     This corresponds to (j1) perturbation
  paw_ij1(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
                                     This corresponds to (j1) perturbation
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawprtvol=control print volume and debugging output for PAW
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij_a(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies for the (j1) perturbation
  pawrhoij_b(natom) <type(pawrhoij_type)>=
    if ipert2> 0: paw rhoij 1st-order occupancies for the (j2) perturbation corrsponding to n1_b^(j2)[r]
    if ipert2<=0: paw rhoij occupancies corresponding to n1_b[r]
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
  xclevel= XC functional level

OUTPUT

  delta_energy(2)= real and imaginary parts of contributions to non-stationary expression for the
              second derivative of the total energy

SIDE EFFECTS

    ==== if paw_an1(:)%has_vxc<2, compute 1st-order XC potentials
      paw_an1(natom)%vxc1(cplex_a*mesh_size,:,nspden) =AE 1st-order XC potential Vxc^(j1)
      paw_an1(natom)%vxct1(cplex_a*mesh_size,:,nspden)=PS 1st-order XC potential tVxc^(j1)
    ==== if paw_ij1(:)%has_dijhartree<2, compute 1st-order Dij_hartree
      paw_ij1(natom)%dijhartree(cplex_a*lmn2_size)=Hartree contribution to Dij^(j1)

PARENTS

      dfpt_nstpaw,newfermie1

CHILDREN

      free_my_atmtab,get_my_atmtab,pawdensities,pawdijhartree,pawxc_dfpt
      pawxcm_dfpt,timab,xmpi_sum

SOURCE

148 subroutine pawdfptenergy(delta_energy,ipert1,ipert2,ixc,my_natom,natom,ntypat,nzlmopt_a,nzlmopt_b,&
149 &                    paw_an0,paw_an1,paw_ij1,pawang,pawprtvol,pawrad,pawrhoij_a,pawrhoij_b,&
150 &                    pawtab,pawxcdev,xclevel, &
151 &                    mpi_atmtab,comm_atom) ! optional arguments (parallelism)
152 
153 
154 !This section has been created automatically by the script Abilint (TD).
155 !Do not modify the following lines by hand.
156 #undef ABI_FUNC
157 #define ABI_FUNC 'pawdfptenergy'
158 !End of the abilint section
159 
160  implicit none
161 
162 !Arguments ---------------------------------------------
163 !scalars
164  integer,intent(in) :: ipert1,ipert2,ixc,my_natom,natom,ntypat,nzlmopt_a,nzlmopt_b
165  integer,intent(in) :: pawprtvol,pawxcdev,xclevel
166  integer,optional,intent(in) :: comm_atom
167  type(pawang_type),intent(in) :: pawang
168 !arrays
169  integer,optional,target,intent(in) :: mpi_atmtab(:)
170  real(dp),intent(out) :: delta_energy(2)
171  type(paw_an_type),intent(in) :: paw_an0(my_natom)
172  type(paw_an_type),intent(inout) :: paw_an1(my_natom)
173  type(paw_ij_type),intent(inout) :: paw_ij1(my_natom)
174  type(pawrad_type),intent(in) :: pawrad(ntypat)
175  type(pawrhoij_type),intent(in) :: pawrhoij_a(my_natom),pawrhoij_b(my_natom)
176  type(pawtab_type),intent(in) :: pawtab(ntypat)
177 
178 !Local variables ---------------------------------------
179 !scalars
180  integer :: cplex_a,cplex_b,cplex_dijh1,cplex_diju1,cplex_vxc1,iatom,iatom_tot,ierr,irhoij,ispden,itypat,jrhoij
181  integer :: klmn,lm_size_a,lm_size_b,nlmn2_dij1,mesh_size,my_comm_atom,nspden,nspdiag
182  integer :: opt_compch,optexc,optvxc,usecore,usepawu,usetcore,usexcnhat
183  logical :: my_atmtab_allocated,paral_atom
184  real(dp) :: compch,eexc,dij_im,eexc_im,ro_im
185  character(len=500) :: msg
186 !arrays
187  integer,pointer :: my_atmtab(:)
188  logical,allocatable :: lmselect_a(:),lmselect_b(:),lmselect_tmp(:)
189  real(dp) :: dij(2),delta_energy_h(2),delta_energy_u(2),delta_energy_xc(2),ro(2),tsec(2)
190  real(dp),allocatable :: diju_im(:,:),kxc_dum(:,:,:),nhat1(:,:,:),rho1(:,:,:),trho1(:,:,:)
191 
192 ! *************************************************************************
193 
194  DBG_ENTER("COLL")
195 
196  call timab(567,1,tsec)
197 
198  usepawu=maxval(pawtab(1:ntypat)%usepawu)
199  if (.not.(ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11 &
200 & .or.ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11)) then
201    if((abs(nzlmopt_a)/=1.and.nzlmopt_a/=0).or.(abs(nzlmopt_b)/=1.and.nzlmopt_b/=0)) then
202      msg='invalid value for nzlmopt!'
203      MSG_BUG(msg)
204    end if
205    if (my_natom>0) then
206      if(paw_ij1(1)%has_dijhartree==0) then
207        msg='dijhartree must be allocated!'
208        MSG_BUG(msg)
209      end if
210      if (usepawu==1.or.usepawu==2.or.usepawu==5.or.usepawu==6) then
211        if(paw_ij1(1)%has_dijU==0) then
212          msg='dijU must be allocated!'
213          MSG_BUG(msg)
214        end if
215        if ((usepawu==1.or.usepawu==5).and.(ipert1==0.or.ipert2==0)) then
216          msg='If usepawu=1 or 5, pawdfptenergy is not implemented when ipert1=0 or ipert2=0!'
217          MSG_BUG(msg)
218        end if
219      end if
220      if(paw_an1(1)%has_vxc==0) then
221        msg='vxc1 and vxct1 must be allocated!'
222        MSG_BUG(msg)
223      end if
224      if(paw_an0(1)%has_kxc==0) then
225        msg='kxc1 must be allocated!'
226        MSG_BUG(msg)
227      end if
228      if ((ipert1<=natom.or.ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11).and.paw_an0(1)%has_kxc/=2) then
229        msg='XC kernels for ground state must be in memory!'
230        MSG_BUG(msg)
231      end if
232      if (paw_ij1(1)%cplex_rf/=paw_an1(1)%cplex) then
233        msg='paw_ij1()%cplex and paw_an1()%cplex must be equal!'
234        MSG_BUG(msg)
235      end if
236      if (pawrhoij_a(1)%cplex<paw_an1(1)%cplex.or.pawrhoij_b(1)%cplex<paw_an1(1)%cplex) then
237        msg='pawrhoij()%cplex must be >=paw_an1()%cplex!'
238        MSG_BUG(msg)
239      end if
240      if (paw_ij1(1)%nspden/=paw_an1(1)%nspden) then
241        msg='paw_ij1()%nspden and paw_an1()%nspden must be equal!'
242        MSG_BUG(msg)
243      end if
244      if (pawrhoij_a(1)%nspden/=pawrhoij_b(1)%nspden) then
245        msg='pawrhoij_a()%nspden must =pawrhoij_b()%nspden !'
246        MSG_BUG(msg)
247      end if
248    end if
249  end if
250 
251 !Set up parallelism over atoms
252  paral_atom=(present(comm_atom).and.(my_natom/=natom))
253  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
254  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
255  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
256 
257 !Init contribution to 1st-order (or 2nd-order) energy
258  delta_energy(1:2)=zero
259 
260 !For some perturbations, nothing else to do
261  if (ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11 .or. &
262 & ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11) return
263 
264 !Various inits
265  opt_compch=0;optvxc=1;optexc=3
266  usecore=0;usetcore=0  ! This is true for phonons and Efield pert.
267  usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
268  delta_energy_xc(1:2)=zero;delta_energy_h(1:2)=zero;delta_energy_u(1:2)=zero
269  dij(1:2)=zero;ro(1:2)=zero
270 
271 
272 !================ Loop on atomic sites =======================
273  do iatom=1,my_natom
274    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
275 
276    itypat=pawrhoij_a(iatom)%itypat
277    mesh_size=pawtab(itypat)%mesh_size
278    nspden=paw_an1(iatom)%nspden
279    cplex_a=pawrhoij_a(iatom)%cplex
280    cplex_b=pawrhoij_b(iatom)%cplex
281    cplex_dijh1=paw_ij1(iatom)%cplex_rf
282    cplex_diju1=paw_ij1(iatom)%cplex_rf
283    cplex_vxc1=paw_an1(iatom)%cplex
284    nlmn2_dij1=paw_ij1(iatom)%lmn2_size
285    lm_size_a=paw_an1(iatom)%lm_size
286    if (ipert2<=0) lm_size_b=paw_an0(iatom)%lm_size
287    if (ipert2> 0) lm_size_b=paw_an1(iatom)%lm_size
288 
289 !  If Vxc potentials are not in memory, compute them
290    if (paw_an1(iatom)%has_vxc/=2) then
291      ABI_ALLOCATE(rho1 ,(cplex_a*mesh_size,lm_size_a,nspden))
292      ABI_ALLOCATE(trho1,(cplex_a*mesh_size,lm_size_a,nspden))
293      ABI_ALLOCATE(nhat1,(cplex_a*mesh_size,lm_size_a,nspden*usexcnhat))
294      ABI_ALLOCATE(lmselect_a,(lm_size_a))
295      lmselect_a(:)=paw_an1(iatom)%lmselect(:)
296      ABI_ALLOCATE(lmselect_tmp,(lm_size_a))
297      lmselect_tmp(:)=.true.
298      if (nzlmopt_a==1) lmselect_tmp(:)=lmselect_a(:)
299 !    Compute on-site 1st-order densities
300      call pawdensities(compch,cplex_a,iatom_tot,lmselect_tmp,lmselect_a,&
301 &     lm_size_a,nhat1,nspden,nzlmopt_a,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
302 &     pawrad(itypat),pawrhoij_a(iatom),pawtab(itypat),rho1,trho1)
303      ABI_DEALLOCATE(lmselect_tmp)
304 !    Compute on-site 1st-order xc potentials
305      if (pawxcdev/=0) then
306        call pawxcm_dfpt(pawtab(itypat)%coredens,cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxc1,&
307 &       lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,&
308 &       pawang,pawrad(itypat),rho1,usecore,0,&
309 &       paw_an1(iatom)%vxc1,xclevel)
310        call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),&
311 &       cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxct1,&
312 &       lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,&
313 &       pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,&
314 &       paw_an1(iatom)%vxct1,xclevel)
315      else
316        call pawxc_dfpt(pawtab(itypat)%coredens,cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxc1,&
317 &       lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,&
318 &       pawang,pawrad(itypat),rho1,usecore,0,&
319 &       paw_an0(iatom)%vxc1,paw_an1(iatom)%vxc1,xclevel)
320        call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),&
321 &       cplex_a,cplex_vxc1,eexc,ixc,paw_an0(iatom)%kxct1,&
322 &       lm_size_a,lmselect_a,nhat1,paw_an0(iatom)%nkxc1,mesh_size,nspden,optvxc,&
323 &       pawang,pawrad(itypat),trho1,usetcore,2*usexcnhat,&
324 &       paw_an0(iatom)%vxct1,paw_an1(iatom)%vxct1,xclevel)
325      end if
326 
327      paw_an1(iatom)%has_vxc=2
328      ABI_DEALLOCATE(lmselect_a)
329      ABI_DEALLOCATE(rho1)
330      ABI_DEALLOCATE(trho1)
331      ABI_DEALLOCATE(nhat1)
332    end if ! has_vxc
333 
334 !  If Dij_hartree are not in memory, compute them
335    if (paw_ij1(iatom)%has_dijhartree/=2) then
336      call pawdijhartree(cplex_dijh1,paw_ij1(iatom)%dijhartree,paw_ij1(iatom)%nspden,&
337 &     pawrhoij_a(iatom),pawtab(itypat))
338      paw_ij1(iatom)%has_dijhartree=2
339    end if
340 
341 !  Compute contribution to 1st-order (or 2nd-order) energy from 1st-order XC potential
342    ABI_ALLOCATE(rho1 ,(cplex_b*mesh_size,lm_size_b,nspden))
343    ABI_ALLOCATE(trho1,(cplex_b*mesh_size,lm_size_b,nspden))
344    ABI_ALLOCATE(nhat1,(cplex_b*mesh_size,lm_size_b,nspden*usexcnhat))
345    ABI_ALLOCATE(lmselect_b,(lm_size_b))
346    if (ipert2<=0) lmselect_b(:)=paw_an0(iatom)%lmselect(:)
347    if (ipert2> 0) lmselect_b(:)=paw_an1(iatom)%lmselect(:)
348    ABI_ALLOCATE(lmselect_tmp,(lm_size_b))
349    lmselect_tmp(:)=.true.
350    if (nzlmopt_b==1) lmselect_tmp(:)=lmselect_b(:)
351 !  Compute on-site 1st-order densities
352    call pawdensities(compch,cplex_b,iatom_tot,lmselect_tmp,lmselect_b,&
353 &   lm_size_b,nhat1,nspden,nzlmopt_b,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
354 &   pawrad(itypat),pawrhoij_b(iatom),pawtab(itypat),rho1,trho1)
355    ABI_DEALLOCATE(lmselect_tmp)
356 !  Compute contributions to 1st-order (or 2nd-order) energy
357    if (pawxcdev/=0) then
358      ABI_ALLOCATE(kxc_dum,(mesh_size,pawang%angl_size,0))
359      call pawxcm_dfpt(pawtab(itypat)%coredens,cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
360 &     lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
361 &     rho1,usecore,0,paw_an1(iatom)%vxc1,xclevel,d2enxc_im=eexc_im)
362 
363      delta_energy_xc(1)=delta_energy_xc(1)+eexc
364      delta_energy_xc(2)=delta_energy_xc(2)+eexc_im
365      call pawxcm_dfpt(pawtab(itypat)%tcoredens(:,1),&
366 &     cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
367 &     lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
368 &     trho1,usetcore,2*usexcnhat,paw_an1(iatom)%vxct1,xclevel,&
369 &     d2enxc_im=eexc_im)
370      ABI_DEALLOCATE(kxc_dum)
371      delta_energy_xc(1)=delta_energy_xc(1)-eexc
372      delta_energy_xc(2)=delta_energy_xc(2)-eexc_im
373    else
374      ABI_ALLOCATE(kxc_dum,(mesh_size,lm_size_b,0))
375      call pawxc_dfpt(pawtab(itypat)%coredens,cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
376 &     lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
377 &     rho1,usecore,0,paw_an0(iatom)%vxc1,paw_an1(iatom)%vxc1,xclevel,d2enxc_im=eexc_im)
378      delta_energy_xc(1)=delta_energy_xc(1)+eexc
379      delta_energy_xc(2)=delta_energy_xc(2)+eexc_im
380      call pawxc_dfpt(pawtab(itypat)%tcoredens(:,1),&
381 &     cplex_b,cplex_vxc1,eexc,ixc,kxc_dum,&
382 &     lm_size_b,lmselect_b,nhat1,0,mesh_size,nspden,optexc,pawang,pawrad(itypat),&
383 &     trho1,usetcore,2*usexcnhat,paw_an0(iatom)%vxct1,paw_an1(iatom)%vxct1,xclevel,&
384 &     d2enxc_im=eexc_im)
385      ABI_DEALLOCATE(kxc_dum)
386      delta_energy_xc(1)=delta_energy_xc(1)-eexc
387      delta_energy_xc(2)=delta_energy_xc(2)-eexc_im
388    end if
389    ABI_DEALLOCATE(lmselect_b)
390    ABI_DEALLOCATE(rho1)
391    ABI_DEALLOCATE(trho1)
392    ABI_DEALLOCATE(nhat1)
393 
394 !  Compute contribution to 1st-order(or 2nd-order) energy from 1st-order Hartree potential
395    nspdiag=1;if (nspden==2) nspdiag=2
396    do ispden=1,nspdiag
397      if (cplex_dijh1==1) then
398        jrhoij=1
399        do irhoij=1,pawrhoij_b(iatom)%nrhoijsel
400          klmn=pawrhoij_b(iatom)%rhoijselect(irhoij)
401          dij(1)=paw_ij1(iatom)%dijhartree(klmn)
402          ro(1)=pawrhoij_b(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn)
403          delta_energy_h(1)=delta_energy_h(1)+ro(1)*dij(1)
404          if (cplex_b==2) then
405            ro(2)=pawrhoij_b(iatom)%rhoijp(jrhoij+1,ispden)*pawtab(itypat)%dltij(klmn)
406            delta_energy_h(2)=delta_energy_h(2)+ro(2)*dij(1)
407          end if
408          jrhoij=jrhoij+cplex_b
409        end do
410      else ! cplex_dijh1==2
411        jrhoij=1
412        do irhoij=1,pawrhoij_b(iatom)%nrhoijsel
413          klmn=pawrhoij_b(iatom)%rhoijselect(irhoij)
414          dij(1)=paw_ij1(iatom)%dijhartree(klmn)
415          dij(2)=paw_ij1(iatom)%dijhartree(klmn+nlmn2_dij1)
416          ro(1)=pawrhoij_b(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn)
417          delta_energy_h(1)=delta_energy_h(1)+ro(1)*dij(1)
418          delta_energy_h(2)=delta_energy_h(2)-ro(1)*dij(2)
419          if (cplex_b==2) then
420            ro(2)=pawrhoij_b(iatom)%rhoijp(jrhoij+1,ispden)*pawtab(itypat)%dltij(klmn)
421            delta_energy_h(1)=delta_energy_h(1)+ro(2)*dij(2)
422            delta_energy_h(2)=delta_energy_h(2)+ro(2)*dij(1)
423          end if
424          jrhoij=jrhoij+cplex_b
425        end do
426      end if
427    end do
428 
429 !  Compute contribution to 1st-order(or 2nd-order) energy from 1st-order PAW+U potential
430    if (usepawu==5.or.usepawu==6) then
431      ABI_CHECK(paw_ij1(iatom)%cplex_dij==1,'cplex_dij not yet available!')
432      ABI_CHECK(paw_ij1(iatom)%ndij/=4,'ndij not yet available!')
433 !    If DijU are not in memory, compute them
434      if (paw_ij1(iatom)%has_dijU/=2.or.cplex_diju1/=1) then ! We force the recomputation of dijU in when cplex=2 to get diju_im
435        if (cplex_diju1==1) then
436          call pawdiju_euijkl(cplex_diju1,paw_ij1(iatom)%cplex_dij,paw_ij1(iatom)%dijU,&
437 &             paw_ij1(iatom)%ndij,pawrhoij_a(iatom),pawtab(itypat))
438        else
439          ABI_ALLOCATE(diju_im,(pawtab(itypat)%lmn2_size,paw_ij1(iatom)%ndij))
440          call pawdiju_euijkl(cplex_diju1,paw_ij1(iatom)%cplex_dij,paw_ij1(iatom)%dijU,&
441 &             paw_ij1(iatom)%ndij,pawrhoij_a(iatom),pawtab(itypat),diju_im=diju_im)
442        end if
443        paw_ij1(iatom)%has_dijU=2
444      end if
445      if (cplex_diju1==1) then
446        do ispden=1,paw_ij1(iatom)%ndij
447          jrhoij=1
448          do irhoij=1,pawrhoij_b(iatom)%nrhoijsel
449            klmn=pawrhoij_b(iatom)%rhoijselect(irhoij)
450            dij(1)=paw_ij1(iatom)%dijU(klmn,ispden)
451            ro(1)=pawrhoij_b(iatom)%rhoijp(jrhoij,ispden)*pawtab(itypat)%dltij(klmn)
452            delta_energy_u(1)=delta_energy_u(1)+ro(1)*dij(1)
453 !          The imaginary part is not implemented (because not used)
454 !          if (cplex_b==2) then
455 !            ...
456 !          end if
457            jrhoij=jrhoij+cplex_b
458          end do
459        end do
460      else ! cplex_diju1==2
461        do ispden=1,paw_ij1(iatom)%ndij
462          jrhoij=1
463          do irhoij=1,pawrhoij_b(iatom)%nrhoijsel
464            klmn=pawrhoij_b(iatom)%rhoijselect(irhoij)
465            dij(1)=paw_ij1(iatom)%dijU(klmn,ispden)
466            dij(2)=paw_ij1(iatom)%dijU(klmn+nlmn2_dij1,ispden)
467            ro(1)=pawrhoij_b(iatom)%rhoijp(jrhoij,ispden)
468            if (cplex_b==1) then
469              delta_energy_u(1)=delta_energy_u(1)+ro(1)*dij(1)*pawtab(itypat)%dltij(klmn)
470            else
471              ro_im = pawrhoij_b(iatom)%rhoijim(klmn,ispden)
472              if (abs(pawtab(itypat)%dltij(klmn)-1)<tol14) then ! case i=j
473                delta_energy_u(1)=delta_energy_u(1)+ro(1)*dij(1)+(ro(2)+ro_im)*dij(2)
474              else ! case i/=j (factor of two for taking into account dltij)
475                dij_im = diju_im(klmn,ispden)
476                delta_energy_u(1)=delta_energy_u(1)+two*ro(1)*dij(1)           ! re  * re
477                delta_energy_u(1)=delta_energy_u(1)+two*ro_im*dij_im           ! ima * ima
478                delta_energy_u(1)=delta_energy_u(1)+two*ro(2)*(dij(2)-dij_im)  ! imb * imb
479              end if
480            end if
481            jrhoij=jrhoij+cplex_b
482          end do
483        end do
484        ABI_DEALLOCATE(diju_im)
485      end if
486    end if
487 
488 !  ================ End loop on atomic sites =======================
489  end do
490 
491 !Final building of 1st-order (or 2nd-order) energy
492  delta_energy(1:2)=delta_energy_xc(1:2)+delta_energy_h(1:2)+delta_energy_u(1:2)
493 
494 !Reduction in case of parallelism
495  if (paral_atom) then
496    call xmpi_sum(delta_energy,my_comm_atom,ierr)
497  end if
498 
499 !Destroy atom table used for parallelism
500  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
501 
502  call timab(567,2,tsec)
503 
504  DBG_EXIT("COLL")
505 
506 end subroutine pawdfptenergy

m_paw_dfpt/pawgrnl [ Functions ]

[ Top ] [ m_paw_dfpt ] [ Functions ]

NAME

 pawgrnl

FUNCTION

 PAW: Add to GRadients of total energy due to non-local term of Hamiltonian
      the contribution due to Dij derivatives
 In particular, compute contribution to forces, stresses, dyn. matrix
 Remember: Vnl=Sum_ij[|p_i>Dij<p_j|]

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  dimnhat=second dimension of array nhat (0 or # of spin components)
  distribfft<type(distribfft_type)>=--optional-- contains all the informations related
                                    to the FFT parallelism and plane sharing
  dyfr_cplex=1 if dyfrnl is real, 2 if it is complex
  gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere
  mgfft=maximum size of 1D FFTs
  me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  comm_fft=--optional-- MPI communicator over FFT components (=mpi_comm_grid is not present)
  mpi_comm_grid=--optional-- MPI communicator over real space grid components (=comm_fft is not present)
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  nattyp(ntypat)=array describing how many atoms of each type in cell
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhat(nfft,dimnhat)=compensation charge density on rectangular grid in real space
  nspden=number of spin-density components
  nsym=number of symmetries in space group
  ntypat=number of types of atoms
  optgr= 1 if gradients with respect to atomic position(s) have to be computed
  optgr2= 1 if 2nd gradients with respect to atomic position(s) have to be computed
  optstr= 1 if gradients with respect to strain(s) have to be computed
  optstr2= 1 if 2nd gradients with respect to strain(s) have to be computed
  paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  qphon(3)=wavevector of the phonon
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates
  typat(natom)=types of atoms
  ucvol=unit cell volume
  vtrial(nfft,nspden)= total local potential
  vxc(nfft,nspden)=XC potential
  xred(3,natom)=reduced dimensionless atomic coordinates

SIDE EFFECTS

  At input, this terms contain contribution from non-local projectors derivatives
  At output, they are updated with the contribution of Dij derivatives
  ==== if optgr=1 ====
   grnl(3*natom) =gradients of NL energy wrt atomic coordinates
  ==== if optstr=1 ====
   nlstr(6) =gradients of NL energy wrt strains
  ==== if optgr2=1 ====
   dyfrnl(dyfr_cplex,3,3,natom,natom) =2nd gradients of NL energy wrt atomic coordinates
  ==== if optstr=2 ====
    eltfrnl(6+3*natom,6)=non-symmetrized non-local contribution to the elastic tensor

NOTES

   In the case of parallelisation over atoms and calculation of dynamical matrix (optgr2=1)
   several data are gathered and no more distributed inside this routine.

PARENTS

      d2frnl,etotfor,forstr

CHILDREN

      pawgrnl_convert,destroy_distribfft,dfpt_atm2fft,free_my_atmtab
      get_my_atmtab,init_distribfft_seq,metric,pawexpiqr,pawfgrtab_free
      pawfgrtab_gather,pawfgrtab_nullify,pawgylm,pawrfgd_fft,pawrhoij_free
      pawrhoij_gather,pawrhoij_nullify,stresssym,xmpi_sum

SOURCE

 590 subroutine pawgrnl(atindx1,dimnhat,dyfrnl,dyfr_cplex,eltfrnl,grnl,gsqcut,mgfft,my_natom,natom,&
 591 &          nattyp,nfft,ngfft,nhat,nlstr,nspden,nsym,ntypat,optgr,optgr2,optstr,optstr2,&
 592 &          pawang,pawfgrtab,pawrhoij,pawtab,ph1d,psps,qphon,rprimd,symrec,typat,ucvol,vtrial,vxc,xred,&
 593 &          mpi_atmtab,comm_atom,comm_fft,mpi_comm_grid,me_g0,paral_kgb,distribfft) ! optional arguments (parallelism)
 594 
 595 
 596 !This section has been created automatically by the script Abilint (TD).
 597 !Do not modify the following lines by hand.
 598 #undef ABI_FUNC
 599 #define ABI_FUNC 'pawgrnl'
 600 !End of the abilint section
 601 
 602  implicit none
 603 
 604 !Arguments ------------------------------------
 605 !scalars
 606  integer,intent(in) :: dimnhat,dyfr_cplex,mgfft,my_natom,natom,nfft,nspden,nsym,ntypat
 607  integer,intent(in) :: optgr,optgr2,optstr,optstr2
 608  integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,mpi_comm_grid,paral_kgb
 609  real(dp),intent(in) :: gsqcut,ucvol
 610  type(distribfft_type),optional,target,intent(in) :: distribfft
 611  type(pawang_type),intent(in) :: pawang
 612  type(pseudopotential_type),intent(in) :: psps
 613 !arrays
 614  integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18)
 615  integer,intent(in) :: symrec(3,3,nsym),typat(natom)
 616  integer,optional,target,intent(in) :: mpi_atmtab(:)
 617  real(dp),intent(in) :: nhat(nfft,dimnhat),ph1d(2,3*(2*mgfft+1)*natom),qphon(3)
 618  real(dp),intent(in) :: rprimd(3,3),vxc(nfft,nspden),xred(3,natom)
 619  real(dp),intent(in),target :: vtrial(nfft,nspden)
 620  real(dp),intent(inout) :: dyfrnl(dyfr_cplex,3,3,natom,natom*optgr2)
 621  real(dp),intent(inout) :: eltfrnl(6+3*natom,6),grnl(3*natom*optgr)
 622  real(dp),intent(inout) :: nlstr(6*optstr)
 623  type(pawfgrtab_type),target,intent(inout) :: pawfgrtab(:)
 624  type(pawrhoij_type),target,intent(inout) ::  pawrhoij(:)
 625  type(pawtab_type),intent(in) :: pawtab(ntypat)
 626 
 627 !Local variables-------------------------------
 628 !scalars
 629  integer :: bufind,bufsiz,cplex,dimvtrial,eps_alpha,eps_beta,eps_gamma,eps_delta,iatm,iatom
 630  integer :: iatom_pawfgrtab,iatom_pawrhoij,iatom_tot,iatshft,ic,idiag,idir,ier,ilm,indx,irhoij
 631  integer :: isel,ishift_grhoij,ishift_gr,ishift2_gr,ishift_gr2,ishift_str,ishift_str2,ishift_str2is,ispden
 632  integer :: ispvtr,itypat,jatom,jatom_tot,jatm,jc,jrhoij,jtypat,klm,klmn,klmn1,ll,lm_size
 633  integer :: lm_sizej,lmax,lmin,lmn2_size,me_fft,mu,mua,mub,mushift,my_me_g0,my_comm_atom,my_comm_fft
 634  integer :: my_comm_grid,my_paral_kgb,n1,n2,n3,nfftot,nfgd,nfgd_jatom
 635  integer :: ngrad,ngrad_nondiag,ngradp,ngradp_nondiag,ngrhat,nsploop
 636  integer :: opt1,opt2,opt3,qne0,usexcnhat
 637  logical,parameter :: save_memory=.true.
 638  logical :: has_phase,my_atmtab_allocated
 639  logical :: paral_atom,paral_atom_pawfgrtab,paral_atom_pawrhoij,paral_grid
 640  real(dp) :: dlt_tmp,fact_ucvol,grhat_x,hatstr_diag,rcut_jatom,ro,ro_d,ucvol_
 641  character(len=500) :: msg
 642  type(distribfft_type),pointer :: my_distribfft
 643  type(pawfgrtab_type),pointer :: pawfgrtab_iatom,pawfgrtab_jatom
 644  type(pawrhoij_type),pointer :: pawrhoij_iatom,pawrhoij_jatom
 645 !arrays
 646  integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
 647  integer,parameter :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/)
 648  integer,parameter :: mu9(9)=(/1,2,3,4,5,6,4,5,6/)
 649  integer,allocatable :: atindx(:),atm_indx(:),mu4(:)
 650  integer,allocatable,target :: ifftsph_tmp(:)
 651  integer,ABI_CONTIGUOUS pointer :: ffti3_local(:),fftn3_distrib(:),ifft_jatom(:)
 652  integer, pointer :: my_atmtab(:)
 653  real(dp) :: gmet(3,3),gprimd(3,3),hatstr(6),rdum(1),rdum2(1),rmet(3,3),tmp(12)
 654  real(dp) :: work1(dyfr_cplex,3,3),work2(dyfr_cplex,3,3)
 655  real(dp),allocatable :: buf(:,:),buf1(:),dyfr(:,:,:,:,:),eltfr(:,:)
 656  real(dp),allocatable :: grhat_tmp(:,:),grhat_tmp2(:,:),hatgr(:)
 657  real(dp),allocatable :: prod(:,:),prodp(:,:),vloc(:),vpsp1_gr(:,:),vpsp1_str(:,:)
 658  real(dp),allocatable,target :: rfgd_tmp(:,:)
 659  real(dp),ABI_CONTIGUOUS pointer :: gylm_jatom(:,:),gylmgr_jatom(:,:,:),gylmgr2_jatom(:,:,:),expiqr_jatom(:,:)
 660  real(dp),ABI_CONTIGUOUS pointer :: rfgd_jatom(:,:),vtrial_(:,:)
 661  type(coeff2_type),allocatable :: prod_nondiag(:),prodp_nondiag(:)
 662  type(pawfgrtab_type),pointer :: pawfgrtab_(:),pawfgrtab_tot(:)
 663  type(pawrhoij_type),pointer :: pawrhoij_(:),pawrhoij_tot(:)
 664 
 665 ! *************************************************************************
 666 
 667  DBG_ENTER("COLL")
 668 
 669 !Compatibility tests
 670  qne0=0;if (qphon(1)**2+qphon(2)**2+qphon(3)**2>=1.d-15) qne0=1
 671  if (my_natom>0) then
 672    if ((optgr2==1.or.optstr2==1).and.pawrhoij(1)%ngrhoij==0) then
 673      msg='Inconsistency between variables optgr2/optstr2 and ngrhoij!'
 674      MSG_BUG(msg)
 675    end if
 676    if (pawfgrtab(1)%rfgd_allocated==0) then
 677      if ((optgr2==1.and.qne0==1).or.optstr2==1) then
 678        MSG_BUG('pawfgrtab()%rfgd array must be allocated!')
 679      end if
 680    end if
 681  end if
 682 
 683 !----------------------------------------------------------------------
 684 !Parallelism setup
 685 
 686 !Set up parallelism over atoms
 687  paral_atom=(present(comm_atom).and.(my_natom/=natom))
 688  paral_atom_pawfgrtab=(size(pawfgrtab)/=natom)
 689  paral_atom_pawrhoij=(size(pawrhoij)/=natom)
 690  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
 691  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
 692  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
 693  if (paral_atom) then
 694    ABI_ALLOCATE(atm_indx,(natom))
 695    atm_indx=-1
 696    do iatom=1,my_natom
 697      atm_indx(my_atmtab(iatom))=iatom
 698    end do
 699  end if
 700 
 701 !Set up parallelism over real space grid and/or FFT
 702  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nfftot=n1*n2*n3
 703  my_comm_grid=xmpi_comm_self;my_comm_fft=xmpi_comm_self;me_fft=0
 704  my_me_g0=1;my_paral_kgb=0;paral_grid=.false.;nullify(my_distribfft)
 705  if (present(mpi_comm_grid).or.present(comm_fft)) then
 706    if (present(mpi_comm_grid)) my_comm_grid=mpi_comm_grid
 707    if (present(comm_fft)) my_comm_fft=comm_fft
 708    if (.not.present(mpi_comm_grid)) my_comm_grid=comm_fft
 709    if (.not.present(comm_fft)) my_comm_fft=mpi_comm_grid
 710    paral_grid=(xmpi_comm_size(my_comm_grid)>1)
 711    me_fft=xmpi_comm_rank(my_comm_fft)
 712  end if
 713  if (optgr2==1.or.optstr2==1) then
 714    if (present(comm_fft)) then
 715      if ((.not.present(paral_kgb)).or.(.not.present(me_g0)).or.(.not.present(distribfft))) then
 716        MSG_BUG(' Need paral_kgb, me_g0 and distribfft with comm_fft !')
 717      end if
 718      my_me_g0=me_g0;my_paral_kgb=paral_kgb
 719      my_distribfft => distribfft
 720    else
 721      ABI_DATATYPE_ALLOCATE(my_distribfft,)
 722      call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp')
 723    end if
 724    if (n2 == my_distribfft%n2_coarse) then
 725      fftn3_distrib => my_distribfft%tab_fftdp3_distrib
 726      ffti3_local => my_distribfft%tab_fftdp3_local
 727    else
 728      fftn3_distrib => my_distribfft%tab_fftdp3dg_distrib
 729      ffti3_local => my_distribfft%tab_fftdp3dg_local
 730    end if
 731  else
 732    nullify(my_distribfft,fftn3_distrib,ffti3_local)
 733  end if
 734 
 735 !----------------------------------------------------------------------
 736 !Initializations
 737 
 738 !Compute different geometric tensors
 739 !ucvol is not computed here but provided as input arg
 740  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol_)
 741  fact_ucvol=ucvol/dble(nfftot)
 742 
 743 !Retrieve local potential according to the use of nhat in XC
 744  usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
 745  if (usexcnhat==0) then
 746    ABI_ALLOCATE(vtrial_,(nfft,1))
 747    dimvtrial=1
 748 !$OMP PARALLEL DO PRIVATE(ic) SHARED(nfft,vtrial,vtrial_,vxc)
 749    do ic=1,nfft
 750      vtrial_(ic,1)=vtrial(ic,1)-vxc(ic,1)
 751    end do
 752  else
 753    dimvtrial=nspden
 754    vtrial_ => vtrial
 755  end if
 756 
 757 !Initializations and allocations
 758  ngrhat=0;ngrad=0;ngradp=0;ngrad_nondiag=0;ngradp_nondiag=0
 759  ishift_grhoij=0;ishift_gr=0;ishift_gr2=0;ishift_str=0;ishift_str2=0;ishift_str2is=0;ishift2_gr=0
 760  cplex=1;if (qne0==1) cplex=2
 761  if (optgr==1) then
 762    ABI_ALLOCATE(hatgr,(3*natom))
 763    hatgr=zero
 764    ngrad=ngrad+3
 765    ngrhat=ngrhat+3
 766    ishift_gr2=ishift_gr2+3
 767  end if
 768  if (optgr2==1) then
 769    mu=min(dyfr_cplex,cplex)
 770    ngrad =ngrad +9
 771    ngradp=ngradp+3
 772    ngrad_nondiag = ngrad_nondiag +9*mu
 773    ngradp_nondiag= ngradp_nondiag+3*mu
 774    ngrhat= ngrhat+9*mu
 775  end if
 776  if (optstr==1) then
 777    hatstr=zero
 778    ngrad=ngrad+6
 779    ngrhat=ngrhat+6
 780    ishift_gr=ishift_gr+6
 781    ishift_gr2=ishift_gr2+6
 782    ishift_str2=ishift_str2+6
 783    ishift_str2is = ishift_str2is+6
 784  end if
 785  if (optstr2==1) then
 786    ngrad =ngrad+6*(6+3)
 787    ngradp=ngradp+(6+3)
 788    ngrad_nondiag =ngrad_nondiag+6*(6+3)
 789    ngradp_nondiag=ngradp_nondiag+3
 790    ishift2_gr=ishift2_gr+3
 791    ngrhat=ngrhat+6*(6+3)
 792    ishift_gr=ishift_gr+(6+3)
 793    ishift_gr2=ishift_gr2+6*(6+3)
 794    ishift_str2is=ishift_str2is+36
 795    ishift_grhoij = 6
 796  end if
 797 
 798  nsploop=nspden;if (dimvtrial<nspden) nsploop=2
 799  if (optgr2/=1.and.optstr2/=1) then
 800    ABI_ALLOCATE(grhat_tmp,(ngrhat,1))
 801  else
 802    ABI_ALLOCATE(grhat_tmp,(ngrhat,natom))
 803    grhat_tmp=zero
 804    ABI_DATATYPE_ALLOCATE(prod_nondiag,(natom))
 805    ABI_DATATYPE_ALLOCATE(prodp_nondiag,(natom))
 806    ABI_ALLOCATE(atindx,(natom))
 807    if(optgr2==1.or.optstr2==1)then
 808      ABI_ALLOCATE(vpsp1_gr,(cplex*nfft,3))
 809      vpsp1_gr(:,:)= zero
 810    end if
 811    if (optgr2==1) then
 812      ABI_ALLOCATE(dyfr,(dyfr_cplex,3,3,natom,natom))
 813      dyfr=zero
 814    end if
 815    if (optstr2==1) then
 816      ABI_ALLOCATE(vpsp1_str,(cplex*nfft,6))
 817      ABI_ALLOCATE(grhat_tmp2,(18,natom))
 818      ABI_ALLOCATE(eltfr,(6+3*natom,6))
 819      eltfr=zero
 820    end if
 821    ABI_ALLOCATE(mu4,(4))
 822    atindx(:)=0
 823    do iatom=1,natom
 824      iatm=0
 825      do while (atindx(iatom)==0.and.iatm<natom)
 826        iatm=iatm+1;if (atindx1(iatm)==iatom) atindx(iatom)=iatm
 827      end do
 828    end do
 829  end if
 830 
 831 !The computation of dynamical matrix and elastic tensor requires the knowledge of
 832 !g_l(r-R).Y_lm(r-R) and derivatives for all atoms
 833 !Compute them here, except memory saving is activated
 834  if ((.not.save_memory).and.(optgr2==1.or.optstr2==1)) then
 835    do jatom=1,size(pawfgrtab)
 836      jatom_tot=jatom;if (paral_atom_pawfgrtab) jatom_tot=my_atmtab(jatom)
 837      pawfgrtab_jatom => pawfgrtab(jatom)
 838      lm_sizej=pawfgrtab_jatom%l_size**2
 839      opt1=0;opt2=0;opt3=0
 840      if (pawfgrtab_jatom%gylm_allocated==0) then
 841        if (allocated(pawfgrtab_jatom%gylm))  then
 842          ABI_DEALLOCATE(pawfgrtab_jatom%gylm)
 843        end if
 844        ABI_ALLOCATE(pawfgrtab_jatom%gylm,(pawfgrtab_jatom%nfgd,lm_sizej))
 845        pawfgrtab_jatom%gylm_allocated=2;opt1=1
 846      end if
 847      if (pawfgrtab_jatom%gylmgr_allocated==0) then
 848        if (allocated(pawfgrtab_jatom%gylmgr))  then
 849          ABI_DEALLOCATE(pawfgrtab_jatom%gylmgr)
 850        end if
 851        ABI_ALLOCATE(pawfgrtab_jatom%gylmgr,(3,pawfgrtab_jatom%nfgd,lm_sizej))
 852        pawfgrtab_jatom%gylmgr_allocated=2;opt2=1
 853      end if
 854      if (opt1+opt2+opt3>0) then
 855        call pawgylm(pawfgrtab_jatom%gylm,pawfgrtab_jatom%gylmgr,&
 856 &       pawfgrtab_jatom%gylmgr2,lm_sizej,pawfgrtab_jatom%nfgd,&
 857 &       opt1,opt2,opt3,pawtab(typat(jatom_tot)),pawfgrtab_jatom%rfgd)
 858      end if
 859      if (optgr2==1.and.qne0==1) then
 860        if (pawfgrtab_jatom%expiqr_allocated==0) then
 861          if (allocated(pawfgrtab_jatom%expiqr))  then
 862            ABI_DEALLOCATE(pawfgrtab_jatom%expiqr)
 863          end if
 864          pawfgrtab_jatom%expiqr_allocated=2
 865          ABI_ALLOCATE(pawfgrtab_jatom%expiqr,(2,nfgd))
 866          call pawexpiqr(pawfgrtab_jatom%expiqr,gprimd,pawfgrtab_jatom%nfgd,&
 867 &         qphon,pawfgrtab_jatom%rfgd,xred(:,jatom_tot))
 868        end if
 869      end if
 870    end do
 871  end if
 872 
 873 !The computation of dynamical matrix and elastic tensor might require some communications
 874  if ((optgr2==1.or.optstr2==1).and.paral_atom.and.paral_atom_pawfgrtab.and.(.not.save_memory)) then
 875    ABI_DATATYPE_ALLOCATE(pawfgrtab_tot,(natom))
 876    call pawfgrtab_nullify(pawfgrtab_tot)
 877    call pawfgrtab_gather(pawfgrtab,pawfgrtab_tot,my_comm_atom,ier,mpi_atmtab=my_atmtab)
 878  else
 879    pawfgrtab_tot => pawfgrtab
 880  end if
 881  if ((optgr2==1.or.optstr2==1).and.paral_atom.and.paral_atom_pawrhoij) then
 882    ABI_DATATYPE_ALLOCATE(pawrhoij_tot,(natom))
 883    call pawrhoij_nullify(pawrhoij_tot)
 884    call pawrhoij_gather(pawrhoij,pawrhoij_tot,-1,my_comm_atom, &
 885 &   with_rhoijres=.false.,with_rhoij_=.false.,with_lmnmix=.false.)
 886  else
 887    pawrhoij_tot => pawrhoij
 888  end if
 889 
 890  if (save_memory) then
 891    pawfgrtab_ => pawfgrtab
 892    pawrhoij_  => pawrhoij
 893  else
 894    pawfgrtab_ => pawfgrtab_tot
 895    pawrhoij_  => pawrhoij_tot
 896  end if
 897 
 898 !----------------------------------------------------------------------
 899 !Loops over types and atoms
 900 
 901  iatshft=0
 902  do itypat=1,ntypat
 903 
 904    lmn2_size=pawtab(itypat)%lmn2_size
 905    lm_size=pawtab(itypat)%lcut_size**2
 906 
 907    do iatm=iatshft+1,iatshft+nattyp(itypat)
 908 
 909      iatom_tot=atindx1(iatm)
 910      iatom=iatom_tot
 911      if (paral_atom) then
 912        if (save_memory.or.(optgr2/=1.and.optstr2/=1)) iatom=atm_indx(iatom_tot)
 913      end if
 914 
 915      if (iatom==-1) cycle
 916      iatom_pawfgrtab=iatom_tot;if (paral_atom_pawfgrtab) iatom_pawfgrtab=iatom
 917      iatom_pawrhoij =iatom_tot;if (paral_atom_pawrhoij)  iatom_pawrhoij =iatom
 918      pawfgrtab_iatom => pawfgrtab_(iatom_pawfgrtab)
 919      pawrhoij_iatom  => pawrhoij_(iatom_pawrhoij)
 920 
 921      idiag=1;if (optgr2==1.or.optstr2==1) idiag=iatm
 922      nfgd=pawfgrtab_iatom%nfgd
 923 
 924      ABI_ALLOCATE(vloc,(nfgd))
 925      if (ngrad>0)  then
 926        ABI_ALLOCATE(prod,(ngrad,lm_size))
 927      end if
 928      if (ngradp>0)  then
 929        ABI_ALLOCATE(prodp,(ngradp,lm_size))
 930      end if
 931      if (ngrad_nondiag>0.and.ngradp_nondiag>0) then
 932        do jatm=1,natom
 933          jtypat=typat(atindx1(jatm))
 934          lm_sizej=pawtab(jtypat)%lcut_size**2
 935          ABI_ALLOCATE(prod_nondiag(jatm)%value,(ngrad_nondiag,lm_sizej))
 936          ABI_ALLOCATE(prodp_nondiag(jatm)%value,(ngradp_nondiag,lm_sizej))
 937        end do
 938      end if
 939 
 940      grhat_tmp=zero
 941      if(optstr2==1) grhat_tmp2=zero
 942 
 943 !    ------------------------------------------------------------------
 944 !    Compute some useful data
 945 
 946 !    Eventually compute g_l(r).Y_lm(r) derivatives for the current atom (if not already done)
 947      if ((optgr==1.or.optstr==1).and.(optgr2/=1).and.(optstr2/=1)) then
 948        if (pawfgrtab_iatom%gylmgr_allocated==0) then
 949          if (allocated(pawfgrtab_iatom%gylmgr))  then
 950            ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr)
 951          end if
 952          ABI_ALLOCATE(pawfgrtab_iatom%gylmgr,(3,pawfgrtab_iatom%nfgd,lm_size))
 953          pawfgrtab_iatom%gylmgr_allocated=2
 954          call pawgylm(rdum,pawfgrtab_iatom%gylmgr,rdum2,lm_size,pawfgrtab_iatom%nfgd,&
 955 &         0,1,0,pawtab(itypat),pawfgrtab_iatom%rfgd)
 956        end if
 957 
 958      end if
 959      if (optgr2==1.or.optstr2==1) then
 960        opt1=0;opt2=0;opt3=0
 961        if (pawfgrtab_iatom%gylm_allocated==0) then
 962          if (allocated(pawfgrtab_iatom%gylm))  then
 963            ABI_DEALLOCATE(pawfgrtab_iatom%gylm)
 964          end if
 965          ABI_ALLOCATE(pawfgrtab_iatom%gylm,(pawfgrtab_iatom%nfgd,lm_size))
 966          pawfgrtab_iatom%gylm_allocated=2;opt1=1
 967        end if
 968        if (pawfgrtab_iatom%gylmgr_allocated==0) then
 969          if (allocated(pawfgrtab_iatom%gylmgr))  then
 970            ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr)
 971          end if
 972          ABI_ALLOCATE(pawfgrtab_iatom%gylmgr,(3,pawfgrtab_iatom%nfgd,lm_size))
 973          pawfgrtab_iatom%gylmgr_allocated=2;opt2=1
 974        end if
 975        if (pawfgrtab_iatom%gylmgr2_allocated==0) then
 976          if (allocated(pawfgrtab_iatom%gylmgr2))  then
 977            ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr2)
 978          end if
 979          ABI_ALLOCATE(pawfgrtab_iatom%gylmgr2,(6,pawfgrtab_iatom%nfgd,lm_size))
 980          pawfgrtab_iatom%gylmgr2_allocated=2;opt3=1
 981        end if
 982        if (opt1+opt2+opt3>0) then
 983          call pawgylm(pawfgrtab_iatom%gylm,pawfgrtab_iatom%gylmgr,&
 984 &         pawfgrtab_iatom%gylmgr2,lm_size,pawfgrtab_iatom%nfgd,&
 985 &         opt1,opt2,opt3,pawtab(itypat),pawfgrtab_iatom%rfgd)
 986        end if
 987      end if
 988 
 989 !    Eventually compute exp(-i.q.r) factors for the current atom (if not already done)
 990      if (optgr2==1.and.qne0==1.and.(pawfgrtab_iatom%expiqr_allocated==0)) then
 991        if (allocated(pawfgrtab_iatom%expiqr))  then
 992          ABI_DEALLOCATE(pawfgrtab_iatom%expiqr)
 993        end if
 994        ABI_ALLOCATE(pawfgrtab_iatom%expiqr,(2,nfgd))
 995        call pawexpiqr(pawfgrtab_iatom%expiqr,gprimd,nfgd,qphon,&
 996 &       pawfgrtab_iatom%rfgd,xred(:,iatom))
 997        pawfgrtab_iatom%expiqr_allocated=2
 998      end if
 999      has_phase=(optgr2==1.and.pawfgrtab_iatom%expiqr_allocated/=0)
1000 
1001 !    Eventually compute 1st-order potential
1002      if (optgr2==1.or.optstr2==1) then
1003        call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,iatom_tot,&
1004 &       mgfft,psps%mqgrid_vl,natom,3,nfft,ngfft,ntypat,ph1d,&
1005 &       psps%qgrid_vl,qphon,typat,ucvol,psps%usepaw,xred,psps,pawtab,atmvlocr1=vpsp1_gr,&
1006 &       vspl=psps%vlspl,comm_fft=my_comm_fft,me_g0=my_me_g0,&
1007 &       paral_kgb=my_paral_kgb,distribfft=my_distribfft)
1008        if (cplex==1) then
1009          do ic=1,nfft
1010            tmp(1:3)=vpsp1_gr(ic,1:3)
1011            do mu=1,3
1012              vpsp1_gr(ic,mu)=-(gprimd(mu,1)*tmp(1)+gprimd(mu,2)*tmp(2)+gprimd(mu,3)*tmp(3))
1013            end do
1014          end do
1015        else ! cplex=2
1016          do ic=1,nfft
1017            jc=2*ic;tmp(1:3)=vpsp1_gr(jc-1,1:3);tmp(4:6)=vpsp1_gr(jc,1:3)
1018            do mu=1,3
1019              vpsp1_gr(jc-1,mu)=-(gprimd(mu,1)*tmp(1)+gprimd(mu,2)*tmp(2)+gprimd(mu,3)*tmp(3))
1020              vpsp1_gr(jc  ,mu)=-(gprimd(mu,1)*tmp(4)+gprimd(mu,2)*tmp(5)+gprimd(mu,3)*tmp(6))
1021            end do
1022          end do
1023        end if
1024      end if
1025      if (optstr2==1) then
1026        vpsp1_str(:,:) = zero
1027        call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,natom+3,&
1028 &       mgfft,psps%mqgrid_vl,natom,6,nfft,ngfft,ntypat,&
1029 &       ph1d,psps%qgrid_vl,qphon,typat,ucvol,psps%usepaw,xred,psps,pawtab,atmvlocr1=vpsp1_str,&
1030 &       vspl=psps%vlspl,comm_fft=my_comm_fft,me_g0=my_me_g0,&
1031 &       paral_kgb=my_paral_kgb,distribfft=my_distribfft)
1032      end if
1033 
1034 !    ------------------------------------------------------------------
1035 !    Loop over spin components
1036 
1037      do ispden=1,nsploop
1038 
1039 !      ----- Retrieve potential (subtle if nspden=4 ;-)
1040        if (nspden/=4) then
1041          ispvtr=min(dimvtrial,ispden)
1042          do ic=1,nfgd
1043            jc = pawfgrtab_iatom%ifftsph(ic)
1044            vloc(ic)=vtrial_(jc,ispvtr)
1045          end do
1046        else
1047          if (ispden==1) then
1048            ispvtr=min(dimvtrial,2)
1049            do ic=1,nfgd
1050              jc=pawfgrtab_iatom%ifftsph(ic)
1051              vloc(ic)=half*(vtrial_(jc,1)+vtrial_(jc,ispvtr))
1052            end do
1053          else if (ispden==4) then
1054            ispvtr=min(dimvtrial,2)
1055            do ic=1,nfgd
1056              jc=pawfgrtab_iatom%ifftsph(ic)
1057              vloc(ic)=half*(vtrial_(jc,1)-vtrial_(jc,ispvtr))
1058            end do
1059          else if (ispden==2) then
1060            ispvtr=min(dimvtrial,3)
1061            do ic=1,nfgd
1062              jc=pawfgrtab_iatom%ifftsph(ic)
1063              vloc(ic)=vtrial_(jc,ispvtr)
1064            end do
1065          else ! ispden=3
1066            ispvtr=min(dimvtrial,4)
1067            do ic=1,nfgd
1068              jc=pawfgrtab_iatom%ifftsph(ic)
1069              vloc(ic)=-vtrial_(jc,ispvtr)
1070            end do
1071          end if
1072        end if
1073 
1074 !      -----------------------------------------------------------------------
1075 !      ----- Compute projected scalars (integrals of vloc and Q_ij^hat) ------
1076 !      ----- and/or their derivatives ----------------------------------------
1077 
1078        if (ngrad>0) prod=zero
1079        if (ngradp>0) prodp=zero
1080 
1081 !      ==== Contribution to forces ====
1082        if (optgr==1) then
1083          do ilm=1,lm_size
1084            do ic=1,pawfgrtab_iatom%nfgd
1085              do mu=1,3
1086                prod(mu+ishift_gr,ilm)=prod(mu+ishift_gr,ilm)-&
1087 &               vloc(ic)*pawfgrtab_iatom%gylmgr(mu,ic,ilm)
1088              end do
1089            end do
1090          end do
1091        end if ! optgr
1092 
1093 !      ==== Contribution to stresses ====
1094        if (optstr==1) then
1095          do ilm=1,lm_size
1096            do ic=1,pawfgrtab_iatom%nfgd
1097              jc=pawfgrtab_iatom%ifftsph(ic)
1098              do mu=1,6
1099                mua=alpha(mu);mub=beta(mu)
1100                prod(mu+ishift_str,ilm)=prod(mu+ishift_str,ilm) &
1101 &               +half*vloc(ic)&
1102 &               *(pawfgrtab_iatom%gylmgr(mua,ic,ilm)*pawfgrtab_iatom%rfgd(mub,ic)&
1103 &               +pawfgrtab_iatom%gylmgr(mub,ic,ilm)*pawfgrtab_iatom%rfgd(mua,ic))
1104              end do
1105            end do
1106          end do
1107        end if ! optstr
1108 
1109 !      ==== Diagonal contribution to frozen wf part of dyn. matrix ====
1110        if (optgr2==1) then
1111 !        Diagonal contribution
1112          do ilm=1,lm_size
1113            do ic=1,pawfgrtab_iatom%nfgd
1114              do mu=1,9
1115                prod(ishift_gr2+mu,ilm)=prod(ishift_gr2+mu,ilm) &
1116 &               +half*vloc(ic)*pawfgrtab_iatom%gylmgr2(mu9(mu),ic,ilm)
1117              end do
1118              do mu=1,3
1119                prodp(ishift_gr+mu,ilm)=prodp(ishift_gr+mu,ilm) &
1120 &               -vloc(ic)*pawfgrtab_iatom%gylmgr(mu,ic,ilm)
1121              end do
1122            end do
1123          end do
1124        end if ! optgr2
1125 
1126 !      ==== Diagonal contribution to elastic tensor ====
1127        if (optstr2==1) then
1128          do ilm=1,lm_size
1129            do ic=1,pawfgrtab_iatom%nfgd
1130              mu=1
1131              jc=pawfgrtab_iatom%ifftsph(ic)
1132              do mua=1,6
1133                eps_alpha=eps1(mua);eps_beta=eps2(mua);
1134                do mub=1,6
1135                  eps_gamma=eps1(mub);eps_delta=eps2(mub);
1136                  mu4 = zero
1137                  call pawgrnl_convert(mu4,eps_alpha,eps_beta,eps_gamma,eps_delta)
1138 !                v_loc*d2glylm
1139                  prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) + half*half*vloc(ic)*( &
1140 &                 pawfgrtab_iatom%rfgd(eps_beta,ic)*pawfgrtab_iatom%rfgd(eps_gamma,ic)*&
1141 &                 pawfgrtab_iatom%gylmgr2(mu4(2),ic,ilm)&
1142 &                 +pawfgrtab_iatom%rfgd(eps_alpha,ic)*pawfgrtab_iatom%rfgd(eps_gamma,ic)*&
1143                  pawfgrtab_iatom%gylmgr2(mu4(4),ic,ilm)&
1144 &                 +pawfgrtab_iatom%rfgd(eps_beta,ic) *pawfgrtab_iatom%rfgd(eps_delta,ic)*&
1145                  pawfgrtab_iatom%gylmgr2(mu4(1),ic,ilm)&
1146 &                 +pawfgrtab_iatom%rfgd(eps_alpha,ic)*pawfgrtab_iatom%rfgd(eps_delta,ic)*&
1147                  pawfgrtab_iatom%gylmgr2(mu4(3),ic,ilm))
1148                  if(eps_gamma==eps_beta)then
1149                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1150 &                   +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic))
1151                  end if
1152                  if(eps_gamma==eps_alpha)then
1153                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1154 &                   +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1155                  end if
1156                  if(eps_delta==eps_beta)then
1157                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1158 &                   +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic))
1159                  end if
1160                  if(eps_delta==eps_alpha)then
1161                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1162 &                   +half*half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1163                  end if
1164 !                d(vloc)/d(eps_gammadelta) * d(gylm)/d(eps_alphabeta)
1165                  prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm)&
1166 &                 +vpsp1_str(jc,mub)*half*(&
1167 &                 (pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic)&
1168 &                 +pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm) *pawfgrtab_iatom%rfgd(eps_alpha,ic)))
1169 !                d(vloc)/d(eps_alphabeta)  * d(gylm)/d(eps_gammadelta)
1170                  prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm)&
1171 &                 +vpsp1_str(jc,mua)*half*(&
1172 &                 (pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_delta,ic)&
1173 &                 +pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_gamma,ic)))
1174 !                delta_alphabeta * dv_loc/depsgammadelta * (gylm)
1175                  if (mua<=3) then
1176                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1177 &                   +vpsp1_str(jc,mub)*pawfgrtab_iatom%gylm(ic,ilm)
1178                  end if
1179 !                delta_gammadelta * dv_loc/depsalphabeta * (gylm)
1180                  if (mub<=3) then
1181                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1182 &                   +vpsp1_str(jc,mua) * pawfgrtab_iatom%gylm(ic,ilm)
1183                  end if
1184 !                delta_gammadelta * v_loc * d(gylm)/d(eps_alphabeta)
1185                  if (mub<=3) then
1186                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1187 &                   +half*vloc(ic)&
1188 &                   *(pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)&
1189 &                   + pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1190                  end if
1191 !                delta_alphabeta * v_loc * d(gylm)/d(eps_gammadelta)
1192                  if (mua<=3) then
1193                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1194 &                   +half*vloc(ic)&
1195 &                   *(pawfgrtab_iatom%gylmgr(eps_gamma,ic,ilm)*pawfgrtab_iatom%rfgd(eps_delta,ic)&
1196 &                   + pawfgrtab_iatom%gylmgr(eps_delta,ic,ilm)*pawfgrtab_iatom%rfgd(eps_gamma,ic))
1197                  end if
1198 !                delta_gammadelta delta_alphabeta * v_loc * (gylm)
1199                  if (mua<=3.and.mub<=3) then
1200                    prod(ishift_str2+mu,ilm)=prod(ishift_str2+mu,ilm) &
1201 &                   +vloc(ic)*pawfgrtab_iatom%gylm(ic,ilm)
1202                  end if
1203                  mu=mu+1
1204                end do !end loop mub
1205              end do !end loop mua
1206 !            vloc * d(gylm)/d(eps_alphabeta)
1207              do mu=1,6
1208                mua=alpha(mu);mub=beta(mu)
1209                prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)&
1210 &               +half*vloc(ic)*&
1211 &               (pawfgrtab_iatom%gylmgr(mua,ic,ilm)*pawfgrtab_iatom%rfgd(mub,ic)&
1212 &               +pawfgrtab_iatom%gylmgr(mub,ic,ilm)*pawfgrtab_iatom%rfgd(mua,ic))
1213 !              d(vloc)/d(eps_alphabeta or gammadelta) * gylm
1214                prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)&
1215 &               +vpsp1_str(jc,mu)*pawfgrtab_iatom%gylm(ic,ilm)
1216 !              delta_alphabeta * vloc * gylm
1217                if (mu<=3) then
1218                  prodp(ishift_str2+mu,ilm)=prodp(ishift_str2+mu,ilm)&
1219 &                 +vloc(ic)*pawfgrtab_iatom%gylm(ic,ilm)
1220                end if
1221 
1222 !              INTERNAL STRAIN CONTRIBUTION:
1223                do idir=1,3
1224 !                v_loc*d2glylm/dR contribution:
1225                  eps_alpha=alpha(mu);eps_beta=beta(mu);
1226                  call pawgrnl_convert(mu4,eps_alpha,eps_beta,idir,idir)
1227                  prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)&
1228 &                 -half*vloc(ic)&
1229 &                 *(pawfgrtab_iatom%gylmgr2(mu4(3),ic,ilm)*pawfgrtab_iatom%rfgd(eps_alpha,ic)&
1230 &                 +pawfgrtab_iatom%gylmgr2(mu4(1),ic,ilm)*pawfgrtab_iatom%rfgd(eps_beta,ic))
1231                  if (idir==eps_beta)then
1232                    prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)&
1233 &                   -half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_alpha,ic,ilm))
1234                  end if
1235                  if (idir==eps_alpha)then
1236                    prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)&
1237 &                   -half*vloc(ic)*(pawfgrtab_iatom%gylmgr(eps_beta,ic,ilm))
1238                  end if
1239 !                delta_gammadelta * v_loc * d(gylm)/dR
1240                  if (mu<=3) then
1241                    prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)-&
1242                    vloc(ic)*pawfgrtab_iatom%gylmgr(idir,ic,ilm)
1243                  end if
1244 !                dv_loc/deps_alph_beta * d(gylm)/dR
1245                  prod(ishift_str2is+(mu-1)*3+idir,ilm)=prod(ishift_str2is+(mu-1)*3+idir,ilm)-&
1246                  vpsp1_str(jc,mu)*pawfgrtab_iatom%gylmgr(idir,ic,ilm)
1247                end do
1248              end do
1249              do idir=1,3
1250 !              v_loc * d(gylm)/dR
1251                prodp(6+idir,ilm) = prodp(6+idir,ilm)-vloc(ic)*pawfgrtab_iatom%gylmgr(idir,ic,ilm)
1252              end do !end loop idir
1253 !            END INTERNAL STRAIN CONTRIBUTION
1254 
1255            end do
1256          end do
1257        end if !optstr2
1258 
1259 !      Off-diagonal contributions
1260        if (optgr2==1.or.optstr2==1) then
1261          do jatm=1,natom
1262            jatom_tot=atindx1(jatm);jtypat=typat(jatom_tot)
1263            jatom=jatom_tot;if (paral_atom.and.save_memory) jatom=atm_indx(jatom_tot)
1264            lm_sizej=pawtab(jtypat)%lcut_size**2
1265 
1266 !          Retrieve data for the atom j
1267            if (save_memory.and.jatom/=iatom) then
1268              rcut_jatom=pawtab(jtypat)%rshp
1269              call pawrfgd_fft(ifftsph_tmp,gmet,n1,n2,n3,nfgd_jatom,rcut_jatom,rfgd_tmp,rprimd,&
1270 &             ucvol,xred(:,jatom_tot),fft_distrib=fftn3_distrib,fft_index=ffti3_local,me_fft=me_fft)
1271              ifft_jatom => ifftsph_tmp ; rfgd_jatom => rfgd_tmp
1272              ABI_ALLOCATE(gylm_jatom,(nfgd_jatom,lm_sizej))
1273              ABI_ALLOCATE(gylmgr_jatom,(3,nfgd_jatom,lm_sizej))
1274              opt1=1;opt2=1;opt3=0;gylmgr2_jatom=>gylmgr_jatom
1275              call pawgylm(gylm_jatom,gylmgr_jatom,gylmgr2_jatom,lm_sizej,nfgd_jatom,&
1276 &             opt1,opt2,opt3,pawtab(typat(jatom_tot)),rfgd_jatom)
1277              if (optgr2==1.and.qne0==1) then
1278                ABI_ALLOCATE(expiqr_jatom,(2,nfgd_jatom))
1279                call pawexpiqr(expiqr_jatom,gprimd,nfgd_jatom,qphon,rfgd_jatom,xred(:,jatom_tot))
1280              end if
1281            else
1282              pawfgrtab_jatom => pawfgrtab_tot(jatom)
1283              nfgd_jatom      =  pawfgrtab_jatom%nfgd
1284              ifft_jatom      => pawfgrtab_jatom%ifftsph
1285              rfgd_jatom      => pawfgrtab_jatom%rfgd
1286              gylm_jatom      => pawfgrtab_jatom%gylm
1287              gylmgr_jatom    => pawfgrtab_jatom%gylmgr
1288              gylmgr2_jatom   => pawfgrtab_jatom%gylmgr2
1289              expiqr_jatom    => pawfgrtab_jatom%expiqr
1290            end if
1291 
1292 !          ==== Off-diagonal contribution to frozen wf part of dyn. matrix ====
1293            if (optgr2==1) then
1294              mu = min(dyfr_cplex,cplex)
1295              prod_nondiag(jatm)%value(ishift_gr2+1:ishift_gr2+(9*mu),:) = zero
1296              prodp_nondiag(jatm)%value(ishift2_gr+1:ishift2_gr +(3*mu),:) = zero
1297              if (has_phase.or.cplex==2) then
1298                if (dyfr_cplex==1.or.cplex==1) then
1299                  do ilm=1,lm_sizej
1300                    do ic=1,nfgd_jatom
1301                      jc=2*ifft_jatom(ic)
1302                      tmp(1:3)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(1,ic) &
1303 &                     -vpsp1_gr(jc  ,1:3)*expiqr_jatom(2,ic)
1304                      do mu=1,9
1305                        mua=alpha(mu);mub=beta(mu)
1306                        prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) &
1307 &                       +tmp(mua)*gylmgr_jatom(mub,ic,ilm)
1308                      end do
1309                      do mu=1,3
1310                        prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) &
1311 &                       -tmp(mu)*gylm_jatom(ic,ilm)
1312                      end do
1313                    end do
1314                  end do
1315                else
1316                  do ilm=1,lm_sizej
1317                    do ic=1,nfgd_jatom
1318                      jc=2*ifft_jatom(ic)
1319                      tmp(1:3)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(1,ic) &
1320 &                     -vpsp1_gr(jc  ,1:3)*expiqr_jatom(2,ic)
1321                      tmp(4:6)=vpsp1_gr(jc-1,1:3)*expiqr_jatom(2,ic) &
1322 &                     +vpsp1_gr(jc  ,1:3)*expiqr_jatom(1,ic)
1323                      do mu=1,9
1324                        mua=alpha(mu);mub=beta(mu)
1325                        prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) &
1326 &                       +tmp(mua  )*gylmgr_jatom(mub,ic,ilm)
1327                        prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm) &
1328 &                       +tmp(3+mua)*gylmgr_jatom(mub,ic,ilm)
1329                      end do
1330                      do mu=1,3
1331                        prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) &
1332 &                       -tmp(  mu)*gylm_jatom(ic,ilm)
1333                        prodp_nondiag(jatm)%value(ishift2_gr+3+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+3+mu,ilm) &
1334 &                       -tmp(3+mu)*gylm_jatom(ic,ilm)
1335                      end do
1336                    end do
1337                  end do
1338                end if
1339              else ! no phase
1340                do ilm=1,lm_sizej
1341                  do ic=1,nfgd_jatom
1342                    jc=ifft_jatom(ic)
1343                    do mu=1,9
1344                      mua=alpha(mu);mub=beta(mu)
1345                      prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)=prod_nondiag(jatm)%value(ishift_gr2+mu,ilm) &
1346 &                     +vpsp1_gr(jc,mua)*gylmgr_jatom(mub,ic,ilm)
1347                    end do
1348                    do mu=1,3
1349                      prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm)=prodp_nondiag(jatm)%value(ishift2_gr+mu,ilm) &
1350 &                     -vpsp1_gr(jc,mu)*gylm_jatom(ic,ilm)
1351                    end do
1352                  end do
1353                end do
1354              end if
1355            end if ! optgr2
1356 
1357 !          ==== Off-diagonal contribution to elastic tensor ====
1358            if (optstr2==1) then
1359              prod_nondiag(jatm)%value(ishift_str2is+1:ishift_str2is+18,:)=zero
1360              prodp_nondiag(jatm)%value(1:3,:)=zero
1361              do ilm=1,lm_sizej
1362                do ic=1,nfgd_jatom
1363                  mu=1;jc=ifft_jatom(ic)
1364 !                INTERNAL STRAIN CONTRIBUTION:
1365                  do mua=1,6
1366                    eps_alpha=eps1(mua);eps_beta=eps2(mua);
1367 !                  d(-vloc)/dR * d(gylm)/d(eps_gamma_delta)
1368                    do idir=1,3
1369                      prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)=&
1370 &                     prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)&
1371 &                     -vpsp1_gr(jc,idir)*half&
1372 &                     *(gylmgr_jatom(eps_alpha,ic,ilm) * rfgd_jatom(eps_beta ,ic)&
1373 &                     +gylmgr_jatom(eps_beta ,ic,ilm) * rfgd_jatom(eps_alpha,ic))
1374 !                    delta_alphabeta * d(-v_loc/dr) * gylm
1375                      if (mua<=3) then
1376                        prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)=&
1377 &                       prod_nondiag(jatm)%value(ishift_str2is+(mua-1)*3+idir,ilm)&
1378 &                       -vpsp1_gr(jc,idir)*gylm_jatom(ic,ilm)
1379                      end if
1380                    end do ! dir
1381                  end do ! mua
1382                  do idir=1,3
1383 !                  d(-v_loc/dr) * gylm
1384                    prodp_nondiag(jatm)%value(idir,ilm) = prodp_nondiag(jatm)%value(idir,ilm)&
1385 &                   -vpsp1_gr(jc,idir)*gylm_jatom(ic,ilm)
1386                  end do !end loop idir
1387 !                END INTERNAL STRAIN CONTRIBUTION
1388                end do
1389              end do
1390            end if ! optstr2
1391 
1392 !          Release temp memory allocated for atom j
1393            if (save_memory.and.jatom/=iatom) then
1394              ABI_DEALLOCATE(ifftsph_tmp)
1395              ABI_DEALLOCATE(rfgd_tmp)
1396              ABI_DEALLOCATE(gylm_jatom)
1397              ABI_DEALLOCATE(gylmgr_jatom)
1398              if (optgr2==1.and.qne0==1) then
1399                ABI_DEALLOCATE(expiqr_jatom)
1400              end if
1401            end if
1402 
1403          end do ! loop on atoms j
1404        end if ! optgr2 or optstr2
1405 
1406 !      --- Apply scaling factor on integrals ---
1407        if (ngrad >0) prod (:,:)=prod (:,:)*fact_ucvol
1408        if (ngradp>0) prodp(:,:)=prodp(:,:)*fact_ucvol
1409        if (ngrad_nondiag>0) then
1410          do jatm=1,natom
1411            prod_nondiag(jatm)%value(:,:)=prod_nondiag(jatm)%value(:,:)*fact_ucvol
1412          end do
1413        end if
1414        if (ngradp_nondiag>0) then
1415          do jatm=1,natom
1416            prodp_nondiag(jatm)%value(:,:)=prodp_nondiag(jatm)%value(:,:)*fact_ucvol
1417          end do
1418        end if
1419 
1420 !      --- Reduction in case of parallelization ---
1421        if (paral_grid) then
1422          if (ngrad>0) then
1423            call xmpi_sum(prod,my_comm_grid,ier)
1424          end if
1425          if (ngradp>0) then
1426            call xmpi_sum(prodp,my_comm_grid,ier)
1427          end if
1428          if (ngrad_nondiag>0.or.ngradp_nondiag>0) then
1429            bufsiz=0;bufind=0
1430            do jatm=1,natom
1431              jtypat=typat(atindx1(jatm))
1432              bufsiz=bufsiz+pawtab(jtypat)%lcut_size**2
1433            end do
1434            ABI_ALLOCATE(buf,(ngrad_nondiag+ngradp_nondiag,bufsiz))
1435            do jatm=1,natom
1436              jtypat=typat(atindx1(jatm))
1437              lm_sizej=pawtab(jtypat)%lcut_size**2
1438              if (ngrad_nondiag> 0) buf(1:ngrad_nondiag,bufind+1:bufind+lm_sizej)= &
1439 &             prod_nondiag(jatm)%value(:,:)
1440              if (ngradp_nondiag>0) buf(ngrad_nondiag+1:ngrad_nondiag+ngradp_nondiag, &
1441 &             bufind+1:bufind+lm_sizej)=prodp_nondiag(jatm)%value(:,:)
1442              bufind=bufind+lm_sizej*(ngrad_nondiag+ngradp_nondiag)
1443            end do
1444            call xmpi_sum(buf,my_comm_grid,ier)
1445            bufind=0
1446            do jatm=1,natom
1447              jtypat=typat(atindx1(jatm))
1448              lm_sizej=pawtab(jtypat)%lcut_size**2
1449              if (ngrad> 0) prod_nondiag(jatm)%value(:,:)= &
1450 &             buf(1:ngrad_nondiag,bufind+1:bufind+lm_sizej)
1451              if (ngradp>0) prodp_nondiag(jatm)%value(:,:)= &
1452 &             buf(ngrad_nondiag+1:ngrad_nondiag+ngradp_nondiag,bufind+1:bufind+lm_sizej)
1453              bufind=bufind+lm_sizej*(ngrad_nondiag+ngradp_nondiag)
1454            end do
1455            ABI_DEALLOCATE(buf)
1456          end if
1457        end if
1458 
1459 !      ----------------------------------------------------------------
1460 !      Compute final sums (i.e. derivatives of Sum_ij[rho_ij.Intg{Qij.Vloc}]
1461 
1462 !      ---- Compute terms common to all gradients
1463        jrhoij=1
1464        do irhoij=1,pawrhoij_iatom%nrhoijsel
1465          klmn=pawrhoij_iatom%rhoijselect(irhoij)
1466          klm =pawtab(itypat)%indklmn(1,klmn)
1467          lmin=pawtab(itypat)%indklmn(3,klmn)
1468          lmax=pawtab(itypat)%indklmn(4,klmn)
1469          ro =pawrhoij_iatom%rhoijp(jrhoij,ispden)
1470          ro_d=ro*pawtab(itypat)%dltij(klmn)
1471          do ll=lmin,lmax,2
1472            do ilm=ll**2+1,(ll+1)**2
1473              isel=pawang%gntselect(ilm,klm)
1474              if (isel>0) then
1475                grhat_x=ro_d*pawtab(itypat)%qijl(ilm,klmn)
1476                do mu=1,ngrad
1477                  grhat_tmp(mu,idiag)=grhat_tmp(mu,idiag)+grhat_x*prod(mu,ilm)
1478                end do
1479              end if
1480            end do
1481          end do
1482          jrhoij=jrhoij+pawrhoij_iatom%cplex
1483        end do
1484 
1485 !      ---- Add additional (diagonal) terms for dynamical matrix
1486 !      ---- Terms including rhoij derivatives
1487        if (optgr2==1) then
1488          klmn1=1
1489          do klmn=1,lmn2_size
1490            klm =pawtab(itypat)%indklmn(1,klmn)
1491            lmin=pawtab(itypat)%indklmn(3,klmn)
1492            lmax=pawtab(itypat)%indklmn(4,klmn)
1493            dlt_tmp=pawtab(itypat)%dltij(klmn)
1494            do ll=lmin,lmax,2
1495              do ilm=ll**2+1,(ll+1)**2
1496                isel=pawang%gntselect(ilm,klm)
1497                if (isel>0) then
1498                  ro_d= dlt_tmp*pawtab(itypat)%qijl(ilm,klmn)
1499                  do mu=1,9
1500                    mua=alpha(mu);mub=beta(mu)
1501                    grhat_tmp(ishift_gr2+mu,idiag)=grhat_tmp(ishift_gr2+mu,idiag)&
1502 &                   +ro_d*pawrhoij_iatom%grhoij(ishift_grhoij+mua,klmn1,ispden)*prodp(mub+ishift_gr,ilm)
1503                  end do
1504                end if
1505              end do
1506            end do
1507            klmn1=klmn1+pawrhoij_iatom%cplex
1508          end do ! klmn
1509        end if ! optgr2
1510 
1511 !      ---- Add additional (diagonal) terms for elastic tensor
1512 !      ---- Terms including rhoij derivatives
1513        if (optstr2==1)then
1514          klmn1=1
1515          do klmn=1,lmn2_size
1516            klm =pawtab(itypat)%indklmn(1,klmn)
1517            lmin=pawtab(itypat)%indklmn(3,klmn)
1518            lmax=pawtab(itypat)%indklmn(4,klmn)
1519            dlt_tmp=pawtab(itypat)%dltij(klmn)
1520            do ll=lmin,lmax,2
1521              do ilm=ll**2+1,(ll+1)**2
1522                isel=pawang%gntselect(ilm,klm)
1523                if (isel>0) then
1524                  ro_d=dlt_tmp*pawtab(itypat)%qijl(ilm,klmn)
1525                  mu=1
1526                  do mua=1,6
1527                    do mub=1,6
1528                      grhat_tmp(ishift_str2+mu,iatm)= grhat_tmp(ishift_str2+mu,iatm)&
1529 &                     +ro_d*pawrhoij_iatom%grhoij(mub,klmn1,ispden)*prodp(mua,ilm)
1530                      grhat_tmp(ishift_str2+mu,iatm)= grhat_tmp(ishift_str2+mu,iatm)&
1531 &                     +ro_d*pawrhoij_iatom%grhoij(mua,klmn1,ispden)*prodp(mub,ilm)
1532                      mu=mu+1
1533                    end do
1534 !                  INTERNAL STRAIN CONTRIBUTION
1535                    do idir=1,3
1536                      grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm) = grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm)&
1537 &                     +ro_d*pawrhoij_iatom%grhoij(ishift_grhoij+idir,klmn1,ispden)*prodp(mua,ilm)
1538                      grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm) = grhat_tmp(ishift_str2is+(mua-1)*3+idir,iatm)&
1539 &                     +ro_d*pawrhoij_iatom%grhoij(mua,klmn1,ispden)*prodp(6+idir,ilm)
1540                    end do
1541                  end do
1542                end if
1543              end do
1544            end do
1545            klmn1=klmn1+pawrhoij_iatom%cplex
1546          end do
1547        end if ! optstr2
1548 
1549 !      ---- Add off-diagonal additional contributions for second gradients
1550        if (optgr2==1.or.optstr2==1) then
1551          do jatm=1,natom
1552            jatom_tot=atindx1(jatm);jtypat=typat(jatom_tot)
1553            pawrhoij_jatom => pawrhoij_tot(jatom_tot)
1554 
1555 !          ---- Dynamical matrix
1556            if (optgr2==1) then
1557 
1558 !            Off-diagonal term including rhoij
1559              if (dyfr_cplex==1.or.cplex==1) then
1560                jrhoij=1
1561                do irhoij=1,pawrhoij_jatom%nrhoijsel
1562                  klmn=pawrhoij_jatom%rhoijselect(irhoij)
1563                  klm =pawtab(jtypat)%indklmn(1,klmn)
1564                  lmin=pawtab(jtypat)%indklmn(3,klmn)
1565                  lmax=pawtab(jtypat)%indklmn(4,klmn)
1566                  ro  =pawrhoij_jatom%rhoijp(jrhoij,ispden)
1567                  ro_d=ro*pawtab(jtypat)%dltij(klmn)
1568                  do ll=lmin,lmax,2
1569                    do ilm=ll**2+1,(ll+1)**2
1570                      isel=pawang%gntselect(ilm,klm)
1571                      if (isel>0) then
1572                        grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn)
1573                        do mu=1,9
1574                          grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) &
1575 &                         +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)
1576                        end do
1577                      end if
1578                    end do
1579                  end do
1580                  jrhoij=jrhoij+pawrhoij_jatom%cplex
1581                end do
1582              else
1583                jrhoij=1;mushift=ishift_gr2+9
1584                do irhoij=1,pawrhoij_jatom%nrhoijsel
1585                  klmn=pawrhoij_jatom%rhoijselect(irhoij)
1586                  klm =pawtab(jtypat)%indklmn(1,klmn)
1587                  lmin=pawtab(jtypat)%indklmn(3,klmn)
1588                  lmax=pawtab(jtypat)%indklmn(4,klmn)
1589                  ro  =pawrhoij_jatom%rhoijp(jrhoij,ispden)
1590                  ro_d=ro*pawtab(jtypat)%dltij(klmn)
1591                  do ll=lmin,lmax,2
1592                    do ilm=ll**2+1,(ll+1)**2
1593                      isel=pawang%gntselect(ilm,klm)
1594                      if (isel>0) then
1595                        grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn)
1596                        do mu=1,9
1597                          grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm)&
1598 &                         +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+mu,ilm)
1599                          grhat_tmp(mushift+mu,jatm)=grhat_tmp(mushift+mu,jatm)&
1600 &                         +grhat_x*prod_nondiag(jatm)%value(ishift_gr2+9+mu,ilm)
1601                        end do
1602                      end if
1603                    end do
1604                  end do
1605                  jrhoij=jrhoij+pawrhoij_jatom%cplex
1606                end do
1607              end if
1608 
1609 !            Off-diagonal term including rhoij derivative
1610              if (dyfr_cplex==1.or.cplex==1) then
1611                klmn1=1
1612                do klmn=1,pawrhoij_jatom%lmn2_size
1613                  klm =pawtab(jtypat)%indklmn(1,klmn)
1614                  lmin=pawtab(jtypat)%indklmn(3,klmn)
1615                  lmax=pawtab(jtypat)%indklmn(4,klmn)
1616                  dlt_tmp=pawtab(jtypat)%dltij(klmn)
1617                  do ll=lmin,lmax,2
1618                    do ilm=ll**2+1,(ll+1)**2
1619                      isel=pawang%gntselect(ilm,klm)
1620                      if (isel>0) then
1621                        ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn)
1622                        do mu=1,9
1623                          mua=alpha(mu);mub=beta(mu)
1624                          grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) &
1625 &                         +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) &
1626 &                         *prodp_nondiag(jatm)%value(ishift2_gr+mub,ilm)
1627                        end do
1628                      end if
1629                    end do
1630                  end do
1631                  klmn1=klmn1+pawrhoij_jatom%cplex
1632                end do ! klmn
1633              else ! ngradp_nondiag>=6
1634                klmn1=1;mushift=ishift_gr2+9
1635                do klmn=1,pawrhoij_jatom%lmn2_size
1636                  klm =pawtab(jtypat)%indklmn(1,klmn)
1637                  lmin=pawtab(jtypat)%indklmn(3,klmn)
1638                  lmax=pawtab(jtypat)%indklmn(4,klmn)
1639                  dlt_tmp=pawtab(jtypat)%dltij(klmn)
1640                  do ll=lmin,lmax,2
1641                    do ilm=ll**2+1,(ll+1)**2
1642                      isel=pawang%gntselect(ilm,klm)
1643                      if (isel>0) then
1644                        ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn)
1645                        do mu=1,9
1646                          mua=alpha(mu);mub=beta(mu)
1647                          grhat_tmp(ishift_gr2+mu,jatm)=grhat_tmp(ishift_gr2+mu,jatm) &
1648 &                         +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) &
1649 &                         *prodp_nondiag(jatm)%value(ishift2_gr+mub,ilm)
1650                          grhat_tmp(mushift+mu,jatm)=grhat_tmp(mushift+mu,jatm) &
1651 &                         +ro_d*pawrhoij_jatom%grhoij(ishift_grhoij+mua,klmn1,ispden) &
1652 &                         *prodp_nondiag(jatm)%value(ishift2_gr+3+mub,ilm)
1653                        end do
1654                      end if
1655                    end do
1656                  end do
1657                  klmn1=klmn1+pawrhoij_jatom%cplex
1658                end do
1659              end if
1660            end if ! optgr2
1661 
1662 !          ---- Elastic tensor
1663            if (optstr2==1)then
1664 
1665 !            Off-diagonal term including rhoij
1666              jrhoij=1;
1667              do irhoij=1,pawrhoij_jatom%nrhoijsel
1668                klmn=pawrhoij_jatom%rhoijselect(irhoij)
1669                klm =pawtab(jtypat)%indklmn(1,klmn)
1670                lmin=pawtab(jtypat)%indklmn(3,klmn)
1671                lmax=pawtab(jtypat)%indklmn(4,klmn)
1672                ro  =pawrhoij_jatom%rhoijp(jrhoij,ispden)
1673                ro_d=ro*pawtab(jtypat)%dltij(klmn)
1674                do ll=lmin,lmax,2
1675                  do ilm=ll**2+1,(ll+1)**2
1676                    isel=pawang%gntselect(ilm,klm)
1677                    if (isel>0) then
1678                      grhat_x=ro_d*pawtab(jtypat)%qijl(ilm,klmn)
1679                      do mu=1,18
1680                        grhat_tmp2(mu,jatm)=grhat_tmp2(mu,jatm) &
1681 &                       +grhat_x*prod_nondiag(jatm)%value(ishift_str2is+mu,ilm)
1682                      end do
1683                    end if
1684                  end do
1685                end do
1686                jrhoij=jrhoij+pawrhoij_jatom%cplex
1687              end do
1688 !            Off-diagonal term including rhoij derivative
1689              klmn1=1
1690              do klmn=1,pawrhoij_jatom%lmn2_size
1691                klm =pawtab(jtypat)%indklmn(1,klmn)
1692                lmin=pawtab(jtypat)%indklmn(3,klmn)
1693                lmax=pawtab(jtypat)%indklmn(4,klmn)
1694                dlt_tmp=pawtab(jtypat)%dltij(klmn)
1695                do ll=lmin,lmax,2
1696                  do ilm=ll**2+1,(ll+1)**2
1697                    isel=pawang%gntselect(ilm,klm)
1698                    if (isel>0) then
1699                      ro_d=dlt_tmp*pawtab(jtypat)%qijl(ilm,klmn)
1700                      mu=1
1701                      do mua=1,6
1702                        do idir=1,3
1703                          grhat_tmp2((mua-1)*3+idir,jatm) = grhat_tmp2((mua-1)*3+idir,jatm) &
1704 &                         +ro_d*pawrhoij_jatom%grhoij(mua,klmn1,ispden) &
1705 &                         *prodp_nondiag(jatm)%value(idir,ilm)
1706                        end do
1707                      end do
1708                    end if
1709                  end do
1710                end do
1711                klmn1=klmn1+pawrhoij_jatom%cplex
1712              end do
1713            end if ! optstr2
1714 
1715          end do ! jatm
1716        end if ! optgr2 or optstr2
1717 
1718 !    ----------------------------------------------------------------
1719 !    End of loop over spin components
1720 
1721      end do ! ispden
1722 
1723 !    Eventually free temporary space for g_l(r).Y_lm(r) factors
1724      if (pawfgrtab_iatom%gylm_allocated==2) then
1725        ABI_DEALLOCATE(pawfgrtab_iatom%gylm)
1726        ABI_ALLOCATE(pawfgrtab_iatom%gylm,(0,0))
1727        pawfgrtab_iatom%gylm_allocated=0
1728      end if
1729      if (pawfgrtab_iatom%gylmgr_allocated==2) then
1730        ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr)
1731        ABI_ALLOCATE(pawfgrtab_iatom%gylmgr,(0,0,0))
1732        pawfgrtab_iatom%gylmgr_allocated=0
1733      end if
1734      if (pawfgrtab_iatom%gylmgr2_allocated==2) then
1735        ABI_DEALLOCATE(pawfgrtab_iatom%gylmgr2)
1736        ABI_ALLOCATE(pawfgrtab_iatom%gylmgr2,(0,0,0))
1737        pawfgrtab_iatom%gylmgr2_allocated=0
1738      end if
1739      if (pawfgrtab_iatom%expiqr_allocated==2) then
1740        ABI_DEALLOCATE(pawfgrtab_iatom%expiqr)
1741        ABI_ALLOCATE(pawfgrtab_iatom%expiqr,(0,0))
1742        pawfgrtab_iatom%expiqr_allocated=0
1743      end if
1744 
1745 !    ----------------------------------------------------------------
1746 !    Copy results in corresponding arrays
1747 
1748 !    ==== Forces ====
1749 !    Convert from cartesian to reduced coordinates
1750      if (optgr==1) then
1751        mushift=3*(iatm-1)
1752        tmp(1:3)=grhat_tmp(ishift_gr+1:ishift_gr+3,idiag)
1753        do mu=1,3
1754          hatgr(mu+mushift)=rprimd(1,mu)*tmp(1)+rprimd(2,mu)*tmp(2)+rprimd(3,mu)*tmp(3)
1755        end do
1756      end if
1757 
1758 !    ==== Stresses ====
1759      if (optstr==1) then
1760        hatstr(1:6)=hatstr(1:6)+grhat_tmp(ishift_str+1:ishift_str+6,idiag)
1761      end if
1762 
1763 !    ==== Frozen wf part of dyn. matrix ====
1764      if (optgr2==1) then
1765        do jatm=1,natom
1766          do mu=1,9
1767            mua=alpha(mu);mub=beta(mu)
1768            dyfr(1,mub,mua,jatm,iatm)=grhat_tmp(ishift_gr2+mu,jatm)
1769          end do
1770          if (dyfr_cplex==2.and.cplex==2) then
1771            mushift=ishift_gr2+9
1772            do mu=1,9
1773              mua=alpha(mu);mub=beta(mu)
1774              dyfr(2,mub,mua,jatm,iatm)=grhat_tmp(mushift+mu,jatm)
1775            end do
1776          end if
1777        end do
1778      end if
1779 
1780 !    ==== Elastic tensor ====
1781      if (optstr2==1) then
1782        eltfr(1:6,1:6)=eltfr(1:6,1:6)+ &
1783 &       reshape(grhat_tmp(ishift_str2+1:ishift_str2+36,iatm),(/6,6/))
1784 !      Convert internal Strain in reduced coordinates
1785        do mua = 1,6
1786          tmp(1:3)=grhat_tmp(ishift_str2is+(mua-1)*3+1:ishift_str2is+(mua-1)*3+3,iatm)
1787          do idir=1,3
1788            eltfr(6+(iatm-1)*3+idir,mua)=eltfr(6+(iatm-1)*3+idir,mua)+ &
1789 &           (rprimd(1,idir)*tmp(1)+rprimd(2,idir)*tmp(2)+rprimd(3,idir)*tmp(3))
1790          end do
1791          do jatm=1,natom
1792            tmp(1:3)=grhat_tmp2((mua-1)*3+1:(mua-1)*3+3,jatm)
1793            do idir=1,3
1794              eltfr(6+(iatm-1)*3+idir,mua)=eltfr(6+(iatm-1)*3+idir,mua)+ &
1795 &             (rprimd(1,idir)*tmp(1)+rprimd(2,idir)*tmp(2)+rprimd(3,idir)*tmp(3))
1796            end do
1797          end do
1798        end do
1799      end if
1800 
1801 !    ----------------------------------------------------------------
1802 !    End loops on types and atoms
1803 
1804      ABI_DEALLOCATE(vloc)
1805      if (ngrad>0)  then
1806        ABI_DEALLOCATE(prod)
1807      end if
1808      if (ngradp>0)  then
1809        ABI_DEALLOCATE(prodp)
1810      end if
1811      if (optgr2==1.or.optstr2==1) then
1812        do jatm=1,natom
1813          ABI_DEALLOCATE(prod_nondiag(jatm)%value)
1814          ABI_DEALLOCATE(prodp_nondiag(jatm)%value)
1815        end do
1816      end if
1817    end do
1818    iatshft=iatshft+nattyp(itypat)
1819  end do
1820 
1821 !Reduction in case of parallelisation over atoms
1822  if (paral_atom) then
1823    bufsiz=3*natom*optgr+6*optstr
1824    if (save_memory) bufsiz=bufsiz+9*dyfr_cplex*natom**2*optgr2+6*(6+3*natom)*optstr2
1825    if (bufsiz>0) then
1826      ABI_ALLOCATE(buf1,(bufsiz))
1827      if (optgr==1) buf1(1:3*natom)=hatgr(1:3*natom)
1828      indx=optgr*3*natom
1829      if (optstr==1) buf1(indx+1:indx+6)=hatstr(1:6)
1830      indx=indx+optstr*6
1831      if (save_memory) then
1832        if (optgr2==1) then
1833          buf1(indx+1:indx+9*dyfr_cplex*natom**2)= &
1834 &         reshape(dyfr,(/9*dyfr_cplex*natom**2/))
1835          indx=indx+9*dyfr_cplex*natom**2
1836        end if
1837        if (optstr2==1) then
1838          buf1(indx+1:indx+6*(6+3*natom))= &
1839 &         reshape(eltfr,(/6*(6+3*natom)/))
1840          indx=indx+6*(6+3*natom)
1841        end if
1842      end if
1843      call xmpi_sum(buf1,my_comm_atom,ier)
1844      if (optgr==1) hatgr(1:3*natom)=buf1(1:3*natom)
1845      indx=optgr*3*natom
1846      if (optstr==1) hatstr(1:6)=buf1(indx+1:indx+6)
1847      indx=indx+optstr*6
1848      if (save_memory) then
1849        if (optgr2==1) then
1850          dyfr(1:dyfr_cplex,1:3,1:3,1:natom,1:natom)= &
1851 &         reshape(buf1(indx+1:indx+9*dyfr_cplex*natom**2),(/dyfr_cplex,3,3,natom,natom/))
1852          indx=indx+9*dyfr_cplex*natom**2
1853        end if
1854        if (optstr2==1) then
1855          eltfr(1:6+3*natom,1:6)= &
1856 &         reshape(buf1(indx+1:indx+6*(6+3*natom)),(/6+3*natom,6/))
1857          indx=indx+6*(6+3*natom)
1858        end if
1859      end if
1860      ABI_DEALLOCATE(buf1)
1861    end if
1862  end if
1863 
1864 !Deallocate additional memory
1865  ABI_DEALLOCATE(grhat_tmp)
1866  if (optgr2==1.or.optstr2==1) then
1867    ABI_DEALLOCATE(mu4)
1868    ABI_DEALLOCATE(atindx)
1869    if (optgr2==1.or.optstr2==1) then
1870      ABI_DEALLOCATE(vpsp1_gr)
1871    end if
1872    if (optstr2==1) then
1873      ABI_DEALLOCATE(grhat_tmp2)
1874      ABI_DEALLOCATE(vpsp1_str)
1875    end if
1876    ABI_DATATYPE_DEALLOCATE(prod_nondiag)
1877    ABI_DATATYPE_DEALLOCATE(prodp_nondiag)
1878    if (.not.save_memory) then
1879      do jatom=1,size(pawfgrtab)
1880        pawfgrtab_jatom => pawfgrtab(jatom)
1881        if (pawfgrtab(jatom)%gylm_allocated==2) then
1882          ABI_DEALLOCATE(pawfgrtab(jatom)%gylm)
1883          ABI_ALLOCATE(pawfgrtab(jatom)%gylm,(0,0))
1884          pawfgrtab(jatom)%gylm_allocated=0
1885        end if
1886        if (pawfgrtab(jatom)%gylmgr_allocated==2) then
1887          ABI_DEALLOCATE(pawfgrtab(jatom)%gylmgr)
1888          ABI_ALLOCATE(pawfgrtab(jatom)%gylmgr,(0,0,0))
1889          pawfgrtab(jatom)%gylmgr_allocated=0
1890        end if
1891        if (pawfgrtab(jatom)%gylmgr2_allocated==2) then
1892          ABI_DEALLOCATE(pawfgrtab(jatom)%gylmgr2)
1893          ABI_ALLOCATE(pawfgrtab(jatom)%gylmgr2,(0,0,0))
1894          pawfgrtab(jatom)%gylmgr2_allocated=0
1895        end if
1896        if (pawfgrtab(jatom)%expiqr_allocated==2) then
1897          ABI_DEALLOCATE(pawfgrtab(jatom)%expiqr)
1898          ABI_ALLOCATE(pawfgrtab(jatom)%expiqr,(0,0))
1899          pawfgrtab(jatom)%expiqr_allocated=0
1900        end if
1901      end do
1902    end if
1903    if (paral_atom) then
1904      if ((.not.save_memory).and.paral_atom_pawfgrtab) then
1905        call pawfgrtab_free(pawfgrtab_tot)
1906        ABI_DATATYPE_DEALLOCATE(pawfgrtab_tot)
1907      end if
1908      if (paral_atom_pawrhoij) then
1909        call pawrhoij_free(pawrhoij_tot)
1910        ABI_DATATYPE_DEALLOCATE(pawrhoij_tot)
1911      end if
1912    end if
1913  end if
1914 
1915 !----------------------------------------------------------------------
1916 !Update non-local gardients
1917 
1918 !===== Update forces =====
1919  if (optgr==1) then
1920    grnl(1:3*natom)=grnl(1:3*natom)+hatgr(1:3*natom)
1921    ABI_DEALLOCATE(hatgr)
1922  end if
1923 
1924 !===== Convert stresses (add diag and off-diag contributions) =====
1925  if (optstr==1) then
1926 
1927 !  Has to compute int[nhat*vtrial]
1928    hatstr_diag=zero
1929    if (nspden==1.or.dimvtrial==1) then
1930      do ic=1,nfft
1931        hatstr_diag=hatstr_diag+vtrial_(ic,1)*nhat(ic,1)
1932      end do
1933    else if (nspden==2) then
1934      do ic=1,nfft
1935        hatstr_diag=hatstr_diag+vtrial_(ic,1)*nhat(ic,2)+vtrial_(ic,2)*(nhat(ic,1)-nhat(ic,2))
1936      end do
1937    else if (nspden==4) then
1938      do ic=1,nfft
1939        hatstr_diag=hatstr_diag+half*(vtrial_(ic,1)*(nhat(ic,1)+nhat(ic,4)) &
1940 &       +vtrial_(ic,2)*(nhat(ic,1)-nhat(ic,4))) &
1941 &       +vtrial_(ic,3)*nhat(ic,2)+vtrial_(ic,4)*nhat(ic,3)
1942      end do
1943    end if
1944    hatstr_diag=hatstr_diag*fact_ucvol
1945    if (paral_grid) then
1946      call xmpi_sum(hatstr_diag,my_comm_grid,ier)
1947    end if
1948 
1949 !  Convert hat contribution
1950 
1951    hatstr(1:3)=(hatstr(1:3)+hatstr_diag)/ucvol
1952    hatstr(4:6)= hatstr(4:6)/ucvol
1953 
1954 !  Add to already computed NL contrib
1955    nlstr(1:6)=nlstr(1:6)+hatstr(1:6)
1956 
1957 !  Apply symmetries
1958    call stresssym(gprimd,nsym,nlstr,symrec)
1959  end if
1960 
1961 !===== Convert dynamical matrix (from cartesian to reduced coordinates) =====
1962  if (optgr2==1) then
1963    do iatm=1,natom
1964      do jatm=1,natom
1965        do mua=1,3
1966          do mub=1,3
1967            work1(1,mua,mub)=dyfr(1,mub,mua,jatm,iatm)+dyfr(1,mua,mub,iatm,jatm)
1968          end do
1969        end do
1970        if (dyfr_cplex==2) then
1971          do mua=1,3
1972            do mub=1,3
1973              work1(2,mua,mub)=dyfr(2,mub,mua,jatm,iatm)-dyfr(2,mua,mub,iatm,jatm)
1974            end do
1975          end do
1976        end if
1977        do mu=1,3
1978          work2(:,:,mu)=rprimd(1,mu)*work1(:,:,1)+rprimd(2,mu)*work1(:,:,2)+rprimd(3,mu)*work1(:,:,3)
1979        end do
1980        do mub=1,3
1981          do mua=1,3
1982            dyfrnl(:,mua,mub,jatm,iatm)=dyfrnl(:,mua,mub,jatm,iatm) &   ! Already contains NL projectors contribution
1983 &          +rprimd(1,mua)*work2(:,1,mub) &
1984 &           +rprimd(2,mua)*work2(:,2,mub) &
1985 &           +rprimd(3,mua)*work2(:,3,mub)
1986          end do
1987        end do
1988      end do
1989    end do
1990    ABI_DEALLOCATE(dyfr)
1991  end if
1992 
1993 !===== Update elastic tensor =====
1994  if (optstr2==1) then
1995    eltfrnl(1:6+3*natom,1:6)=eltfrnl(1:6+3*natom,1:6)+eltfr(1:6+3*natom,1:6)
1996    ABI_DEALLOCATE(eltfr)
1997  end if
1998 
1999 !----------------------------------------------------------------------
2000 !End
2001 
2002 !Destroy temporary space
2003  if (usexcnhat==0)  then
2004    ABI_DEALLOCATE(vtrial_)
2005  end if
2006 
2007 !Destroy atom tables used for parallelism
2008  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2009  if (paral_atom) then
2010    ABI_DEALLOCATE(atm_indx)
2011  end if
2012 
2013 !Destroy FFT tables used for parallelism
2014  if ((optgr2==1.or.optstr2==1).and.(.not.present(comm_fft))) then
2015    call destroy_distribfft(my_distribfft)
2016    ABI_DATATYPE_DEALLOCATE(my_distribfft)
2017  end if
2018 
2019  DBG_ENTER("COLL")
2020 
2021  CONTAINS

pawgrnl/pawgrnl_convert [ Functions ]

[ Top ] [ pawgrnl ] [ Functions ]

NAME

  pawgrnl_convert

FUNCTION

  notation: Convert index of the elastic tensor:
    - voigt notation       => 32
    - normal notation      => 3 3 2 2
    - notation for gylmgr2 => 32 32 32 32 => 4 4 4

INPUTS

  eps_alpha, eps_beta, eps_delta, eps_gamma

OUTPUT

  mu4(4) = array with index for the second derivative of gylm

SIDE EFFECTS

  mu4(4) = input : array with index for the second derivative of gylm
           output: the 4 indexes for the calculation of the second derivative of gylm

SOURCE

2047 subroutine pawgrnl_convert(mu4,eps_alpha,eps_beta,eps_gamma,eps_delta)
2048 
2049 
2050 !This section has been created automatically by the script Abilint (TD).
2051 !Do not modify the following lines by hand.
2052 #undef ABI_FUNC
2053 #define ABI_FUNC 'pawgrnl_convert'
2054 !End of the abilint section
2055 
2056  implicit none
2057 
2058 !Arguments ------------------------------------
2059  !scalar
2060  integer,intent(in)  :: eps_alpha,eps_beta
2061  integer,optional,intent(in)  :: eps_gamma,eps_delta
2062  !array
2063  integer,intent(inout) :: mu4(4)
2064 
2065 !Local variables-------------------------------
2066  integer :: eps1,eps2,i,j,k
2067  integer,allocatable :: mu_temp(:)
2068 
2069 ! *************************************************************************
2070 
2071  ABI_ALLOCATE(mu_temp,(4))
2072  if (present(eps_gamma).and.present(eps_delta)) then
2073    mu_temp(1)=eps_alpha
2074    mu_temp(2)=eps_beta
2075    mu_temp(3)=eps_gamma
2076    mu_temp(4)=eps_delta
2077  else
2078    mu_temp(1)=eps_alpha
2079    mu_temp(2)=eps_beta
2080    mu_temp(3)= 0
2081    mu_temp(4)= 0
2082  end if
2083  k=1
2084  do i=1,2
2085    eps1=mu_temp(i)
2086    do j=1,2
2087      eps2=mu_temp(2+j)
2088      if(eps1==eps2) then
2089        if(eps1==1) mu4(k)=1;
2090        if(eps1==2) mu4(k)=2;
2091        if(eps1==3) mu4(k)=3;
2092      else
2093        if((eps1==3.and.eps2==2).or.(eps1==2.and.eps2==3)) mu4(k)=4;
2094        if((eps1==3.and.eps2==1).or.(eps1==1.and.eps2==3)) mu4(k)=5;
2095        if((eps1==1.and.eps2==2).or.(eps1==2.and.eps2==1)) mu4(k)=6;
2096      end if
2097      k=k+1
2098    end do
2099  end do
2100  ABI_DEALLOCATE(mu_temp)
2101 
2102 end subroutine pawgrnl_convert
2103 ! ------------------------------------------------
2104 
2105 end subroutine pawgrnl