TABLE OF CONTENTS


ABINIT/m_pead_nl_loop [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pead_nl_loop

FUNCTION

COPYRIGHT

  Copyright (C) 2002-2016 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 .

PARENTS

CHILDREN

SOURCE

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

PARENTS

      nonlinear

CHILDREN

      appdig,dfpt_mkcore,dfpt_mkvxc,dfpt_vlocal,pead_nl_mv,pead_nl_resp
      dotprod_vn,fourdp,getph,hartre,initylmg,inwffil,read_rhor,status,timab
      wffclose,wrtout

SOURCE

154 subroutine pead_nl_loop(blkflg,cg,cgindex,dtfil,dtset,d3lo,&
155 & gmet,gprimd,gsqcut, &
156 & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,mband,mgfft,mkmem,mkmem_max,mk1mem,&
157 & mpert,mpi_enreg,mpw,mvwtk,natom,nfft,nkpt,nkpt3,nkxc,nk3xc,nneigh,nspinor,nsppol,&
158 & npwarr,occ,psps,pwind,&
159 & rfpert,rprimd,ucvol,xred)
160 
161 
162 !This section has been created automatically by the script Abilint (TD).
163 !Do not modify the following lines by hand.
164 #undef ABI_FUNC
165 #define ABI_FUNC 'pead_nl_loop'
166 !End of the abilint section
167 
168  implicit none
169 
170 !Arguments ------------------------------------
171 !scalars
172  integer,intent(in) :: mband,mgfft,mk1mem,mkmem,mkmem_max,mpert,mpw,natom,nfft
173  integer,intent(in) :: nk3xc,nkpt,nkpt3,nkxc,nneigh,nspinor,nsppol
174  real(dp),intent(in) :: gsqcut,ucvol
175  type(MPI_type),intent(inout) :: mpi_enreg
176  type(datafiles_type),intent(in) :: dtfil
177  type(dataset_type),intent(inout) :: dtset
178  type(hdr_type),intent(inout) :: hdr
179  type(pseudopotential_type),intent(in) :: psps
180 !arrays
181  integer,intent(in) :: cgindex(nkpt,nsppol),kg(3,mk1mem*mpw),kneigh(30,nkpt)
182  integer,intent(in) :: kg_neigh(30,nkpt,3)
183  integer,intent(in) :: kptindex(2,nkpt3),npwarr(nkpt),pwind(mpw,nneigh,mkmem)
184  integer,intent(in) :: rfpert(3,mpert,3,mpert,3,mpert)
185  integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) !vz_i
186  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),gmet(3,3)
187  real(dp),intent(in) :: gprimd(3,3),k3xc(nfft,nk3xc),kpt3(3,nkpt3)
188  real(dp),intent(in) :: kxc(nfft,nkxc),mvwtk(30,nkpt),rprimd(3,3)
189  real(dp),intent(in) :: xred(3,natom)
190  real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
191  real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert) !vz_i
192 
193 !Local variables-------------------------------
194 !scalars
195  integer,parameter :: level=51
196  integer :: ask_accurate,counter,cplex,formeig,i1dir
197  integer :: i1pert,i2dir,i2pert,i3dir,i3pert,iatom,ierr,iexit,ifft,index,ir
198  integer :: ireadwf,itypat,mcg,mpsang,n1,n2,n3,n3xccc,nfftot,nspden,option,optorth
199  integer :: pert1case,pert2case,pert3case,rdwrpaw,timrev,comm_cell
200  real(dp) :: ecut_eff,exc3,valuei
201  character(len=500) :: message
202  character(len=fnlen) :: fiden1i,fiwf1i,fiwf3i
203  type(wffile_type) :: wff1,wff2,wfft1,wfft2
204  type(wvl_data) :: wvl
205  type(hdr_type) :: hdr_den
206 !arrays
207  integer,allocatable :: atindx(:),atindx1(:),nattyp(:)
208  real(dp) :: d3_berry(2,3),rho_dum(1),tsec(2),ylmgr_dum(1)
209  real(dp),allocatable :: cg1(:,:),cg3(:,:),eigen1(:),ph1d(:,:),rho1r1(:,:)
210  real(dp),allocatable :: rho2g1(:,:),rho2r1(:,:),rho3r1(:,:),vhartr1(:)
211  real(dp),allocatable :: vpsp1(:),vtrial1(:,:),vxc1(:,:),work(:),xc_tmp(:,:)
212  real(dp),allocatable :: xccc3d1(:),xccc3d2(:),xccc3d3(:),ylm(:,:,:)
213  type(pawrhoij_type),allocatable :: rhoij_dum(:)
214 
215 ! ***********************************************************************
216 
217  call timab(502,1,tsec)
218  call status(0,dtfil%filstat,iexit,level,'enter         ')
219 
220  comm_cell = mpi_enreg%comm_cell
221 
222  timrev = 1
223  cplex = 2 - timrev
224  nspden = dtset%nspden
225  ecut_eff = (dtset%ecut)*(dtset%dilatmx)**2
226  mpsang = psps%mpsang
227  optorth=1;if (psps%usepaw==1) optorth=0
228 
229  ABI_ALLOCATE(cg1,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
230  ABI_ALLOCATE(cg3,(2,dtset%mpw*dtset%nspinor*mband*dtset%mk1mem*dtset%nsppol))
231  ABI_ALLOCATE(eigen1,(2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol))
232  ABI_ALLOCATE(rho1r1,(cplex*nfft,dtset%nspden))
233  ABI_ALLOCATE(rho2r1,(cplex*nfft,dtset%nspden))
234  ABI_ALLOCATE(rho2g1,(2,nfft))
235  ABI_ALLOCATE(rho3r1,(cplex*nfft,dtset%nspden))
236  ABI_ALLOCATE(ylm,(2,dtset%mpw*dtset%mkmem,mpsang*mpsang*psps%useylm))
237 
238  ask_accurate=1 ; formeig = 1 ; ireadwf = 1
239  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
240  nfftot=n1*n2*n3
241 
242 !Generate an index table of atoms, in order for them to be used
243 !type after type.
244  ABI_ALLOCATE(atindx,(natom))
245  ABI_ALLOCATE(atindx1,(natom))
246  ABI_ALLOCATE(nattyp,(psps%ntypat))
247  index=1
248  do itypat=1,psps%ntypat
249    nattyp(itypat)=0
250    do iatom=1,natom
251      if(dtset%typat(iatom)==itypat)then
252        atindx(iatom)=index
253        atindx1(index)=iatom
254        index=index+1
255        nattyp(itypat)=nattyp(itypat)+1
256      end if
257    end do
258  end do
259 
260 !Generate the 1-dimensional phases
261  ABI_ALLOCATE(ph1d,(2,3*(2*mgfft+1)*natom))
262  call getph(atindx,natom,n1,n2,n3,ph1d,xred)
263 
264 !Set up the Ylm for each k point
265  if (psps%useylm==1) then
266    call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,&
267 &   dtset%mpw,dtset%nband,dtset%nkpt,&
268 &   npwarr,dtset%nsppol,0,rprimd,ylm,ylmgr_dum)
269  end if
270 
271  ABI_ALLOCATE(vpsp1,(cplex*nfft))
272  ABI_ALLOCATE(xccc3d1,(cplex*nfft))
273  ABI_ALLOCATE(xccc3d2,(cplex*nfft))
274  ABI_ALLOCATE(xccc3d3,(cplex*nfft))
275  ABI_ALLOCATE(vhartr1,(cplex*nfft))
276  ABI_ALLOCATE(vxc1,(cplex*nfft,dtset%nspden))
277  ABI_ALLOCATE(vtrial1,(cplex*nfft,dtset%nspden))
278 
279 !Loop over the perturbations j1, j2, j3
280 
281  pert1case = 0 ; pert2case = 0 ; pert3case = 0
282 
283  do i1pert = 1, mpert
284    do i1dir = 1, 3
285 
286      if ((maxval(rfpert(i1dir,i1pert,:,:,:,:))==1)) then
287 
288        pert1case = i1dir + (i1pert-1)*3
289        counter = pert1case
290        call appdig(pert1case,dtfil%fnamewff1,fiwf1i)
291 
292        call status(counter,dtfil%filstat,iexit,level,'call inwffil  ')
293 
294        mcg=mpw*nspinor*mband*mkmem*nsppol
295        call inwffil(ask_accurate,cg1,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
296 &       formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
297 &       dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
298 &       dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
299 &       dtset%nsppol,dtset%nsym,&
300 &       occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
301 &       dtfil%unkg1,wff1,wfft1,dtfil%unwff1,fiwf1i,wvl)
302 
303        if (ireadwf==1) then
304          call WffClose (wff1,ierr)
305        end if
306 
307        rho1r1(:,:) = 0._dp
308        if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
309          rdwrpaw=0
310          call appdig(pert1case,dtfil%fildens1in,fiden1i)
311          call status(counter,dtfil%filstat,iexit,level,'call ioarr    ')
312 
313          call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho1r1, &
314          hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
315          call hdr_free(hdr_den)
316        end if
317 
318        xccc3d1(:) = 0._dp
319        if ((psps%n1xccc/=0).and.(i1pert <= natom)) then
320          call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkcore   ')
321          call dfpt_mkcore(cplex,i1dir,i1pert,natom,psps%ntypat,n1,psps%n1xccc,&
322 &         n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
323 &         psps%xcccrc,psps%xccc1d,xccc3d1,xred)
324        end if ! psps%n1xccc/=0
325 
326        do i3pert = 1, mpert
327          do i3dir = 1, 3
328 
329            if ((maxval(rfpert(i1dir,i1pert,:,:,i3dir,i3pert))==1)) then
330 
331              pert3case = i3dir + (i3pert-1)*3
332              counter = 100*pert3case + pert1case
333              call appdig(pert3case,dtfil%fnamewff1,fiwf3i)
334 
335              call status(counter,dtfil%filstat,iexit,level,'call inwffil  ')
336              mcg=mpw*nspinor*mband*mkmem*nsppol
337              call inwffil(ask_accurate,cg3,dtset,dtset%ecut,ecut_eff,eigen1,dtset%exchn2n3d,&
338 &             formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,dtset%localrdwf,&
339 &             dtset%mband,mcg,dtset%mk1mem,mpi_enreg,mpw,&
340 &             dtset%nband,dtset%ngfft,dtset%nkpt,npwarr,&
341 &             dtset%nsppol,dtset%nsym,&
342 &             occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,&
343 &             dtfil%unkg1,wff2,wfft2,dtfil%unwff2,&
344 &             fiwf3i,wvl)
345              if (ireadwf==1) then
346                call WffClose (wff2,ierr)
347              end if
348 
349              rho3r1(:,:) = 0._dp
350              if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
351                rdwrpaw=0
352                call appdig(pert3case,dtfil%fildens1in,fiden1i)
353                call status(counter,dtfil%filstat,iexit,level,'call ioarr    ')
354 
355                call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho3r1, &
356                hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
357                call hdr_free(hdr_den)
358              end if
359 
360              xccc3d3(:) = 0._dp
361              if ((psps%n1xccc/=0).and.(i3pert <= natom)) then
362                call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkcore   ')
363                call dfpt_mkcore(cplex,i3dir,i3pert,natom,psps%ntypat,n1,psps%n1xccc,&
364 &               n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
365 &               psps%xcccrc,psps%xccc1d,xccc3d3,xred)
366              end if ! psps%n1xccc/=0
367 
368              do i2pert = 1, mpert
369 
370 !              In case of electric field perturbation, evaluate the ddk
371 !              using the finite difference expression of
372 !              Marzari and Vanderbilt PRB 56, 12847 (1997) [[cite:Marzari1997]].
373 
374                d3_berry(:,:) = 0._dp
375 
376                if ((i2pert==dtset%natom+2).and.&
377 &               (maxval(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert)) == 1)) then
378 
379                  call timab(511,1,tsec)
380                  call status(counter,dtfil%filstat,iexit,level,'call pead_nl_mv  ')
381                  call pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,&
382 &                 i1pert,i3pert,i1dir,i3dir,&
383 &                 kneigh,kg_neigh,kptindex,kpt3,mband,mkmem,mkmem_max,mk1mem,&
384 &                 mpi_enreg,mpw,mvwtk,natom,nkpt,nkpt3,nneigh,npwarr,nspinor,nsppol,pwind)
385                  call timab(511,2,tsec)
386 
387                end if
388 
389                if (mpi_enreg%me == 0) then
390 
391                  if(sum(rfpert(i1dir,i1pert,:,i2pert,i3dir,i3pert))>0)then
392                    write(message,'(a,a,a,a,a,a)')ch10,ch10,&
393 &                   ' Decomposition of the third-order energy for the set of perturbations',ch10
394                    call wrtout(std_out,message,'COLL')
395                    call wrtout(ab_out,message,'COLL')
396                    if (i1pert < natom + 1) then
397                      write(message,'(a,i3,a,i3)') &
398 &                     ' j1 : displacement of atom ',i1pert,' along direction ', i1dir
399                    end if
400                    if (i1pert == dtset%natom + 2) then
401                      write(message,'(a,i4)')' j1 : homogeneous electric field along direction ',i1dir
402                    end if
403                    call wrtout(std_out,message,'COLL')
404                    call wrtout(ab_out,message,'COLL')
405                    if (i3pert < natom + 1) then
406                      write(message,'(a,i3,a,i3,a)') &
407 &                     ' j3 : displacement of atom ',i3pert,' along direction ', i3dir,ch10
408                    end if
409                    if (i3pert == dtset%natom + 2) then
410                      write(message,'(a,i4,a)')' j3 : homogeneous electric field along direction ',i3dir,ch10
411                    end if
412                    call wrtout(std_out,message,'COLL')
413                    call wrtout(ab_out,message,'COLL')
414                  end if
415 
416                end if ! mpi_enreg%me == 0
417 
418                do i2dir = 1, 3
419 
420                  if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
421                    pert2case = i2dir + (i2pert-1)*3
422 
423                    blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
424 
425 !                  Read the first-order densities from disk-files
426                    rho2r1(:,:) = 0._dp ; rho2g1(:,:) = 0._dp
427 
428                    if (dtset%get1den /= 0 .or. dtset%ird1den /= 0) then
429                      rdwrpaw=0
430                      call appdig(pert2case,dtfil%fildens1in,fiden1i)
431                      call status(counter,dtfil%filstat,iexit,level,'call ioarr    ')
432 
433                      call read_rhor(fiden1i, cplex, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, mpi_enreg, rho2r1, &
434                      hdr_den, rhoij_dum, comm_cell, check_hdr=hdr)
435                      call hdr_free(hdr_den)
436 
437 !                    Compute up+down rho1(G) by fft
438                      ABI_ALLOCATE(work,(cplex*nfft))
439                      work(:)=rho2r1(:,1)
440                      call status(counter,dtfil%filstat,iexit,level,'call fourdp   ')
441                      call fourdp(cplex,rho2g1,work,-1,mpi_enreg,nfft,dtset%ngfft,dtset%paral_kgb,0)
442                      ABI_DEALLOCATE(work)
443 
444                    end if
445 
446 !                  Compute first-order local potentials
447 !                  (hartree, xc and pseudopotential)
448 
449                    n3xccc=0; if(psps%n1xccc/=0)n3xccc=nfft
450                    xccc3d2(:)=0._dp ; vpsp1(:)=0._dp
451 
452                    if (i2pert <= natom) then
453 
454                      call status(counter,dtfil%filstat,iexit,level,'call dfpt_vlocal   ')
455                      call dfpt_vlocal(atindx,cplex,gmet,gsqcut,i2dir,i2pert,mpi_enreg,psps%mqgrid_vl,natom,&
456 &                     nattyp,nfft,dtset%ngfft,psps%ntypat,n1,n2,n3,dtset%paral_kgb,ph1d,psps%qgrid_vl,&
457 &                     dtset%qptn,ucvol,psps%vlspl,vpsp1,xred)
458 
459                      if (psps%n1xccc/=0) then
460                        call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkcore   ')
461                        call dfpt_mkcore(cplex,i2dir,i2pert,natom,psps%ntypat,n1,psps%n1xccc,&
462 &                       n2,n3,dtset%qptn,rprimd,dtset%typat,ucvol,&
463 &                       psps%xcccrc,psps%xccc1d,xccc3d2,xred)
464                      end if ! psps%n1xccc/=0
465 
466                    end if  ! i2pert <= natom
467 
468                    call status(counter,dtfil%filstat,iexit,level,'call hartre   ')
469                    call hartre(cplex,gsqcut,0,mpi_enreg,nfft,dtset%ngfft,dtset%paral_kgb,rho2g1,rprimd,vhartr1)
470                    option=1
471                    call status(counter,dtfil%filstat,iexit,level,'call dfpt_mkvxc   ')
472                    call dfpt_mkvxc(cplex,dtset%ixc,kxc,mpi_enreg,nfft,dtset%ngfft,&
473 &                   rho_dum,0,rho_dum,0,nkxc,dtset%nspden,n3xccc,option,&
474 &                   dtset%paral_kgb,dtset%qptn,rho2r1,rprimd,0,vxc1,xccc3d2)
475 
476                    if(dtset%nsppol==1)then
477                      if(cplex==1)then
478                        do ir=1,nfft
479                          vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
480                        end do
481                      else
482                        do ir=1,nfft
483                          vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
484                          vtrial1(2*ir  ,1)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,1)
485                        end do
486                      end if
487                    else
488                      if(cplex==1)then
489                        do ir=1,nfft
490                          vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
491                          vtrial1(ir,2)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,2)
492                        end do
493                      else
494 !                      fab: I think there was an error in the definition of  vtrial1(2*ir-1,2); I have corrected it...
495                        do ir=1,nfft
496                          vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
497                          vtrial1(2*ir  ,1)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,1)
498                          vtrial1(2*ir-1,2)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1  ,2)
499                          vtrial1(2*ir  ,2)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,2)
500                        end do
501                      end if
502                    end if
503 
504 !                  Compute the third-order xc energy
505 !                  take into account the contribution of the term
506 !$
507 !                  \frac{d}{d \lambda}
508 !                  \frac{\delta^2 E_{Hxc}}{\delta n(r) \delta n(r\prim)}
509 !$
510 !                  (seventh term of Eq. (110) of X. Gonze, PRA 52, 1096 (1995) [[cite:Gonze1995]]).
511 
512 !                  the following are essentially the 4th and the 3rd terms of PRB 71,125107 [[cite:Veithen2005]], but the
513 !                  multiplication for rho1 will be done by dotprod_vn later
514 
515 !                  in the non spin polarized case xc_tmp has only 1 component
516                    if (nspden==1)then
517 
518                      ABI_ALLOCATE(xc_tmp,(cplex*nfft,1))
519 
520                      if (cplex==1) then
521 !                      This, and the next lines, have to be changed in case cplex=2
522                        do ifft=1,nfft
523                          xc_tmp(ifft,1)= k3xc(ifft,1)*(rho2r1(ifft,1)+3*xccc3d2(ifft))*rho3r1(ifft,1)
524                        end do
525                      else
526                        do ifft=1,nfft   ! 2*ifft-1 denotes the real part, 2*ifft the imaginary part
527                          xc_tmp(2*ifft-1,1)= k3xc(ifft,1)*( (rho2r1(2*ifft-1,1)+3*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,1) &
528 &                         -( rho2r1(2*ifft,1)+3*xccc3d2(2*ifft))*rho3r1(2*ifft,1))
529 
530                          xc_tmp(2*ifft,1)= k3xc(ifft,1)*( (rho2r1(2*ifft-1,1)+3*xccc3d2(2*ifft-1))*rho3r1(2*ifft,1) &
531 &                         +( rho2r1(2*ifft,1)+3*xccc3d2(2*ifft))*rho3r1(2*ifft-1,1))
532                        end do
533 
534                      end if
535 
536                    end if
537 
538 !                  fab: modifications for the spin polarized raman part:
539 !                  in the spin polarized case xc_tmp has 2 components
540 !                  note that now the non linear core correction is divided by 2
541                    if (nspden==2) then
542 
543                      ABI_ALLOCATE(xc_tmp,(cplex*nfft,2))
544 
545                      if (cplex==1) then
546                        do ifft=1,nfft
547                          xc_tmp(ifft,1)= k3xc(ifft,1)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ &
548 &                         k3xc(ifft,2)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2))+ &
549 &                         k3xc(ifft,2)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ &
550 &                         k3xc(ifft,3)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2))
551                          xc_tmp(ifft,2)= k3xc(ifft,2)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ &
552 &                         k3xc(ifft,3)*(rho2r1(ifft,2)+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2))+ &
553 &                         k3xc(ifft,3)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*rho3r1(ifft,2)+ &
554 &                         k3xc(ifft,4)*((rho2r1(ifft,1)-rho2r1(ifft,2))+(3._dp/2._dp)*xccc3d2(ifft))*(rho3r1(ifft,1)-rho3r1(ifft,2))
555                        end do
556 
557                      else
558                        do ifft=1,nfft
559 !                        These sections should be rewritten, to be easier to read ... (defining intermediate scalars)
560                          xc_tmp(2*ifft-1,1)= k3xc(ifft,1)*&
561 &                         ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,2)- &
562 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+   &
563 &                         k3xc(ifft,2)*&
564 &                         ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*(rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- &
565 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft,1)-rho3r1(2*ifft,2)))+ &
566 &                         k3xc(ifft,2)*&
567 &                         ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,2)- &
568 &                         ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+ &
569 &                         k3xc(ifft,3)*&
570 &                         ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
571 &                         (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- &
572 &                         ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*&
573 &                         (rho3r1(2*ifft,1)-rho3r1(2*ifft,2)))
574                          xc_tmp(2*ifft,1)=k3xc(ifft,1)*&
575 &                         ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft,2)+ &
576 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+   &
577 &                         k3xc(ifft,2)*&
578 &                         ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*(rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ &
579 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2)))+ &
580 &                         k3xc(ifft,2)*&
581 &                         ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft,2)+ &
582 &                         ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+ &
583 &                         k3xc(ifft,3)*&
584 &                         ( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
585 &                         (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ &
586 &                         ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*&
587 &                         (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2)))
588 !                        fab: now the spin down component
589                          xc_tmp(2*ifft-1,2)= k3xc(ifft,2)*&
590 &                         ( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*rho3r1(2*ifft-1,2)- &
591 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+   &
592 &                         k3xc(ifft,3)*( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
593 &                         (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- &
594 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft,1)-rho3r1(2*ifft,2)))+ &
595 &                         k3xc(ifft,3)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
596 &                         rho3r1(2*ifft-1,2)- &
597 &                         ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft,2))+ &
598 &                         k3xc(ifft,4)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
599 &                         (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2))- &
600                          ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*&
601 &                         (rho3r1(2*ifft,1)-rho3r1(2*ifft,2)))
602                          xc_tmp(2*ifft,2)=k3xc(ifft,1)*( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
603 &                         rho3r1(2*ifft,2)+ &
604 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+   &
605 &                         k3xc(ifft,3)*( (rho2r1(2*ifft-1,2)+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
606 &                         (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ &
607 &                         (rho2r1(2*ifft,2)+(3._dp/2._dp)*xccc3d2(2*ifft))*(rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2)))+ &
608 &                         k3xc(ifft,3)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
609 &                         rho3r1(2*ifft,2)+ &
610 &                         ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*rho3r1(2*ifft-1,2))+ &
611 &                         k3xc(ifft,4)*( ((rho2r1(2*ifft-1,1)-rho2r1(2*ifft-1,2))+(3._dp/2._dp)*xccc3d2(2*ifft-1))*&
612 &                         (rho3r1(2*ifft,1)-rho3r1(2*ifft,2))+ &
613 &                         ((rho2r1(2*ifft,1)-rho2r1(2*ifft,2))+(3._dp/2._dp)*xccc3d2(2*ifft))*&
614 &                         (rho3r1(2*ifft-1,1)-rho3r1(2*ifft-1,2)))
615                        end do
616 
617 !                      fab: this is the end if over cplex
618                      end if
619 !                    fab: this is the enf if over nspden
620                    end if
621 
622                    call dotprod_vn(1,rho1r1,exc3,valuei,nfft,nfftot,nspden,1,xc_tmp,ucvol,mpi_comm_sphgrid=mpi_enreg%comm_fft)
623                    ABI_DEALLOCATE(xc_tmp)
624 
625 !                  Perform DFPT part of the 3dte calculation
626 
627                    call timab(512,1,tsec)
628                    call status(counter,dtfil%filstat,iexit,level,'call pead_nl_resp ')
629                    call pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,&
630 &                   kg,mband,mgfft,mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,&
631 &                   nspinor,nsppol,npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm)
632                    call timab(512,2,tsec)
633 
634                    call status(counter,dtfil%filstat,iexit,level,'after pead_nl_resp')
635 
636 !                  Describe the perturbation and write out the result
637                    if (mpi_enreg%me == 0) then
638                      if (i2pert < natom + 1) then
639                        write(message,'(a,i3,a,i3)') &
640 &                       ' j2 : displacement of atom ',i2pert,&
641 &                       ' along direction ', i2dir
642                      end if
643                      if (i2pert == dtset%natom + 2) then
644                        write(message,'(a,i4)') &
645 &                       ' j2 : homogeneous electric field along direction ',&
646 &                       i2dir
647                      end if
648                      call wrtout(std_out,message,'COLL')
649                      call wrtout(ab_out,message,'COLL')
650                      write(ab_out,'(20x,a,13x,a)')'real part','imaginary part'
651                      write(ab_out,'(5x,a2,1x,f22.10,3x,f22.10)')'xc',exc3*sixth,zero
652                      if (i2pert == natom + 2) then
653                        write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',&
654 &                       d3_berry(1,i2dir),d3_berry(2,i2dir)
655                      end if
656                      write(ab_out,'(5x,a3,f22.10,3x,f22.10)')'dft',&
657 &                     d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
658 &                     d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
659                      write(ab_out,*)
660                      write(std_out,'(18x,a,11x,a)')'real part','imaginary part'
661                      write(std_out,'(5x,a2,1x,f20.10,3x,f20.10)')'xc',exc3*sixth,zero
662                      write(std_out,'(5x,a3,f22.10,3x,f22.10)')'ddk',&
663 &                     d3_berry(1,i2dir),d3_berry(2,i2dir)
664                      write(std_out,'(5x,a3,f22.10,3x,f22.10)')'dft',&
665 &                     d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert),&
666 &                     d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
667                      write(std_out,*)
668                    end if  ! mpi_enreg%me == 0
669 
670                    d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = &
671 &                   d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + exc3*sixth
672                    d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = &
673 &                   d3lo(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + d3_berry(:,i2dir)
674 
675                  end if   !rfpert
676                end do    !i2dir
677              end do     ! i2pert
678 
679            end if   ! rfpert
680          end do    ! i3dir
681        end do     ! i3pert
682 
683      end if   ! rfpert
684    end do    ! i1dir
685  end do     ! i1pert
686 
687  call status(0,dtfil%filstat,iexit,level,'exit          ')
688 
689  ABI_DEALLOCATE(cg1)
690  ABI_DEALLOCATE(cg3)
691  ABI_DEALLOCATE(eigen1)
692  ABI_DEALLOCATE(rho1r1)
693  ABI_DEALLOCATE(rho2r1)
694  ABI_DEALLOCATE(rho2g1)
695  ABI_DEALLOCATE(rho3r1)
696  ABI_DEALLOCATE(atindx1)
697  ABI_DEALLOCATE(atindx)
698  ABI_DEALLOCATE(nattyp)
699  ABI_DEALLOCATE(ph1d)
700  ABI_DEALLOCATE(ylm)
701  ABI_DEALLOCATE(vtrial1)
702  ABI_DEALLOCATE(vxc1)
703  ABI_DEALLOCATE(vhartr1)
704  ABI_DEALLOCATE(vpsp1)
705  ABI_DEALLOCATE(xccc3d1)
706  ABI_DEALLOCATE(xccc3d2)
707  ABI_DEALLOCATE(xccc3d3)
708 
709  call timab(502,2,tsec)
710 
711 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

PARENTS

      pead_nl_loop

CHILDREN

      dzgedi,dzgefa,mpi_recv,mpi_send,status,wrtout,xmpi_sum

SOURCE

1104 subroutine pead_nl_mv(cg,cgindex,cg1,cg3,dtset,dtfil,d3_berry,gmet,&
1105 &                   i1pert,i3pert,i1dir,i3dir,kneigh,kg_neigh,kptindex,&
1106 &                   kpt3,mband,mkmem,mkmem_max,mk1mem,mpi_enreg,&
1107 &                   mpw,mvwtk,natom,nkpt2,nkpt3,nneigh,npwarr,nspinor,&
1108 &                   nsppol,pwind)
1109 
1110  use m_hide_lapack, only : dzgedi, dzgefa
1111 
1112 !This section has been created automatically by the script Abilint (TD).
1113 !Do not modify the following lines by hand.
1114 #undef ABI_FUNC
1115 #define ABI_FUNC 'pead_nl_mv'
1116 !End of the abilint section
1117 
1118  implicit none
1119 
1120 !Arguments ------------------------------------
1121 !
1122 !---  Arguments : integer scalars
1123  integer, intent(in) :: i1dir,i1pert,i3dir,i3pert,mband,mk1mem
1124  integer, intent(in) :: mkmem,mkmem_max,mpw,natom
1125  integer, intent(in) :: nkpt2,nkpt3,nneigh,nspinor,nsppol
1126 !
1127 !---  Arguments : integer arrays
1128  integer, intent(in) :: cgindex(nkpt2,nsppol)
1129  integer, intent(in) :: kneigh(30,nkpt2),kg_neigh(30,nkpt2,3),kptindex(2,nkpt3)
1130  integer, intent(in) :: npwarr(nkpt2),pwind(mpw,nneigh,mkmem)
1131 !
1132 !---  Arguments : real(dp) scalars
1133 !
1134 !---  Arguments : real(dp) arrays
1135  real(dp), intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
1136  real(dp), intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
1137  real(dp), intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
1138  real(dp), intent(in) :: gmet(3,3),kpt3(3,nkpt3)
1139  real(dp), intent(in) :: mvwtk(30,nkpt2)
1140  real(dp), intent(out) :: d3_berry(2,3)
1141 !
1142 !---  Arguments : structured datatypes
1143  type(MPI_type), intent(in) :: mpi_enreg
1144  type(datafiles_type), intent(in) :: dtfil
1145  type(dataset_type), intent(in) :: dtset
1146 
1147 !Local variables-------------------------------
1148 !
1149 !---- Local variables : integer scalars
1150  integer :: count,counter,count1,iband,icg
1151  integer :: ierr,iexit,ii,ikpt,ikpt_loc,ikpt2
1152  integer :: ikpt_rbz,ineigh,info,ipw,isppol,jband,jcg,jj,jkpt,job,jpw, jkpt2, jkpt_rbz
1153  integer :: lband,lpband,nband_occ,npw_k,npw_k1,my_source,his_source,dest,tag
1154  integer :: spaceComm
1155  integer,parameter :: level=52
1156  integer :: bdtot_index
1157 !
1158 !---- Local variables : integer arrays
1159  integer,allocatable :: ipvt(:)
1160  integer, allocatable :: bd_index(:,:)
1161 !
1162 !---- Local variables : real(dp) scalars
1163  real(dp) :: dotnegi,dotnegr,dotposi,dotposr
1164 ! real(dp) :: c1,c2 ! appear commented out below
1165 !
1166 !---- Local variables : real(dp) arrays
1167  real(dp) :: d3_aux(2,3),det(2,2),dk(3),dk_(3)
1168  real(dp) :: z1(2),z2(2)
1169  real(dp),allocatable :: buffer(:,:),cgq(:,:),cg1q(:,:),cg3q(:,:)
1170  real(dp),allocatable :: qmat(:,:,:),s13mat(:,:,:),s1mat(:,:,:),s3mat(:,:,:)
1171  real(dp),allocatable :: smat(:,:,:),zgwork(:,:)
1172 !
1173 !---- Local variables : character variables
1174  character(len=500) :: message
1175 !
1176 !---- Local variables : structured datatypes
1177 
1178 #if defined HAVE_MPI
1179 integer :: status1(MPI_STATUS_SIZE)
1180 spaceComm=mpi_enreg%comm_cell
1181 #endif
1182 
1183 
1184 ! ***********************************************************************
1185 
1186  call status(0,dtfil%filstat,iexit,level,'enter         ')
1187 
1188  write(message,'(8a)') ch10,&
1189 & ' pead_nl_mv : finite difference expression of the k-point derivative',ch10,&
1190 & '           is performed using the PEAD formulation of ',&
1191 & 'the third-order energy',ch10,&
1192 & '           (see Nunes and Gonze PRB 63, 155107 (2001) [[cite:Nunes2001]] Eq. 102)',ch10
1193 !call wrtout(ab_out,message,'COLL')
1194  call wrtout(std_out,  message,'COLL')
1195 
1196 
1197 !fab: I think that the following restriction must be eliminated:
1198 !isppol = 1
1199 
1200  ikpt_loc = 0
1201  d3_aux(:,:) = 0_dp
1202 
1203  ABI_ALLOCATE(s13mat,(2,mband,mband))
1204  ABI_ALLOCATE(smat,(2,mband,mband))
1205  ABI_ALLOCATE(s1mat,(2,mband,mband))
1206  ABI_ALLOCATE(qmat,(2,mband,mband))
1207  ABI_ALLOCATE(ipvt,(mband))
1208  ABI_ALLOCATE(s3mat,(2,mband,mband))
1209  ABI_ALLOCATE(zgwork,(2,mband))
1210  ABI_ALLOCATE(bd_index, (nkpt2, nsppol))
1211 
1212  bdtot_index = 0
1213  do isppol = 1, nsppol
1214    do ikpt_rbz = 1, nkpt2
1215      bd_index(ikpt_rbz,isppol) = bdtot_index
1216      bdtot_index = bdtot_index + dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1217    end do
1218  end do
1219 
1220 !fab: I think here I have to add the loop over spin
1221 
1222  do isppol = 1, nsppol
1223 
1224 !  Loop over k-points
1225 !  COMMENT: Every processor has to make mkmem_max iterations
1226 !  even if mkmem < mkemem_max. This is due to the fact
1227 !  that it still has to communicate its wavefunctions
1228 !  to other processors even if it has no more overlap
1229 !  matrices to compute.
1230 
1231    ikpt_loc = 0 ; ikpt = 0
1232 
1233    do while (ikpt_loc < mkmem_max)
1234 
1235      if (ikpt_loc < mkmem) ikpt = ikpt + 1
1236 
1237      if (xmpi_paral == 1) then
1238 !      if ((minval(abs(mpi_enreg%proc_distrb(ikpt,1:mband,1:dtset%nsppol) &
1239 !      &       - mpi_enreg%me)) /= 0).and.(ikpt_loc < mkmem)) cycle
1240        if(ikpt>nkpt2)then
1241          ikpt_loc=mkmem_max
1242          cycle
1243        end if
1244        if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me)) then
1245          if(ikpt==nkpt2) ikpt_loc=mkmem_max
1246          cycle
1247        end if
1248      end if
1249 
1250      ikpt_loc = ikpt_loc + 1
1251      npw_k = npwarr(ikpt)
1252      counter = 100*ikpt
1253      call status(counter,dtfil%filstat,iexit,level,'loop over k   ')
1254 
1255      ii = cgindex(ikpt,isppol)
1256 
1257 !    Loop on the  neighbours
1258 
1259      do ineigh = 1,nneigh
1260 
1261        s13mat(:,:,:) = zero
1262        smat(:,:,:) = zero
1263        s1mat(:,:,:) = zero
1264        s3mat(:,:,:) = zero
1265        qmat(:,:,:) = zero
1266 
1267        ikpt2  = kneigh(ineigh,ikpt)
1268        ikpt_rbz = kptindex(1,ikpt2)   ! index of the k-point in the reduced BZ
1269        jj = cgindex(ikpt_rbz,isppol)
1270        ! previous fixed value for nband_k now called nband_occ:
1271        !nband_occ = dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1272        ! TODO: check if all these bands are occupied in nsppol = 2 case
1273        nband_occ = 0
1274        do iband = 1, dtset%nband(ikpt_rbz+nkpt2*(isppol-1))
1275          !Note, only one image is allowed here (or occ_orig should be the same or all images)
1276          if (dtset%occ_orig(bd_index(ikpt_rbz,isppol) + iband,1) > tol10) nband_occ = nband_occ + 1
1277        end do
1278        npw_k1 = npwarr(ikpt_rbz)
1279        dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt)
1280        dk(:)  = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp)
1281 
1282        count = nspinor*mband*npw_k1
1283        ABI_ALLOCATE(cgq,(2,count))
1284        ABI_ALLOCATE(cg1q,(2,count))
1285        ABI_ALLOCATE(cg3q,(2,count))
1286 
1287 #if defined HAVE_MPI
1288 
1289        my_source = mpi_enreg%proc_distrb(ikpt_rbz,1,1)
1290 
1291 !      do dest = 0, mpi_enreg%nproc-1
1292        do dest = 0, maxval(mpi_enreg%proc_distrb(1:nkpt2,1:mband,1:dtset%nsppol))
1293 
1294          if ((dest==mpi_enreg%me).and.(ikpt_loc <= mkmem)) then
1295 !          I am dest and have something to do
1296 
1297            if (my_source == mpi_enreg%me) then
1298 !            I am destination and source
1299              jcg = cgindex(ikpt_rbz,isppol)
1300 
1301              cgq(:,1:count)  = cg(:,jcg+1:jcg+count)
1302              cg1q(:,1:count) = cg1(:,jcg+1:jcg+count)
1303              cg3q(:,1:count) = cg3(:,jcg+1:jcg+count)
1304 
1305            else
1306 !            I am the destination but not the source -> receive
1307 
1308              tag = ikpt_rbz
1309 
1310              ABI_ALLOCATE(buffer,(2,3*count))
1311 
1312              call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,my_source,tag,spaceComm,status1,ierr)
1313 
1314              cgq(:,1:count)  = buffer(:,1:count)
1315              cg1q(:,1:count) = buffer(:,count+1:2*count)
1316              cg3q(:,1:count) = buffer(:,2*count+1:3*count)
1317              ABI_DEALLOCATE(buffer)
1318 
1319            end if
1320 
1321          else if (ikpt_loc <= mpi_enreg%mkmem(dest)) then  ! dest != me and the dest has a k-point to treat
1322 
1323            jkpt=mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,1)
1324            jkpt2  = kneigh(ineigh,jkpt)
1325            jkpt_rbz = kptindex(1,jkpt2)   ! index of the k-point in the reduced BZ
1326 
1327            his_source = mpi_enreg%proc_distrb(jkpt_rbz,1,1)
1328 
1329            if (his_source == mpi_enreg%me) then
1330 
1331              jcg = cgindex(jkpt_rbz,isppol)
1332 
1333              tag = jkpt_rbz
1334              count1 = npwarr(jkpt_rbz)*mband*nspinor
1335              ABI_ALLOCATE(buffer,(2,3*count1))
1336              buffer(:,1:count1)            = cg(:,jcg+1:jcg+count1)
1337              buffer(:,count1+1:2*count1)   = cg1(:,jcg+1:jcg+count1)
1338              buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1)
1339 
1340              call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,dest,tag,spaceComm,ierr)
1341 
1342              ABI_DEALLOCATE(buffer)
1343 
1344            end if
1345 
1346          end if
1347 
1348        end do
1349 !
1350 !      do jkpt = 1, nkpt2
1351 !
1352 !      if ((jkpt == ikpt_rbz).and.(source /= mpi_enreg%me).and.&
1353 !      &         (ikpt_loc <= mkmem)) then
1354 !
1355 !      tag = jkpt
1356 !
1357 !      allocate(buffer(2,3*count))
1358 !      call MPI_RECV(buffer,2*3*count,MPI_DOUBLE_PRECISION,&
1359 !      source,tag,spaceComm,status1,ierr)
1360 !
1361 !      cgq(:,1:count)  = buffer(:,1:count)
1362 !      cg1q(:,1:count) = buffer(:,count+1:2*count)
1363 !      cg3q(:,1:count) = buffer(:,2*count+1:3*count)
1364 !      deallocate(buffer)
1365 !
1366 !      end if
1367 !
1368 !      !        ----------------------------------------------------------------------------
1369 !      !        --------------- Here: send the WF to all the cpus that need it -------------
1370 !      !        ----------------------------------------------------------------------------
1371 !
1372 !      do dest = 1, mpi_enreg%nproc
1373 !
1374 !      if ((minval(abs(mpi_enreg%proc_distrb(jkpt,1:mband,1:dtset%nsppol) &
1375 !      &           - mpi_enreg%me)) == 0).and.&
1376 !      &           (mpi_enreg%kptdstrb(dest,ineigh,ikpt_loc) == jkpt)) then
1377 !
1378 !
1379 !
1380 !      jcg = cgindex(jkpt,isppol)
1381 !
1382 !      if (((dest-1) == mpi_enreg%me)) then
1383 !
1384 !      cgq(:,1:count)  = cg(:,jcg+1:jcg+count)
1385 !      cg1q(:,1:count) = cg1(:,jcg+1:jcg+count)
1386 !      cg3q(:,1:count) = cg3(:,jcg+1:jcg+count)
1387 !
1388 !      else
1389 !
1390 !      tag = jkpt
1391 !      count1 = npwarr(jkpt)*mband*nspinor
1392 !      allocate(buffer(2,3*count1))
1393 !      buffer(:,1:count1)            = cg(:,jcg+1:jcg+count1)
1394 !      buffer(:,count1+1:2*count1)   = cg1(:,jcg+1:jcg+count1)
1395 !      buffer(:,2*count1+1:3*count1) = cg3(:,jcg+1:jcg+count1)
1396 !
1397 !      call MPI_SEND(buffer,2*3*count1,MPI_DOUBLE_PRECISION,(dest-1),tag,spaceComm,ierr)
1398 !
1399 !      deallocate(buffer)
1400 !
1401 !      end if
1402 !
1403 !      end if
1404 !
1405 !      end do          ! loop over dest
1406 !
1407 !      end do          ! loop over jkpt
1408 
1409        if (ikpt_loc > mkmem) then
1410          ABI_DEALLOCATE(cgq)
1411          ABI_DEALLOCATE(cg1q)
1412          ABI_DEALLOCATE(cg3q)
1413          cycle
1414        end if
1415 
1416 #else
1417 !      no // over k-points
1418 
1419        cgq(:,1:count)  = cg(:,jj+1:jj+count)
1420        cg1q(:,1:count) = cg1(:,jj+1:jj+count)
1421        cg3q(:,1:count) = cg3(:,jj+1:jj+count)
1422 
1423 #endif
1424 
1425 !      Compute overlap matrices
1426 
1427        if (kptindex(2,ikpt2) == 0) then  ! no time-reversal symmetry
1428 
1429          do ipw = 1, npw_k
1430 
1431            jpw = pwind(ipw,ineigh,ikpt_loc)
1432            if (jpw /= 0) then
1433 
1434              do iband = 1, nband_occ
1435                do jband = 1, nband_occ
1436 
1437                  icg = ii + (iband-1)*npw_k + ipw
1438                  jcg = (jband-1)*npw_k1 + jpw
1439 
1440                  smat(1,iband,jband) = smat(1,iband,jband) + &
1441 &                 cg(1,icg)*cgq(1,jcg) + cg(2,icg)*cgq(2,jcg)
1442                  smat(2,iband,jband) = smat(2,iband,jband) + &
1443 &                 cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg)
1444 
1445                  s13mat(1,iband,jband) = s13mat(1,iband,jband) + &
1446 &                 cg1(1,icg)*cg3q(1,jcg) + cg1(2,icg)*cg3q(2,jcg)
1447                  s13mat(2,iband,jband) = s13mat(2,iband,jband) + &
1448 &                 cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg)
1449 
1450                  s1mat(1,iband,jband) = s1mat(1,iband,jband) + &
1451 &                 cg1(1,icg)*cgq(1,jcg) + cg1(2,icg)*cgq(2,jcg) + &
1452 &                 cg(1,icg)*cg1q(1,jcg) + cg(2,icg)*cg1q(2,jcg)
1453                  s1mat(2,iband,jband) = s1mat(2,iband,jband) + &
1454 &                 cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) + &
1455 &                 cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg)
1456 
1457                  s3mat(1,iband,jband) = s3mat(1,iband,jband) + &
1458 &                 cg3(1,icg)*cgq(1,jcg) + cg3(2,icg)*cgq(2,jcg) + &
1459 &                 cg(1,icg)*cg3q(1,jcg) + cg(2,icg)*cg3q(2,jcg)
1460                  s3mat(2,iband,jband) = s3mat(2,iband,jband) + &
1461 &                 cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) + &
1462 &                 cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg)
1463 
1464                end do
1465              end do
1466 
1467            end if
1468 
1469          end do   ! ipw
1470 
1471        else                              ! use time-reversal symmetry
1472 
1473          do ipw = 1,npw_k
1474 
1475            jpw = pwind(ipw,ineigh,ikpt_loc)
1476            if (jpw /= 0) then
1477 
1478              do iband = 1, nband_occ
1479                do jband = 1, nband_occ
1480 
1481                  icg = ii + (iband-1)*npw_k + ipw
1482                  jcg = (jband-1)*npw_k1 + jpw
1483 
1484                  smat(1,iband,jband) = smat(1,iband,jband) + &
1485 &                 cg(1,icg)*cgq(1,jcg) - cg(2,icg)*cgq(2,jcg)
1486                  smat(2,iband,jband) = smat(2,iband,jband) - &
1487 &                 cg(1,icg)*cgq(2,jcg) - cg(2,icg)*cgq(1,jcg)
1488 
1489                  s13mat(1,iband,jband) = s13mat(1,iband,jband) + &
1490 &                 cg1(1,icg)*cg3q(1,jcg) - cg1(2,icg)*cg3q(2,jcg)
1491                  s13mat(2,iband,jband) = s13mat(2,iband,jband) - &
1492 &                 cg1(1,icg)*cg3q(2,jcg) - cg1(2,icg)*cg3q(1,jcg)
1493 
1494                  s1mat(1,iband,jband) = s1mat(1,iband,jband) + &
1495 &                 cg1(1,icg)*cgq(1,jcg) - cg1(2,icg)*cgq(2,jcg) + &
1496 &                 cg(1,icg)*cg1q(1,jcg) - cg(2,icg)*cg1q(2,jcg)
1497                  s1mat(2,iband,jband) = s1mat(2,iband,jband) - &
1498 &                 cg1(1,icg)*cgq(2,jcg) - cg1(2,icg)*cgq(1,jcg) - &
1499 &                 cg(1,icg)*cg1q(2,jcg) - cg(2,icg)*cg1q(1,jcg)
1500 
1501                  s3mat(1,iband,jband) = s3mat(1,iband,jband) + &
1502 &                 cg3(1,icg)*cgq(1,jcg) - cg3(2,icg)*cgq(2,jcg) + &
1503 &                 cg(1,icg)*cg3q(1,jcg) - cg(2,icg)*cg3q(2,jcg)
1504                  s3mat(2,iband,jband) = s3mat(2,iband,jband) - &
1505 &                 cg3(1,icg)*cgq(2,jcg) - cg3(2,icg)*cgq(1,jcg) - &
1506 &                 cg(1,icg)*cg3q(2,jcg) - cg(2,icg)*cg3q(1,jcg)
1507 
1508                end do
1509              end do
1510 
1511            end if
1512 
1513          end do   ! ipw
1514 
1515        end if
1516 
1517        ABI_DEALLOCATE(cgq)
1518        ABI_DEALLOCATE(cg1q)
1519        ABI_DEALLOCATE(cg3q)
1520 
1521 !      Compute qmat, the inverse of smat
1522 
1523        job = 1  ! compute inverse only
1524        qmat(:,:,:) = smat(:,:,:)
1525 
1526        call dzgefa(qmat,mband,nband_occ,ipvt,info)
1527        call dzgedi(qmat,mband,nband_occ,ipvt,det,zgwork,job)
1528 
1529 !      DEBUG
1530 !      write(100,*)
1531 !      write(100,*)'ikpt = ',ikpt,'ineigh = ',ineigh
1532 !      do iband = 1,nband_occ
1533 !      do jband = 1,nband_occ
1534 !      c1 = 0_dp ; c2 = 0_dp
1535 !      do lband = 1,nband_occ
1536 !      c1 = c1 + smat(1,iband,lband)*qmat(1,lband,jband) - &
1537 !      &           smat(2,iband,lband)*qmat(2,lband,jband)
1538 !      c2 = c2 + smat(1,iband,lband)*qmat(2,lband,jband) + &
1539 !      &           smat(2,iband,lband)*qmat(1,lband,jband)
1540 !      end do
1541 !      write(100,'(2(2x,i2),2(2x,f16.9))')iband,jband,&
1542 !      & c1,c2
1543 !      end do
1544 !      end do
1545 !      ENDDEBUG
1546 
1547 
1548 
1549 !      Accumulate sum over bands
1550 
1551        dotposr = 0_dp ; dotposi = 0_dp
1552        dotnegr = 0_dp ; dotnegi = 0_dp
1553        do iband = 1, nband_occ
1554          do jband = 1, nband_occ
1555 
1556            dotposr = dotposr + &
1557 &           s13mat(1,iband,jband)*qmat(1,jband,iband) - &
1558 &           s13mat(2,iband,jband)*qmat(2,jband,iband)
1559            dotposi = dotposi + &
1560 &           s13mat(1,iband,jband)*qmat(2,jband,iband) + &
1561 &           s13mat(2,iband,jband)*qmat(1,jband,iband)
1562 
1563 
1564            do lband = 1, nband_occ
1565              do lpband= 1, nband_occ
1566 
1567                z1(1) = s1mat(1,iband,jband)*qmat(1,jband,lband) - &
1568 &               s1mat(2,iband,jband)*qmat(2,jband,lband)
1569                z1(2) = s1mat(1,iband,jband)*qmat(2,jband,lband) + &
1570 &               s1mat(2,iband,jband)*qmat(1,jband,lband)
1571 
1572                z2(1) = s3mat(1,lband,lpband)*qmat(1,lpband,iband) - &
1573 &               s3mat(2,lband,lpband)*qmat(2,lpband,iband)
1574                z2(2) = s3mat(1,lband,lpband)*qmat(2,lpband,iband) + &
1575 &               s3mat(2,lband,lpband)*qmat(1,lpband,iband)
1576 
1577                dotnegr = dotnegr + &
1578 &               z1(1)*z2(1) - z1(2)*z2(2)
1579                dotnegi = dotnegi + &
1580 &               z1(1)*z2(2) + z1(2)*z2(1)
1581 
1582              end do   ! lpband
1583            end do   ! lband
1584 
1585          end do   ! jband
1586        end do   ! iband
1587 
1588        d3_aux(1,:) = d3_aux(1,:) + &
1589 &       dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposr-dotnegr)
1590        d3_aux(2,:) = d3_aux(2,:) + &
1591 &       dk(:)*mvwtk(ineigh,ikpt)*dtset%wtk(ikpt)*(2_dp*dotposi-dotnegi)
1592 
1593      end do        ! End loop over neighbours
1594 
1595 
1596    end do      ! End loop over k-points
1597 
1598  end do  ! fab: end loop over spin
1599 
1600 
1601 
1602 
1603  call xmpi_sum(d3_aux,spaceComm,ierr)
1604 
1605 
1606  ABI_DEALLOCATE(s13mat)
1607  ABI_DEALLOCATE(smat)
1608  ABI_DEALLOCATE(s1mat)
1609  ABI_DEALLOCATE(qmat)
1610  ABI_DEALLOCATE(ipvt)
1611  ABI_DEALLOCATE(s3mat)
1612  ABI_DEALLOCATE(zgwork)
1613  ABI_DEALLOCATE(bd_index)
1614 
1615 
1616 !fab: I think that in the following we have to make a distinction:
1617 !for the spin unpolarized case we leave the PEAD expression as it is, while
1618 !in the spin polarized case we have simply to divide by 2
1619 !(see eq.19 di PRB 63,155107 [[cite:Nunes2001]], eq. 7 di PRB 71,125107 [[cite:Veithen2005]]
1620 ! and eq 13 di PRB 71, 125107 [[cite:Veithen2005]] ...
1621 !in this latter equation the 2 must be simply replaced by the sum over the spin components...
1622 !and indeed we have inserted the loop over the spin,
1623 !but there was a factor 2 already present in the routine due to spin degenracy that had to be removed)
1624 
1625 
1626  if (nsppol==1) then
1627 
1628 !  Take minus the imaginary part
1629 
1630    d3_berry(1,:) = -1_dp*d3_aux(2,:)
1631    d3_berry(2,:) = d3_aux(1,:)
1632 
1633    d3_berry(2,:) = 0_dp
1634 
1635  else
1636 
1637    d3_berry(1,:) = -1_dp*d3_aux(2,:)/2._dp
1638    d3_berry(2,:) = d3_aux(1,:)/2._dp
1639 
1640    d3_berry(2,:) = 0_dp/2._dp
1641 
1642  end if
1643 
1644 !DEBUG
1645 !write(100,*)'pead_nl_mv.f : d3_berry'
1646 !write(100,*)'Perturbation',i1dir,i3dir
1647 !write(100,*)
1648 !write(100,*)'before transformation'
1649 !write(100,*)'real part'
1650 !write(100,'(3(2x,f20.9))')d3_berry(1,:)
1651 !write(100,*)
1652 !write(100,*)'imaginary part'
1653 !write(100,'(3(2x,f20.9))')d3_berry(2,:)
1654 !write(100,*)
1655 !write(100,*)'after transformation'
1656 !ENDDEBUG
1657 
1658 !Compute the projection on the basis vectors of
1659 !reciprocal space
1660 
1661  d3_aux(:,:) = 0_dp
1662  do ii = 1,3
1663    do jj = 1,3
1664      d3_aux(:,ii) = d3_aux(:,ii) + gmet(ii,jj)*d3_berry(:,jj)
1665    end do
1666  end do
1667  d3_berry(:,:) = d3_aux(:,:)
1668 
1669 !Write out the berryphase part of the third order energy
1670 
1671  if (mpi_enreg%me == 0) then
1672 
1673    write(message,'(a,a,a)')ch10,' Berryphase part of the third-order energy:',ch10
1674    call wrtout(std_out,  message,'COLL')
1675 
1676    if (i1pert < natom + 1) then
1677      write(message,'(a,i3,a,i3)')&
1678 &     '            j1: Displacement of atom ',i1pert,&
1679 &     ' along direction ',i1dir
1680    else if (i1pert == natom + 2) then
1681      write(message,'(a,i3)')&
1682 &     '            j1: homogenous electric field along direction ',i1dir
1683    end if
1684    call wrtout(std_out,  message,'COLL')
1685 
1686    write(message,'(a)')&
1687 &   '            j2: k-point derivative along direction i2dir '
1688    call wrtout(std_out,  message,'COLL')
1689 
1690    if (i3pert < natom + 1) then
1691      write(message,'(a,i3,a,i3,a)')&
1692 &     '            j3: Displacement of atom ',i3pert,&
1693 &     ' along direction ',i3dir,ch10
1694    else if (i3pert == natom + 2) then
1695      write(message,'(a,i3,a)')&
1696 &     '            j3: homogenous electric field along direction ',i3dir,ch10
1697    end if
1698    call wrtout(std_out,  message,'COLL')
1699 
1700 !  write(ab_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part'
1701    write(std_out,'(5x,a5,8x,a9,5x,a14)')'i2dir','real part','imaginary part'
1702    do ii = 1,3
1703      write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,&
1704 &     d3_berry(1,ii),d3_berry(2,ii)
1705      write(std_out,'(7x,i1,3x,f16.9,3x,f16.9)')ii,&
1706 &     d3_berry(1,ii),d3_berry(2,ii)
1707    end do
1708 
1709  end if    ! mpi_enreg%me == 0
1710 
1711 !DEBUG
1712 !write(100,*)'real part'
1713 !write(100,'(3(2x,f20.9))')d3_berry(1,:)
1714 !write(100,*)
1715 !write(100,*)'imaginary part'
1716 !write(100,'(3(2x,f20.9))')d3_berry(2,:)
1717 !ENDDEBUG
1718 
1719  call status(0,dtfil%filstat,iexit,level,'exit          ')
1720 
1721 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

PARENTS

      pead_nl_loop

CHILDREN

      destroy_hamiltonian,dotprod_g,fftpac,fourwf,init_hamiltonian
      load_k_hamiltonian,mkffnl,mkkpg,nonlop,status,xmpi_sum

SOURCE

 770 subroutine pead_nl_resp(cg,cg1,cg3,cplex,dtfil,dtset,d3lo,&
 771 & i1dir,i2dir,i3dir,i1pert,i2pert,i3pert,&
 772 & kg,mband,mgfft,mkmem,mk1mem,&
 773 & mpert,mpi_enreg,mpsang,mpw,natom,nfft,nkpt,nspden,nspinor,nsppol,&
 774 & npwarr,occ,ph1d,psps,rprimd,vtrial1,xred,ylm)
 775 
 776 
 777 !This section has been created automatically by the script Abilint (TD).
 778 !Do not modify the following lines by hand.
 779 #undef ABI_FUNC
 780 #define ABI_FUNC 'pead_nl_resp'
 781 !End of the abilint section
 782 
 783  implicit none
 784 
 785 !Arguments ------------------------------------
 786 !scalars
 787  integer,intent(in) :: cplex,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,mband,mgfft
 788  integer,intent(in) :: mk1mem,mkmem,mpert,mpsang,mpw,natom,nfft,nkpt,nspden
 789  integer,intent(in) :: nspinor,nsppol
 790  type(MPI_type),intent(in) :: mpi_enreg
 791  type(datafiles_type),intent(in) :: dtfil
 792  type(dataset_type),intent(in) :: dtset
 793  type(pseudopotential_type),intent(in) :: psps
 794 !arrays
 795  integer,intent(in) :: kg(3,mpw*mkmem),npwarr(nkpt)
 796  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
 797  real(dp),intent(in) :: cg1(2,mpw*nspinor*mband*mk1mem*nsppol)
 798  real(dp),intent(in) :: cg3(2,mpw*nspinor*mband*mk1mem*nsppol)
 799  real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom),rprimd(3,3)
 800  real(dp),intent(in) :: xred(3,natom),ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
 801  real(dp),intent(inout) :: vtrial1(cplex*nfft,nspden)
 802  real(dp),intent(inout) :: d3lo(2,3,mpert,3,mpert,3,mpert)
 803 
 804 !Local variables-------------------------------
 805 !scalars
 806  integer,parameter :: level=52
 807  integer :: bantot,choice,counter,cpopt,dimffnl,iband,icg0,ider,ierr,iexit
 808  integer :: ii,ikg,ikpt,ilm,ipw,isppol,istwf_k,jband,jj
 809  integer :: me,n1,n2,n3,n4,n5,n6,nband_k,nkpg,nnlout,npw_k
 810  integer :: option,paw_opt,signs,spaceComm,tim_fourwf,tim_nonlop
 811  real(dp) :: dot1i,dot1r,dot2i,dot2r,doti,dotr,lagi,lagr,sumi,sumr,weight
 812  type(gs_hamiltonian_type) :: gs_hamk
 813 !arrays
 814  integer,allocatable :: kg_k(:,:)
 815  real(dp) :: buffer(2),enlout(3),kpq(3),kpt(3)
 816  real(dp) :: dum_svectout(1,1),dum(1),rmet(3,3),ylmgr_dum(1,1,1)
 817  real(dp),allocatable :: cwave0(:,:),cwavef3(:,:),ffnlk(:,:,:,:)
 818  real(dp),allocatable :: gh0(:,:),gh1(:,:),gvnl(:,:),kpg_k(:,:)
 819  real(dp),allocatable :: vlocal1(:,:,:),wfraug(:,:,:,:),ylm_k(:,:)
 820  type(pawcprj_type) :: cprj_dum(1,1)
 821  type(pawtab_type) :: pawtab_dum(0)
 822 
 823 !***********************************************************************
 824 
 825  me = mpi_enreg%me
 826  spaceComm=mpi_enreg%comm_cell
 827 
 828  call status(0,dtfil%filstat,iexit,level,'enter         ')
 829 
 830  bantot = 0
 831  icg0 = 0
 832  option = 2
 833  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
 834  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
 835 
 836  ABI_ALLOCATE(vlocal1,(cplex*n4,n5,n6))
 837  ABI_ALLOCATE(wfraug,(2,n4,n5,n6))
 838 
 839 !Initialize Hamiltonian (k-independent terms) - NCPP only
 840  call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,&
 841 & dtset%typat,xred,nfft,mgfft,dtset%ngfft,rprimd,dtset%nloalg,ph1d=ph1d,&
 842 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
 843 !& paw_ij=paw_ij)
 844  rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
 845 
 846  sumr = zero ; sumi = zero
 847 
 848 !Loop over spins
 849 
 850  do isppol = 1, nsppol
 851 
 852    call status(0,dtfil%filstat,iexit,level,'call fftpac   ')
 853    call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,vtrial1,vlocal1,option)
 854 
 855 !  Loop over k-points
 856 
 857    ikg = 0
 858    do ikpt = 1, nkpt
 859 
 860      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me))cycle
 861 
 862      counter = 100*ikpt
 863 
 864      nband_k = dtset%nband(ikpt+(isppol-1)*nkpt)
 865      npw_k = npwarr(ikpt)
 866      istwf_k = dtset%istwfk(ikpt)
 867 
 868      kpt(:) = dtset%kptns(:,ikpt)
 869      kpq(:) = dtset%kptns(:,ikpt) ! In case of non zero q, kpt = kpt + q
 870 
 871      ABI_ALLOCATE(cwave0,(2,npw_k*dtset%nspinor))
 872      ABI_ALLOCATE(cwavef3,(2,npw_k*dtset%nspinor))
 873      ABI_ALLOCATE(gh0,(2,npw_k*dtset%nspinor))
 874      ABI_ALLOCATE(gvnl,(2,npw_k*dtset%nspinor))
 875      ABI_ALLOCATE(gh1,(2,npw_k*dtset%nspinor))
 876 
 877      ABI_ALLOCATE(kg_k,(3,npw_k))
 878      ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
 879      kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg)
 880      if (psps%useylm==1) then
 881        do ilm=1,mpsang*mpsang
 882          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
 883        end do
 884      end if
 885 
 886 !    Compute (k+G) and (k+q+G) vectors (only if useylm=1)
 887      nkpg=0;if (i2pert<natom+1) nkpg=3*dtset%nloalg(3)
 888      ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
 889      if (nkpg>0) then
 890        call mkkpg(kg_k,kpg_k,kpt,nkpg,npw_k)
 891      end if
 892 
 893 !    Compute nonlocal form factors ffnl at (k+G), for all atoms
 894      dimffnl=1
 895      ABI_ALLOCATE(ffnlk,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
 896      if (i2pert<natom+1) then
 897        ider=0
 898        call status(counter,dtfil%filstat,iexit,level,'call mkffnl  ')
 899        call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnlk,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
 900 &       ider,ider,psps%indlmn,kg_k,kpg_k,kpt,psps%lmnmax,psps%lnmax,psps%mpsang,&
 901 &       psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,&
 902 &       psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
 903      end if
 904 
 905 !    Load k-dependent part in the Hamiltonian datastructure
 906      call load_k_hamiltonian(gs_hamk,kpt_k=kpt,npw_k=npw_k,istwf_k=istwf_k,&
 907 &     kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnlk,compute_gbound=.true.)
 908 !    Load k+q-dependent part in the Hamiltonian datastructure
 909 !    call load_kprime_hamiltonian...  !! To be activated when q/=0
 910 
 911 !    Loop over bands
 912 
 913      do iband = 1,nband_k
 914 
 915        cwave0(:,:)=cg(:,1+(iband - 1)*npw_k*dtset%nspinor+icg0:&
 916 &       iband*npw_k*dtset%nspinor+icg0)
 917        cwavef3(:,:)=cg3(:,1+(iband-1)*npw_k*dtset%nspinor+icg0:&
 918 &       iband*npw_k*dtset%nspinor+icg0)
 919 
 920 !      Compute vtrial1 | cwafef3 >
 921        tim_fourwf = 0 ; weight = one
 922        call status(counter,dtfil%filstat,iexit,level,'call fourwf  ')
 923        call fourwf(cplex,vlocal1,cwavef3,gh1,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
 924 &       istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,&
 925 &       dtset%paral_kgb,tim_fourwf,weight,weight,&
 926 &       use_gpu_cuda=dtset%use_gpu_cuda)
 927 
 928 !      In case i2pert = phonon-type perturbation
 929 !      add first-order change in the nonlocal potential
 930        if (i2pert<natom+1) then
 931          signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1
 932          call status(counter,dtfil%filstat,iexit,level,'call nonlop  ')
 933          call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,&
 934 &         signs,dum_svectout,tim_nonlop,cwavef3,gvnl,iatom_only=i2pert)
 935          gh1(:,:) = gh1(:,:) + gvnl(:,:)
 936        end if
 937 
 938        ii = (iband-1)*npw_k*dtset%nspinor + icg0
 939        call dotprod_g(dotr,doti,istwf_k,npw_k,2,cg1(:,ii+1:ii+npw_k),gh1,mpi_enreg%me_g0,xmpi_comm_self)
 940 
 941 !      Compute vtrial1 | cwave0 >
 942        tim_fourwf = 0 ; weight = one
 943        call status(counter,dtfil%filstat,iexit,level,'call fourwf  ')
 944        call fourwf(cplex,vlocal1,cwave0,gh0,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
 945 &       istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,dtset%ngfft,npw_k,npw_k,n4,n5,n6,option,&
 946 &       dtset%paral_kgb,tim_fourwf,weight,weight,use_gpu_cuda=dtset%use_gpu_cuda)
 947 
 948 !      In case i2pert = phonon-type perturbation
 949 !      add first-order change in the nonlocal potential
 950        if (i2pert<natom+1) then
 951          signs=2 ; choice=2 ; nnlout=3 ; tim_nonlop = 0 ; paw_opt=0 ; cpopt=-1
 952          call status(counter,dtfil%filstat,iexit,level,'call nonlop  ')
 953          call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk,i2dir,dum,mpi_enreg,1,nnlout,paw_opt,&
 954 &         signs,dum_svectout,tim_nonlop,cwave0,gvnl,iatom_only=i2pert)
 955          gh0(:,:) = gh0(:,:) + gvnl(:,:)
 956        end if
 957 
 958 !      Compute the dft contribution to the Lagrange multiplier
 959 !      cwavef3 and cwave0 have been transferred to gh1 and gh0
 960 !      these vectors will be used to store the wavefunctions of band iband
 961 !      cg1 and gh0 contain the wavefunctions of band jband
 962 
 963        lagr = zero ; lagi = zero
 964        do jband = 1, nband_k
 965 
 966          ii = (jband - 1)*npw_k*dtset%nspinor + icg0
 967          jj = (iband - 1)*npw_k*dtset%nspinor + icg0
 968 
 969 !        dot1r and dot1i contain < u_mk | v^(1) | u_nk >
 970 !        dot2r and dot2i contain < u_nk^(1) | u_mk^(1) >
 971 !        m -> jband and n -> iband
 972 
 973          dot1r = zero ; dot1i = zero
 974          dot2r = zero ; dot2i = zero
 975          do ipw = 1, npw_k
 976            ii = ii + 1 ; jj = jj + 1
 977            dot1r = dot1r + cg(1,ii)*gh0(1,ipw) + cg(2,ii)*gh0(2,ipw)
 978            dot1i = dot1i + cg(1,ii)*gh0(2,ipw) - cg(2,ii)*gh0(1,ipw)
 979            dot2r = dot2r + cg1(1,jj)*cg3(1,ii) + &
 980 &           cg1(2,jj)*cg3(2,ii)
 981            dot2i = dot2i + cg1(1,jj)*cg3(2,ii) - &
 982 &           cg1(2,jj)*cg3(1,ii)
 983          end do  !  ipw
 984 
 985          lagr = lagr + dot1r*dot2r - dot1i*dot2i
 986          lagi = lagi + dot1r*dot2i + dot1i*dot2r
 987 
 988        end do    ! jband
 989 
 990        sumr = sumr + &
 991 &       dtset%wtk(ikpt)*occ(bantot+iband)*(dotr-lagr)
 992        sumi = sumi + &
 993 &       dtset%wtk(ikpt)*occ(bantot+iband)*(doti-lagi)
 994 
 995      end do   ! end loop over bands
 996 
 997      bantot = bantot + nband_k
 998      icg0 = icg0 + npw_k*dtset%nspinor*nband_k
 999      ikg = ikg + npw_k
