TABLE OF CONTENTS


ABINIT/m_rf2 [ Modules ]

[ Top ] [ Modules ]

NAME

 m_rf2

FUNCTION

 This module defines structures and provides procedures used to compute the 2nd order Sternheimer
 equation.

COPYRIGHT

  Copyright (C) 2015-2018 ABINIT group (LB,MT)
  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

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_rf2
29 
30  use defs_basis
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_hamiltonian
35  use m_cgtools
36 
37  use m_getgh1c,     only : getgh1c
38  use m_pawcprj,     only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_free, pawcprj_output
39  use m_getghc,      only : getghc
40  use m_nonlop,      only : nonlop
41  use m_getgh2c,     only : getgh2c
42 
43  implicit none
44 
45  private

m_rf2/rf2_accumulate_bands [ Functions ]

[ Top ] [ m_rf2 ] [ Functions ]

NAME

 rf2_init

FUNCTION

 Compute the scalar product < vi | v1j > and add it to rf2%lambda_mn
 If necessary, compute < vi | v2j > and add it to rf2%amn.

INPUTS

  rf2 : rf2_t object containing all rf2 data
  choice : 4 possible values, see "select case (choice)" below
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  mpi_enreg=information about MPI parallelization
  iband : band index for vi
  idir1  (used only if debug_mode/=0)  : direction of the 1st perturbation
  idir2  (used only if debug_mode/=0)  : direction of the 2nd perturbation
  ipert1 (used only if debug_mode/=0)  : 1st perturbation
  ipert2 (used only if debug_mode/=0)  : 2nd perturbation
  jband : band index for v1j and v2j
  debug_mode : if /=0 : all < vi | v1j > and < vi | v2j > are printed in std_out
  vi,v1j,v2j : input vectors

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

295 subroutine rf2_accumulate_bands(rf2,choice,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,&
296                                  jband,debug_mode,vi,v1j,v2j)
297 
298 
299 !This section has been created automatically by the script Abilint (TD).
300 !Do not modify the following lines by hand.
301 #undef ABI_FUNC
302 #define ABI_FUNC 'rf2_accumulate_bands'
303 !End of the abilint section
304 
305  implicit none
306 
307 !Arguments ---------------------------------------------
308 !scalars
309  integer,intent(in) :: choice,iband,idir1,idir2,ipert1,ipert2,jband,debug_mode
310  type(rf2_t),intent(inout) :: rf2
311  type(gs_hamiltonian_type),intent(in) :: gs_hamkq
312  type(MPI_type),intent(in) :: mpi_enreg
313 
314 !arrays
315  real(dp),intent(in) :: vi(2,rf2%size_wf),v1j(2,rf2%size_wf),v2j(2,rf2%size_wf)
316 
317 !Local variables ---------------------------------------
318 !scalars
319  integer :: nband_k,size_wf
320  real(dp) :: dotr,dot2r,doti,dot2i,factor
321  character(len=500) :: msg
322  character(len=15) :: bra_i,ket_j,op1,op2
323  character(len=2) :: pert1,pert2
324 
325 ! *************************************************************************
326 
327  DBG_ENTER("COLL")
328 
329  nband_k = rf2%nband_k
330  size_wf = rf2%size_wf
331  factor = one
332  if(rf2%ndir==1 .and. choice /= 3) factor = two ! in order to not compute same terms twice
333 
334  call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,vi,v1j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
335 
336  if(debug_mode/=0) then
337    if (ipert1 == gs_hamkq%natom+1) then
338      pert1 = "dk"
339    else
340      pert1 = "dE"
341    end if
342    if (ipert2 == gs_hamkq%natom+1) then
343      pert2 = "dk"
344    else
345      pert2 = "dE"
346    end if
347    select case (choice)
348     case(1) ! < u^(pert2) | H^(0) | u^(pert1) >
349       write(bra_i,'(2a,2(i1,a))') ' < du/',pert2,idir2,'(',iband,') | '
350       write(ket_j,'(2a,2(i1,a))') ' | du/',pert1,idir1,'(',jband,') > '
351       write(op1,'(a)')            '     H^(0)     '
352       write(op2,'(a)')            '     S^(0)     '
353     case(2) ! < u^(pert2) | H^(pert1) | u^(0) >
354       write(bra_i,'(2a,2(i1,a))') ' < du/',pert2,idir2,'(',iband,') | '
355       write(ket_j,'(a,i1,a)')     ' | u^(0) (',jband,') > '
356       write(op1,'(2a,i1,a)')      '     dH/',pert1,idir1,'    '
357       write(op2,'(2a,i1,a)')      '     dS/',pert1,idir1,'    '
358     case(3) ! < u^(0) | H^(pert1pert2) | u^(0) >
359       write(bra_i,'(a,i1,a)')     ' < u^(0) (',iband,') | '
360       write(ket_j,'(a,i1,a)')     ' | u^(0) (',jband,') > '
361       write(op1,'(2(2a,i1),a)')   'd^2H/(',pert1,idir1,' ',pert2,idir2,')'
362       write(op2,'(2(2a,i1),a)')   'd^2S/(',pert1,idir1,' ',pert2,idir2,')'
363     case(4) ! < u^(0) | H^(pert2) | u^(pert1) >
364       write(bra_i,'(a,i1,a)')     ' < u^(0) (',iband,') | '
365       write(ket_j,'(2a,2(i1,a))') ' | du/',pert1,idir1,'(',jband,') > '
366       write(op1,'(2a,i1,a)')      '     dH/',pert2,idir2,'    '
367       write(op2,'(2a,i1,a)')      '     dS/',pert2,idir2,'    '
368    end select
369    dot2r = dotr ; if (abs(dot2r)<tol9) dot2r = zero ! in order to hide the numerical noise
370    dot2i = doti ; if (abs(dot2i)<tol9) dot2i = zero ! in order to hide the numerical noise
371    write(msg,'(3a,2(a,es17.8E3))') bra_i,op1,ket_j,' = ',dot2r,',',dot2i
372    call wrtout(std_out,msg)
373  end if
374 
375  rf2%lambda_mn(:,iband+(jband-1)*nband_k) = factor*(/dotr,doti/) + rf2%lambda_mn(:,iband+(jband-1)*nband_k)
376 
377  if (choice == 1 .or. gs_hamkq%usepaw==1) then
378    call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,vi,v2j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
379 
380    if(debug_mode/=0) then
381      dot2r = dotr ; if (abs(dot2r)<tol9) dot2r = zero ! in order to hide the numerical noise
382      dot2i = doti ; if (abs(dot2i)<tol9) dot2i = zero ! in order to hide the numerical noise
383      write(msg,'(3a,2(a,es17.8E3))') bra_i,op2,ket_j,' = ',dot2r,',',dot2i
384      call wrtout(std_out,msg)
385    end if
386 
387    rf2%amn(:,iband+(jband-1)*nband_k) = factor*(/dotr,doti/) + rf2%amn(:,iband+(jband-1)*nband_k)
388 
389  end if ! end choice
390 
391  DBG_EXIT("COLL")
392 
393 end subroutine rf2_accumulate_bands

