TABLE OF CONTENTS


ABINIT/dfpt_mkvxc [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_mkvxc

FUNCTION

 Compute the first-order change of exchange-correlation potential
 due to atomic displacement: assemble the first-order
 density change with the frozen-core density change, then use
 the exchange-correlation kernel.

INPUTS

  cplex= if 1, real space 1-order functions on FFT grid are REAL,
         if 2, COMPLEX
  ixc= choice of exchange-correlation scheme
  kxc(nfft,nkxc)=exchange and correlation kernel (see below)
  mpi_enreg=information about MPI parallelization
  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
  nhat1(cplex*nfft,2nspden*nhat1dim)= -PAW only- 1st-order compensation density
  nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise
  nhat1gr(cplex*nfft,nspden,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density
  nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise
  nkxc=second dimension of the kxc array
  nspden=number of spin-density components
  n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used, otherwise, nfft
  option=if 0, work only with the XC core-correction,
         if 1, treat both density change and XC core correction
         if 2, treat only density change
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  rhor1(cplex*nfft,nspden)=array for electron density in electrons/bohr**3.
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
  xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc

OUTPUT

  vxc1(cplex*nfft,nspden)=change in exchange-correlation potential (including
   core-correction, if applicable)

NOTES

  Content of Kxc array:
   ===== if LDA
    if nspden==1: kxc(:,1)= d2Exc/drho2
                 (kxc(:,2)= d2Exc/drho_up drho_dn)
    if nspden>=2: kxc(:,1)= d2Exc/drho_up drho_up
                  kxc(:,2)= d2Exc/drho_up drho_dn
                  kxc(:,3)= d2Exc/drho_dn drho_dn
   ===== if GGA
    if nspden==1:
       kxc(:,1)= d2Exc/drho2
       kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
       kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
       kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
       kxc(:,5)= gradx(rho)
       kxc(:,6)= grady(rho)
       kxc(:,7)= gradz(rho)
    if nspden>=2:
       kxc(:,1)= d2Exc/drho_up drho_up
       kxc(:,2)= d2Exc/drho_up drho_dn
       kxc(:,3)= d2Exc/drho_dn drho_dn
       kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
       kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
       kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
       kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
       kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
       kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
       kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
       kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
       kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
       kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
       kxc(:,14)=gradx(rho_up)
       kxc(:,15)=gradx(rho_dn)
       kxc(:,16)=grady(rho_up)
       kxc(:,17)=grady(rho_dn)
       kxc(:,18)=gradz(rho_up)
       kxc(:,19)=gradz(rho_dn)

PARENTS

      dfpt_dyxc1,dfpt_mkvxc_noncoll,dfpt_nstdy,dfpt_nstpaw,dfpt_rhotov
      dfptnl_loop,m_kxc,nres2vres

CHILDREN

      dfpt_mkvxcgga,matr3inv,timab

SOURCE

137 subroutine dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1,nhat1dim,nhat1gr,nhat1grdim,&
138 &          nkxc,nspden,n3xccc,option,paral_kgb,qphon,rhor1,rprimd,usexcnhat,vxc1,xccc3d1)
139 
140 
141 !This section has been created automatically by the script Abilint (TD).
142 !Do not modify the following lines by hand.
143 #undef ABI_FUNC
144 #define ABI_FUNC 'dfpt_mkvxc'
145 !End of the abilint section
146 
147  implicit none
148 
149 !Arguments ------------------------------------
150 !scalars
151  integer,intent(in) :: cplex,ixc,n3xccc,nfft,nhat1dim,nhat1grdim
152  integer,intent(in) :: nkxc,nspden,option,paral_kgb,usexcnhat
153  type(MPI_type),intent(in) :: mpi_enreg
154 !arrays
155  integer,intent(in) :: ngfft(18)
156  real(dp),intent(in),target :: nhat1(cplex*nfft,nspden*nhat1dim)
157  real(dp),intent(in),target :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim)
158  real(dp),intent(in) :: kxc(nfft,nkxc),qphon(3)
159  real(dp),intent(in),target :: rhor1(cplex*nfft,nspden)
160  real(dp),intent(in) :: rprimd(3,3),xccc3d1(cplex*n3xccc)
161  real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
162 
163 !Local variables-------------------------------
164 !scalars
165  integer :: ii,ir,ispden,nhat1dim_,nhat1rgdim_
166  real(dp) :: rho1_dn,rho1_up,rho1im_dn,rho1im_up,rho1re_dn,rho1re_up
167  real(dp) :: spin_scale
168 !arrays
169  real(dp) :: gprimd(3,3),tsec(2)
170  real(dp), ABI_CONTIGUOUS pointer :: nhat1_(:,:),nhat1gr_(:,:,:),rhor1_(:,:)
171 
172 ! *************************************************************************
173 
174  DBG_ENTER("COLL")
175 
176  call timab(181,1,tsec)
177 
178  if(nspden/=1 .and. nspden/=2) then
179    MSG_BUG('For nspden==4 please use dfpt_mkvxc_noncoll!')
180  end if
181 
182 !Special case: no XC applied
183  if (ixc==0.or.nkxc==0) then
184    MSG_WARNING('Note that no xc is applied (ixc=0)')
185    vxc1=zero
186    return
187  end if
188 
189 !Treat first LDA
190  if(nkxc==1.or.nkxc==3)then
191 
192 !  PAW: eventually substract compensation density
193    if (option/=0) then
194      if (usexcnhat==0.and.nhat1dim==1) then
195        ABI_ALLOCATE(rhor1_,(cplex*nfft,nspden))
196        rhor1_(:,:)=rhor1(:,:)-nhat1(:,:)
197      else
198        rhor1_ => rhor1
199      end if
200    end if
201 
202 !  Case without non-linear core correction
203    if(n3xccc==0 .or. option==2)then
204 
205      if(option==0)then  ! No straight XC to compute
206 
207        vxc1(:,:)=zero
208 
209      else               ! XC, without non-linear XC correction
210 
211 !      Non-spin-polarized
212        if(nspden==1)then
213          if(cplex==1)then
214            do ir=1,nfft
215              vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,1)
216            end do
217          else
218            do ir=1,nfft
219              vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,1)
220              vxc1(2*ir  ,1)=kxc(ir,1)*rhor1_(2*ir  ,1)
221            end do
222          end if ! cplex==1
223 
224 !        Spin-polarized
225        else
226          if(cplex==1)then
227            do ir=1,nfft
228              rho1_dn=rhor1_(ir,1)-rhor1_(ir,2)
229              vxc1(ir,1)=kxc(ir,1)*rhor1_(ir,2)+kxc(ir,2)*rho1_dn
230              vxc1(ir,2)=kxc(ir,2)*rhor1_(ir,2)+kxc(ir,3)*rho1_dn
231            end do
232          else
233            do ir=1,nfft
234              rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2)
235              rho1im_dn=rhor1_(2*ir  ,1)-rhor1_(2*ir  ,2)
236              vxc1(2*ir-1,1)=kxc(ir,1)*rhor1_(2*ir-1,2)+kxc(ir,2)*rho1re_dn
237              vxc1(2*ir  ,1)=kxc(ir,1)*rhor1_(2*ir  ,2)+kxc(ir,2)*rho1im_dn
238              vxc1(2*ir-1,2)=kxc(ir,2)*rhor1_(2*ir-1,2)+kxc(ir,3)*rho1re_dn
239              vxc1(2*ir  ,2)=kxc(ir,2)*rhor1_(2*ir  ,2)+kxc(ir,3)*rho1im_dn
240            end do
241          end if ! cplex==1
242        end if ! nspden==1
243 
244      end if ! option==0
245 
246 !    Treat case with non-linear core correction
247    else
248 
249      if(option==0)then
250 
251        if(nspden==1)then
252          if(cplex==1)then
253            do ir=1,nfft
254              vxc1(ir,1)=kxc(ir,1)*xccc3d1(ir)
255            end do
256          else
257            do ir=1,nfft
258              vxc1(2*ir-1,1)=kxc(ir,1)*xccc3d1(2*ir-1)
259              vxc1(2*ir  ,1)=kxc(ir,1)*xccc3d1(2*ir  )
260            end do
261          end if ! cplex==1
262        else
263          if(cplex==1)then
264            do ir=1,nfft
265              vxc1(ir,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(ir)*half
266              vxc1(ir,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(ir)*half
267            end do
268          else
269            do ir=1,nfft
270              vxc1(2*ir-1,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir-1)*half
271              vxc1(2*ir  ,1)=(kxc(ir,1)+kxc(ir,2))*xccc3d1(2*ir  )*half
272              vxc1(2*ir-1,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir-1)*half
273              vxc1(2*ir  ,2)=(kxc(ir,2)+kxc(ir,3))*xccc3d1(2*ir  )*half
274            end do
275          end if ! cplex==1
276        end if ! nspden==1
277 
278      else ! option/=0
279 
280        if(nspden==1)then
281          if(cplex==1)then
282            do ir=1,nfft
283              vxc1(ir,1)=kxc(ir,1)*(rhor1_(ir,1)+xccc3d1(ir))
284            end do
285          else
286            do ir=1,nfft
287              vxc1(2*ir-1,1)=kxc(ir,1)*(rhor1_(2*ir-1,1)+xccc3d1(2*ir-1))
288              vxc1(2*ir  ,1)=kxc(ir,1)*(rhor1_(2*ir  ,1)+xccc3d1(2*ir  ))
289            end do
290          end if ! cplex==1
291        else
292          if(cplex==1)then
293            do ir=1,nfft
294              rho1_dn=rhor1_(ir,1)-rhor1_(ir,2) + xccc3d1(ir)*half
295              rho1_up=rhor1_(ir,2)             + xccc3d1(ir)*half
296              vxc1(ir,1)=kxc(ir,1)*rho1_up+kxc(ir,2)*rho1_dn
297              vxc1(ir,2)=kxc(ir,2)*rho1_up+kxc(ir,3)*rho1_dn
298            end do
299          else
300            do ir=1,nfft
301              rho1re_dn=rhor1_(2*ir-1,1)-rhor1_(2*ir-1,2) + xccc3d1(2*ir-1)*half
302              rho1im_dn=rhor1_(2*ir  ,1)-rhor1_(2*ir  ,2) + xccc3d1(2*ir  )*half
303              rho1re_up=rhor1_(2*ir-1,2)                 + xccc3d1(2*ir-1)*half
304              rho1im_up=rhor1_(2*ir  ,2)                 + xccc3d1(2*ir  )*half
305              vxc1(2*ir-1,1)=kxc(ir,1)*rho1re_up+kxc(ir,2)*rho1re_dn
306              vxc1(2*ir  ,1)=kxc(ir,1)*rho1im_up+kxc(ir,2)*rho1im_dn
307              vxc1(2*ir-1,2)=kxc(ir,2)*rho1re_up+kxc(ir,3)*rho1re_dn
308              vxc1(2*ir  ,2)=kxc(ir,2)*rho1im_up+kxc(ir,3)*rho1im_dn
309            end do
310          end if ! cplex==1
311        end if ! nspden==1
312 
313      end if ! option==0
314 
315    end if ! n3xccc==0
316 
317    if (option/=0.and.usexcnhat==0.and.nhat1dim==1) then
318      ABI_DEALLOCATE(rhor1_)
319    end if
320 
321 !  Treat GGA
322  else if (nkxc==7.or.nkxc==19) then
323 
324 ! Transfer the data to spin-polarized storage
325 
326 ! Treat the density change
327    ABI_ALLOCATE(rhor1_,(cplex*nfft,nspden))
328    if (option==1 .or. option==2) then
329      if (nspden==1) then
330        do ir=1,cplex*nfft
331          rhor1_(ir,1)=rhor1(ir,1)
332        end do
333      else
334        do ir=1,cplex*nfft
335          rho1_dn=rhor1(ir,1)-rhor1(ir,2)
336          rhor1_(ir,1)=rhor1(ir,2)
337          rhor1_(ir,2)=rho1_dn
338        end do
339      end if
340    else
341      do ispden=1,nspden
342        do ir=1,cplex*nfft
343          rhor1_(ir,ispden)=zero
344        end do
345      end do
346    end if
347    if( (option==0 .or. option==1) .and. n3xccc/=0)then
348      spin_scale=one;if (nspden==2) spin_scale=half
349      do ispden=1,nspden
350        do ir=1,cplex*nfft
351          rhor1_(ir,ispden)=rhor1_(ir,ispden)+xccc3d1(ir)*spin_scale
352        end do
353      end do
354    end if
355 
356 !  PAW: treat also compensation density (and gradients)
357    nhat1dim_=nhat1dim ; nhat1rgdim_=nhat1grdim
358    if (option/=0.and.nhat1dim==1.and.nspden==2) then
359      ABI_ALLOCATE(nhat1_,(cplex*nfft,nspden))
360      do ir=1,cplex*nfft
361        rho1_dn=nhat1(ir,1)-nhat1(ir,2)
362        nhat1_(ir,1)=nhat1(ir,2)
363        nhat1_(ir,2)=rho1_dn
364      end do
365    else if (option==0) then
366      ABI_ALLOCATE(nhat1_,(0,0))
367      nhat1dim_=0
368    else
369      nhat1_ => nhat1
370    end if
371    if (option/=0.and.nhat1grdim==1.and.nspden==2) then
372      ABI_ALLOCATE(nhat1gr_,(cplex*nfft,nspden,3))
373      do ii=1,3
374        do ir=1,cplex*nfft
375          rho1_dn=nhat1gr(ir,1,ii)-nhat1gr(ir,2,ii)
376          nhat1gr_(ir,1,ii)=nhat1gr(ir,2,ii)
377          nhat1gr_(ir,2,ii)=rho1_dn
378        end do
379      end do
380    else if (option==0) then
381      ABI_ALLOCATE(nhat1gr_,(0,0,0))
382      nhat1rgdim_=0
383    else
384      nhat1gr_ => nhat1gr
385    end if
386 
387    call matr3inv(rprimd,gprimd)
388 
389    call dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat1_,nhat1dim_,&
390 &   nhat1gr_,nhat1rgdim_,nkxc,nspden,paral_kgb,qphon,rhor1_,usexcnhat,vxc1)
391 
392    ABI_DEALLOCATE(rhor1_)
393    if ((option==0).or.(nhat1dim==1.and.nspden==2)) then
394      ABI_DEALLOCATE(nhat1_)
395    end if
396    if ((option==0).or.(nhat1grdim==1.and.nspden==2)) then
397      ABI_DEALLOCATE(nhat1gr_)
398    end if
399 
400  else
401    MSG_BUG('Invalid nkxc!')
402 
403  end if ! LDA or GGA
404 
405  call timab(181,2,tsec)
406 
407  DBG_EXIT("COLL")
408 
409 end subroutine dfpt_mkvxc

ABINIT/dfpt_mkvxc_noncoll [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_mkvxc_noncoll

FUNCTION

 Compute the first-order change of exchange-correlation potential
 due to atomic displacement for non-collinear spins: assemble the first-order
 density change with the frozen-core density change, then use
 the exchange-correlation kernel.

INPUTS

  cplex= if 1, real space 1-order functions on FFT grid are REAL,
         if 2, COMPLEX
  ixc= choice of exchange-correlation scheme
  kxc(nfft,nkxc)=exchange and correlation kernel (see rhotoxc.F90)
  mpi_enreg=information about MPI parallelization
  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
  nhat(nfft,nspden*nhatdim)= -PAW only- GS compensation density
  nhatdim= -PAW only- 1 if nhat array is used ; 0 otherwise
  nhat1(cplex*nfft,nspden*nhat1dim)= -PAW only- 1st-order compensation density
  nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise
  nhat1gr(cplex*nfft,nspden,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density
  nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise
  nkxc=second dimension of the kxc array
  nspden=number of spin-density components
  n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used, otherwise, nfft
  optnc=option for non-collinear magnetism (nspden=4):
       1: the whole 2x2 Vres matrix is computed
       2: only Vres^{11} and Vres^{22} are computed
  option=if 0, work only with the XC core-correction,
         if 1, treat both density change and XC core correction
         if 2, treat only density change
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  rhor(nfft,nspden)=GS electron density in real space
  rhor1(cplex*nfft,nspden)=1st-order electron density in real space
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
  vxc(nfft,nspden)=GS XC potential

OUTPUT

  vxc1(cplex*nfft,nspden)=change in exchange-correlation potential (including
   core-correction, if applicable)

PARENTS

      dfpt_dyxc1,dfpt_nstdy,dfpt_nstpaw,dfpt_rhotov,nres2vres

CHILDREN

      dfpt_mkvxc,rotate_back_mag,rotate_back_mag_dfpt,rotate_mag,timab

SOURCE

772 subroutine dfpt_mkvxc_noncoll(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat,nhatdim,nhat1,nhat1dim,&
773 &          nhat1gr,nhat1grdim,nkxc,nspden,n3xccc,optnc,option,paral_kgb,qphon,rhor,rhor1,&
774 &          rprimd,usexcnhat,vxc,vxc1,xccc3d1,ixcrot)
775 
776 
777 !This section has been created automatically by the script Abilint (TD).
778 !Do not modify the following lines by hand.
779 #undef ABI_FUNC
780 #define ABI_FUNC 'dfpt_mkvxc_noncoll'
781 !End of the abilint section
782 
783  implicit none
784 
785 !Arguments ------------------------------------
786 !scalars
787  integer,intent(in) :: cplex,ixc,n3xccc,nfft,nhatdim,nhat1dim,nhat1grdim,optnc
788  integer,intent(in) :: nkxc,nspden,option,paral_kgb,usexcnhat
789  type(MPI_type),intent(in) :: mpi_enreg
790 !arrays
791  integer,intent(in) :: ngfft(18)
792  real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim)
793  real(dp),intent(in) :: kxc(nfft,nkxc)
794  real(dp),intent(in) :: vxc(nfft,nspden)
795  real(dp),intent(in) :: nhat(nfft,nspden*nhatdim),nhat1(cplex*nfft,nspden*nhat1dim)
796  real(dp),intent(in),target :: rhor(nfft,nspden),rhor1(cplex*nfft,nspden)
797  real(dp),intent(in) :: qphon(3),rprimd(3,3),xccc3d1(cplex*n3xccc)
798  real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
799  integer,optional,intent(in) :: ixcrot
800 !Local variables-------------------------------
801 !scalars
802 !arrays
803  real(dp) :: nhat1_zero(0,0),nhat1gr_zero(0,0,0),tsec(2)
804  real(dp),allocatable :: m_norm(:),rhor1_diag(:,:),vxc1_diag(:,:)
805  real(dp), ABI_CONTIGUOUS pointer :: mag(:,:),rhor_(:,:),rhor1_(:,:)
806 ! *************************************************************************
807 
808 !  Non-collinear magnetism
809 !  Has to locally "rotate" rho(r)^(1) (according to magnetization),
810 !  Compute Vxc(r)^(1) in the spin frame aligned with \vec{m} and rotate it back
811 
812  DBG_ENTER("COLL")
813  ABI_UNUSED(nhat1gr)
814 
815  call timab(181,1,tsec)
816 
817  if(nspden/=4) then
818    MSG_BUG('only for nspden=4!')
819  end if
820 
821  if(nkxc/=2*min(nspden,2)-1) then
822    MSG_BUG('nspden=4 works only with LSDA.')
823  end if
824 
825 !Special case: no XC applied
826  if (ixc==0.or.nkxc==0) then
827    MSG_WARNING('Note that no xc is applied (ixc=0)')
828    vxc1(:,:)=zero
829    return
830  end if
831 
832 
833 
834 !Treat first LDA
835  if(nkxc==1.or.nkxc==3)then
836 
837    vxc1(:,:)=zero
838 
839 !  PAW: possibly substract compensation density
840    if (usexcnhat==0.and.nhatdim==1) then
841      ABI_ALLOCATE(rhor_,(nfft,nspden))
842      rhor_(:,:) =rhor(:,:)-nhat(:,:)
843    else
844      rhor_ => rhor
845    end if
846    if (usexcnhat==0.and.nhat1dim==1) then
847      ABI_ALLOCATE(rhor1_,(cplex*nfft,nspden))
848      rhor1_(:,:)=rhor1(:,:)-nhat1(:,:)
849    else
850      rhor1_ => rhor1
851    end if
852 
853 !  Magnetization
854    mag => rhor_(:,2:4)
855    ABI_ALLOCATE(rhor1_diag,(cplex*nfft,2))
856    ABI_ALLOCATE(vxc1_diag,(cplex*nfft,2))
857    ABI_ALLOCATE(m_norm,(nfft))
858 
859 !  -- Rotate rho(r)^(1)
860 !  SPr: for option=0 the rhor is not used, only core density xccc3d1
861 !       rotate_mag is only to compute the m_norm
862    call rotate_mag(rhor1_,rhor1_diag,mag,nfft,cplex,mag_norm_out=m_norm,&
863 &   rho_out_format=2)
864 
865 !  -- Compute Vxc(r)^(1)=Kxc(r).rho(r)^(1)_rotated
866 !  Note for PAW: nhat has already been substracted; don't use it in dfpt_mkvxc
867 !                 (put all nhat options to zero).
868 !  The collinear routine dfpt_mkvxc wants a general density built as (tr[rho],rho_upup)
869    call dfpt_mkvxc(cplex,ixc,kxc,mpi_enreg,nfft,ngfft,nhat1_zero,0,nhat1gr_zero,0,&
870 &   nkxc,2,n3xccc,option,paral_kgb,qphon,rhor1_diag,rprimd,0,vxc1_diag,xccc3d1)
871 
872    !call test_rotations(0,1)
873 
874 !  -- Rotate back Vxc(r)^(1)
875    if (optnc==1) then
876      if(present(ixcrot)) then
877        call rotate_back_mag_dfpt(option,vxc1_diag,vxc1,vxc,kxc,rhor1_,mag,nfft,cplex,&
878 &       mag_norm_in=m_norm,rot_method=ixcrot)
879      else
880        call rotate_back_mag_dfpt(option,vxc1_diag,vxc1,vxc,kxc,rhor1_,mag,nfft,cplex,&
881 &       mag_norm_in=m_norm)
882      end if
883    else
884      call rotate_back_mag(vxc1_diag,vxc1,mag,nfft,mag_norm_in=m_norm)
885      vxc1(:,3:4)=zero
886    end if
887 
888    ABI_DEALLOCATE(rhor1_diag)
889    ABI_DEALLOCATE(vxc1_diag)
890    ABI_DEALLOCATE(m_norm)
891    if (usexcnhat==0.and.nhatdim==1) then
892      ABI_DEALLOCATE(rhor_)
893    end if
894    if (usexcnhat==0.and.nhat1dim==1) then
895      ABI_DEALLOCATE(rhor1_)
896    end if
897 
898  end if ! nkxc=1 or nkxc=3
899 
900  call timab(181,2,tsec)
901 
902  DBG_EXIT("COLL")
903 
904 end subroutine dfpt_mkvxc_noncoll

ABINIT/dfpt_mkvxcgga [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_mkvxcgga

FUNCTION

 Compute the first-order change of exchange-correlation potential
 in case of GGA functionals
 Use the exchange-correlation kernel.

INPUTS

  cplex= if 1, real space 1-order functions on FFT grid are REAL,
    if 2, COMPLEX
  gmet(3,3)=metrix tensor in G space in Bohr**-2.
  gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
  gsqcut=cutoff value on G**2 for sphere inside fft box.
  kxc(nfft,nkxc)=exchange and correlation kernel (see below)
  mpi_enreg=informations about MPI parallelization
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT
  nhat1(cplex*nfft,2*nhat1dim)= -PAW only- 1st-order compensation density
  nhat1dim= -PAW only- 1 if nhat1 array is used ; 0 otherwise
  nhat1gr(cplex*nfft,2,3*nhat1grdim)= -PAW only- gradients of 1st-order compensation density
  nhat1grdim= -PAW only- 1 if nhat1gr array is used ; 0 otherwise
  nkxc=second dimension of the kxc array
  nspden=number of spin-density components
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  rhor1tmp(cplex*nfft,2)=array for first-order electron spin-density
   in electrons/bohr**3 (second index corresponds to spin-up and spin-down)
  usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc

OUTPUT

  vxc1(cplex*nfft,nspden)=change in exchange-correlation potential

NOTES

  For the time being, a rather crude coding, to be optimized ...
  Content of Kxc array:
   ===== if GGA
    if nspden==1:
       kxc(:,1)= d2Exc/drho2
       kxc(:,2)= 1/|grad(rho)| dExc/d|grad(rho)|
       kxc(:,3)= 1/|grad(rho)| d2Exc/d|grad(rho)| drho
       kxc(:,4)= 1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dExc/d|grad(rho)| )
       kxc(:,5)= gradx(rho)
       kxc(:,6)= grady(rho)
       kxc(:,7)= gradz(rho)
    if nspden>=2:
       kxc(:,1)= d2Exc/drho_up drho_up
       kxc(:,2)= d2Exc/drho_up drho_dn
       kxc(:,3)= d2Exc/drho_dn drho_dn
       kxc(:,4)= 1/|grad(rho_up)| dEx/d|grad(rho_up)|
       kxc(:,5)= 1/|grad(rho_dn)| dEx/d|grad(rho_dn)|
       kxc(:,6)= 1/|grad(rho_up)| d2Ex/d|grad(rho_up)| drho_up
       kxc(:,7)= 1/|grad(rho_dn)| d2Ex/d|grad(rho_dn)| drho_dn
       kxc(:,8)= 1/|grad(rho_up)| * d/d|grad(rho_up)| ( 1/|grad(rho_up)| dEx/d|grad(rho_up)| )
       kxc(:,9)= 1/|grad(rho_dn)| * d/d|grad(rho_dn)| ( 1/|grad(rho_dn)| dEx/d|grad(rho_dn)| )
       kxc(:,10)=1/|grad(rho)| dEc/d|grad(rho)|
       kxc(:,11)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_up
       kxc(:,12)=1/|grad(rho)| d2Ec/d|grad(rho)| drho_dn
       kxc(:,13)=1/|grad(rho)| * d/d|grad(rho)| ( 1/|grad(rho)| dEc/d|grad(rho)| )
       kxc(:,14)=gradx(rho_up)
       kxc(:,15)=gradx(rho_dn)
       kxc(:,16)=grady(rho_up)
       kxc(:,17)=grady(rho_dn)
       kxc(:,18)=gradz(rho_up)
       kxc(:,19)=gradz(rho_dn)

PARENTS

      dfpt_mkvxc

CHILDREN

      xcden,xcpot

SOURCE

486 subroutine dfpt_mkvxcgga(cplex,gprimd,kxc,mpi_enreg,nfft,ngfft,&
487 &                    nhat1,nhat1dim,nhat1gr,nhat1grdim,nkxc,&
488 &                    nspden,paral_kgb,qphon,rhor1,usexcnhat,vxc1)
489 
490 
491 !This section has been created automatically by the script Abilint (TD).
492 !Do not modify the following lines by hand.
493 #undef ABI_FUNC
494 #define ABI_FUNC 'dfpt_mkvxcgga'
495 !End of the abilint section
496 
497  implicit none
498 
499 !Arguments ------------------------------------
500 !scalars
501  integer,intent(in) :: cplex,nfft,nhat1dim,nhat1grdim,nkxc,nspden,paral_kgb,usexcnhat
502  type(MPI_type),intent(in) :: mpi_enreg
503 !arrays
504  integer,intent(in) :: ngfft(18)
505  real(dp),intent(in) :: gprimd(3,3)
506  real(dp),intent(in) :: kxc(nfft,nkxc)
507  real(dp),intent(in) :: nhat1(cplex*nfft,nspden*nhat1dim)
508  real(dp),intent(in) :: nhat1gr(cplex*nfft,nspden,3*nhat1grdim)
509  real(dp),intent(in) :: qphon(3)
510  real(dp),intent(in),target :: rhor1(cplex*nfft,nspden)
511  real(dp),intent(out) :: vxc1(cplex*nfft,nspden)
512 
513 !Local variables-------------------------------
514 !scalars
515  integer :: ii,ir,ishift,mgga,ngrad,nspgrad
516  logical :: test_nhat
517  real(dp) :: coeff_grho,coeff_grho_corr,coeff_grho_dn,coeff_grho_up
518  real(dp) :: coeffim_grho,coeffim_grho_corr,coeffim_grho_dn,coeffim_grho_up
519  real(dp) :: gradrho_gradrho1,gradrho_gradrho1_dn,gradrho_gradrho1_up
520  real(dp) :: gradrho_gradrho1im,gradrho_gradrho1im_dn,gradrho_gradrho1im_up
521  character(len=500) :: msg
522 !arrays
523  real(dp) :: r0(3),r0_dn(3),r0_up(3),r1(3),r1_dn(3),r1_up(3)
524  real(dp) :: r1im(3),r1im_dn(3),r1im_up(3)
525  real(dp),allocatable :: dnexcdn(:,:),rho1now(:,:,:)
526  real(dp),ABI_CONTIGUOUS pointer :: rhor1_ptr(:,:)
527 
528 ! *************************************************************************
529 
530  DBG_ENTER("COLL")
531 
532  if (nkxc/=12*min(nspden,2)-5) then
533    msg='Wrong nkxc value for GGA!'
534    MSG_BUG(msg)
535  end if
536 
537 !metaGGA contributions are not taken into account here
538  mgga=0
539 
540 !PAW: substract 1st-order compensation density from 1st-order density
541  test_nhat=((nhat1dim==1).and.(usexcnhat==0.or.nhat1grdim==1))
542  if (test_nhat) then
543    ABI_ALLOCATE(rhor1_ptr,(cplex*nfft,nspden))
544    rhor1_ptr(:,:)=rhor1(:,:)-nhat1(:,:)
545  else
546    rhor1_ptr => rhor1
547  end if
548 
549 !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,2,qphon,rhor1_ptr)
550 
551 !Compute the gradients of the first-order density
552 !rho1now(:,:,1) contains the first-order density, and
553 !rho1now(:,:,2:4) contains the gradients of the first-order density
554  ishift=0 ; ngrad=2
555  ABI_ALLOCATE(rho1now,(cplex*nfft,nspden,ngrad*ngrad))
556  call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,paral_kgb,qphon,rhor1_ptr,rho1now)
557 
558 !PAW: add "exact" gradients of compensation density
559  if (test_nhat.and.usexcnhat==1) then
560    rho1now(:,1:nspden,1)=rho1now(:,1:nspden,1)+nhat1(:,1:nspden)
561  end if
562  if (nhat1grdim==1) then
563    do ii=1,3
564      rho1now(:,1:nspden,ii+1)=rho1now(:,1:nspden,ii+1)+nhat1gr(:,1:nspden,ii)
565    end do
566  end if
567  if (test_nhat) then
568    ABI_DEALLOCATE(rhor1_ptr)
569  end if
570 
571 !Apply the XC kernel
572  nspgrad=2; if (nspden==2) nspgrad=5
573  ABI_ALLOCATE(dnexcdn,(cplex*nfft,nspgrad))
574 
575  if (cplex==1) then  ! Treat real case first
576    if (nspden==1) then
577      do ir=1,nfft
578        r0(:)=kxc(ir,5:7) ; r1(:)=rho1now(ir,1,2:4)
579        gradrho_gradrho1=dot_product(r0,r1)
580        dnexcdn(ir,1)=kxc(ir,1)*rho1now(ir,1,1) + kxc(ir,3)*gradrho_gradrho1
581        coeff_grho=kxc(ir,3)*rho1now(ir,1,1) + kxc(ir,4)*gradrho_gradrho1
582   !    Reuse the storage in rho1now
583        rho1now(ir,1,2:4)=r1(:)*kxc(ir,2)+r0(:)*coeff_grho
584      end do
585    else
586      do ir=1,nfft
587        do ii=1,3  ! grad of spin-up ans spin_dn GS rho
588          r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii)
589        end do
590        r0(:)=r0_up(:)+r0_dn(:)      ! grad of GS rho
591        r1_up(:)=rho1now(ir,1,2:4)   ! grad of spin-up rho1
592        r1_dn(:)=rho1now(ir,2,2:4)   ! grad of spin-down rho1
593        r1(:)=r1_up(:)+r1_dn(:)      ! grad of GS rho1
594        gradrho_gradrho1_up=dot_product(r0_up,r1_up)
595        gradrho_gradrho1_dn=dot_product(r0_dn,r1_dn)
596        gradrho_gradrho1   =dot_product(r0,r1)
597        dnexcdn(ir,1)=kxc(ir, 1)*rho1now(ir,1,1)     &
598 &       +kxc(ir, 2)*rho1now(ir,2,1)     &
599 &       +kxc(ir, 6)*gradrho_gradrho1_up &
600 &       +kxc(ir,11)*gradrho_gradrho1
601        dnexcdn(ir,2)=kxc(ir, 3)*rho1now(ir,2,1)     &
602 &       +kxc(ir, 2)*rho1now(ir,1,1)     &
603 &       +kxc(ir, 7)*gradrho_gradrho1_dn &
604 &       +kxc(ir,12)*gradrho_gradrho1
605        coeff_grho_corr=kxc(ir,11)*rho1now(ir,1,1) &
606 &       +kxc(ir,12)*rho1now(ir,2,1) &
607 &       +kxc(ir,13)*gradrho_gradrho1
608        coeff_grho_up=kxc(ir,6)*rho1now(ir,1,1)+kxc(ir,8)*gradrho_gradrho1_up
609        coeff_grho_dn=kxc(ir,7)*rho1now(ir,2,1)+kxc(ir,9)*gradrho_gradrho1_dn
610   !    Reuse the storage in rho1now
611        rho1now(ir,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) &
612 &       +kxc(ir,10)            *r1_dn(:) &
613 &       +coeff_grho_up         *r0_up(:) &
614 &       +coeff_grho_corr       *r0(:)
615        rho1now(ir,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) &
616 &       +kxc(ir,10)            *r1_up(:) &
617 &       +coeff_grho_dn         *r0_dn(:) &
618 &       +coeff_grho_corr       *r0(:)
619      end do
620    end if ! nspden
621 
622  else ! if cplex==2
623    if (nspden==1) then
624      do ir=1,nfft
625        r0(:)=kxc(ir,5:7)
626        r1(:)  =rho1now(2*ir-1,1,2:4)
627        r1im(:)=rho1now(2*ir  ,1,2:4)
628        gradrho_gradrho1  =dot_product(r0,r1)
629        gradrho_gradrho1im=dot_product(r0,r1im)
630        dnexcdn(2*ir-1,1)=kxc(ir,1)*rho1now(2*ir-1,1,1) + kxc(ir,3)*gradrho_gradrho1
631        dnexcdn(2*ir  ,1)=kxc(ir,1)*rho1now(2*ir  ,1,1) + kxc(ir,3)*gradrho_gradrho1im
632        coeff_grho  =kxc(ir,3)*rho1now(2*ir-1,1,1) + kxc(ir,4)*gradrho_gradrho1
633        coeffim_grho=kxc(ir,3)*rho1now(2*ir  ,1,1) + kxc(ir,4)*gradrho_gradrho1im
634   !    Reuse the storage in rho1now
635        rho1now(2*ir-1,1,2:4)=r1(:)  *kxc(ir,2)+r0(:)*coeff_grho
636        rho1now(2*ir  ,1,2:4)=r1im(:)*kxc(ir,2)+r0(:)*coeffim_grho
637      end do
638    else
639      do ir=1,nfft
640        do ii=1,3  ! grad of spin-up ans spin_dn GS rho
641          r0_up(ii)=kxc(ir,13+2*ii);r0_dn(ii)=kxc(ir,12+2*ii)-kxc(ir,13+2*ii)
642        end do
643        r0(:)=r0_up(:)+r0_dn(:)          ! grad of GS rho
644        r1_up(:)=rho1now(2*ir-1,1,2:4)   ! grad of spin-up rho1
645        r1im_up(:)=rho1now(2*ir,1,2:4)   ! grad of spin-up rho1 , im part
646        r1_dn(:)=rho1now(2*ir-1,2,2:4)   ! grad of spin-down rho1
647        r1im_dn(:)=rho1now(2*ir,2,2:4)   ! grad of spin-down rho1 , im part
648        r1(:)=r1_up(:)+r1_dn(:)      ! grad of GS rho1
649        r1im(:)=r1im_up(:)+r1im_dn(:)      ! grad of GS rho1, im part
650        gradrho_gradrho1_up  =dot_product(r0_up,r1_up)
651        gradrho_gradrho1_dn  =dot_product(r0_dn,r1_dn)
652        gradrho_gradrho1     =dot_product(r0,r1)
653        gradrho_gradrho1im_up=dot_product(r0_up,r1im_up)
654        gradrho_gradrho1im_dn=dot_product(r0_dn,r1im_dn)
655        gradrho_gradrho1im   =dot_product(r0,r1im)
656        dnexcdn(2*ir-1,1)=kxc(ir, 1)*rho1now(2*ir-1,1,1) &
657 &       +kxc(ir, 2)*rho1now(2*ir-1,2,1) &
658 &       +kxc(ir, 6)*gradrho_gradrho1_up &
659 &       +kxc(ir,11)*gradrho_gradrho1
660        dnexcdn(2*ir-1,2)=kxc(ir, 3)*rho1now(2*ir-1,2,1) &
661 &       +kxc(ir, 2)*rho1now(2*ir-1,1,1) &
662 &       +kxc(ir, 7)*gradrho_gradrho1_dn &
663 &       +kxc(ir,12)*gradrho_gradrho1
664        dnexcdn(2*ir  ,1)=kxc(ir, 1)*rho1now(2*ir  ,1,1) &
665 &       +kxc(ir, 2)*rho1now(2*ir  ,2,1) &
666 &       +kxc(ir, 6)*gradrho_gradrho1im_up &
667 &       +kxc(ir,11)*gradrho_gradrho1im
668        dnexcdn(2*ir  ,2)=kxc(ir, 3)*rho1now(2*ir  ,2,1) &
669 &       +kxc(ir, 2)*rho1now(2*ir  ,1,1) &
670 &       +kxc(ir, 7)*gradrho_gradrho1im_dn &
671 &       +kxc(ir,12)*gradrho_gradrho1im
672        coeff_grho_corr  =kxc(ir,11)*rho1now(2*ir-1,1,1) &
673 &       +kxc(ir,12)*rho1now(2*ir-1,2,1) &
674 &       +kxc(ir,13)*gradrho_gradrho1
675        coeffim_grho_corr=kxc(ir,11)*rho1now(2*ir  ,1,1) &
676 &       +kxc(ir,12)*rho1now(2*ir  ,2,1) &
677 &       +kxc(ir,13)*gradrho_gradrho1im
678        coeff_grho_up  =kxc(ir,6)*rho1now(2*ir-1,1,1)+kxc(ir,8)*gradrho_gradrho1_up
679        coeff_grho_dn  =kxc(ir,7)*rho1now(2*ir-1,2,1)+kxc(ir,9)*gradrho_gradrho1_dn
680        coeffim_grho_up=kxc(ir,6)*rho1now(2*ir  ,1,1)+kxc(ir,8)*gradrho_gradrho1im_up
681        coeffim_grho_dn=kxc(ir,7)*rho1now(2*ir  ,2,1)+kxc(ir,9)*gradrho_gradrho1im_dn
682 !      Reuse the storage in rho1now
683        rho1now(2*ir-1,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1_up(:) &
684 &       +kxc(ir,10)            *r1_dn(:) &
685 &       +coeff_grho_up         *r0_up(:) &
686 &       +coeff_grho_corr*r0(:)
687        rho1now(2*ir-1,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1_dn(:) &
688 &       +kxc(ir,10)            *r1_up(:) &
689 &       +coeff_grho_dn         *r0_dn(:) &
690 &       +coeff_grho_corr*r0(:)
691        rho1now(2*ir  ,1,2:4)=(kxc(ir,4)+kxc(ir,10))*r1im_up(:) &
692 &       +kxc(ir,10)            *r1im_dn(:) &
693 &       +coeffim_grho_up       *r0_up(:)   &
694 &       +coeffim_grho_corr     *r0(:)
695        rho1now(2*ir  ,2,2:4)=(kxc(ir,5)+kxc(ir,10))*r1im_dn(:) &
696 &       +kxc(ir,10)            *r1im_up(:) &
697 &       +coeffim_grho_dn       *r0_dn(:)   &
698 &       +coeffim_grho_corr     *r0(:)
699      end do
700    end if ! nspden
701 
702  end if
703 
704  vxc1(:,:)=zero
705  call xcpot(cplex,dnexcdn,gprimd,ishift,mgga,mpi_enreg,nfft,ngfft,ngrad,nspden,&
706 & nspgrad,paral_kgb,qphon,rho1now,vxc1)
707 
708 !call filterpot(paral_kgb,cplex,gmet,gsqcut,nfft,ngfft,nspden,qphon,vxc1)
709 
710  ABI_DEALLOCATE(dnexcdn)
711  ABI_DEALLOCATE(rho1now)
712 
713  DBG_EXIT("COLL")
714 
715 end subroutine dfpt_mkvxcgga

ABINIT/m_dfpt_mkvxc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_mkvxc

FUNCTION

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_dfpt_mkvxc
28 
29  use defs_basis
30  use defs_abitypes
31  use m_errors
32  use m_abicore
33  use m_xc_noncoll
34 
35  use m_time,     only : timab
36  use m_symtk,    only : matr3inv
37  use m_xctk,     only : xcden, xcpot
38 
39  implicit none
40 
41  private