1000 
1001      ABI_DEALLOCATE(cwave0)
1002      ABI_DEALLOCATE(cwavef3)
1003      ABI_DEALLOCATE(gh0)
1004      ABI_DEALLOCATE(gh1)
1005      ABI_DEALLOCATE(gvnl)
1006      ABI_DEALLOCATE(kg_k)
1007      ABI_DEALLOCATE(ylm_k)
1008      ABI_DEALLOCATE(ffnlk)
1009      ABI_DEALLOCATE(kpg_k)
1010 
1011    end do   ! end loop over k-points
1012 
1013  end do   ! end loop over spins
1014 
1015  if (xmpi_paral == 1) then
1016    buffer(1) = sumr ; buffer(2) = sumi
1017    call xmpi_sum(buffer,spaceComm,ierr)
1018    sumr = buffer(1) ; sumi = buffer(2)
1019  end if
1020 
1021 
1022  d3lo(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
1023 !d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
1024 
1025 !In some cases, the imaginary part is /= 0 because of the
1026 !use of time reversal symmetry
1027  d3lo(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero
1028 
1029  call destroy_hamiltonian(gs_hamk)
1030 
1031  ABI_DEALLOCATE(vlocal1)
1032  ABI_DEALLOCATE(wfraug)
1033 
1034  call status(0,dtfil%filstat,iexit,level,'exit          ')
1035 
1036 end subroutine pead_nl_resp