m_rf2/rf2_apply_hamiltonian [ Functions ]

[ Top ] [ m_rf2 ] [ Functions ]

NAME

 rf2_apply_hamiltonian

FUNCTION

 Apply the KS Hamiltonian (or derivative) to an input wave function.
 If asked, it also does some checks.

INPUTS

  rf2 : rf2_t object containing all rf2 data
  cg_jband (used only if debug_mode/=0) : array containing |u^(0)(jband)> for all bands
  cprj_jband(natom,nspinor*usecprj)= u^(0) wave functions for all bands
              projected with non-local projectors: cprj_jband=<p_i|u^(0)(jband)>
  cwave(2,size_wf) : input wave function |u>
  cwaveprj(natom,nspinor) : input wave function |u>
              projected with non-local projectors <p_i|u>
  eig0 : 0-order eigenvalue for the present wavefunction at k
  [enl]=optional (if not present, use hamk%ekb); non-local coeffs connecting projectors
  eig1_k_jband : first-order lagrange multipliers for the band j (Lambda^(1)_ji for all i)
  [ffnl1]=nonlocal form factors (needed for tests of ipert>=natom+10)
  [ffnl1_test]=nonlocal form factors used for tests (needed for tests of ipert>=natom+10)
  jband : band index of the input wavefunction
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  gvnl1(2,npw1*nspinor)=  part of <G|K^(1)+Vnl^(1)|C> not depending on VHxc^(1)              (sij_opt/=-1)
                       or part of <G|K^(1)+Vnl^(1)-lambda.S^(1)|C> not depending on VHxc^(1) (sij_opt==-1)
  idir=direction of the perturbation
  ipert=type of the perturbation of the Hamiltonian :
     ipert = 0           : GS calculation, call of getghc
     ipert = natom+1,2   : 1st order calculation, call of getgh1c
     ipert = natom+10,11 : 2nd order calculation, call of getgh2c
  ikpt=number of the k-point
  isppol=1 index of current spin component
  mkmem =number of k points trated by this node (GS data).
  mpi_enreg=information about MPI parallelization
  nband_k=number of bands at this k point for that spin polarization
  nsppol=1 for unpolarized, 2 for spin-polarized
  debug_mode : if /=0 : some tests are done (see NOTES below). Wrong results are printed in std_out
  prtvol=control print volume and debugging output (for getghc)
  rf_hamk_idir <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q (here q=0)
  size_cprj=size of a cprj array (=gs_hamkq%nspinor)
  size_wf=size of a wavefunction array (=gs_hamkq%npw_k*gs_hamkq%nspinor)

