TABLE OF CONTENTS


ABINIT/dfpt_nselt [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_nselt

FUNCTION

 This routine compute the non-stationary expression for the
 second derivative of the total energy, wrt strain for a whole row of
 mixed strain derivatives.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DRH, XG, DCA, GMR, MM, AR, MV, 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
  cplex: if 1, real space 1-order functions on FFT grid are REAL,
    if 2, COMPLEX
  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha)
  effmass_free=effective mass for electrons (1. in common case)
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  idir=direction of the perturbation
  ipert=type of the perturbation
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the
     storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points in the reduced BZ
  kxc(nfft,nkxc)=exchange and correlation kernel
  mband=maximum number of bands
  mgfft=maximum size of 1D FFTs
  mkmem =number of k points treated by this node.
  mk1mem =number of k points treated by this node  (RF data).
  mpert =maximum number of ipert
  mpi_enreg=information about MPI parallelization
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  maximum dimension for q points in grids for nonlocal form factors
  natom=number of atoms in cell.
  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT,
    see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt_rbz=number of k points in the reduced BZ for this perturbation
  nkxc=second dimension of the kxc array. If /=0,
   the exchange-correlation kernel must be computed.
  nloalg(3)=governs the choice of the algorithm for non-local operator.
  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
  nspden=number of spin-density components
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym1=number of symmetry elements in space group consistent with
    perturbation
  ntypat=number of types of atoms in unit cell.
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band
   and k in the reduced Brillouin zone (usually =2)
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  prtbbb=if 1, band-by-band decomposition (also dim of d2bbb)
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  qphon(3)=reduced coordinates for the phonon wavelength
  rhog(2,nfft)=array for Fourier transform of GS electron density
  rhor(nfft,nspden)=GS electron density in electrons/bohr**3.
  rhor1(cplex*nfft,nspden)=RF electron density in electrons/bohr**3.
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
  typat(natom)=type integer for each atom in cell
  ucvol=unit cell volume in bohr**3.
  wtk_rbz(nkpt_rbz)=weight assigned to each k point in the reduced BZ
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang)= real spherical harmonics for each G and k point
  ylm1(mpw1*mk1mem,mpsang*mpsang)= real spherical harmonics for each G and k+q point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k point
  ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical for each G and k+g point

OUTPUT

  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
  d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
       second order derivatives
  d2lo(2,3,mpert,3,mpert)=local contributions to the 2DTEs
  d2nl(2,3,mpert,3,mpert)=non-local contributions to the 2DTEs

PARENTS

      dfpt_scfcv

CHILDREN

      destroy_hamiltonian,dfpt_mkcore,dfpt_mkvxcstr,dfpt_nsteltwf,dotprod_vn
      hartrestr,init_hamiltonian,stresssym,vlocalstr,wrtout,xmpi_barrier

SOURCE

