TABLE OF CONTENTS


ABINIT/m_paw_dfptnl [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_dfptnl

FUNCTION

  This module contains several routines used to compute PAW contributions to a 3rd-order energy
   or 2nd-order PAW occupancies.

COPYRIGHT

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

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 MODULE m_paw_dfptnl
25 
26  use defs_basis
27  use defs_abitypes
28  use m_abicore
29  use m_errors
30  use m_xmpi, only : xmpi_comm_self,xmpi_sum
31 
32  use m_pawang,     only : pawang_type
33  use m_pawrad,     only : pawrad_type,simp_gen
34  use m_pawtab,     only : pawtab_type
35  use m_paw_an,     only : paw_an_type
36  use m_pawrhoij,   only : pawrhoij_type
37  use m_pawcprj,    only : pawcprj_type
38  use m_paw_denpot, only : pawdensities
39  use m_paral_atom, only : get_my_atmtab, free_my_atmtab
40 
41  implicit none
42 
43  private
44 
45 !public procedures.
46  public :: paw_dfptnl_energy   ! Compute the XC PAW on-site contributions to a 3rd-order energy
47  public :: paw_dfptnl_xc       ! Compute a contribution of the 3rd-derivative of XC energy of ONE PAW sphere
48  public :: paw_dfptnl_accrhoij ! Accumulate the 2nd order PAW quantities rhoij^(2)
49 
50 CONTAINS  !========================================================================================

ABINIT/paw_dfptnl_accrhoij [ Functions ]

[ Top ] [ Functions ]

NAME

 paw_dfptnl_accrhoij

FUNCTION

 Accumulate the 2nd order PAW quantities rhoij^(2) (augmentation occupancies)
 This routine is similar to pawaccrhoij.F90 but is implemented independently
 in order to not overload the original routine.

INPUTS

  atindx(natom)=index table for atoms (sorted-->random), inverse of atindx.
  cplex: if 1, WFs (or 1st-order WFs) are REAL, if 2, COMPLEX
  cwaveprj0_pert1(natom,nspinor) = wave function at given n,k projected with non-local projectors:
                                  cwaveprj0%cp    =<p_i|
                                  cwaveprj0%dcp(1)=<p_i^(pert1)|Cnk>
  cwaveprj0_pert2(natom,nspinor) = wave function at given n,k projected with non-local projectors:
                                  cwaveprj0%cp    =<p_i|Cnk>
                                  cwaveprj0%dcp(1)=<p_i^(pert1)|Cnk>
  cwaveprj1_pert12(natom,nspinor)= 1st order wave function at given n,k projected with non-local projectors:
                                  cwaveprj1%cp    =<p_i|Cnk^(pert1)>
                                  cwaveprj1%dcp(1)=<p_i^(pert2)|Cnk^(pert1)>
  cwaveprj1_pert21(natom,nspinor)= 1st order wave function at given n,k projected with non-local projectors:
                                  cwaveprj1%cp    =<p_i|Cnk^(pert2)>
                                  cwaveprj1%dcp(1)=<p_i^(pert1)|Cnk^(pert2)>
  ipert1=index of the first perturbation
  ipert2=index of the second perturbation
  isppol=index of current spin component
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell
  nspinor=number of spinorial components (on current proc)
  occ_k=occupation number for current band n,k
  wtk_k=weight assigned to current k-point

SIDE EFFECTS

  pawrhoij(natom) <type(pawrhoij_type)>= 2-nd order paw rhoij occupancies and related data
  On output, has been updated with the contribution of current n,k
        pawrhoij(:)%rhoij_(lmn2_size,nspden) (non symetrized)

PARENTS

      paw_dfptnl_pert

CHILDREN

      free_my_atmtab,get_my_atmtab

SOURCE

536  subroutine paw_dfptnl_accrhoij(atindx,cplex,cwaveprj0_pert1,cwaveprj0_pert2,&
537 &                       cwaveprj1_pert12,cwaveprj1_pert21,ipert1,ipert2,isppol,my_natom,natom,&
538 &                       nspinor,occ_k,pawrhoij,wtk_k,&
539 &                       comm_atom,mpi_atmtab ) ! optional (parallelism)
540 
541 
542 !This section has been created automatically by the script Abilint (TD).
543 !Do not modify the following lines by hand.
544 #undef ABI_FUNC
545 #define ABI_FUNC 'paw_dfptnl_accrhoij'
546 !End of the abilint section
547 
548  implicit none
549 
550 !Arguments ---------------------------------------------
551 !scalars
552  integer,intent(in) :: cplex,ipert1,ipert2,isppol,my_natom,natom,nspinor
553  integer,optional,intent(in) :: comm_atom
554  real(dp),intent(in) :: occ_k,wtk_k
555 !arrays
556  integer,intent(in) :: atindx(natom)
557  integer,optional,target,intent(in) :: mpi_atmtab(:)
558  type(pawcprj_type),intent(in) :: cwaveprj0_pert1(natom,nspinor),cwaveprj0_pert2(natom,nspinor)
559  type(pawcprj_type),intent(in) :: cwaveprj1_pert12(natom,nspinor),cwaveprj1_pert21(natom,nspinor)
560  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom)
561 
562 !Local variables ---------------------------------------
563 !scalars
564  integer :: cplex_rhoij,iatm,iatom,iatom1,ilmn,iplex,j0lmn,jlmn,klmn,klmn_im,klmn_re
565  integer :: my_comm_atom,ncpgr
566  logical :: compute_impart,compute_impart_cplex
567  logical :: my_atmtab_allocated,paral_atom
568  real(dp) :: ro11_im,ro11_re,weight
569  character(len=500) :: message
570 !arrays
571  integer,pointer :: my_atmtab(:)
572  real(dp) :: cpi0(2,nspinor),d1cpi0(2,nspinor),d2cpi0(2,nspinor)
573  real(dp) :: cpj0(2,nspinor),d1cpj0(2,nspinor),d2cpj0(2,nspinor)
574  real(dp) :: cpi1(2,nspinor),d2cpi1(2,nspinor)
575  real(dp) :: cpj1(2,nspinor),d2cpj1(2,nspinor)
576  real(dp) :: cpi2(2,nspinor),d1cpi2(2,nspinor)
577  real(dp) :: cpj2(2,nspinor),d1cpj2(2,nspinor)
578 
579 ! ***********************************************************************
580 
581  DBG_ENTER("COLL")
582 
583  if (my_natom==0) return
584 
585  ncpgr=1
586  if (ipert1<0.or.ipert1>natom+2.or.ipert2<0.or.ipert2>natom+2) then
587    message = 'Wrong inputs. Necessary conditions on ipert1 or ipert2 : 0 <= ipert <= natom+2'
588    MSG_BUG(message)
589  end if
590 
591 !Set up parallelism over atoms
592  paral_atom=(present(comm_atom).and.(my_natom/=natom))
593  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
594  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
595  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,&
596 & my_natom_ref=my_natom)
597 
598  weight=wtk_k*occ_k
599  if (pawrhoij(1)%nspden==2.and.pawrhoij(1)%nsppol==1.and.nspinor==1) weight=half*weight
600 
601 !!  ==================================================================
602 !!  === Accumulate (n,k) contribution to partial 2nd-order rhoij   ===
603 !!  ==================================================================
604 
605  compute_impart=(pawrhoij(1)%cplex==2)
606  compute_impart_cplex=((pawrhoij(1)%cplex==2).and.(cplex==2))
607 
608 ! NOT USED FOR PAWRHO21! => only for PAWRHO2 (full second derivative)
609 !!Accumulate :   < Psi^(pert1) | p_i^(0) > < p_j^(0) | Psi^(pert2) >
610 !!             + < Psi^(pert2) | p_i^(0) > < p_j^(0) | Psi^(pert1) >
611 ! if (nspinor==1) then
612 !   do iatom=1,my_natom
613 !     iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
614 !     iatm=atindx(iatom1)
615 !     cplex_rhoij=pawrhoij(iatom)%cplex
616 !     do jlmn=1,pawrhoij(iatom)%lmn_size
617 !       j0lmn=jlmn*(jlmn-1)/2
618 !       cpj1(1:2,1)=cwaveprj1_pert12(iatm,1)%cp(1:2,jlmn)   ! < p_j^(0) | Psi^(pert1) >
619 !       cpj2(1:2,1)=cwaveprj1_pert21(iatm,1)%cp(1:2,jlmn)   ! < p_j^(0) | Psi^(pert2) >
620 !       do ilmn=1,jlmn
621 !         klmn=j0lmn+ilmn
622 !         klmn_re=cplex_rhoij*(klmn-1)+1
623 !         cpi1(1:2,1)=cwaveprj1_pert12(iatm,1)%cp(1:2,ilmn) ! < p_i^(0) | Psi^(pert1) >
624 !         cpi2(1:2,1)=cwaveprj1_pert21(iatm,1)%cp(1:2,ilmn) ! < p_i^(0) | Psi^(pert2) >
625 !         ro11_re=zero
626 !         do iplex=1,cplex
627 !           ro11_re=ro11_re+cpi1(iplex,1)*cpj2(iplex,1)+cpj1(iplex,1)*cpi2(iplex,1)
628 !         end do
629 !         pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
630 !         if (compute_impart_cplex) then
631 !           klmn_im=klmn_re+1
632 !           ro11_im=        cpi1(1,1)*cpj2(2,1)-cpi1(2,1)*cpj2(1,1)
633 !           ro11_im=ro11_im+cpj1(1,1)*cpi2(2,1)-cpj1(2,1)*cpi2(1,1)
634 !           pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
635 !         end if
636 !       end do
637 !     end do
638 !   end do
639 ! else ! nspinor=2
640 !   MSG_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2")
641 ! end if
642 
643 !Accumulate :   < Psi^(pert1) | p_i^(pert2) > < p_j^(0)     | Psi^(0)     >
644 !             + < Psi^(pert1) | p_i^(0)     > < p_j^(pert2) | Psi^(0)     >
645 !             + < Psi^(0)     | p_i^(pert2) > < p_j^(0)     | Psi^(pert1) >
646 !             + < Psi^(0)     | p_i^(0)     > < p_j^(pert2) | Psi^(pert1) >
647  if (ipert2>0.and.ipert2<=natom) then
648    if (nspinor==1) then
649      do iatom=1,my_natom
650        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
651        iatm=atindx(iatom1)
652        if (iatom/=ipert2) cycle ! To move atom "ipert2" does not change projectors of other atoms
653        cplex_rhoij=pawrhoij(iatom)%cplex
654        do jlmn=1,pawrhoij(iatom)%lmn_size
655          j0lmn=jlmn*(jlmn-1)/2
656          cpj0(1:2,1)  =cwaveprj0_pert2 (iatm,1)% cp(1:2,  jlmn)   ! < p_j^(0)     | Psi^(0)     >
657          d2cpj0(1:2,1)=cwaveprj0_pert2 (iatm,1)%dcp(1:2,1,jlmn)   ! < p_j^(pert2) | Psi^(0)     >
658          cpj1(1:2,1)  =cwaveprj1_pert12(iatm,1)% cp(1:2,  jlmn)   ! < p_j^(0)     | Psi^(pert1) >
659          d2cpj1(1:2,1)=cwaveprj1_pert12(iatm,1)%dcp(1:2,1,jlmn)   ! < p_j^(pert2) | Psi^(pert1) >
660          do ilmn=1,jlmn
661            klmn=j0lmn+ilmn
662            klmn_re=cplex_rhoij*(klmn-1)+1
663            cpi0(1:2,1)  =cwaveprj0_pert2 (iatm,1)% cp(1:2,  ilmn) ! < p_i^(0)     | Psi^(0)     >
664            d2cpi0(1:2,1)=cwaveprj0_pert2 (iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(0)     >
665            cpi1(1:2,1)  =cwaveprj1_pert12(iatm,1)% cp(1:2,  ilmn) ! < p_i^(0)     | Psi^(pert1) >
666            d2cpi1(1:2,1)=cwaveprj1_pert12(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(pert1) >
667            ro11_re=zero
668            do iplex=1,cplex
669              ro11_re=ro11_re+d2cpi1(iplex,1)*  cpj0(iplex,1)
670              ro11_re=ro11_re+  cpi1(iplex,1)*d2cpj0(iplex,1)
671              ro11_re=ro11_re+d2cpi0(iplex,1)*  cpj1(iplex,1)
672              ro11_re=ro11_re+  cpi0(iplex,1)*d2cpj1(iplex,1)
673            end do
674            pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
675            if (compute_impart_cplex) then
676              klmn_im=klmn_re+1
677              ro11_im=        d2cpi1(1,1)*  cpj0(2,1)-d2cpi1(2,1)*  cpj0(1,1)
678              ro11_im=ro11_im+  cpi1(1,1)*d2cpj0(2,1)-  cpi1(2,1)*d2cpj0(1,1)
679              ro11_im=ro11_im+d2cpi0(1,1)*  cpj1(2,1)-d2cpi0(2,1)*  cpj1(1,1)
680              ro11_im=ro11_im+  cpi0(1,1)*d2cpj1(2,1)-  cpi0(2,1)*d2cpj1(1,1)
681              pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
682            end if
683          end do
684        end do
685      end do
686    else ! nspinor=2
687      MSG_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2")
688    end if
689  end if
690 
691 !Accumulate : < Psi^(pert2) | p_i^(pert1) > < p_j^(0)     | Psi^(0)     >
692 !           + < Psi^(pert2) | p_i^(0)     > < p_j^(pert1) | Psi^(0)     >
693 !           + < Psi^(0)     | p_i^(pert1) > < p_j^(0)     | Psi^(pert2) >
694 !           + < Psi^(0)     | p_i^(0) >     < p_j^(pert1) | Psi^(pert2) >
695  if (ipert1>0.and.ipert1<=natom) then
696    if (nspinor==1) then
697      do iatom=1,my_natom
698        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
699        iatm=atindx(iatom1)
700        cplex_rhoij=pawrhoij(iatom)%cplex
701        if (iatom/=ipert1) cycle ! To move atom "ipert1" does not change projectors of other atoms
702        do jlmn=1,pawrhoij(iatom)%lmn_size
703          j0lmn=jlmn*(jlmn-1)/2
704          cpj0(1:2,1)  =cwaveprj0_pert1 (iatm,1)% cp(1:2,  jlmn)   ! < p_j^(0)     | Psi^(0)     >
705          d1cpj0(1:2,1)=cwaveprj0_pert1 (iatm,1)%dcp(1:2,1,jlmn)   ! < p_j^(pert1) | Psi^(0)     >
706          cpj2(1:2,1)  =cwaveprj1_pert21(iatm,1)% cp(1:2,  jlmn)   ! < p_j^(0)     | Psi^(pert2) >
707          d1cpj2(1:2,1)=cwaveprj1_pert21(iatm,1)%dcp(1:2,1,jlmn)   ! < p_j^(pert1) | Psi^(pert2) >
708          do ilmn=1,jlmn
709            klmn=j0lmn+ilmn
710            klmn_re=cplex_rhoij*(klmn-1)+1
711            cpi0(1:2,1)  =cwaveprj0_pert1 (iatm,1)% cp(1:2,  ilmn) ! < p_i^(0)     | Psi^(0)     >
712            d1cpi0(1:2,1)=cwaveprj0_pert1 (iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(0)     >
713            cpi2(1:2,1)  =cwaveprj1_pert21(iatm,1)% cp(1:2,  ilmn) ! < p_i^(0)     | Psi^(pert2) >
714            d1cpi2(1:2,1)=cwaveprj1_pert21(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(pert2) >
715            ro11_re=zero
716            do iplex=1,cplex
717              ro11_re=ro11_re+d1cpi2(iplex,1)*  cpj0(iplex,1)
718              ro11_re=ro11_re+  cpi2(iplex,1)*d1cpj0(iplex,1)
719              ro11_re=ro11_re+d1cpi0(iplex,1)*  cpj2(iplex,1)
720              ro11_re=ro11_re+  cpi0(iplex,1)*d1cpj2(iplex,1)
721            end do
722            pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
723            if (compute_impart_cplex) then
724              klmn_im=klmn_re+1
725              ro11_im=        d1cpi2(1,1)*  cpj0(2,1)-d1cpi2(2,1)*  cpj0(1,1)
726              ro11_im=ro11_im+  cpi2(1,1)*d1cpj0(2,1)-  cpi2(2,1)*d1cpj0(1,1)
727              ro11_im=ro11_im+d1cpi0(1,1)*  cpj2(2,1)-d1cpi0(2,1)*  cpj2(1,1)
728              ro11_im=ro11_im+  cpi0(1,1)*d1cpj2(2,1)-  cpi0(2,1)*d1cpj2(1,1)
729              pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
730            end if
731          end do
732        end do
733      end do
734    else ! nspinor=2
735      MSG_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2")
736    end if
737  end if
738 !  End
739 
740 !Accumulate :   < Psi^(0) | p_i^(pert1) > < p_j^(pert2) | Psi^(0) >
741 !             + < Psi^(0) | p_i^(pert2) > < p_j^(pert1) | Psi^(0) >
742  if (ipert1>0.and.ipert1<=natom.and.ipert2>0.and.ipert2<=natom) then
743    if (nspinor==1) then
744      do iatom=1,my_natom
745        iatom1=iatom;if (paral_atom) iatom1=my_atmtab(iatom)
746        iatm=atindx(iatom1)
747        if (iatom/=ipert1.or.iatom/=ipert2) cycle ! To move atom "ipert" does not change projectors of other atoms
748        cplex_rhoij=pawrhoij(iatom)%cplex
749        do jlmn=1,pawrhoij(iatom)%lmn_size
750          j0lmn=jlmn*(jlmn-1)/2
751          d1cpj0(1:2,1)=cwaveprj0_pert1(iatm,1)%dcp(1:2,1,jlmn)   ! < p_j^(pert1) | Psi^(0) >
752          d2cpj0(1:2,1)=cwaveprj0_pert2(iatm,1)%dcp(1:2,1,jlmn)   ! < p_j^(pert2) | Psi^(0) >
753          do ilmn=1,jlmn
754            klmn=j0lmn+ilmn
755            klmn_re=cplex_rhoij*(klmn-1)+1
756            d1cpi0(1:2,1)=cwaveprj0_pert1(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert1) | Psi^(0) >
757            d2cpi0(1:2,1)=cwaveprj0_pert2(iatm,1)%dcp(1:2,1,ilmn) ! < p_i^(pert2) | Psi^(0) >
758            ro11_re=zero
759            do iplex=1,cplex
760              ro11_re=ro11_re+d1cpi0(iplex,1)*d2cpj0(iplex,1)+d2cpi0(iplex,1)*d1cpj0(iplex,1)
761            end do
762            pawrhoij(iatom)%rhoij_(klmn_re,isppol)=pawrhoij(iatom)%rhoij_(klmn_re,isppol)+weight*ro11_re
763            if (compute_impart_cplex) then
764              klmn_im=klmn_re+1
765              ro11_im=        d1cpi0(1,1)*d2cpj0(2,1)-d1cpi0(2,1)*d2cpj0(1,1)
766              ro11_im=ro11_im+d2cpi0(1,1)*d1cpj0(2,1)-d2cpi0(2,1)*d1cpj0(1,1)
767              pawrhoij(iatom)%rhoij_(klmn_im,isppol)=pawrhoij(iatom)%rhoij_(klmn_im,isppol)+weight*ro11_im
768            end if
769          end do
770        end do
771      end do
772    else ! nspinor=2
773      MSG_BUG("paw_dfptnl_accrhoij is not implemented for nspinor=2")
774    end if
775  end if
776 
777 !Destroy atom table used for parallelism
778  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
779 
780  DBG_EXIT("COLL")
781 
782 end subroutine paw_dfptnl_accrhoij

m_paw_dfptnl/paw_dfptnl_energy [ Functions ]

[ Top ] [ m_paw_dfptnl ] [ Functions ]

NAME

 paw_dfptnl_energy

FUNCTION

 Compute the XC PAW on-site contributions to a 3rd-order energy.
 It is equal to:
    E_onsite= \sum_at [ E_at(kxc,rho1,rho2,rho3) - E_at(tkxc,trho1,trho2,trho3) ]
 where E_at(...) is computed in paw_dfptnl_xc.F90.
 The atomic densities are computed from pawrhoij1,pawrhoij2 and pawrhoij3.
 This routine is similar to pawdfptenergy.F90 but is implemented independently
 in order to not overload the original routine.
 LDA ONLY - USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)

INPUTS

  ixc= choice of exchange-correlation scheme
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  ntypat=number of types of atoms in unit cell.
  paw_an0(natom) <type(paw_an_type)>=paw arrays for 0th-order quantities given on angular mesh
  paw_an1(natom) <type(paw_an_type)>=paw arrays for 1st-order quantities given on angular mesh
                                     This corresponds to (j1) perturbation
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawprtvol=control print volume and debugging output for PAW
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij_1-2-3(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)

OUTPUT

  d3exc= real and imaginary parts of the contribution to the third derivative of the total energy

SIDE EFFECTS

PARENTS

      dfptnl_pert

CHILDREN

      free_my_atmtab,get_my_atmtab,pawdensities,pawxc_dfpt,pawxcm3
      timab,xmpi_sum

SOURCE

100 subroutine paw_dfptnl_energy(d3exc,ixc,my_natom,natom,ntypat,&
101 &                    paw_an0,pawang,pawprtvol,pawrad,&
102 &                    pawrhoij_1,pawrhoij_2,pawrhoij_3,&
103 &                    pawtab,pawxcdev,&
104 &                    mpi_atmtab,comm_atom) ! optional arguments (parallelism)
105 
106 
107 !This section has been created automatically by the script Abilint (TD).
108 !Do not modify the following lines by hand.
109 #undef ABI_FUNC
110 #define ABI_FUNC 'paw_dfptnl_energy'
111 !End of the abilint section
112 
113  implicit none
114 
115 !Arguments ---------------------------------------------
116 !scalars
117  integer,intent(in) :: ixc,my_natom,natom,ntypat
118  integer,intent(in) :: pawprtvol,pawxcdev
119  integer,optional,intent(in) :: comm_atom
120  type(pawang_type),intent(in) :: pawang
121 !arrays
122  integer,optional,target,intent(in) :: mpi_atmtab(:)
123  real(dp),intent(out) :: d3exc(2)
124  type(paw_an_type),intent(in) :: paw_an0(my_natom)
125  type(pawrad_type),intent(in) :: pawrad(ntypat)
126  type(pawrhoij_type),intent(in) :: pawrhoij_1(natom)
127  type(pawrhoij_type),intent(in) :: pawrhoij_2(natom)
128  type(pawrhoij_type),intent(in) :: pawrhoij_3(natom)
129  type(pawtab_type),intent(in) :: pawtab(ntypat)
130 
131 !Local variables ---------------------------------------
132 !scalars
133  integer :: cplex_1,cplex_2,cplex_3,iatom,iatom_tot,itypat
134  integer :: lm_size_all,mesh_size,my_comm_atom,npts,nspden,nzlmopt
135  integer :: opt_compch,usecore,usetcore,usexcnhat
136  logical :: my_atmtab_allocated,paral_atom
137  real(dp) :: compch,d3exc1_iat(2)
138 ! character(len=500) :: msg
139 !arrays
140  integer,pointer :: my_atmtab(:)
141  logical,allocatable :: lmselect_1(:),lmselect_2(:),lmselect_3(:),lmselect_tmp(:)
142  real(dp),allocatable :: nhat1_1(:,:,:),rho1_1(:,:,:),trho1_1(:,:,:)
143  real(dp),allocatable :: nhat1_2(:,:,:),rho1_2(:,:,:),trho1_2(:,:,:)
144  real(dp),allocatable :: nhat1_3(:,:,:),rho1_3(:,:,:),trho1_3(:,:,:)
145 
146 ! *************************************************************************
147 
148  DBG_ENTER("COLL")
149 
150  nzlmopt = 0 ! compute all LM-moments of the density and use all LM-moments
151 
152  if (pawxcdev/=0) then
153    MSG_BUG("paw_dfptnl_energy is not implemented for pawxcdev /=0")
154  end if
155 
156 !Set up parallelism over atoms
157  paral_atom=(present(comm_atom).and.(my_natom/=natom))
158  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
159  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
160  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
161 
162 !!Various inits
163  opt_compch=0; !optvxc=1;optexc=3
164  usecore=0;usetcore=0  ! This is true for phonons and Efield pert.
165  usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat)
166 
167  npts=pawang%angl_size
168 
169  d3exc = zero
170 
171 !================ Loop on atomic sites =======================
172  do iatom=1,my_natom
173    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
174 
175    itypat=pawrhoij_1(iatom)%itypat
176    mesh_size=pawtab(itypat)%mesh_size
177    nspden=pawrhoij_1(iatom)%nspden
178    cplex_1=pawrhoij_1(iatom)%cplex
179    cplex_2=pawrhoij_2(iatom)%cplex
180    cplex_3=pawrhoij_3(iatom)%cplex
181    lm_size_all=paw_an0(iatom)%lm_size
182 
183    ABI_ALLOCATE(lmselect_tmp,(lm_size_all))
184    lmselect_tmp(:)=.true.
185 
186 !  Compute on-site 1st-order densities (pert1)
187    ABI_ALLOCATE(lmselect_1,(lm_size_all))
188    lmselect_1(:)=paw_an0(iatom)%lmselect(:)
189    ABI_ALLOCATE(rho1_1,(cplex_1*mesh_size,lm_size_all,nspden))
190    ABI_ALLOCATE(trho1_1,(cplex_1*mesh_size,lm_size_all,nspden))
191    ABI_ALLOCATE(nhat1_1,(cplex_1*mesh_size,lm_size_all,nspden*usexcnhat))
192    call pawdensities(compch,cplex_1,iatom_tot,lmselect_tmp,lmselect_1,&
193 &   lm_size_all,nhat1_1,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
194 &   pawrad(itypat),pawrhoij_1(iatom),pawtab(itypat),rho1_1,trho1_1)
195 !  Compute on-site 1st-order densities (pert2)
196    ABI_ALLOCATE(lmselect_2,(lm_size_all))
197    lmselect_2(:)=paw_an0(iatom)%lmselect(:)
198    ABI_ALLOCATE(rho1_2,(cplex_2*mesh_size,lm_size_all,nspden))
199    ABI_ALLOCATE(trho1_2,(cplex_2*mesh_size,lm_size_all,nspden))
200    ABI_ALLOCATE(nhat1_2,(cplex_2*mesh_size,lm_size_all,nspden*usexcnhat))
201    call pawdensities(compch,cplex_2,iatom_tot,lmselect_tmp,lmselect_2,&
202 &   lm_size_all,nhat1_2,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
203 &   pawrad(itypat),pawrhoij_2(iatom),pawtab(itypat),rho1_2,trho1_2)
204 !  Compute on-site 1st-order densities (pert3)
205    ABI_ALLOCATE(lmselect_3,(lm_size_all))
206    lmselect_3(:)=paw_an0(iatom)%lmselect(:)
207    ABI_ALLOCATE(rho1_3,(cplex_3*mesh_size,lm_size_all,nspden))
208    ABI_ALLOCATE(trho1_3,(cplex_3*mesh_size,lm_size_all,nspden))
209    ABI_ALLOCATE(nhat1_3,(cplex_3*mesh_size,lm_size_all,nspden*usexcnhat))
210    call pawdensities(compch,cplex_3,iatom_tot,lmselect_tmp,lmselect_3,&
211 &   lm_size_all,nhat1_3,nspden,nzlmopt,opt_compch,1-usexcnhat,-1,0,pawang,pawprtvol,&
212 &   pawrad(itypat),pawrhoij_3(iatom),pawtab(itypat),rho1_3,trho1_3)
213    ABI_DEALLOCATE(lmselect_tmp)
214 
215    call paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,paw_an0(iatom)%k3xc1,lm_size_all,&
216 &   lmselect_1,lmselect_2,lmselect_3,nhat1_1,nhat1_2,nhat1_3,&
217 &   paw_an0(iatom)%nk3xc1,mesh_size,nspden,pawang,pawrad(itypat),&
218 &   rho1_1,rho1_2,rho1_3,0)
219    d3exc = d3exc + d3exc1_iat
220 
221    call paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,paw_an0(iatom)%k3xct1,lm_size_all,&
222 &   lmselect_1,lmselect_2,lmselect_3,nhat1_1,nhat1_2,nhat1_3,&
223 &   paw_an0(iatom)%nk3xc1,mesh_size,nspden,pawang,pawrad(itypat),&
224 &   trho1_1,trho1_2,trho1_3,usexcnhat)
225    d3exc = d3exc - d3exc1_iat
226 
227    ABI_DEALLOCATE(lmselect_1)
228    ABI_DEALLOCATE(lmselect_2)
229    ABI_DEALLOCATE(lmselect_3)
230    ABI_DEALLOCATE(nhat1_1)
231    ABI_DEALLOCATE(nhat1_2)
232    ABI_DEALLOCATE(nhat1_3)
233    ABI_DEALLOCATE(rho1_1)
234    ABI_DEALLOCATE(rho1_2)
235    ABI_DEALLOCATE(rho1_3)
236    ABI_DEALLOCATE(trho1_1)
237    ABI_DEALLOCATE(trho1_2)
238    ABI_DEALLOCATE(trho1_3)
239 
240 !  ================ End loop oon atomic sites =======================
241  end do
242 
243 !!Reduction in case of parallelism
244 ! if (paral_atom) then
245 !   call xmpi_sum(delta_energy,my_comm_atom,ierr)
246 ! end if
247 
248 !!Destroy atom table used for parallelism
249 ! call free_my_atmtab(my_atmtab,my_atmtab_allocated)
250 
251 ! call timab(567,2,tsec)
252 
253  DBG_EXIT("COLL")
254 
255 end subroutine paw_dfptnl_energy