OUTPUT

  h_cwave(2,size_wf) : array containing H^(ipert)|cwave>
  s_cwave(2,size_wf) : array containing S^(ipert)|cwave>

NOTES

  * Tests are done if debug_mode/=0. In that case : cg_jband(:,jband,1) = |u^(0)(jband)>
    For ipert==natom+2 (electric field) : cg_jband(:,jband,2) = i * |du/dk_idir(jband)>
    According to ipert, we check that :
    -- ipert = 0 :
      < u^(0)(jband)| H^(0) | u^(0)(jband)> = eig0(jband)
    -- ipert = natom+1 or 2 :
      < u^(0)(jband)| ( H^(1) - (eig0(jband)+eig0(iband))/2 * S^(1) ) | u^(0)(iband)> = Lambda^(1)(jband,iband)
    --ipert >= natom+10 :
      < u^(0) | H^(2) | u^(0) > from getgh2c is equal to nonlop with signs=1
  * Use of cprj :
    The cprj array is computed in dfpt_looppert.F90. In the context of rf2 calculation, it contains
    <Proj_i^(0)|u^(0)> and <Proj_i^(1)|u^(0)> (dir=1,2 and 3) for all wavefunctions u^(0) or
    <Proj_i^(0)|u^(1)> and <Proj_i^(1)|u^(1)> (dir=1,2 and 3) for all 1st-order wavefunctions u^(1).
    Note that <Proj_i^(2)|u^(0)> is always computed on the fly.

PARENTS

      rf2_init

CHILDREN

SOURCE

