TABLE OF CONTENTS


ABINIT/newfermie1 [ Functions ]

[ Top ] [ Functions ]

NAME

 newfermie1

FUNCTION

 This routine computes the derivative of the fermi energy wrt
 the active perturbation for use in evaluating the edocc term
 and active subspace contribution to the first-order wavefunctions
 in the case of metals. This is presently used only for the
 strain and magnetic field perturbations, and only for Q = 0.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DRH)
 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

  cplex: if 1, real space 1-order functions on FFT grid are REAL,
    if 2, COMPLEX
  fe1fixed=fixed contribution to the first-order Fermi energy
  ipert=index of perturbation
  istep=index of the number of steps in the routine scfcv
  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=number of atoms
  nfft=(effective) number of FFT grid points (for this processor)
  nfftot= total number of FFT grid points
  nhatfermi(nfft,nspden)=fermi-level compensation charge density (PAW only)
  nspden=number of spin-density components
  ntypat=number of atom types
  occopt=option for occupancies
  paw_an(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
  paw_ij1(natom) <type(paw_ij_type)>=(1st-order) paw arrays given on (i,j) channels
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawnzlm=-- PAW only -- option for the computation of non-zero
          lm moments of the on-sites densities
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij1(natom) <type(pawrhoij_type)>= paw rhoij 1st-order occupancies
  pawrhoijfermi(natom) <type(pawrhoij_type)>=paw rhoij occupancies at Fermi level
  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)
  xclevel= XC functional level
  prtvol=control print volume and debugging output
  rhorfermi(nfft,nspden)=fermi-level electronic density
  ucvol=unit cell volume in bohr**3
  usepaw=1 if PAW is activated
  usexcnhat= -PAW only- flag controling use of compensation density in Vxc
  vtrial1(cplex*nfft,nspden)=1-st order potential
  vxc1(cplex*nfft,nspden)=1-st order XC potential

OUTPUT

  (see side effects)

SIDE EFFECTS

  fermie1=derivative of fermi energy wrt perturbation
   at input  : old value
   at output : updated value

PARENTS

      dfpt_scfcv

CHILDREN

      dotprod_vn,free_my_atmtab,get_my_atmtab,pawdfptenergy,wrtout

SOURCE

 74 #if defined HAVE_CONFIG_H
 75 #include "config.h"
 76 #endif
 77 
 78 #include "abi_common.h"
 79 
 80 
 81 subroutine newfermie1(cplex,fermie1,fe1fixed,ipert,istep,ixc,my_natom,natom,nfft,nfftot,&
 82 &                     nhatfermi,nspden,ntypat,occopt,paw_an,paw_an1,paw_ij1,pawang,pawnzlm,pawrad,&
 83 &                     pawrhoij1,pawrhoijfermi,pawtab,pawxcdev,prtvol,rhorfermi,&
 84 &                     ucvol,usepaw,usexcnhat,vtrial1,vxc1,xclevel,&
 85 &                     mpi_atmtab,comm_atom) ! optional arguments (parallelism)
 86 
 87  use defs_basis
 88  use m_profiling_abi
 89  use m_errors
 90  use m_xmpi
 91 
 92  use m_pawang,     only : pawang_type
 93  use m_pawrad,     only : pawrad_type
 94  use m_pawtab,     only : pawtab_type
 95  use m_paw_an,     only : paw_an_type
 96  use m_paw_ij,     only : paw_ij_type
 97  use m_pawrhoij,   only : pawrhoij_type
 98  use m_paral_atom, only : get_my_atmtab, free_my_atmtab
 99  use m_cgtools,    only : dotprod_vn