m_paw_dfptnl/paw_dfptnl_xc [ Functions ]

[ Top ] [ m_paw_dfptnl ] [ Functions ]

NAME

 paw_dfptnl_xc

FUNCTION

 Compute a contribution of the third derivative of XC energy of ONE PAW sphere Om_a.
 It is equal to:
   E_at(kxc,rho1,rho2,rho3) = Int_{Om_a} dr kxc(r) * rho1(r) * rho2(r) * rho3(r)
 where kxc,rho1,rho2 and rho3 are inputs.
 This routine is similar to m_pawxc.F90:pawxc_dfpt(...) but is implemented independently
 in order to not overload the original routine.
 LDA ONLY - USE THE DENSITY OVER A WHOLE SPHERICAL GRID (r,theta,phi)

INPUTS

  cplex_1-2-3= if 1, 1st-order densities are REAL, if 2, COMPLEX
  d3exc1_iat=third-order derivative to compute
  ixc= choice of exchange-correlation scheme
  kxc(nrad,pawang%angl_size,nkxc)=GS xc kernel
  lm_size=size of density array rhor (see below)
  lmselect1-2-3(lm_size)=select the non-zero LM-moments of input density rhor1-2-3
  nhat1-2-3(cplex_den*nrad,lm_size,nspden)=first-order change of compensation density
                                        (total in 1st half and spin-up in 2nd half if nspden=2)
  nkxc=second dimension of the kxc array
  nrad=size of radial mesh for densities/potentials (might be different from pawrad%mesh_size)
  nspden=number of spin-density components
  option=0  compute both 2nd-order XC energy and 1st-order potential
         1  compute only 1st-order XC potential
         2  compute only 2nd-order XC energy, XC potential is temporary computed here
         3  compute only 2nd-order XC energy, XC potential is input in vxc1(:)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad <type(pawrad_type)>=paw radial mesh and related data
  rhor1-2-3(cplex_den*nrad,lm_size,nspden)=first-order change of density
  usexcnhat= 0 if compensation density does not have to be used
             1 if compensation density has to be used in d2Exc only