467 subroutine rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave,cwaveprj,h_cwave,s_cwave,eig0,eig1_k_jband,&
468 &                                jband,gs_hamkq,gvnl1,idir,ipert,ikpt,isppol,mkmem,mpi_enreg,nband_k,nsppol,&
469 &                                debug_mode,prtvol,rf_hamk_idir,size_cprj,size_wf,&
470 &                                conj,enl,ffnl1,ffnl1_test) ! optional
471 
472 
473 !This section has been created automatically by the script Abilint (TD).
474 !Do not modify the following lines by hand.
475 #undef ABI_FUNC
476 #define ABI_FUNC 'rf2_apply_hamiltonian'
477 !End of the abilint section
478 
479  implicit none
480 
481 !Arguments ---------------------------------------------
482 !scalars
483  logical,intent(in),optional :: conj
484  integer,intent(in) :: idir,ipert,ikpt,isppol,jband,mkmem,nband_k,nsppol,debug_mode,prtvol,size_wf,size_cprj
485  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
486 ! type(rf2_t),intent(in) :: rf2
487  type(rf_hamiltonian_type),intent(inout),target :: rf_hamk_idir
488  type(MPI_type),intent(in) :: mpi_enreg
489 
490 !arrays
491  real(dp),intent(in),target :: cg_jband(2,size_wf*debug_mode*nband_k,2)
492  real(dp),intent(in),optional,target :: enl(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,gs_hamkq%dimekbq)
493  real(dp),intent(in),optional :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
494  real(dp),intent(in) :: eig0(nband_k),eig1_k_jband(2*nband_k)
495  real(dp),intent(inout) :: gvnl1(2,size_wf)
496  real(dp),intent(inout) :: cwave(2,size_wf),h_cwave(2,size_wf),s_cwave(2,size_wf)
497  type(pawcprj_type),intent(inout),target :: cprj_jband(:,:),cwaveprj(:,:)
498 
499 !Local variables ---------------------------------------
500 !scalars
501  integer,parameter :: berryopt=0,tim_getghc=0,tim_getgh1c=0,tim_getgh2c=0,tim_nonlop=0
502  integer,parameter :: alpha(9)=(/1,2,3,2,1,1,3,3,2/),beta(9)=(/1,2,3,3,3,2,2,1,1/)
503  integer :: choice,cpopt,iatom,iband,idirc,idir1,idir2,idir_dum,natom,nnlout,paw_opt,signs,sij_opt
504  integer :: opt_gvnl1,opt_gvnl2,optlocal,optnl,usevnl
505  logical :: compute_conjugate,has_cprj_jband,has_cwaveprj,pert_phon_elfd
506  real(dp) :: dotr,doti,dotr2,doti2,enlout(18*gs_hamkq%natom),tol_test
507  real(dp), pointer :: enl_ptr(:,:,:,:)
508  real(dp),allocatable,target :: enl_temp(:,:,:,:)
509  character(len=500) :: msg
510 
511 !arrays
512  real(dp),target :: cwave_empty(0,0),svectout_dum(0,0)
513  real(dp),allocatable :: gvnlc(:,:),iddk(:,:)
514  real(dp), ABI_CONTIGUOUS pointer :: cwave_i(:,:),cwave_j(:,:)
515  type(pawcprj_type),target :: cprj_empty(0,0)
516  type(pawcprj_type),pointer :: cprj_j(:,:)
517 ! *********************************************************************
518 
519  DBG_ENTER("COLL")
520 
521  compute_conjugate = .false.
522  if(present(conj)) compute_conjugate = conj
523 
524  ABI_UNUSED(ikpt)
525  ABI_UNUSED(isppol)
526  ABI_UNUSED(mkmem)
527  ABI_UNUSED(nsppol)
528 
529 !Check sizes
530  if (size(cprj_jband)/=0) then
531    if (size(cprj_jband)/=nband_k*gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj) then
532      write(msg,'(2(a,i10))') &
533      'Wrong cprj size! actual size = ',size(cprj_jband),&
534      ' good size = ',nband_k*gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj
535      MSG_BUG(msg)
536    end if
537  end if
538  if (size(cwaveprj)/=0) then
539    if (size(cwaveprj)/=gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj) then
540      write(msg,'(2(a,i10))') &
541      'Wrong cwaveprj size! actual size = ',size(cwaveprj),&
542      ' good size = ',gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj
543      MSG_BUG(msg)
544    end if
545  end if
546 
547  natom = gs_hamkq%natom
548 
549  pert_phon_elfd = .false.
550  if (ipert>natom+11.and.ipert<=2*natom+11) pert_phon_elfd = .true.
551  usevnl     = 0
552  opt_gvnl1  = 0
553  opt_gvnl2  = 0
554  optnl      = 2
555  sij_opt=1;if (gs_hamkq%usepaw==0) sij_opt=0
556  if (ipert==natom+2.or.(ipert==natom+11.and.gs_hamkq%usepaw==1).or.pert_phon_elfd) then
557    usevnl = 1
558    opt_gvnl1 = 2
559    opt_gvnl2 = 1
560  end if
561  optlocal = 1
562  if (pert_phon_elfd) then
563    optlocal = 0
564    sij_opt = 0
565    optnl = 1
566  end if
567  tol_test=tol8
568 
569 ! In the PAW case : manage cprj_jband, cwaveprj
570  has_cprj_jband=(gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.size(cprj_jband)/=0)
571  has_cwaveprj=(gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.size(cwaveprj)/=0)
572  cprj_j => cprj_empty
573 
574  if (debug_mode/=0) then
575    if (ipert/=0) then
576      write(msg,'(5(a,i4))') 'RF2 TEST rf2_apply_hamiltonian ipert = ',ipert,' idir =',idir,' isppol =',isppol,&
577      & ' ikpt = ',ikpt,' jband = ',jband
578    else
579      write(msg,'(3(a,i4))') 'RF2 TEST rf2_apply_hamiltonian ipert =    0 idir =    0 isppol =',isppol,&
580      & ' ikpt = ',ikpt,' jband = ',jband
581    end if
582    call wrtout(std_out,msg)
583  end if
584 
585 ! *******************************************************************************************
586 ! apply H^(0)
587 ! *******************************************************************************************
588  if (ipert == 0) then
589 
590    ABI_ALLOCATE(gvnlc,(2,size_wf))
591    gvnlc(:,:) = zero
592 
593 !  Test if < u^(0) | H^(0) | u^(0) > = eig0(jband)
594    if(debug_mode/=0) then
595      cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1)
596      if (has_cprj_jband) cprj_j => cprj_jband(:,1+(jband-1)*size_cprj:jband*size_cprj)
597      cpopt = -1+3*gs_hamkq%usecprj*gs_hamkq%usepaw
598      call getghc(cpopt,cwave_j,cprj_j,h_cwave,s_cwave,gs_hamkq,gvnlc,zero,mpi_enreg,1,prtvol,&
599                  sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
600 
601      call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
602      dotr = dotr - eig0(jband)
603      dotr = sqrt(dotr**2+doti**2)
604      if (dotr > tol_test) then
605        write(msg,'(a,es17.8E3)') 'RF2 TEST GETGHC : NOT PASSED dotr = ',dotr
606        call wrtout(ab_out,msg)
607        call wrtout(std_out,msg)
608      end if
609    end if ! end tests
610 
611    cpopt=-1;if (has_cwaveprj) cpopt=2
612    call getghc(cpopt,cwave,cwaveprj,h_cwave,s_cwave,gs_hamkq,gvnlc,zero,mpi_enreg,1,prtvol,&
613                sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
614    ABI_DEALLOCATE(gvnlc)
615 
616 ! *******************************************************************************************
617 ! apply H^(1)
618 ! *******************************************************************************************
619  else if (ipert<=natom+2) then
620 
621 !  Test if < u^(0) | ( H^(1) - eps^(0) S^(1) ) | u^(0) > = eig^(1)
622    if(debug_mode/=0) then
623      ABI_ALLOCATE(iddk,(2,size_wf))
624      cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1)
625      if (has_cprj_jband) cprj_j => cprj_jband(:,1+(jband-1)*size_cprj:jband*size_cprj)
626      iddk(:,:) = zero;if (ipert==natom+2) iddk(:,:)=cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,2)
627      call getgh1c(berryopt,cwave_j,cprj_j,h_cwave,cwave_empty,s_cwave,gs_hamkq,iddk,idir,ipert,zero,&
628                   mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamk_idir,sij_opt,tim_getgh1c,usevnl,conj=compute_conjugate)
629      do iband=1,nband_k
630        cwave_i => cg_jband(:,1+(iband-1)*size_wf:iband*size_wf,1)
631        call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_i,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
632        if (gs_hamkq%usepaw==1.and.ipert/=natom+2) then ! S^(1) is zero for ipert=natom+2
633          call dotprod_g(dotr2,doti2,gs_hamkq%istwf_k,size_wf,2,cwave_i,s_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
634          dotr = dotr - (eig0(iband)+eig0(jband))*dotr2/two
635          doti = doti - (eig0(iband)+eig0(jband))*doti2/two
636        end if
637        dotr = dotr - eig1_k_jband(1+2*(iband-1))
638        doti = doti - eig1_k_jband(2+2*(iband-1))
639        dotr = sqrt(dotr**2+doti**2)
640        if (dotr > tol_test) then
641          write(msg,'(4(a,i2),a,es17.8E3)') 'RF2 TEST GETGH1 : ipert=',ipert,' idir=',idir,&
642                                             ' jband=',jband,' iband=',iband,' NOT PASSED dotr = ',dotr
643          call wrtout(ab_out,msg)
644          call wrtout(std_out,msg)
645        end if
646      end do ! end iband
647      ABI_DEALLOCATE(iddk)
648    end if ! end tests
649 
650    call getgh1c(berryopt,cwave,cwaveprj,h_cwave,cwave_empty,s_cwave,gs_hamkq,gvnl1,idir,ipert,zero,&
651                 mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamk_idir,sij_opt,tim_getgh1c,usevnl,conj=compute_conjugate)
652 
653 ! *******************************************************************************************
654 ! apply H^(2)
655 ! *******************************************************************************************
656  else if (ipert==natom+10.or.ipert==natom+11.or.pert_phon_elfd) then
657 
658 ! *******************************************************************************************
659 !  Test if < u^(0) | H^(2) | u^(0) > from getgh2c is equal to nonlop with signs=1
660    if(debug_mode/=0.and.present(ffnl1).and.present(ffnl1_test)) then
661 
662      cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1)
663 
664      ABI_ALLOCATE(iddk,(2,size_wf))
665      iddk(:,:) = cwave_j(:,:)
666 
667      call getgh2c(cwave_j,cprj_empty,h_cwave,s_cwave,gs_hamkq,iddk,idir,ipert,zero,&
668                   mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamk_idir,sij_opt,tim_getgh2c,usevnl,&
669                   conj=compute_conjugate,optkin=0,enl=enl)
670      ABI_DEALLOCATE(iddk)
671 
672      call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
673 
674      idir1=alpha(idir);idir2=beta(idir)
675      idirc=3*(idir1-1)+idir2
676 
677      enlout = zero
678 
679 !    Change the pointer ffnl_k to ffnl1_test (idir_ffnl=0, for nonlop with signs=1)
680      call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1_test)
681 
682      signs=1
683      dotr2 = 0
684      doti2 = 0
685 
686 !    ************************
687 !    IPERT == NATOM+10
688 !    ************************
689      if (ipert==natom+10) then
690 
691        cpopt=-1; choice=8; paw_opt=1; if (gs_hamkq%usepaw==0) paw_opt=0 ; nnlout = 6
692 
693        call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
694 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
695        dotr2 = enlout(idir)
696 
697      end if ! IPERT == NATOM+10
698 
699 !    ************************
700 !    IPERT == NATOM+11
701 !    ************************
702      if (ipert==natom+11) then
703 
704        if (present(enl)) then
705          enl_ptr => enl
706        else if (associated(rf_hamk_idir%e1kbfr).and.associated(rf_hamk_idir%e1kbsc).and.optnl==2) then
707          ABI_ALLOCATE(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,rf_hamk_idir%cplex))
708          enl_temp(:,:,:,:) = rf_hamk_idir%e1kbfr(:,:,:,:) + rf_hamk_idir%e1kbsc(:,:,:,:)
709          enl_ptr => enl_temp
710        else if (associated(rf_hamk_idir%e1kbfr)) then
711          enl_ptr => rf_hamk_idir%e1kbfr
712        end if
713 
714 !      Compute application of dS/dk1 to cwavef
715 !      sum_{i,j} s_ij  < psi^(0) | d(|p_i><p_j|)/dk(idir1) | psi^(0) >
716        cpopt=-1 ; choice=5 ; paw_opt=3 ; nnlout=3
717        call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
718 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
719        dotr2 = enlout(idir1)
720 
721 !      Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
722 !      sum_{i,j} chi_ij(idir2) < psi^(0) | d(|p_i><p_j|)/dk(idir1) | psi^(0) >
723        cpopt=-1 ; choice=5 ; paw_opt=1 ; nnlout=3
724        call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
725 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum,enl=enl_ptr)
726        dotr2 = dotr2 + enlout(idir1)
727 
728 !      Compute derivatives due to projectors |d^2[p_i]/dk1dk2>,|d[p_i]/dk1>,|d[p_i]/dk2>
729 !      i * sum_{i,j} < psi^(0) | (d(|p_i><dp_j/dk(idir2)|)/dk(idir1) | psi^(0) >
730        cpopt=-1 ; choice=81 ; paw_opt=3 ; nnlout=18
731        call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
732 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
733        dotr2 = dotr2 - enlout(2*idirc  )
734        doti2 = doti2 + enlout(2*idirc-1)
735 
736      end if ! IPERT == NATOM+11
737 
738 !    ************************
739 !    PERT_PHON_ELFD
740 !    ************************
741      if (pert_phon_elfd) then
742 
743        if (present(enl)) then
744          enl_ptr => enl
745        else if (associated(rf_hamk_idir%e1kbfr).and.associated(rf_hamk_idir%e1kbsc).and.optnl==2) then
746          ABI_ALLOCATE(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,rf_hamk_idir%cplex))
747          enl_temp(:,:,:,:) = rf_hamk_idir%e1kbfr(:,:,:,:) + rf_hamk_idir%e1kbsc(:,:,:,:)
748          enl_ptr => enl_temp
749        else if (associated(rf_hamk_idir%e1kbfr)) then
750          enl_ptr => rf_hamk_idir%e1kbfr
751        end if
752 
753        iatom = gs_hamkq%atindx(ipert - natom - 11)  ! Atoms are type-sorted in enlout()
754 
755 !      Compute application of dS/dtau1 to i*d[cwavef]/dk2
756 !      sum_{i,j} s_ij < psi^(0) | d(|p_i><p_j|)/dtau(idir1) | i*psi^(k(idir2)) >
757        cpopt=-1 ; choice=2 ; paw_opt=3 ; nnlout=3*natom
758        call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
759 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
760        dotr2 = enlout(3*(iatom-1)+idir1)
761 
762 !      Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
763 !      sum_{i,j} chi_ij(idir2) < psi^(0) | d(|p_i><p_j|)/dtau(idir1) | psi^(0) >
764        cpopt=-1 ; choice=2 ; paw_opt=1 ; nnlout=3*natom
765        call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
766 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum,enl=enl_ptr)
767        dotr2 = dotr2 + enlout(3*(iatom-1)+idir1)
768 
769 !      Compute derivatives due to projectors |d^2[p_i]/dtau1dk2>,|d[p_i]/dtau1>,|d[p_i]/dk2>
770 !      i * sum_{i,j} < psi^(0) | (d(|p_i><dp_j/dk(idir2)|)/dtau(idir1) | psi^(0) >
771        cpopt=-1 ; choice=54 ; paw_opt=3 ; nnlout=18*natom
772        call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,&
773 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum)
774        dotr2 = dotr2 - enlout(18*(iatom-1)+2*idirc  )
775        doti2 = doti2 + enlout(18*(iatom-1)+2*idirc-1)
776 
777      end if ! PERT_PHON_ELFD
778 
779      dotr = dotr - dotr2
780      doti = doti - doti2
781      dotr = sqrt(dotr**2+doti**2)
782      if (dotr > 10*tol_test) then
783        write(msg,'(a,es17.8E3)') 'RF2 TEST GETGH2C : NOT PASSED dotr = ',dotr
784        call wrtout(ab_out,msg)
785        call wrtout(std_out,msg)
786      end if
787 
788 !    Change the pointer ffnl_k back to ffnl1 (idir_ffnl=4, for nonlop with signs=2)
789      call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1)
790 
791      if (associated(enl_ptr)) then
792        nullify(enl_ptr)
793      end if
794      if (allocated(enl_temp)) then
795        ABI_DEALLOCATE(enl_temp)
796      end if
797 
798    end if ! end tests
799 ! *******************************************************************************************
800 
801    call getgh2c(cwave,cwaveprj,h_cwave,s_cwave,gs_hamkq,gvnl1,idir,ipert,zero,&
802                 mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamk_idir,sij_opt,tim_getgh2c,usevnl,&
803                 conj=compute_conjugate,enl=enl)
804 
805  else
806 
807    h_cwave = zero
808    MSG_ERROR(" rf2_apply_hamiltonian can be used only for : 0<=ipert<=natom+2 and natom+10<=ipert<=2*natom+11.")
809    return
810 
811  end if
812 
813  DBG_EXIT("COLL")
814 
815 end subroutine rf2_apply_hamiltonian

