TABLE OF CONTENTS


ABINIT/m_pead_nl_loop [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pead_nl_loop

FUNCTION

COPYRIGHT

  Copyright (C) 2002-2022 ABINIT group (MVeithen,MB)
  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

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_pead_nl_loop
23 
24  use defs_basis
25  use defs_wvltypes
26  use m_wffile
27  use m_abicore
28  use m_xmpi
29  use m_hdr
30  use m_dtset
31  use m_dtfil
32 #if defined HAVE_MPI2
33  use mpi
34 #endif
35 
36  use defs_datatypes, only : pseudopotential_type
37  use defs_abitypes, only : MPI_type
38  use m_time,     only : timab
39  use m_kg,       only : getph, mkkpg
40  use m_cgtools,  only : dotprod_vn, dotprod_g
41  use m_fft,      only : fourdp, fftpac, fourwf
42  use m_ioarr,    only : read_rhor
43  use m_pawtab,   only : pawtab_type
44  use m_pawrhoij, only : pawrhoij_type
45  use m_pawcprj,    only : pawcprj_type
46  use m_inwffil,  only : inwffil
47  use m_spacepar, only : hartre
48  use m_initylmg, only : initylmg
49  use m_dfpt_mkvxc, only : dfpt_mkvxc
50  use m_mkcore,     only : dfpt_mkcore
51  use m_mklocl,     only : dfpt_vlocal
52  use m_hamiltonian,only : init_hamiltonian, gs_hamiltonian_type
53  use m_mkffnl,     only : mkffnl
54  use m_mpinfo,     only : proc_distrb_cycle
55  use m_nonlop,     only : nonlop
56  use m_dfptnl_pert, only : dfptnl_exc3
57 
58  implicit none
59 
60  private
61 
62 #if defined HAVE_MPI1
63  include 'mpif.h'
64 #endif

ABINIT/pead_nl_loop [ Functions ]

[ Top ] [ Functions ]

NAME

 pead_nl_loop

FUNCTION

 Loop over the perturbations j1, j2 and j3

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions
  cgindex(nkpt,nsppol) = for each k-point, cgindex tores the location
                         of the WF in the cg array
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
  gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double
   that of the basis sphere--appropriate for charge density rho(G),
   Hartree potential, and pseudopotentials
  kg(3,mpw*mkmem)=reduced planewave coordinates
  kneigh(30,nkpt) = index of the neighbours of each k-point
  kg_neigh(30,nkpt,3) = necessary to construct the vector joining a k-point
                         to its nearest neighbour in case of a single k-point,
                         a line of k-points or a plane of k-points.
  kptindex(2,nkpt3)= index of the k-points in the reduced BZ
                     related to a k-point in the full BZ
  kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
  kxc(nfft,nkxc)=exchange-correlation kernel
  k3xc(nfft,nk3xc)=third-order exchange-correlation kernel
  mband = maximum number of bands
  mgfft = maximum single fft dimension
  mkmem = Number of k points treated by this node.
  mkmem_max = maximal number of k-points on each processor (MPI //)
  mk1mem = Number of k points for first-order WF treated by this node.
  mpert =maximum number of ipert
  mpi_enreg=MPI-parallelisation information
  mpw   = maximum number of planewaves in basis sphere (large number)
  mvwtk(30,nkpt) = weights to compute the finite difference ddk
  natom = number of atoms in unit cell
  nfft  = (effective) number of FFT grid points (for this processor)
  nkpt  = number of k points
  nkpt3 = number of k-points in the full BZ
  nkxc=second dimension of the array kxc, see rhohxc.f for a description
  nneigh  = total number of neighbours required to evaluate the finite
          difference formula
  nspinor = number of spinorial components of the wavefunctions
  nsppol = number of channels for spin-polarization (1 or 2)
  npwarr(nkpt) = array holding npw for each k point
  occ(mband*nkpt*nsppol) = occupation number for each band and k
  psps <type(pseudopotential_type)> = variables related to pseudopotentials
  pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat
                           between k-points
  rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
       that have to be computed
       1   ->   element has to be computed explicitely
      -1   ->   use symmetry operations to obtain the corresponding element
  rprimd(3,3)=dimensional primitive translations (bohr)
  ucvol = unit cell volume (bohr^3)
  xred(3,natom) = reduced atomic coordinates

OUTPUT

  blkflg(3,mpert,3,mpert) = flags for each element of the 3DTE
                             (=1 if computed)
  d3lo(2,3,mpert,3,mpert,3,mpert) = matrix of the 3DTEs

SIDE EFFECTS

  hdr <type(hdr_type)>=the header of wf, den and pot files

SOURCE

142 subroutine pead_nl_loop(blkflg,cg,cgindex,dtfil,dtset,d3lo,&
143 & gmet,gprimd,gsqcut, &
144 & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,mband,mgfft,mkmem,mkmem_max,mk1mem,&
145 & mpert,mpi_enreg,mpw,mvwtk,natom,nfft,nkpt,nkpt3,nkxc,nk3xc,nneigh,nspinor,nsppol,&
146 & npwarr,occ,psps,pwind,&
147 & rfpert,rprimd,ucvol,xred)
148 
149 !Arguments ------------------------------------
150 !scalars
151  integer,intent(in) :: mband,mgfft,mk1mem,mkmem,mkmem_max,mpert,mpw,natom,nfft
152  integer,intent(in) :: nk3xc,nkpt,nkpt3,nkxc,nneigh,nspinor,nsppol
153  real(dp),intent(in) :: gsqcut,ucvol
154  type(MPI_type),intent(inout) :: mpi_enreg
155  type(datafiles_type),intent(in) :: dtfil
156  type(dataset_type),intent(inout) :: dtset
157  type(hdr_type),intent(inout) :: hdr
158  type(pseudopotential_type),intent(in) :: psps
159 !arrays
160  integer,intent(in) :: cgindex(nkpt,nsppol),kg(3,mk1mem*mpw),kneigh(30,nkpt)
161  integer,intent(in) :: kg_neigh(30,nkpt,3)
162  integer,intent(in) :: kptindex(2,nkpt3),npwarr(nkpt),pwind(mpw,nneigh,mkmem)
163  integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
164  integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i
165  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3)
166  real(dp),intent(in) :: gprimd(3,3),k3xc(nfft,nk3xc),kpt3(3,nkpt3)
167  real(dp),intent(in) :: kxc(nfft,nkxc),mvwtk(30,nkpt),rprimd(3,3)
168  real(dp),intent(in) :: xred(3,natom)
169  real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
170  real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert) !vz_i
171 
172 !Local variables-------------------------------
173 !scalars
174  integer,parameter :: level=51
175  integer :: ask_accurate,counter,cplex,formeig,i1dir
176  integer :: i1pert,i2dir,i2pert,i3dir,i3pert,iatom,ierr,index,ir
177  integer :: ireadwf,itypat,mcg,mpsang,n1,n2,n3,n3xccc,nfftot,nspden,option,optorth
178  integer :: pert1case,pert2case,pert3case,rdwrpaw,timrev,comm_cell
179  logical :: nmxc
180  real(dp) :: ecut_eff,exc3(2)
181  character(len=500) :: message
182  character(len=fnlen) :: fiden1i,fiwf1i,fiwf3i
183  type(wffile_type) :: wff1,wff2,wfft1,wfft2
184  type(wvl_data) :: wvl
185  type(hdr_type) :: hdr_den
186 !arrays
187  integer,allocatable :: atindx(:),atindx1(:),nattyp(:)
188  real(dp) :: d3_berry(2,3),rho_dum(1),tsec(2),ylmgr_dum(1)
189  real(dp),allocatable :: cg1(:,:),cg3(:,:),eigen1(:),ph1d(:,:),rho1r1(:,:)
190  real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1(:)
191  real(dp),allocatable :: vpsp1(:),vtrial1(:,:),vxc1(:,:),work(:)
192  real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:),ylm(:,:,:)
193  type(pawrhoij_type),allocatable :: rhoij_dum(:)
194 
195 ! ***********************************************************************
196 
197  call timab(502,1,tsec)
198 
199  comm_cell = mpi_enreg%comm_cell
200 
201  timrev = 1
202  cplex = 2 - timrev
203  nspden = dtset%nspden
204  ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2
205  mpsang = psps%mpsang
206  optorth=1;if (psps%usepaw==1) optorth=0
207 
208  ABI_MALLOC(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
209  ABI_MALLOC(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
210  ABI_MALLOC(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
211  ABI_MALLOC(rho1r1,(cplex*nfft,dtset%nspden))
212  ABI_MALLOC(rho2r1,(cplex*nfft,dtset%nspden))
213  ABI_MALLOC(rho2g1,(2,nfft))
214  ABI_MALLOC(rho3r1,(cplex*nfft,dtset%nspden))
215  ABI_MALLOC(ylm,(2,dtset%mpw*dtset%mkmem,mpsang*mpsang*psps%useylm))
216 
217  ask_accurate=1 ; formeig = 1 ; ireadwf = 1
218  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
219  nfftot=n1*n2*n3
220 
221 !Generate an index table of atoms, in order for them to be used
222 !type after type.
223  ABI_MALLOC(atindx,(natom))
224  ABI_MALLOC(atindx1,(natom))
225  ABI_MALLOC(nattyp,(psps%ntypat))
226  index=1
227  do itypat=1,psps%ntypat
228    nattyp(itypat)=0
229    do iatom=1,natom
230      if(dtset%typat(iatom)==itypat)then
231        atindx(iatom)=index
232        atindx1(index)=iatom
233        index=index+1
234        nattyp(itypat)=nattyp(itypat)+1
235      end if
236    end do
237  end do
238 
239 !Generate the 1-dimensional phases
240  ABI_MALLOC(ph1d,(2,3*(2*mgfft+1)*natom))
241  call getph(atindx,natom,n1,n2,n3,ph1d,xred)
242 
243 !Set up the Ylm for each k point
244  if (psps%useylm==1) then
245    call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,&
246 &   dtset%mpw,dtset%nband,dtset%nkpt,&
247 &   npwarr,dtset%nsppol,0,rprimd,ylm,ylmgr_dum)
248  end if
249 
250  ABI_MALLOC(vpsp1,(cplex*nfft))
251  ABI_MALLOC(xccc3d1,(cplex*nfft))
252  ABI_MALLOC(xccc3d2,(cplex*nfft))
253  ABI_MALLOC(xccc3d3,(cplex*nfft))
254  ABI_MALLOC(vhartr1,(cplex*nfft))
255  ABI_MALLOC(vxc1,(cplex*nfft,dtset%nspden))
256  ABI_MALLOC(vtrial1,(cplex*nfft,dtset%nspden))
257 
258 !Loop over the perturbations j1, j2, j3
259 
260  pert1case = 0 ; pert2case = 0 ; pert3case = 0
261 
262  do i1pert = 1, mpert
263    do i1dir = 1, 3
264 
265      if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then
266 
267        pert1case = i1dir + (i1pert-1)*3
268        counter = pert1case
269        call appdig(pert1case,dtfil%fnamewff1,fiwf1i)
270 
271        mcg=mpw*nspinor*mband*mkmem*nsppol
272        call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
273 &       formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
274 &       dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
275 &       dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
276 &       dtset%nsppol,dtset%nsym,&
277 &       occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
278 &       dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl)
279 
280        if (ireadwf==1) then
281          call WffClose (wff1,ierr)
282        end if
283 
284        rho1r1(:,:) = 0._dp
285        if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
286          rdwrpaw=0
287          call appdig(pert1case,dtfil%fildens1in,fiden1i)
288 
289          call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho1r1, &
290          hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
291          call hdr_den%free()
292        end if
293 
294        xccc3d1(:) = 0._dp
295        if ((psps%n1xccc/=0).and.(i1pert <= natom)) then
296          call dfpt_mkcore(cplex,i1dir,i1pert,natom,psps%ntypat,n1,psps%n1xccc,&
297 &         n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
298 &         psps%xcccrc,psps%xccc1d,xccc3d1,xred)
299        end if ! psps%n1xccc/=0
300 
301        do i3pert = 1, mpert
302          do i3dir = 1, 3
303 
304            if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then
305 
306              pert3case = i3dir + (i3pert-1)*3
307              counter = 100*pert3case + pert1case
308              call appdig(pert3case,dtfil%fnamewff1,fiwf3i)
309 
310              mcg=mpw*nspinor*mband*mkmem*nsppol
311              call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
312 &             formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
313 &             dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
314 &             dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
315 &             dtset%nsppol,dtset%nsym,&
316 &             occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
317 &             dtfil%unkg1,wff2,wfft2,dtfil%unwff2,&
318 &             fiwf3i,wvl)
319              if (ireadwf==1) then
320                call WffClose (wff2,ierr)
321              end if
322 
323              rho3r1(:,:) = 0._dp
324              if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
325                rdwrpaw=0
326                call appdig(pert3case,dtfil%fildens1in,fiden1i)
327 
328                call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho3r1, &
329                hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
330                call hdr_den%free()
331              end if
332 
333              xccc3d3(:) = 0._dp
334              if ((psps%n1xccc/=0).and.(i3pert <= natom)) then
335                call dfpt_mkcore(cplex,i3dir,i3pert,natom,psps%ntypat,n1,psps%n1xccc,&
336 &               n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
337 &               psps%xcccrc,psps%xccc1d,xccc3d3,xred)
338              end if ! psps%n1xccc/=0
339 
340              do i2pert = 1, mpert
341 
342 !              In case of electric field perturbation, evaluate the ddk
343 !              using the finite difference expression of
344 !              Marzari and Vanderbilt PRB 56, 12847 (1997) [[cite:Marzari1997]].
345 
346                d3_berry(:,:) = 0._dp
347 
348                if ((i2pert==dtset%natom+2).and.&
349 &               (maxval(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert)) == 1)) then
350 
351                  call timab(511,1,tsec)
352                  call pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,&
353 &                 i1pert,i3pert,i1dir,i3dir,&
354 &                 kneigh,kg_neigh,kptindex,kpt3,mband,mkmem,mkmem_max,mk1mem,&
355 &                 mpi_enreg,mpw,mvwtk,natom,nkpt,nkpt3,nneigh,npwarr,nspinor,nsppol,pwind)
356                  call timab(511,2,tsec)
357 
358                end if
359 
360                if (mpi_enreg%me == 0) then
361 
362                  if(sum(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert))>0)then
363                    write(message,'(a,a,a,a,a,a)')ch10,ch10,&
364 &                   ' Decomposition of the third-order energy for the set of perturbations',ch10
365                    call wrtout(std_out,message,'COLL')
366                    call wrtout(ab_out,message,'COLL')
367                    if (i1pert < natom + 1) then
368                      write(message,'(a,i3,a,i3)') &
369 &                     ' j1 : displacement of atom ',i1pert,' along direction ', i1dir
370                    end if
371                    if (i1pert == dtset%natom + 2) then
372                      write(message,'(a,i4)')' j1 : homogeneous electric field along direction ',i1dir
373                    end if
374                    call wrtout(std_out,message,'COLL')
375                    call wrtout(ab_out,message,'COLL')
376                    if (i3pert < natom + 1) then
377                      write(message,'(a,i3,a,i3,a)') &
378 &                     ' j3 : displacement of atom ',i3pert,' along direction ', i3dir,ch10
379                    end if
380                    if (i3pert == dtset%natom + 2) then
381                      write(message,'(a,i4,a)')' j3 : homogeneous electric field along direction ',i3dir,ch10
382                    end if
383                    call wrtout(std_out,message,'COLL')
384                    call wrtout(ab_out,message,'COLL')
385                  end if
386 
387                end if ! mpi_enreg%me == 0
388 
389                do i2dir = 1, 3
390 
391                  if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
392                    pert2case = i2dir + (i2pert-1)*3
393 
394                    blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
395 
396 !                  Read the first-order densities from disk-files
397                    rho2r1(:,:) = 0._dp ; rho2g1(:,:) = 0._dp
398 
399                    if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
400                      rdwrpaw=0
401                      call appdig(pert2case,dtfil%fildens1in,fiden1i)
402 
403                      call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho2r1, &
404                      hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
405                      call hdr_den%free()
406 
407 !                    Compute up+down rho1(G) by fft
408                      ABI_MALLOC(work,(cplex*nfft))
409                      work(:)=rho2r1(:,1)
410                      call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfft,1,dtset%ngfft,0)
411                      ABI_FREE(work)
412 
413                    end if
414 
415 !                  Compute first-order local potentials
416 !                  (hartree, xc and pseudopotential)
417 
418                    n3xccc=0; if(psps%n1xccc/=0)n3xccc=nfft
419                    xccc3d2(:)=0._dp ; vpsp1(:)=0._dp
420 
421                    if (i2pert <= natom) then
422 
423                      call dfpt_vlocal(atindx,cplex,gmet,gsqcut,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,natom,&
424 &                     nattyp,nfft,dtset%ngfft,psps%ntypat,n1,n2,n3,ph1d,psps%qgrid_vl,&
425 &                     dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
426 
427                      if (psps%n1xccc/=0) then
428                        call dfpt_mkcore(cplex,i2dir,i2pert,natom,psps%ntypat,n1,psps%n1xccc,&
429 &                       n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
430 &                       psps%xcccrc,psps%xccc1d,xccc3d2,xred)
431                      end if ! psps%n1xccc/=0
432 
433                    end if  ! i2pert <= natom
434 
435                    call hartre(cplex,gsqcut,3,0,mpi_enreg,nfft,dtset%ngfft,dtset%nkpt,&
436                                &dtset%rcut,rho2g1,rprimd,dtset%vcutgeo,vhartr1)
437                    option=1 ; nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
438                    call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,dtset%ngfft,&
439 &                   rho_dum,0,rho_dum,0,nkxc,nmxc,dtset%nspden,n3xccc,option,&
440 &                   dtset%qptn,rho2r1,rprimd,0,vxc1,xccc3d2)
441 
442                    if(dtset%nsppol==1)then
443                      if(cplex==1)then
444                        do ir=1,nfft
445                          vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
446                        end do
447                      else
448                        do ir=1,nfft
449                          vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
450                          vtrial1(2*ir  ,1)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,1)
451                        end do
452                      end if
453                    else
454                      if(cplex==1)then
455                        do ir=1,nfft
456                          vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
457                          vtrial1(ir,2)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,2)
458                        end do
459                      else
460 !                      fab: I think there was an error in the definition of  vtrial1(2*ir-1,2); I have corrected it...
461                        do ir=1,nfft
462                          vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
463                          vtrial1(2*ir  ,1)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,1)
464                          vtrial1(2*ir-1,2)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1  ,2)
465                          vtrial1(2*ir  ,2)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,2)
466                        end do
467                      end if
468                    end if
469 
470 !                  Compute the third-order xc energy
471                    call dfptnl_exc3(cplex,exc3,k3xc,mpi_enreg,nk3xc,nfft,nfftot,dtset%nspden,&
472 &                   rho1r1,rho2r1,rho3r1,ucvol,xccc3d1,xccc3d2,xccc3d3)
473 
474 !                  Perform DFPT part of the 3dte calculation
475 
476                    call timab(512,1,tsec)
477                    call pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,&
478 &                   kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,&
479 &                   nspinor,nsppol,npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm)
480                    call timab(512,2,tsec)
481 
482 
483 !                  Describe the perturbation and write out the result
484                    if (mpi_enreg%me == 0) then
485                      if (i2pert < natom + 1) then
486                        write(message,'(a,i3,a,i3)') &
487 &                       ' j2 : displacement of atom ',i2pert,&
488 &                       ' along direction ', i2dir
489                      end if
490                      if (i2pert == dtset%natom + 2) then
491                        write(message,'(a,i4)') &
492 &                       ' j2 : homogeneous electric field along direction ',&
493 &                       i2dir
494                      end if
495                      call wrtout(std_out,message,'COLL')
496                      call wrtout(ab_out,message,'COLL')
497                      write(ab_out,'(20x,a,13x,a)')'real part','imaginary part'
498                      write(ab_out,'(5x,a2,1x,f22.10,3x,f22.10)')'xc',exc3(1)*sixth,zero
499                      if (i2pert == natom + 2) then
500                        write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',&
501 &                       d3_berry(1,i2dir),d3_berry(2,i2dir)
502                      end if
503                      write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'dft',&
504 &                     d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
505 &                     d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
506                      write(ab_out,*)
507                      write(std_out,'(18x,a,11x,a)')'real part','imaginary part'
508                      write(std_out,'(5x,a2,1x,f20.10,3x,f20.10)')'xc',exc3(1)*sixth,zero
509                      write(std_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',&
510 &                     d3_berry(1,i2dir),d3_berry(2,i2dir)
511                      write(std_out,'(5x,a3,f22.10,3x,f22.10)')'dft',&
512 &                     d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
513 &                     d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
514                      write(std_out,*)
515                    end if  ! mpi_enreg%me == 0
516 
517                    d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = &
518 &                   d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + exc3(1)*sixth
519                    d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = &
520 &                   d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + d3_berry(:,i2dir)
521 
522                  end if   !rfpert
523                end do    !i2dir
524              end do     ! i2pert
525 
526            end if   ! rfpert
527          end do    ! i3dir
528        end do     ! i3pert
529 
530      end if   ! rfpert
531    end do    ! i1dir
532  end do     ! i1pert
533 
534 
535  ABI_FREE(cg1)
536  ABI_FREE(cg3)
537  ABI_FREE(eigen1)
538  ABI_FREE(rho1r1)
539  ABI_FREE(rho2r1)
540  ABI_FREE(rho2g1)
541  ABI_FREE(rho3r1)
542  ABI_FREE(atindx1)
543  ABI_FREE(atindx)
544  ABI_FREE(nattyp)
545  ABI_FREE(ph1d)
546  ABI_FREE(ylm)
547  ABI_FREE(vtrial1)
548  ABI_FREE(vxc1)
549  ABI_FREE(vhartr1)
550  ABI_FREE(vpsp1)
551  ABI_FREE(xccc3d1)
552  ABI_FREE(xccc3d2)
553  ABI_FREE(xccc3d3)
554 
555  call timab(502,2,tsec)
556 
557 end subroutine pead_nl_loop

ABINIT/pead_nl_mv [ Functions ]

[ Top ] [ Functions ]

NAME

 pead_nl_mv

FUNCTION

 Compute the finite difference expression of the k-point derivative
 using the PEAD formulation of the third-order energy
 (see Nunes and Gonze PRB 63, 155107 (2001) [[cite:Nunes2001]] Eq. 102)
 and the finite difference formula of Marzari and Vanderbilt
 (see Marzari and Vanderbilt, PRB 56, 12847 (1997) [[cite:Marzari1997]], Appendix B)

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave coefficients of wavefunctions
  cgindex(nkpt2,nsppol) = for each k-point, cgindex stores the location of the WF in the cg array
  cg1 = first-order wavefunction relative to the perturbations i1pert
  cg3 = first-order wavefunction relative to the perturbations i3pert
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  i1pert,i3pert = type of perturbation that has to be computed
  i1dir,i3dir=directions of the corresponding perturbations
  kneigh(30,nkpt2) = index of the neighbours of each k-point
  kg_neigh(30,nkpt2,3) = necessary to construct the vector joining a k-point
                         to its nearest neighbour in case of a single k-point,
                         a line of k-points or a plane of k-points.
                         See getshell.F90 for details
  kptindex(2,nkpt3)= index of the k-points in the reduced BZ
                     related to a k-point in the full BZ
  kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ
  mband = maximum number of bands
  mkmem = maximum number of k points which can fit in core memory
  mkmem_max = maximal number of k-points on each processor (MPI //)
  mk1mem = maximum number of k points for first-order WF which can fit in core memory
  mpi_enreg=MPI-parallelisation information
  mpw   = maximum number of planewaves in basis sphere (large number)
  mvwtk(30,nkpt) = weights to compute the finite difference ddk
  natom = number of atoms in unit cell
  nkpt2 = number of k-points in the reduced part of the BZ
          nkpt2 = nkpt/2 in case of time-reversal symmetry (kptopt = 2)
  nkpt3 = number of k-points in the full BZ
  nneigh = total number of neighbours required to evaluate the finite difference formula
  npwarr(nkpt) = array holding npw for each k point
  nspinor = number of spinorial components of the wavefunctions
  nsppol = number of channels for spin-polarization (1 or 2)
  pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat between k-points

OUTPUT

  d3_berry(2,3) = Berry-phase part of the third-order energy

SIDE EFFECTS

  mpi_enreg=MPI-parallelisation information

NOTES

 For a given set of values of i1pert,i3pert,i1dir and
 i3dir, the routine computes the k-point derivatives for
 12dir = 1,2,3

SOURCE

 920 subroutine pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,&
 921 &                   i1pert,i3pert,i1dir,i3dir,kneigh,kg_neigh,kptindex,&
 922 &                   kpt3,mband,mkmem,mkmem_max,mk1mem,mpi_enreg,&
 923 &                   mpw,mvwtk,natom,nkpt2,nkpt3,nneigh,npwarr,nspinor,&
 924 &                   nsppol,pwind)
 925 
 926  use m_hide_lapack, only : dzgedi, dzgefa
 927 
 928 !Arguments ------------------------------------
 929 !
 930 !---  Arguments : integer scalars
 931  integer, intent(in) :: i1dir,i1pert,i3dir,i3pert,mband,mk1mem
 932  integer, intent(in) :: mkmem,mkmem_max,mpw,natom
 933  integer, intent(in) :: nkpt2,nkpt3,nneigh,nspinor,nsppol
 934 !
 935 !---  Arguments : integer arrays
 936  integer, intent(in) :: cgindex(nkpt2,nsppol)
 937  integer, intent(in) :: kneigh(30,nkpt2),kg_neigh(30,nkpt2,3),kptindex(2,nkpt3)
 938  integer, intent(in) :: npwarr(nkpt2),pwind(mpw,nneigh,mkmem)
 939 !
 940 !---  Arguments : real(dp) scalars
 941 !
 942 !---  Arguments : real(dp) arrays
 943  real(dp), intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
 944  real(dp), intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
 945  real(dp), intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
 946  real(dp), intent(in) :: gmet(3,3),kpt3(3,nkpt3)
 947  real(dp), intent(in) :: mvwtk(30,nkpt2)
 948  real(dp), intent(out) :: d3_berry(2,3)
 949 !
 950 !---  Arguments : structured datatypes
 951  type(MPI_type), intent(in) :: mpi_enreg
 952  type(datafiles_type), intent(in) :: dtfil
 953  type(dataset_type), intent(in) :: dtset
 954 
 955 !Local variables-------------------------------
 956 !
 957 !---- Local variables : integer scalars
 958  integer :: count,counter,count1,iband,icg
 959  integer :: ierr,ii,ikpt,ikpt_loc,ikpt2
 960  integer :: ikpt_rbz,ineigh,info,ipw,isppol,jband,jcg,jj,jkpt,job,jpw, jkpt2, jkpt_rbz
 961  integer :: lband,lpband,nband_occ,npw_k,npw_k1,my_source,his_source,dest,tag
 962  integer :: spaceComm
 963  integer,parameter :: level=52
 964  integer :: bdtot_index
 965 !
 966 !---- Local variables : integer arrays
 967  integer,allocatable :: ipvt(:)
 968  integer, allocatable :: bd_index(:,:)
 969 !
 970 !---- Local variables : real(dp) scalars
 971  real(dp) :: dotnegi,dotnegr,dotposi,dotposr
 972 ! real(dp) :: c1,c2 ! appear commented out below
 973 !
 974 !---- Local variables : real(dp) arrays
 975  real(dp) :: d3_aux(2,3),det(2,2),dk(3),dk_(3)
 976  real(dp) :: z1(2),z2(2)
 977  real(dp),allocatable :: buffer(:,:),cgq(:,:),cg1q(:,:),cg3q(:,:)
 978  real(dp),allocatable :: qmat(:,:,:),s13mat(:,:,:),s1mat(:,:,:),s3mat(:,:,:)
 979  real(dp),allocatable :: smat(:,:,:),zgwork(:,:)
 980 !
 981 !---- Local variables : character variables
 982  character(len=500) :: message
 983 !
 984 !---- Local variables : structured datatypes
 985 
 986 
 987 #if defined HAVE_MPI
 988 integer :: status1(MPI_STATUS_SIZE)
 989 spaceComm=mpi_enreg%comm_cell
 990 #endif
 991 
 992  ABI_UNUSED(dtfil%ireadwf)
 993 
 994 ! ***********************************************************************
 995 
 996  write(message,'(8a)') ch10,&
 997 & ' pead_nl_mv : finite difference expression of the k-point derivative',ch10,&
 998 & '           is performed using the PEAD formulation of ',&
 999 & 'the third-order energy',ch10,&
1000 & '           (see Nunes and Gonze PRB 63, 155107 (2001) [[cite:Nunes2001]] Eq. 102)',ch10
1001 !call wrtout(ab_out,message,'COLL')
1002  call wrtout(std_out,  message,'COLL')
1003 
1004 
1005 !fab: I think that the following restriction must be eliminated:
1006 !isppol = 1
1007 
1008  ikpt_loc = 0
1009  d3_aux(:,:) = 0_dp
1010 
1011  ABI_MALLOC(s13mat,(2,mband,mband))
1012  ABI_MALLOC(smat,(2,mband,mband))
1013  ABI_MALLOC(s1mat,(2,mband,mband))
1014  ABI_MALLOC(qmat,(2,mband,mband))
1015  ABI_MALLOC(ipvt,(mband))
1016  ABI_MALLOC(s3mat,(2,mband,mband))
1017  ABI_MALLOC(zgwork,(2,mband))
1018  ABI_MALLOC(bd_index, (nkpt2, nsppol))
1019 
1020  bdtot_index = 0
1021  do isppol = 1, nsppol
1022    do ikpt_rbz = 1, nkpt2
1023      bd_index(ikpt_rbz,isppol) = bdtot_index
1024      bdtot_index = bdtot_index + dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1025    end do
1026  end do
1027 
1028 !fab: I think here I have to add the loop over spin
1029 
1030  do isppol = 1, nsppol
1031 
1032 !  Loop over k-points
1033 !  COMMENT: Every processor has to make mkmem_max iterations
1034 !  even if mkmem < mkemem_max. This is due to the fact
1035 !  that it still has to communicate its wavefunctions
1036 !  to other processors even if it has no more overlap
1037 !  matrices to compute.
1038 
1039    ikpt_loc = 0 ; ikpt = 0
1040 
1041    do while (ikpt_loc < mkmem_max)
1042 
1043      if (ikpt_loc < mkmem) ikpt = ikpt + 1
1044 
1045      if (xmpi_paral == 1) then
1046 !      if ((minval(abs(mpi_enreg%proc_distrb(ikpt,1:mband,1:dtset%nsppol) &
1047 !      &       - mpi_enreg%me)) /= 0).and.(ikpt_loc < mkmem)) cycle
1048        if(ikpt>nkpt2)then
1049          ikpt_loc=mkmem_max
1050          cycle
1051        end if
1052        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me)) then
1053          if(ikpt==nkpt2) ikpt_loc=mkmem_max
1054          cycle
1055        end if
1056      end if
1057 
1058      ikpt_loc = ikpt_loc + 1
1059      npw_k = npwarr(ikpt)
1060      counter = 100*ikpt
1061 
1062      ii = cgindex(ikpt,isppol)
1063 
1064 !    Loop on the  neighbours
1065 
1066      do ineigh = 1,nneigh
1067 
1068        s13mat(:,:,:) = zero
1069        smat(:,:,:) = zero
1070        s1mat(:,:,:) = zero
1071        s3mat(:,:,:) = zero
1072        qmat(:,:,:) = zero
1073 
1074        ikpt2  = kneigh(ineigh,ikpt)
1075        ikpt_rbz = kptindex(1,ikpt2)   ! index of the k-point in the reduced BZ
1076        jj = cgindex(ikpt_rbz,isppol)
1077        ! previous fixed value for nband_k now called nband_occ:
1078        !nband_occ = dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1079        ! TODO: check if all these bands are occupied in nsppol = 2 case
1080        nband_occ = 0
1081        do iband = 1, dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1082          !Note, only one image is allowed here (or occ_orig should be the same or all images)
1083          if (dtset%occ_orig(bd_index(ikpt_rbz,isppol) + iband,1) > tol10) nband_occ = nband_occ + 1
1084        end do
1085        npw_k1 = npwarr(ikpt_rbz)
1086        dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt)
1087        dk(:)  = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp)
1088 
1089        count = nspinor*mband*npw_k1
1090        ABI_MALLOC(cgq,(2,count))
1091        ABI_MALLOC(cg1q,(2,count))
1092        ABI_MALLOC(cg3q,(2,count))
1093 
1094 #if defined HAVE_MPI
1095 
1096        my_source = mpi_enreg%proc_distrb(ikpt_rbz,1,1)
1097 
1098 !      do dest = 0, mpi_enreg%nproc-1
1099        do dest = 0, maxval(mpi_enreg%proc_distrb(1:nkpt2,1:mband,1:dtset%nsppol))
1100 
1101          if ((dest==mpi_enreg%me).and.(ikpt_loc <= mkmem)) then
1102 !          I am dest and have something to do
1103 
1104            if (my_source == mpi_enreg%me) then
1105 !            I am destination and source
1106              jcg = cgindex(ikpt_rbz,isppol)
1107 
1108              cgq(:,1:count)  = cg(:,jcg+1:jcg+count)
1109              cg1q(:,1:count) = cg1(:,jcg+1:jcg+count)
1110              cg3q(:,1:count) = cg3(:,jcg+1:jcg+count)
1111 
1112            else
1113 !            I am the destination but not the source -> receive
1114 
1115              tag = ikpt_rbz
1116 
1117              ABI_MALLOC(buffer,(2,3*count))
1118 
1119              call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,my_source,tag,spaceComm,status1,ierr)
1120 
1121              cgq(:,1:count)  = buffer(:,1:count)
1122              cg1q(:,1:count) = buffer(:,count+1:2*count)
1123              cg3q(:,1:count) = buffer(:,2*count+1:3*count)
1124              ABI_FREE(buffer)
1125 
1126            end if
1127 
1128          else if (ikpt_loc <= mpi_enreg%mkmem(dest)) then  ! dest != me and the dest has a k-point to treat
1129 
1130            jkpt=mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,1)
1131            jkpt2  = kneigh(ineigh,jkpt)
1132            jkpt_rbz = kptindex(1,jkpt2)   ! index of the k-point in the reduced BZ
1133 
1134            his_source = mpi_enreg%proc_distrb(jkpt_rbz,1,1)
1135 
1136            if (his_source == mpi_enreg%me) then
1137 
1138              jcg = cgindex(jkpt_rbz,isppol)
1139 
1140              tag = jkpt_rbz
1141              count1 = npwarr(jkpt_rbz)*mband*nspinor
1142              ABI_MALLOC(buffer,(2,3*count1))
1143              buffer(:,1:count1)            = cg(:,jcg+1:jcg+count1)
1144              buffer(:,count1+1:2*count1)   = cg1(:,jcg+1:jcg+count1)
1145              buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1)
1146 
1147              call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,dest,tag,spaceComm,ierr)
1148 
1149              ABI_FREE(buffer)
1150 
1151            end if
1152 
1153          end if
1154 
1155        end do
1156 !
1157 !      do jkpt = 1, nkpt2
1158 !
1159 !      if ((jkpt == ikpt_rbz).and.(source /= mpi_enreg%me).and.&
1160 !      &         (ikpt_loc <= mkmem)) then
1161 !
1162 !      tag = jkpt
1163 !
1164 !      allocate(buffer(2,3*count))
1165 !      call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,&
1166 !      source,tag,spaceComm,status1,ierr)
1167 !
1168 !      cgq(:,1:count)  = buffer(:,1:count)
1169 !      cg1q(:,1:count) = buffer(:,count+1:2*count)
1170 !      cg3q(:,1:count) = buffer(:,2*count+1:3*count)
1171 !      deallocate(buffer)
1172 !
1173 !      end if
1174 !
1175 !      !        ----------------------------------------------------------------------------
1176 !      !        --------------- Here: send the WF to all the cpus that need it -------------
1177 !      !        ----------------------------------------------------------------------------
1178 !
1179 !      do dest = 1, mpi_enreg%nproc
1180 !
1181 !      if ((minval(abs(mpi_enreg%proc_distrb(jkpt,1:mband,1:dtset%nsppol) &
1182 !      &           - mpi_enreg%me)) == 0).and.&
1183 !      &           (mpi_enreg%kptdstrb(dest,ineigh,ikpt_loc) == jkpt)) then
1184 !
1185 !
1186 !
1187 !      jcg = cgindex(jkpt,isppol)
1188 !
1189 !      if (((dest-1) == mpi_enreg%me)) then
1190 !
1191 !      cgq(:,1:count)  = cg(:,jcg+1:jcg+count)
1192 !      cg1q(:,1:count) = cg1(:,jcg+1:jcg+count)
1193 !      cg3q(:,1:count) = cg3(:,jcg+1:jcg+count)
1194 !
1195 !      else
1196 !
1197 !      tag = jkpt
1198 !      count1 = npwarr(jkpt)*mband*nspinor
1199 !      allocate(buffer(2,3*count1))
1200 !      buffer(:,1:count1)            = cg(:,jcg+1:jcg+count1)
1201 !      buffer(:,count1+1:2*count1)   = cg1(:,jcg+1:jcg+count1)
1202 !      buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1)
1203 !
1204 !      call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,(dest-1),tag,spaceComm,ierr)
1205 !
1206 !      deallocate(buffer)
1207 !
1208 !      end if
1209 !
1210 !      end if
1211 !
1212 !      end do          ! loop over dest
1213 !
1214 !      end do          ! loop over jkpt
1215 
1216        if (ikpt_loc > mkmem) then
1217          ABI_FREE(cgq)
1218          ABI_FREE(cg1q)
1219          ABI_FREE(cg3q)
1220          cycle
1221        end if
1222 
1223 #else
1224 !      no // over k-points
1225 
1226        cgq(:,1:count)  = cg(:,jj+1:jj+count)
1227        cg1q(:,1:count) = cg1(:,jj+1:jj+count)
1228        cg3q(:,1:count) = cg3(:,jj+1:jj+count)
1229 
1230 #endif
1231 
1232 !      Compute overlap matrices
1233 
1234        if (kptindex(2,ikpt2) == 0) then  ! no time-reversal symmetry
1235 
1236          do ipw = 1, npw_k
1237 
1238            jpw = pwind(ipw,ineigh,ikpt_loc)
1239            if (jpw /= 0) then
1240 
1241              do iband = 1, nband_occ
1242                do jband = 1, nband_occ
1243 
1244                  icg = ii + (iband-1)*npw_k + ipw
1245                  jcg = (jband-1)*npw_k1 + jpw
1246 
1247                  smat(1,iband,jband) = smat(1,iband,jband) + &
1248 &                 cg(1,icg)*cgq(1,jcg) + cg(2,icg)*cgq(2,jcg)
1249                  smat(2,iband,jband) = smat(2,iband,jband) + &
1250 &                 cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg)
1251 
1252                  s13mat(1,iband,jband) = s13mat(1,iband,jband) + &
1253 &                 cg1(1,icg)*cg3q(1,jcg) + cg1(2,icg)*cg3q(2,jcg)
1254                  s13mat(2,iband,jband) = s13mat(2,iband,jband) + &
1255 &                 cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg)
1256 
1257                  s1mat(1,iband,jband) = s1mat(1,iband,jband) + &
1258 &                 cg1(1,icg)*cgq(1,jcg) + cg1(2,icg)*cgq(2,jcg) + &
1259 &                 cg(1,icg)*cg1q(1,jcg) + cg(2,icg)*cg1q(2,jcg)
1260                  s1mat(2,iband,jband) = s1mat(2,iband,jband) + &
1261 &                 cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) + &
1262 &                 cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg)
1263 
1264                  s3mat(1,iband,jband) = s3mat(1,iband,jband) + &
1265 &                 cg3(1,icg)*cgq(1,jcg) + cg3(2,icg)*cgq(2,jcg) + &
1266 &                 cg(1,icg)*cg3q(1,jcg) + cg(2,icg)*cg3q(2,jcg)
1267                  s3mat(2,iband,jband) = s3mat(2,iband,jband) + &
1268 &                 cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) + &
1269 &                 cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg)
1270 
1271                end do
1272              end do
1273 
1274            end if
1275 
1276          end do   ! ipw
1277 
1278        else                              ! use time-reversal symmetry
1279 
1280          do ipw = 1,npw_k
1281 
1282            jpw = pwind(ipw,ineigh,ikpt_loc)
1283            if (jpw /= 0) then
1284 
1285              do iband = 1, nband_occ
1286                do jband = 1, nband_occ
1287 
1288                  icg = ii + (iband-1)*npw_k + ipw
1289                  jcg = (jband-1)*npw_k1 + jpw
1290 
1291                  smat(1,iband,jband) = smat(1,iband,jband) + &
1292 &                 cg(1,icg)*cgq(1,jcg) - cg(2,icg)*cgq(2,jcg)
1293                  smat(2,iband,jband) = smat(2,iband,jband) - &
1294 &                 cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg)
1295 
1296                  s13mat(1,iband,jband) = s13mat(1,iband,jband) + &
1297 &                 cg1(1,icg)*cg3q(1,jcg) - cg1(2,icg)*cg3q(2,jcg)
1298                  s13mat(2,iband,jband) = s13mat(2,iband,jband) - &
1299 &                 cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg)
1300 
1301                  s1mat(1,iband,jband) = s1mat(1,iband,jband) + &
1302 &                 cg1(1,icg)*cgq(1,jcg) - cg1(2,icg)*cgq(2,jcg) + &
1303 &                 cg(1,icg)*cg1q(1,jcg) - cg(2,icg)*cg1q(2,jcg)
1304                  s1mat(2,iband,jband) = s1mat(2,iband,jband) - &
1305 &                 cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) - &
1306 &                 cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg)
1307 
1308                  s3mat(1,iband,jband) = s3mat(1,iband,jband) + &
1309 &                 cg3(1,icg)*cgq(1,jcg) - cg3(2,icg)*cgq(2,jcg) + &
1310 &                 cg(1,icg)*cg3q(1,jcg) - cg(2,icg)*cg3q(2,jcg)
1311                  s3mat(2,iband,jband) = s3mat(2,iband,jband) - &
1312 &                 cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) - &
1313 &                 cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg)
1314 
1315                end do
1316              end do
1317 
1318            end if
1319 
1320          end do   ! ipw
1321 
1322        end if
1323 
1324        ABI_FREE(cgq)
1325        ABI_FREE(cg1q)
1326        ABI_FREE(cg3q)
1327 
1328 !      Compute qmat, the inverse of smat
1329 
1330        job = 1  ! compute inverse only
1331        qmat(:,:,:) = smat(:,:,:)
1332 
1333        call dzgefa(qmat,mband,nband_occ,ipvt,info)
1334        call dzgedi(qmat,mband,nband_occ,ipvt,det,zgwork,job)
1335 
1336 !      DEBUG
1337 !      write(100,*)
1338 !      write(100,*)'ikpt = ',ikpt,'ineigh = ',ineigh
1339 !      do iband = 1,nband_occ
1340 !      do jband = 1,nband_occ
1341 !      c1 = 0_dp ; c2 = 0_dp
1342 !      do lband = 1,nband_occ
1343 !      c1 = c1 + smat(1,iband,lband)*qmat(1,lband,jband) - &
1344 !      &           smat(2,iband,lband)*qmat(2,lband,jband)
1345 !      c2 = c2 + smat(1,iband,lband)*qmat(2,lband,jband) + &
1346 !      &           smat(2,iband,lband)*qmat(1,lband,jband)
1347 !      end do
1348 !      write(100,'(2(2x,i2),2(2x,f16.9))')iband,jband,&
1349 !      & c1,c2
1350 !      end do
1351 !      end do
1352 !      ENDDEBUG
1353 
1354 
1355 
1356 !      Accumulate sum over bands
1357 
1358        dotposr = 0_dp ; dotposi = 0_dp
1359        dotnegr = 0_dp ; dotnegi = 0_dp
1360        do iband = 1, nband_occ
1361          do jband = 1, nband_occ
1362 
1363            dotposr = dotposr + &
1364 &           s13mat(1,iband,jband)*qmat(1,jband,iband) - &
1365 &           s13mat(2,iband,jband)*qmat(2,jband,iband)
1366            dotposi = dotposi + &
1367 &           s13mat(1,iband,jband)*qmat(2,jband,iband) + &
1368 &           s13mat(2,iband,jband)*qmat(1,jband,iband)
1369 
1370 
1371            do lband = 1, nband_occ
1372              do lpband= 1, nband_occ
1373 
1374                z1(1) = s1mat(1,iband,jband)*qmat(1,jband,lband) - &
1375 &               s1mat(2,iband,jband)*qmat(2,jband,lband)
1376                z1(2) = s1mat(1,iband,jband)*qmat(2,jband,lband) + &
1377 &               s1mat(2,iband,jband)*qmat(1,jband,lband)
1378 
1379                z2(1) = s3mat(1,lband,lpband)*qmat(1,lpband,iband) - &
1380 &               s3mat(2,lband,lpband)*qmat(2,lpband,iband)
1381                z2(2) = s3mat(1,lband,lpband)*qmat(2,lpband,iband) + &
1382 &               s3mat(2,lband,lpband)*qmat(1,lpband,iband)
1383 
1384                dotnegr = dotnegr + &
1385 &               z1(1)*z2(1) - z1(2)*z2(2)
1386                dotnegi = dotnegi + &
1387 &               z1(1)*z2(2) + z1(2)*z2(1)
1388 
1389              end do   ! lpband
1390            end do   ! lband
1391 
1392          end do   ! jband
1393        end do   ! iband
1394 
1395        d3_aux(1,:) = d3_aux(1,:) + &
1396 &       dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposr-dotnegr)
1397        d3_aux(2,:) = d3_aux(2,:) + &
1398 &       dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposi-dotnegi)
1399 
1400      end do        ! End loop over neighbours
1401 
1402 
1403    end do      ! End loop over k-points
1404 
1405  end do  ! fab: end loop over spin
1406 
1407 
1408 
1409 
1410  call xmpi_sum(d3_aux,spaceComm,ierr)
1411 
1412 
1413  ABI_FREE(s13mat)
1414  ABI_FREE(smat)
1415  ABI_FREE(s1mat)
1416  ABI_FREE(qmat)
1417  ABI_FREE(ipvt)
1418  ABI_FREE(s3mat)
1419  ABI_FREE(zgwork)
1420  ABI_FREE(bd_index)
1421 
1422 
1423 !fab: I think that in the following we have to make a distinction:
1424 !for the spin unpolarized case we leave the PEAD expression as it is, while
1425 !in the spin polarized case we have simply to divide by 2
1426 !(see eq.19 di PRB 63,155107 [[cite:Nunes2001]], eq. 7 di PRB 71,125107 [[cite:Veithen2005]]
1427 ! and eq 13 di PRB 71, 125107 [[cite:Veithen2005]] ...
1428 !in this latter equation the 2 must be simply replaced by the sum over the spin components...
1429 !and indeed we have inserted the loop over the spin,
1430 !but there was a factor 2 already present in the routine due to spin degenracy that had to be removed)
1431 
1432 
1433  if (nsppol==1) then
1434 
1435 !  Take minus the imaginary part
1436 
1437    d3_berry(1,:) = -1_dp*d3_aux(2,:)
1438    d3_berry(2,:) = d3_aux(1,:)
1439 
1440    d3_berry(2,:) = 0_dp
1441 
1442  else
1443 
1444    d3_berry(1,:) = -1_dp*d3_aux(2,:)/2._dp
1445    d3_berry(2,:) = d3_aux(1,:)/2._dp
1446 
1447    d3_berry(2,:) = 0_dp/2._dp
1448 
1449  end if
1450 
1451 !DEBUG
1452 !write(100,*)'pead_nl_mv.f : d3_berry'
1453 !write(100,*)'Perturbation',i1dir,i3dir
1454 !write(100,*)
1455 !write(100,*)'before transformation'
1456 !write(100,*)'real part'
1457 !write(100,'(3(2x,f20.9))')d3_berry(1,:)
1458 !write(100,*)
1459 !write(100,*)'imaginary part'
1460 !write(100,'(3(2x,f20.9))')d3_berry(2,:)
1461 !write(100,*)
1462 !write(100,*)'after transformation'
1463 !ENDDEBUG
1464 
1465 !Compute the projection on the basis vectors of
1466 !reciprocal space
1467 
1468  d3_aux(:,:) = 0_dp
1469  do ii = 1,3
1470    do jj = 1,3
1471      d3_aux(:,ii) = d3_aux(:,ii) + gmet(ii,jj)*d3_berry(:,jj)
1472    end do
1473  end do
1474  d3_berry(:,:) = d3_aux(:,:)
1475 
1476 !Write out the berryphase part of the third order energy
1477 
1478  if (mpi_enreg%me == 0) then
1479 
1480    write(message,'(a,a,a)')ch10,' Berryphase part of the third-order energy:',ch10
1481    call wrtout(std_out,  message,'COLL')
1482 
1483    if (i1pert < natom + 1) then
1484      write(message,'(a,i3,a,i3)')&
1485 &     '            j1: Displacement of atom ',i1pert,&
1486 &     ' along direction ',i1dir
1487    else if (i1pert == natom + 2) then
1488      write(message,'(a,i3)')&
1489 &     '            j1: homogenous electric field along direction ',i1dir
1490    end if
1491    call wrtout(std_out,  message,'COLL')
1492 
1493    write(message,'(a)')&
1494 &   '            j2: k-point derivative along direction i2dir '
1495    call wrtout(std_out,  message,'COLL')
1496 
1497    if (i3pert < natom + 1) then
1498      write(message,'(a,i3,a,i3,a)')&
1499 &     '            j3: Displacement of atom ',i3pert,&
1500 &     ' along direction ',i3dir,ch10
1501    else if (i3pert == natom + 2) then
1502      write(message,'(a,i3,a)')&
1503 &     '            j3: homogenous electric field along direction ',i3dir,ch10
1504    end if
1505    call wrtout(std_out,  message,'COLL')
1506 
1507 !  write(ab_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part'
1508    write(std_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part'
1509    do ii = 1,3
1510      write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,&
1511 &     d3_berry(1,ii),d3_berry(2,ii)
1512      write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,&
1513 &     d3_berry(1,ii),d3_berry(2,ii)
1514    end do
1515 
1516  end if    ! mpi_enreg%me == 0
1517 
1518 !DEBUG
1519 !write(100,*)'real part'
1520 !write(100,'(3(2x,f20.9))')d3_berry(1,:)
1521 !write(100,*)
1522 !write(100,*)'imaginary part'
1523 !write(100,'(3(2x,f20.9))')d3_berry(2,:)
1524 !ENDDEBUG
1525 
1526 end subroutine pead_nl_mv

ABINIT/pead_nl_resp [ Functions ]

[ Top ] [ Functions ]

NAME

 pead_nl_resp

FUNCTION

 Compute the linear response part to the 3dte

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol) = array for planewave
                                          coefficients of wavefunctions
  cg1 = first-order wavefunction relative to the perturbations i1pert
  cg3 = first-order wavefunction relative to the perturbations i3pert
  cplex= if 1, real space 1-order functions on FFT grid are REAL,
          if 2, COMPLEX
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  i1dir,i2dir,i3dir=directions of the corresponding perturbations
  i1pert,i2pert,i3pert = type of perturbation that has to be computed
  kg(3,mpw*mkmem)=reduced planewave coordinates
  mband = maximum number of bands
  mgfft=maximum size of 1D FFTs
  mkmem = maximum number of k points which can fit in core memory
  mk1mem = maximum number of k points for first-order WF
           which can fit in core memory
  mpert =maximum number of ipert
  mpi_enreg=MPI-parallelisation information
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mpw   = maximum number of planewaves in basis sphere (large number)
  natom = number of atoms in unit cell
  nfft  = (effective) number of FFT grid points (for this processor)
  nkpt  = number of k points
  nspden = number of spin-density components
  nspinor = number of spinorial components of the wavefunctions
  nsppol = number of channels for spin-polarization (1 or 2)
  npwarr(nkpt) = array holding npw for each k point
  occ(mband*nkpt*nsppol) = occupation number for each band and k
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  psps <type(pseudopotential_type)> = variables related to pseudopotentials
  rprimd(3,3) = dimensional primitive translations (bohr)
  vtrial1(cplex*nfft,nspden)=firs-order local potential
  xred(3,natom) = reduced atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= spherical harmonics for
       each G and k point

OUTPUT

  d3lo(2,3,mpert,3,mpert,3,mpert) = matrix of the 3DTEs

SOURCE

609 subroutine pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,&
610 & i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,&
611 & kg,mband,mgfft,mkmem,mk1mem,&
612 & mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,nspinor,nsppol,&
613 & npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm)
614 
615 !Arguments ------------------------------------
616 !scalars
617  integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft
618  integer,intent(in) :: mk1mem,mkmem,mpert,mpsang,mpw,natom,nfft,nkpt,nspden
619  integer,intent(in) :: nspinor,nsppol
620  type(MPI_type),intent(in) :: mpi_enreg
621  type(datafiles_type),intent(in) :: dtfil
622  type(dataset_type),intent(in) :: dtset
623  type(pseudopotential_type),intent(in) :: psps
624 !arrays
625  integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
626  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
627  real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
628  real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
629  real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom),rprimd(3,3)
630  real(dp),intent(in) :: xred(3,natom),ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
631  real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden)
632  real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert)
633 
634 !Local variables-------------------------------
635 !scalars
636  integer,parameter :: level=52
637  integer :: bantot,choice,counter,cpopt,dimffnl,iband,icg0,ider,ierr
638  integer :: ii,ikg,ikpt,ilm,ipw,isppol,istwf_k,jband,jj
639  integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nnlout,npw_k
640  integer :: option,paw_opt,signs,spaceComm,tim_fourwf,tim_nonlop
641  real(dp) :: dot1i,dot1r,dot2i,dot2r,doti,dotr,lagi,lagr,sumi,sumr,weight
642  type(gs_hamiltonian_type) :: gs_hamk
643 !arrays
644  integer,allocatable :: kg_k(:,:)
645  real(dp) :: buffer(2),enlout(3),kpq(3),kpt(3)
646  real(dp) :: dum_svectout(1,1),dum(1),rmet(3,3),ylmgr_dum(1,1,1)
647  real(dp),allocatable :: cwave0(:,:),cwavef3(:,:),ffnlk(:,:,:,:)
648  real(dp),allocatable :: gh0(:,:),gh1(:,:),gvnl(:,:),kpg_k(:,:)
649  real(dp),allocatable :: vlocal1(:,:,:),wfraug(:,:,:,:),ylm_k(:,:)
650  type(pawcprj_type) :: cprj_dum(1,1)
651  type(pawtab_type) :: pawtab_dum(0)
652 
653 !***********************************************************************
654 
655  ABI_UNUSED(dtfil%ireadwf)
656 
657  me = mpi_enreg%me
658  spaceComm=mpi_enreg%comm_cell
659 
660  bantot = 0
661  icg0 = 0
662  option = 2
663  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
664  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
665 
666  ABI_MALLOC(vlocal1,(cplex*n4,n5,n6))
667  ABI_MALLOC(wfraug,(2,n4,n5,n6))
668 
669 !Initialize Hamiltonian (k-independent terms) - NCPP only
670  call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,&
671 & dtset%typat,xred,nfft,mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
672 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
673 !& paw_ij=paw_ij)
674  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
675 
676  sumr = zero ; sumi = zero
677 
678 !Loop over spins
679 
680  do isppol = 1, nsppol
681 
682    call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,vtrial1,vlocal1,option)
683 
684 !  Loop over k-points
685 
686    ikg = 0
687    do ikpt = 1, nkpt
688 
689      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me))cycle
690 
691      counter = 100*ikpt
692 
693      nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
694      npw_k = npwarr(ikpt)
695      istwf_k = dtset%istwfk(ikpt)
696 
697      kpt(:) = dtset%kptns(:,ikpt)
698      kpq(:) = dtset%kptns(:,ikpt) ! In case of non zero q, kpt = kpt + q
699 
700      ABI_MALLOC(cwave0,(2,npw_k*dtset%nspinor))
701      ABI_MALLOC(cwavef3,(2,npw_k*dtset%nspinor))
702      ABI_MALLOC(gh0,(2,npw_k*dtset%nspinor))
703      ABI_MALLOC(gvnl,(2,npw_k*dtset%nspinor))
704      ABI_MALLOC(gh1,(2,npw_k*dtset%nspinor))
705 
706      ABI_MALLOC(kg_k,(3,npw_k))
707      ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
708      kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg)
709      if (psps%useylm==1) then
710        do ilm=1,mpsang*mpsang
711          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
712        end do
713      end if
714 
715 !    Compute (k+G) and (k+q+G) vectors (only if useylm=1)
716      nkpg=0;if (i2pert<natom+1) nkpg=3*dtset%nloalg(3)
717      ABI_MALLOC(kpg_k,(npw_k,nkpg))
718      if (nkpg>0) then
719        call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k)
720      end if
721 
722 !    Compute nonlocal form factors ffnl at (k+G), for all atoms
723      dimffnl=1
724      ABI_MALLOC(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
725      if (i2pert<natom+1) then
726        ider=0
727        call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
728 &       ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,&
729 &       psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,&
730 &       psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
731      end if
732 
733 !    Load k-dependent part in the Hamiltonian datastructure
734      call gs_hamk%load_k(kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,&
735 &     kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk,compute_gbound=.true.)
736 !    Load k+q-dependent part in the Hamiltonian datastructure
737 !    call load_kprime_hamiltonian...  !! To be activated when q/=0
738 
739 !    Loop over bands
740 
741      do iband = 1,nband_k
742 
743        cwave0(:,:)=cg(:,1+(iband - 1)*npw_k*dtset%nspinor+icg0:&
744 &       iband*npw_k*dtset%nspinor+icg0)
745        cwavef3(:,:)=cg3(:,1+(iband-1)*npw_k*dtset%nspinor+icg0:&
746 &       iband*npw_k*dtset%nspinor+icg0)
747 
748 !      Compute vtrial1 | cwafef3 >
749        tim_fourwf = 0 ; weight = one
750        call fourwf(cplex,vlocal1,cwavef3,gh1,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
751 &       istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,&
752 &       tim_fourwf,weight,weight,&
753 &       use_gpu_cuda=dtset%use_gpu_cuda)
754 
755 !      In case i2pert = phonon-type perturbation
756 !      add first-order change in the nonlocal potential
757        if (i2pert<natom+1) then
758          signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1
759          call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,&
760 &         signs,dum_svectout,tim_nonlop,cwavef3,gvnl,iatom_only=i2pert)
761          gh1(:,:) = gh1(:,:) + gvnl(:,:)
762        end if
763 
764        ii = (iband-1)*npw_k*dtset%nspinor + icg0
765        call dotprod_g(dotr,doti,istwf_k,npw_k,2,cg1(:,ii+1:ii+npw_k),gh1,mpi_enreg%me_g0,xmpi_comm_self)
766 
767 !      Compute vtrial1 | cwave0 >
768        tim_fourwf = 0 ; weight = one
769        call fourwf(cplex,vlocal1,cwave0,gh0,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
770 &       istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,&
771 &       tim_fourwf,weight,weight,use_gpu_cuda=dtset%use_gpu_cuda)
772 
773 !      In case i2pert = phonon-type perturbation
774 !      add first-order change in the nonlocal potential
775        if (i2pert<natom+1) then
776          signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1
777          call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,&
778 &         signs,dum_svectout,tim_nonlop,cwave0,gvnl,iatom_only=i2pert)
779          gh0(:,:) = gh0(:,:) + gvnl(:,:)
780        end if
781 
782 !      Compute the dft contribution to the Lagrange multiplier
783 !      cwavef3 and cwave0 have been transferred to gh1 and gh0
784 !      these vectors will be used to store the wavefunctions of band iband
785 !      cg1 and gh0 contain the wavefunctions of band jband
786 
787        lagr = zero ; lagi = zero
788        do jband = 1, nband_k
789 
790          ii = (jband - 1)*npw_k*dtset%nspinor + icg0
791          jj = (iband - 1)*npw_k*dtset%nspinor + icg0
792 
793 !        dot1r and dot1i contain < u_mk | v^(1) | u_nk >
794 !        dot2r and dot2i contain < u_nk^(1) | u_mk^(1) >
795 !        m -> jband and n -> iband
796 
797          dot1r = zero ; dot1i = zero
798          dot2r = zero ; dot2i = zero
799          do ipw = 1, npw_k
800            ii = ii + 1 ; jj = jj + 1
801            dot1r = dot1r + cg(1,ii)*gh0(1,ipw) + cg(2,ii)*gh0(2,ipw)
802            dot1i = dot1i + cg(1,ii)*gh0(2,ipw) - cg(2,ii)*gh0(1,ipw)
803            dot2r = dot2r + cg1(1,jj)*cg3(1,ii) + &
804 &           cg1(2,jj)*cg3(2,ii)
805            dot2i = dot2i + cg1(1,jj)*cg3(2,ii) - &
806 &           cg1(2,jj)*cg3(1,ii)
807          end do  !  ipw
808 
809          lagr = lagr + dot1r*dot2r - dot1i*dot2i
810          lagi = lagi + dot1r*dot2i + dot1i*dot2r
811 
812        end do    ! jband
813 
814        sumr = sumr + &
815 &       dtset%wtk(ikpt)*occ(bantot+iband)*(dotr-lagr)
816        sumi = sumi + &
817 &       dtset%wtk(ikpt)*occ(bantot+iband)*(doti-lagi)
818 
819      end do   ! end loop over bands
820 
821      bantot = bantot + nband_k
822      icg0 = icg0 + npw_k*dtset%nspinor*nband_k
823      ikg = ikg + npw_k
824 
825      ABI_FREE(cwave0)
826      ABI_FREE(cwavef3)
827      ABI_FREE(gh0)
828      ABI_FREE(gh1)
829      ABI_FREE(gvnl)
830      ABI_FREE(kg_k)
831      ABI_FREE(ylm_k)
832      ABI_FREE(ffnlk)
833      ABI_FREE(kpg_k)
834 
835    end do   ! end loop over k-points
836 
837  end do   ! end loop over spins
838 
839  if (xmpi_paral == 1) then
840    buffer(1) = sumr ; buffer(2) = sumi
841    call xmpi_sum(buffer,spaceComm,ierr)
842    sumr = buffer(1) ; sumi = buffer(2)
843  end if
844 
845 
846  d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
847 !d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
848 
849 !In some cases, the imaginary part is /= 0 because of the
850 !use of time reversal symmetry
851  d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero
852 
853  call gs_hamk%free()
854 
855  ABI_FREE(vlocal1)
856  ABI_FREE(wfraug)
857 
858 end subroutine pead_nl_resp