OUTPUT

  d3exc1_iat = E_at(kxc,rho1,rho2,rho3) (see FUNCTION above)

SIDE EFFECTS

PARENTS

      paw_dfptnl_energy

CHILDREN

SOURCE

307 subroutine paw_dfptnl_xc(cplex_1,cplex_2,cplex_3,d3exc1_iat,ixc,kxc,lm_size,lmselect1,lmselect2,lmselect3,&
308 &                 nhat1,nhat2,nhat3,nkxc,nrad,nspden,pawang,pawrad,rhor1,rhor2,rhor3,usexcnhat)
309 
310 
311 !This section has been created automatically by the script Abilint (TD).
312 !Do not modify the following lines by hand.
313 #undef ABI_FUNC
314 #define ABI_FUNC 'paw_dfptnl_xc'
315 !End of the abilint section
316 
317  implicit none
318 
319 !Arguments ------------------------------------
320 !scalars
321  integer,intent(in) :: cplex_1,cplex_2,cplex_3,ixc,lm_size,nkxc,nrad,nspden,usexcnhat
322  type(pawang_type),intent(in) :: pawang
323  type(pawrad_type),intent(in) :: pawrad
324 !arrays
325  logical,intent(in) :: lmselect1(lm_size),lmselect2(lm_size),lmselect3(lm_size)
326  real(dp),intent(out) :: d3exc1_iat(2)
327  real(dp),intent(in) :: kxc(nrad,pawang%angl_size,nkxc)
328  real(dp),intent(in) :: nhat1(cplex_1*nrad,lm_size,nspden*((usexcnhat+1)/2))
329  real(dp),intent(in) :: nhat2(cplex_2*nrad,lm_size,nspden*((usexcnhat+1)/2))
330  real(dp),intent(in) :: nhat3(cplex_3*nrad,lm_size,nspden*((usexcnhat+1)/2))
331  real(dp),intent(in) :: rhor1(cplex_1*nrad,lm_size,nspden)
332  real(dp),intent(in) :: rhor2(cplex_2*nrad,lm_size,nspden)
333  real(dp),intent(in) :: rhor3(cplex_3*nrad,lm_size,nspden)
334 
335 !Local variables-------------------------------
336 !scalars
337  integer :: ii,ilm,ipts,ispden,lm_size_eff,npts
338  real(dp) :: d3exc1_int,rho1u,rho1d,rho2u,rho2d,rho3u,rho3d
339  character(len=500) :: msg
340 !arrays
341 ! real(dp) :: tsec(2)
342  real(dp),allocatable :: ff(:),rho1arr(:,:),rho2arr(:,:),rho3arr(:,:)
343 
344 ! *************************************************************************
345 
346 !----------------------------------------------------------------------
347 !----- Initializations
348 !----------------------------------------------------------------------
349 
350  npts=pawang%angl_size
351  lm_size_eff=min(lm_size,pawang%ylm_size)
352 
353  d3exc1_iat(:) = zero
354 
355 !Special case: no XC applied
356  if (ixc==0) then
357    msg='Note that no xc is applied (ixc=0). Returning'
358    MSG_WARNING(msg)
359    return
360  end if
361 
362  ABI_ALLOCATE(rho1arr,(cplex_1*nrad,nspden))
363  ABI_ALLOCATE(rho2arr,(cplex_2*nrad,nspden))
364  ABI_ALLOCATE(rho3arr,(cplex_3*nrad,nspden))
365 
366 !Restriction : all cplex must be 1
367  if (cplex_1/=1.or.cplex_2/=1.or.cplex_3/=1) then
368    msg='All cplex must be one (for the moment...)'
369    MSG_BUG(msg)
370  end if
371 !Restriction : nspden must be 1
372 ! if (nkxc>1) then
373 !   msg='nkxc must be one (<=> nspden=1) (for the moment...)'
374 !   MSG_BUG(msg)
375 ! end if
376 
377  ABI_ALLOCATE(ff,(nrad))
378 
379 !!----------------------------------------------------------------------
380 !!----- Loop on the angular part and inits
381 !!----------------------------------------------------------------------
382 
383 !Do loop on the angular part (theta,phi)
384  do ipts=1,npts
385 
386 !  Copy the input 1st-order density for this (theta,phi) - PERT1
387    rho1arr(:,:)=zero
388    if (usexcnhat==0) then
389      do ispden=1,nspden
390        do ilm=1,lm_size_eff
391          if (lmselect1(ilm)) rho1arr(:,ispden)=rho1arr(:,ispden) &
392 &         +rhor1(:,ilm,ispden)*pawang%ylmr(ilm,ipts)
393        end do
394      end do
395    else
396      do ispden=1,nspden
397        do ilm=1,lm_size_eff
398          if (lmselect1(ilm)) rho1arr(:,ispden)=rho1arr(:,ispden) &
399 &         +(rhor1(:,ilm,ispden)+nhat1(:,ilm,ispden))*pawang%ylmr(ilm,ipts)
400        end do
401      end do
402    end if
403 
404 !  Copy the input 1st-order density for this (theta,phi) - PERT2
405    rho2arr(:,:)=zero
406    if (usexcnhat==0) then
407      do ispden=1,nspden
408        do ilm=1,lm_size_eff
409          if (lmselect2(ilm)) rho2arr(:,ispden)=rho2arr(:,ispden) &
410 &         +rhor2(:,ilm,ispden)*pawang%ylmr(ilm,ipts)
411        end do
412      end do
413    else
414      do ispden=1,nspden
415        do ilm=1,lm_size_eff
416          if (lmselect2(ilm)) rho2arr(:,ispden)=rho2arr(:,ispden) &
417 &         +(rhor2(:,ilm,ispden)+nhat2(:,ilm,ispden))*pawang%ylmr(ilm,ipts)
418        end do
419      end do
420    end if
421 
422 !  Copy the input 1st-order density for this (theta,phi) - PERT3
423    rho3arr(:,:)=zero
424    if (usexcnhat==0) then
425      do ispden=1,nspden
426        do ilm=1,lm_size_eff
427          if (lmselect3(ilm)) rho3arr(:,ispden)=rho3arr(:,ispden) &
428 &         +rhor3(:,ilm,ispden)*pawang%ylmr(ilm,ipts)
429        end do
430      end do
431    else
432      do ispden=1,nspden
433        do ilm=1,lm_size_eff
434          if (lmselect3(ilm)) rho3arr(:,ispden)=rho3arr(:,ispden) &
435 &         +(rhor3(:,ilm,ispden)+nhat3(:,ilm,ispden))*pawang%ylmr(ilm,ipts)
436        end do
437      end do
438    end if
439 
440 !  ----------------------------------------------------------------------
441 !  ----- Accumulate and store 3nd-order change of XC energy
442 !  ----------------------------------------------------------------------
443 
444    if (cplex_1==1.and.cplex_2==1.and.cplex_3==1) then ! all cplex are 1 :
445      if (nspden==1) then
446        ff(:)=kxc(:,ipts,1)*rho1arr(:,1)*rho2arr(:,1)*rho3arr(:,1)
447      else if (nspden==2) then
448        do ii=1,nrad
449          rho1u=rho1arr(ii,2)
450          rho1d=rho1arr(ii,1)-rho1arr(ii,2)
451          rho2u=rho2arr(ii,2)
452          rho2d=rho2arr(ii,1)-rho2arr(ii,2)
453          rho3u=rho3arr(ii,2)
454          rho3d=rho3arr(ii,1)-rho3arr(ii,2)
455          ff(ii)=&
456 !          uuu                                uud
457 &         kxc(ii,ipts,1)*rho1u*rho2u*rho3u + kxc(ii,ipts,2)*rho1u*rho2u*rho3d + &
458 !          udu                                udd
459 &         kxc(ii,ipts,2)*rho1u*rho2d*rho3u + kxc(ii,ipts,3)*rho1u*rho2d*rho3d + &
460 !          duu                                dud
461 &         kxc(ii,ipts,2)*rho1d*rho2u*rho3u + kxc(ii,ipts,3)*rho1d*rho2u*rho3d + &
462 !          ddu                                ddd
463 &         kxc(ii,ipts,3)*rho1d*rho2d*rho3u + kxc(ii,ipts,4)*rho1d*rho2d*rho3d
464        end do
465      end if
466    end if
467 
468    ff(1:nrad)=ff(1:nrad)*pawrad%rad(1:nrad)**2
469    call simp_gen(d3exc1_int,ff,pawrad)
470    d3exc1_iat(1)=d3exc1_iat(1)+d3exc1_int*pawang%angwgth(ipts)
471 
472 !  ----- End of the loop on npts (angular part)
473  end do
474 
475  d3exc1_iat = d3exc1_iat*four_pi
476 
477  ABI_DEALLOCATE(ff)
478  ABI_DEALLOCATE(rho1arr)
479  ABI_DEALLOCATE(rho2arr)
480  ABI_DEALLOCATE(rho3arr)
481 
482 end subroutine paw_dfptnl_xc