m_rf2/rf2_destroy [ Functions ]

[ Top ] [ m_rf2 ] [ Functions ]

NAME

 rf2_init

FUNCTION

  Free all allocated arrays in a rf2_t object

INPUTS

 rf2 : the rf2_t object to destroy

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

840 subroutine rf2_destroy(rf2)
841 
842 
843 !This section has been created automatically by the script Abilint (TD).
844 !Do not modify the following lines by hand.
845 #undef ABI_FUNC
846 #define ABI_FUNC 'rf2_destroy'
847 !End of the abilint section
848 
849  implicit none
850 
851 !Arguments ---------------------------------------------
852 !scalars
853  type(rf2_t),intent(inout) :: rf2
854 !arrays
855 
856 !Local variables ---------------------------------------
857 !scalars
858 
859 ! *************************************************************************
860 
861  if (associated(rf2%RHS_Stern)) then
862    ABI_DEALLOCATE(rf2%RHS_Stern)
863  end if
864  if (allocated(rf2%dcwavef)) then
865    ABI_DEALLOCATE(rf2%dcwavef)
866  end if
867  if (allocated(rf2%amn)) then
868    ABI_DEALLOCATE(rf2%amn)
869  end if
870  if (allocated(rf2%lambda_mn)) then
871    ABI_DEALLOCATE(rf2%lambda_mn)
872  end if
873 
874 end subroutine rf2_destroy