101 #if defined HAVE_CONFIG_H
102 #include "config.h"
103 #endif
104 
105 #include "abi_common.h"
106 
107 
108 subroutine dfpt_nselt(blkflg,cg,cg1,cplex,&
109 & d2bbb,d2lo,d2nl,ecut,ecutsm,effmass_free,&
110 & gmet,gprimd,gsqcut,idir,&
111 & ipert,istwfk_rbz,kg,kg1,kpt_rbz,kxc,mband,mgfft,&
112 & mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,mpw1,&
113 & natom,nband_rbz,nfft,ngfft,&
114 & nkpt_rbz,nkxc,nloalg,npwarr,npwar1,nspden,nspinor,nsppol,&
115 & nsym1,ntypat,occ_rbz,&
116 & paral_kgb,ph1d,prtbbb,psps,qphon,rhog,&
117 & rhor,rhor1,rmet,rprimd,symrc1,typat,ucvol,&
118 & wtk_rbz,&
119 & xred,ylm,ylm1,ylmgr,ylmgr1)
120 
121  use defs_basis
122  use defs_datatypes
123  use defs_abitypes
124  use m_errors
125  use m_profiling_abi
126  use m_xmpi
127 
128  use m_cgtools,    only : dotprod_vn
129  use m_pawtab,     only : pawtab_type
130  use m_hamiltonian,only : init_hamiltonian,destroy_hamiltonian,gs_hamiltonian_type
131 
132 !This section has been created automatically by the script Abilint (TD).
133 !Do not modify the following lines by hand.
134 #undef ABI_FUNC
135 #define ABI_FUNC 'dfpt_nselt'
136  use interfaces_14_hidewrite
137  use interfaces_32_util
138  use interfaces_41_geometry
139  use interfaces_72_response, except_this_one => dfpt_nselt
140 !End of the abilint section
141 
142  implicit none
143 
144 !Arguments -------------------------------
145 !scalars
146  integer,intent(in) :: cplex,idir,ipert,mband,mgfft,mk1mem
147  integer,intent(in) :: mkmem,mpert,mpsang,mpw,mpw1,natom,nfft,nkpt_rbz
148  integer,intent(in) :: nkxc,nspden,nspinor,nsppol,nsym1,ntypat
149  integer,intent(in) :: paral_kgb,prtbbb
150  real(dp),intent(in) :: ecut,ecutsm,effmass_free,gsqcut,ucvol
151  type(MPI_type),intent(in) :: mpi_enreg
152  type(pseudopotential_type),intent(in) :: psps
153 !arrays
154  integer,intent(in) :: istwfk_rbz(nkpt_rbz)
155  integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
156  integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18)
157  integer,intent(in) :: nloalg(3),npwar1(nkpt_rbz),npwarr(nkpt_rbz)
158  integer,intent(in) :: symrc1(3,3,nsym1),typat(natom)
159  integer,intent(inout) :: blkflg(3,mpert,3,mpert)
160  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol)
161  real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)
162  real(dp),intent(in) :: gmet(3,3)
163  real(dp),intent(in) :: gprimd(3,3),kpt_rbz(3,nkpt_rbz),kxc(nfft,nkxc)
164  real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol)
165  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qphon(3),rhog(2,nfft)
166  real(dp),intent(in) :: rhor(nfft,nspden)
167  real(dp),intent(in) :: rhor1(cplex*nfft,nspden),rmet(3,3),rprimd(3,3)
168  real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom)
169  real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
170  real(dp),intent(in) :: ylm1(mpw1*mk1mem,mpsang*mpsang*psps%useylm)
171  real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm)
172  real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*psps%useylm)
173  real(dp),intent(out) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb)
174  real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert)
175  real(dp),intent(inout) :: d2nl(2,3,mpert,3,mpert)
176 
177 !Local variables-------------------------------
178 !scalars
179  integer :: ban2tot,bantot,bd2tot_index,bdtot_index
180  integer :: icg,icg1,idir1,ifft,ii,ikg,ikg1,ikpt,comm
181  integer :: ilm,ipert1,ispden,isppol,istr1,istwf_k
182  integer :: mbd2kpsp,mbdkpsp,me,n1,n2,n3,n3xccc,n4,n5,n6
183  integer :: nband_k,nfftot,npw1_k,npw_k,option
184  real(dp) :: doti,dotr
185  real(dp) :: wtk_k
186  character(len=500) :: message
187  type(gs_hamiltonian_type) :: gs_hamk
188 !arrays
189  integer :: ikpt_fbz(3)
190  integer,allocatable :: kg1_k(:,:),kg_k(:,:)
191  real(dp) :: kpoint(3),restr(6),dummy(0,0)
192  real(dp),allocatable :: d2bbb_k(:,:,:,:),d2nl_k(:,:,:)
193  real(dp),allocatable :: occ_k(:)
194  real(dp),allocatable :: vhartr01(:),vpsp1(:),vxc1(:,:),xccc3d1(:),ylm1_k(:,:)
195  real(dp),allocatable :: ylm_k(:,:),ylmgr1_k(:,:,:),ylmgr_k(:,:,:)
196  type(pawtab_type) :: pawtab_dum(0)
197 
198 ! *********************************************************************
199 
200 !Init me
201  comm = mpi_enreg%comm_cell
202  me   = mpi_enreg%me_kpt
203 
204 !Unit numbers
205 
206 !Zero only portion of nonlocal matrix to be computed here
207  d2nl(:,:,natom+3:natom+4,idir,ipert)=zero
208  bdtot_index=0
209  bd2tot_index=0
210  icg=0
211  icg1=0
212  mbdkpsp=mband*nkpt_rbz*nsppol
213  mbd2kpsp=2*mband**2*nkpt_rbz*nsppol
214 
215 !Update list of computed matrix elements
216  if((ipert==natom+3) .or. (ipert==natom+4)) then
217 !  Eventually expand when strain coupling to other perturbations is
218 !  implemented
219    do ipert1=natom+3,natom+4
220      do idir1=1,3
221        blkflg(idir1,ipert1,idir,ipert)=1
222      end do
223    end do
224  end if
225 
226 !allocate(enl1nk(mbdkpsp))
227  ABI_ALLOCATE(d2bbb_k,(2,3,mband,mband*prtbbb))
228  ABI_ALLOCATE(d2nl_k,(2,3,mpert))
229 
230 
231  ABI_ALLOCATE(kg_k,(3,mpw))
232  ABI_ALLOCATE(kg1_k,(3,mpw1))
233 
234 
235  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
236  n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
237  nfftot=n1*n2*n3
238 
239 !Initialize Hamiltonian (k-independent terms) - NCPP only
240  call init_hamiltonian(gs_hamk,psps,pawtab_dum,nspinor,nsppol,nspden,natom,&
241 & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,ph1d=ph1d)
242 
243  bantot = 0
244  ban2tot = 0
245 
246 !LOOP OVER SPINS
247  do isppol=1,nsppol
248 
249    if (nsppol/=1) then
250      write(message,*)' ****  In dfpt_nselt for isppol=',isppol
251      call wrtout(std_out,message,'COLL')
252    end if
253 
254    ikg=0
255    ikg1=0
256 
257    ikpt_fbz(1:3)=0
258 
259 !  BIG FAT k POINT LOOP
260    do ikpt=1,nkpt_rbz
261 
262      nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
263      istwf_k=istwfk_rbz(ikpt)
264      npw_k=npwarr(ikpt)
265      npw1_k=npwar1(ikpt)
266      kpoint(:)=kpt_rbz(:,ikpt)
267 
268      bantot = bantot + nband_k
269      ban2tot = ban2tot + 2*nband_k**2
270 
271      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
272        bdtot_index=bdtot_index+nband_k
273        bd2tot_index=bd2tot_index+2*nband_k**2
274 !      Skip the rest of the k-point loop
275        cycle
276      end if
277 
278      ABI_ALLOCATE(occ_k,(nband_k))
279 
280      ABI_ALLOCATE(ylm_k,(npw_k,mpsang*mpsang))
281      ABI_ALLOCATE(ylm1_k,(npw1_k,mpsang*mpsang))
282      if (ipert==natom+3.or.ipert==natom+4) then
283        ABI_ALLOCATE(ylmgr_k,(npw_k,3,mpsang*mpsang))
284        ABI_ALLOCATE(ylmgr1_k,(npw1_k,3,mpsang*mpsang))
285      end if
286 
287 !    enl1_k(:)=zero
288      d2nl_k(:,:,:)=zero
289      if(prtbbb==1)d2bbb_k(:,:,:,:)=zero
290      occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
291 
292      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
293      if (psps%useylm==1) then
294        do ilm=1,mpsang*mpsang
295          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
296        end do
297        if (ipert==natom+3.or.ipert==natom+4) then
298          do ilm=1,mpsang*mpsang
299            do ii=1,3
300              ylmgr_k(1:npw_k,ii,ilm)=ylmgr(1+ikg:npw_k+ikg,ii,ilm)
301            end do
302          end do
303        end if
304      end if
305 
306      wtk_k=wtk_rbz(ikpt)
307 
308      kg1_k(:,:) = 0
309 
310      kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
311      if (psps%useylm==1) then
312        do ilm=1,mpsang*mpsang
313          ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
314        end do
315        if (ipert==natom+3.or.ipert==natom+4) then
316          do ilm=1,mpsang*mpsang
317            do ii=1,3
318              ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm)
319            end do
320          end do
321        end if
322      end if
323 
324 !    Compute the eigenvalues, wavefunction,
325 !    contributions to kinetic energy, nonlocal energy, forces,
326 !    and update of rhor1 to this k-point and this spin polarization.
327 
328 !    Note that dfpt_nsteltwf is called with kpoint, while kpt is used inside dfpt_vtowfk
329      call dfpt_nsteltwf(cg,cg1,d2nl_k,ecut,ecutsm,effmass_free,gs_hamk,icg,icg1,ikpt,isppol,&
330 &     istwf_k,kg_k,kg1_k,kpoint,mband,mkmem,mk1mem,mpert,mpi_enreg,mpw,mpw1,natom,nband_k,&
331 &     npw_k,npw1_k,nspinor,nsppol,ntypat,occ_k,psps,rmet,wtk_k,ylm_k,ylmgr_k)
332      d2nl(:,:,:,idir,ipert)=d2nl(:,:,:,idir,ipert)+d2nl_k(:,:,:)
333      if(prtbbb==1)then
334        d2bbb(:,:,idir,ipert,:,:) = d2bbb(:,:,idir,ipert,:,:) + &
335 &       d2bbb_k(:,:,:,:)
336      end if
337 
338      ABI_DEALLOCATE(occ_k)
339 
340 !    Keep track of total number of bands (all k points so far, even for
341 !    k points not treated by me)
342      bdtot_index=bdtot_index+nband_k
343      bd2tot_index=bd2tot_index+2*nband_k**2
344 
345 !    Shift array memory
346      if (mkmem/=0) then
347        icg=icg+npw_k*nspinor*nband_k
348        ikg=ikg+npw_k
349      end if
350      if (mk1mem/=0) then
351        icg1=icg1+npw1_k*nspinor*nband_k
352        ikg1=ikg1+npw1_k
353      end if
354      ABI_DEALLOCATE(ylm_k)
355      ABI_DEALLOCATE(ylm1_k)
356      if (ipert==natom+3.or.ipert==natom+4)  then
357        ABI_DEALLOCATE(ylmgr_k)
358        ABI_DEALLOCATE(ylmgr1_k)
359      end if
360 
361    end do ! End big k point loop
362  end do ! End loop over spins
363 
364  if(xmpi_paral==1)then
365    call xmpi_barrier(comm)
366    call wrtout(std_out,' dfpt_nselt: loop on k-points and spins done in parallel','COLL')
367  end if
368 
369 !Treat now varying occupation numbers
370 !if(occopt>=3 .and. occopt <=8) then
371 !SUPPRESSED metallic coding of vtorho
372 
373 !Treat fixed occupation numbers
374 !else
375 
376 !Accumulation over parallel processed now carried out for all terms
377 !in dfpt_nstdy.f
378 
379 !End of test on varying or fixed occupation numbers
380 !end if
381 
382 !The imaginary part of d2nl will be must be set to zero here since
383 !time-reversal symmetry will always be true for the strain peturbation.
384 !The symmetry-reduced kpt set will leave a non-zero imaginary part.
385 
386  d2nl(2,:,natom+3:natom+4,idir,ipert)=zero
387 
388 !Symmetrize the non-local contributions,
389 !as was needed for the stresses in a ground-state calculation
390 
391  if (nsym1>1) then
392 !  Pack like symmetric-storage cartesian stress tensor
393    ii=0
394    do ipert1=natom+3,natom+4
395      do idir1=1,3
396        ii=ii+1
397        restr(ii)=d2nl(1,idir1,ipert1,idir,ipert)
398      end do
399    end do
400 !  Do the symmetrization using the ground state routine
401    call stresssym(gprimd,nsym1,restr,symrc1)
402 !  Unpack symmetrized stress tensor
403    ii=0
404    do ipert1=natom+3,natom+4
405      do idir1=1,3
406        ii=ii+1
407        d2nl(1,idir1,ipert1,idir,ipert)=restr(ii)
408      end do
409    end do
410  end if !nsym>1
411 
412 !----------------------------------------------------------------------------
413 !Now, treat the local contribution
414 
415  ABI_ALLOCATE(vpsp1,(cplex*nfft))
416  n3xccc=0
417  if(psps%n1xccc/=0)n3xccc=nfft
418  ABI_ALLOCATE(xccc3d1,(cplex*n3xccc))
419  ABI_ALLOCATE(vxc1,(cplex*nfft,nspden))
420  ABI_ALLOCATE(vhartr01,(nfft))
421  xccc3d1(:)=zero
422 
423 !Double loop over strain perturbations
424  do ipert1=natom+3,natom+4
425    do idir1=1,3
426      if(ipert1==natom+3) then
427        istr1=idir1
428      else
429        istr1=idir1+3
430      end if
431 
432 !    Get first-order local potential.
433      call vlocalstr(gmet,gprimd,gsqcut,istr1,mgfft,mpi_enreg,&
434 &     psps%mqgrid_vl,natom,gs_hamk%nattyp,nfft,ngfft,ntypat,paral_kgb,ph1d,psps%qgrid_vl,&
435 &     ucvol,psps%vlspl,vpsp1)
436 
437 !    Get first-order hartree potential.
438      call hartrestr(gsqcut,idir1,ipert1,mpi_enreg,natom,nfft,ngfft,&
439 &     paral_kgb,rhog,rprimd,vhartr01)
440 
441 !    Get first-order exchange-correlation potential
442      if(psps%n1xccc/=0)then
443        call dfpt_mkcore(cplex,idir1,ipert1,natom,ntypat,n1,psps%n1xccc,&
444 &       n2,n3,qphon,rprimd,typat,ucvol,psps%xcccrc,psps%xccc1d,xccc3d1,xred)
445      end if ! psps%n1xccc/=0
446 
447      option=0
448      call dfpt_mkvxcstr(cplex,idir1,ipert1,kxc,mpi_enreg,natom,nfft,ngfft,&
449 &     dummy,dummy,nkxc,nspden,n3xccc,option,paral_kgb,qphon,rhor,rhor1,rprimd,&
450 &     0,0,vxc1,xccc3d1)
451 
452 !    Combines density j2 with local potential j1
453      do ispden=1,min(nspden,2)
454        do ifft=1,cplex*nfft
455          vxc1(ifft,ispden)=vxc1(ifft,ispden)+vpsp1(ifft)+vhartr01(ifft)
456        end do
457      end do
458      call dotprod_vn(cplex,rhor1,dotr,doti,nfft,nfftot,nspden,2,vxc1,ucvol)
459      write(std_out,*)
460      d2lo(1,idir1,ipert1,idir,ipert)=dotr
461      d2lo(2,idir1,ipert1,idir,ipert)=doti
462    end do ! istr1
463  end do ! ipert1
464 
465  call destroy_hamiltonian(gs_hamk)
466 
467  ABI_DEALLOCATE(vxc1)
468  ABI_DEALLOCATE(xccc3d1)
469  ABI_DEALLOCATE(vhartr01)
470 
471  ABI_DEALLOCATE(d2bbb_k)
472  ABI_DEALLOCATE(d2nl_k)
473  ABI_DEALLOCATE(kg_k)
474  ABI_DEALLOCATE(kg1_k)
475  ABI_DEALLOCATE(vpsp1)
476 
477 end subroutine dfpt_nselt