100 
101 !This section has been created automatically by the script Abilint (TD).
102 !Do not modify the following lines by hand.
103 #undef ABI_FUNC
104 #define ABI_FUNC 'newfermie1'
105  use interfaces_14_hidewrite
106  use interfaces_65_paw
107 !End of the abilint section
108 
109  implicit none
110 
111 !Arguments -------------------------------
112 !scalars
113  integer,intent(in) :: cplex,ipert,istep,ixc,my_natom,natom,nfft,nfftot,nspden,ntypat
114  integer,intent(in) :: occopt,pawnzlm,pawxcdev,prtvol,usepaw,usexcnhat,xclevel
115  integer,optional,intent(in) :: comm_atom
116  real(dp),intent(in) :: fe1fixed,ucvol
117  real(dp),intent(inout) :: fermie1
118  type(pawang_type),intent(in) :: pawang
119 !arrays
120  integer,optional,target,intent(in) :: mpi_atmtab(:)
121  real(dp),intent(in) :: rhorfermi(nfft,nspden),vtrial1(cplex*nfft,nspden)
122  real(dp),intent(in) :: nhatfermi(:,:),vxc1(:,:)
123  type(paw_an_type),intent(in) :: paw_an(my_natom*usepaw)
124  type(paw_an_type),intent(inout) :: paw_an1(my_natom*usepaw)
125  type(paw_ij_type),intent(inout) :: paw_ij1(my_natom*usepaw)
126  type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw)
127  type(pawrhoij_type),intent(in) :: pawrhoij1(my_natom*usepaw),pawrhoijfermi(my_natom*usepaw)
128  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
129 
130 !Local variables-------------------------------
131 !scalars
132  integer :: ipert0,my_comm_atom,nzlmopt,nzlmopt_fermi,option,pawprtvol
133  logical :: my_atmtab_allocated,paral_atom
134  real(dp) :: doti,fe1_scf,fe1_tmp,fermie1_new,fermie1rs
135  character(len=500) :: msg
136 !arrays
137  integer, pointer :: my_atmtab(:)
138  real(dp) :: fe1_paw(2)
139  real(dp), allocatable :: rhor_nonhat(:,:),vtrial1_novxc(:,:)
140 
141 ! *********************************************************************
142 
143 !Tests
144  if (cplex==2) then
145    msg='Not compatible with cplex=2!'
146    MSG_BUG(msg)
147  end if
148  if (usepaw==1.and.usexcnhat==0.and.(size(nhatfermi)<=0.or.size(vxc1)<=0)) then
149    msg='Should have nhatfermi and vxc1 allocated with usexcnhat=0!'
150    MSG_BUG(msg)
151  end if
152 
153 !Set up parallelism over atoms
154  paral_atom=(present(comm_atom).and.(my_natom/=natom))
155  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
156  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
157  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
158 
159  if(occopt>=3 .and. occopt <=8) then
160 
161 !  The product of the current trial potential and the so-called Fermi level
162 !  density is integrated to give the local potential contributions to the
163 !  first-order Fermi level.
164    option=1
165    if (usepaw==1.and.usexcnhat==0) then
166      ABI_ALLOCATE(rhor_nonhat,(nfft,nspden))
167      ABI_ALLOCATE(vtrial1_novxc,(nfft,nspden))
168      rhor_nonhat(1:nfft,1:nspden)=rhorfermi(1:nfft,1:nspden)-nhatfermi(1:nfft,1:nspden)
169      vtrial1_novxc(1:nfft,1:nspden)=vtrial1(1:nfft,1:nspden)-vxc1(1:nfft,1:nspden)
170      call dotprod_vn(cplex,rhor_nonhat,fe1_scf,doti,nfft,nfftot,&
171 &     nspden,option,vtrial1,ucvol)
172      call dotprod_vn(cplex,nhatfermi,fe1_tmp,doti,nfft,nfftot,&
173 &     nspden,option,vtrial1_novxc,ucvol)
174      fe1_scf=fe1_scf+fe1_tmp
175      ABI_DEALLOCATE(rhor_nonhat)
176      ABI_DEALLOCATE(vtrial1_novxc)
177    else
178      call dotprod_vn(cplex,rhorfermi,fe1_scf,doti,nfft,nfftot,&
179 &     nspden,option,vtrial1,ucvol)
180    end if
181 
182    fe1_paw(:)=zero
183 !  PAW on-site contribution (use Fermi level occupation matrix)
184    if (usepaw==1) then
185      ipert0=0;pawprtvol=0
186      nzlmopt=0;if (istep>1) nzlmopt=pawnzlm
187      if (istep==1.and.pawnzlm>0) nzlmopt=-1
188      nzlmopt_fermi=0;if (pawnzlm>0) nzlmopt_fermi=-1
189      call pawdfptenergy(fe1_paw,ipert,ipert0,ixc,my_natom,natom,ntypat,nzlmopt,&
190 &     nzlmopt_fermi,paw_an,paw_an1,paw_ij1,pawang,pawprtvol,pawrad,&
191 &     pawrhoij1,pawrhoijfermi,pawtab,pawxcdev,xclevel,&
192 &     mpi_atmtab=my_atmtab, comm_atom=my_comm_atom)
193    end if
194 
195 !  The fixed contributions consisting of non-local potential and kinetic terms
196 !  are added
197    fermie1_new=fe1fixed+fe1_scf+fe1_paw(1)
198    fermie1rs=(fermie1-fermie1_new)**2
199    fermie1=fermie1_new
200 
201    if(prtvol>=10)then
202      write(msg, '(a,i5,2es18.8)' ) ' fermie1, residual squared',istep,fermie1,fermie1rs
203      call wrtout(std_out,msg,'COLL')
204    end if
205 
206  else
207    fermie1=zero
208  end if
209 
210 !Destroy atom table used for parallelism
211  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
212 
213 end subroutine newfermie1