m_rf2/rf2_getidir [ Functions ]

[ Top ] [ m_rf2 ] [ Functions ]

NAME

 rf2_getidir

FUNCTION

  Get the direction of the 2nd order perturbation according to inputs idir1 and idir2

INPUTS

  idir1 : index of the 1st direction (1<=idir1<=3)
  idir2 : index of the 2nd direction (1<=idir2<=3)

OUTPUT

  idir : integer between 1 and 9

NOTES

PARENTS

      dfpt_looppert,dfpt_scfcv,rf2_init

CHILDREN

SOURCE

158 subroutine rf2_getidir(idir1,idir2,idir)
159 
160 
161 !This section has been created automatically by the script Abilint (TD).
162 !Do not modify the following lines by hand.
163 #undef ABI_FUNC
164 #define ABI_FUNC 'rf2_getidir'
165 !End of the abilint section
166 
167  implicit none
168 
169 !Arguments ---------------------------------------------
170 !scalars
171  integer,intent(in) :: idir1,idir2
172  integer,intent(out) :: idir
173  integer :: dir_from_dirs(9) = (/1,6,5,9,2,4,8,7,3/)
174 
175 ! *************************************************************************
176 
177  idir = dir_from_dirs(3*(idir1-1)+idir2)
178 
179 end subroutine rf2_getidir

m_rf2/rf2_getidirs [ Functions ]

[ Top ] [ m_rf2 ] [ Functions ]

NAME

 rf2_getidirs

FUNCTION

  Get directions of the 1st and 2nd perturbations according to the input idir

INPUTS

  idir : integer between 1 and 9

OUTPUT

  idir1 : index of the 1st direction (1<=idir1<=3)
  idir2 : index of the 2nd direction (1<=idir2<=3)

NOTES

PARENTS

      dfpt_looppert,dfpt_scfcv,rf2_init

CHILDREN

SOURCE

207 subroutine rf2_getidirs(idir,idir1,idir2)
208 
209 
210 !This section has been created automatically by the script Abilint (TD).
211 !Do not modify the following lines by hand.
212 #undef ABI_FUNC
213 #define ABI_FUNC 'rf2_getidirs'
214 !End of the abilint section
215 
216  implicit none
217 
218 !Arguments ---------------------------------------------
219 !scalars
220  integer,intent(in) :: idir
221  integer,intent(out) :: idir1,idir2
222 
223 ! *************************************************************************
224 
225  select case(idir)
226 !  Diagonal terms :
227    case(1)
228      idir1 = 1
229      idir2 = 1
230    case(2)
231      idir1 = 2
232      idir2 = 2
233    case(3)
234      idir1 = 3
235      idir2 = 3
236 !  Upper triangular terms :
237    case(4)
238      idir1 = 2
239      idir2 = 3
240    case(5)
241      idir1 = 1
242      idir2 = 3
243    case(6)
244      idir1 = 1
245      idir2 = 2
246 !  Lower triangular terms :
247    case(7)
248      idir1 = 3
249      idir2 = 2
250    case(8)
251      idir1 = 3
252      idir2 = 1
253    case(9)
254      idir1 = 2
255      idir2 = 1
256  end select
257 
258 end subroutine rf2_getidirs

m_rf2/rf2_t [ Types ]

[ Top ] [ m_rf2 ] [ Types ]

NAME

  rf2_t

FUNCTION

  Datatype gathering all needed data for the computation of the 2nd order Sternheimer equation.

SOURCE

 59  type,public :: rf2_t
 60 
 61 !scalars
 62   integer :: ndir ! number of directions to consider
 63   integer :: nband_k ! number of bands
 64   integer :: size_wf ! number of components in a wavefunction
 65   integer :: size_cprj ! number of components of a cprj variable (i.e. <p_i|WF>)
 66 
 67 !arrays
 68   integer :: iperts(2) ! perturbations to compute
 69    ! ipert = natom + 10 :
 70    !   iperts(1) = iperts(2) = natom+1 (wavevector k)
 71    !
 72    ! ipert = natom + 11 :
 73    !   iperts(1) = natom+1 (wavevector k)
 74    !   iperts(2) = natom+2 (electric field)
 75 
 76   integer :: idirs(2) ! directions of the perturbations (ndir=1 : idirs(1)=idirs(2) , ndir=2 : idirs(1)/=idirs(2))
 77 
 78   real(dp), ABI_CONTIGUOUS pointer :: RHS_Stern(:,:) => null() ! need a pointer because is a ptr target
 79    ! Right-hand side of the 2nd order Sternheimer equation, for every bands.
 80    ! Namely, for a band "n" :
 81    ! |(RHS_Stern)_n> = (H^(2)-epsilon_n S^(2)) |u^(0)_n> + 2(H^(1)-epsilon_n S^(1))|u^(1)_n>
 82    !                   - 2 sum_m ( lambda_mn^(1) ( S^(1) |u^(0)_m> + S^(0) |u^(1)_m> )
 83    !                      ( /!\ : in the sum, the index m runs over occupied bands )
 84    ! where :
 85    !  - epsilon_n = eigenvalue of the GS Hamiltonian
 86    !  - lambda_mn^(1) = <u^(0)_m| (H^(1)-(epsilon_n+epsilon_m)/2 S^(1) |u^(0)_n> (1st order Lagrange multiplier)
 87    !  - X^(2) = d^2X/(dk_dir1 dlambda_dir2)
 88    !  - 2X^(1)Y^(1) = dX/dk_dir1 dY/dlambda_dir2 + dX/dlambda_dir2 dY/dk_dir1
 89    !  lambda_dir2 = k_dir2 (ipert==natom+10) , E_dir2 (ipert==natom+11)
 90    ! Note : P_c^* will be apply to |(RHS_Stern)_n> in dfpt_cgwf.F90
 91    ! **
 92    ! Computed in "rf2_init"
 93 
 94   real(dp),allocatable :: amn(:,:)
 95    ! Scalar needed for the "orhtonormalization condition", see "dcwavef(:,:)" below
 96    ! Namely :
 97    ! A_mn = <u^(0)_m| S^(2) |u^(0)_n> + 2 <u^(1)_m| S^(0) |u^(1)_n>
 98    !    + 2 <u^(1)_m| S^(1) |u^(0)_n> + 2 <u^(0)_m| S^(1) |u^(1)_n>
 99    ! **
100    ! Computed in "rf2_init", stored only for testing purposes
101 
102   real(dp),allocatable :: dcwavef(:,:)
103    ! Vector needed to enforce the "orhtonormalization condition" on 2nd order wave functions.
104    ! Namely :
105    ! |dcwavef_n> = -1/2 sum_m A_mn |u^(0)_m>
106    ! **
107    ! Computed in "rf2_init"
108 
109   real(dp),allocatable :: lambda_mn(:,:)
110    ! 2nd order Lagrange multiplier.
111    ! Namely :
112    ! lambda_mn =  <u^(0)_m| H^(2) |u^(0)_n> + 2 <u^(1)_m| H^(0) |u^(1)_n>
113    !          + 2 <u^(1)_m| H^(1) |u^(0)_n> + 2 <u^(0)_m| H^(1) |u^(1)_n>
114    !          - A_mn * (epsilon_m + epsilon_n) / 2
115    ! **
116    ! Computed in "rf2_init"
117 
118  end type rf2_t
119 
120  public :: rf2_getidir
121  public :: rf2_getidirs
122  public :: rf2_accumulate_bands
123  public :: rf2_apply_hamiltonian
124  public :: rf2_destroy