TABLE OF CONTENTS


ABINIT/bare_vqg [ Functions ]

[ Top ] [ Functions ]

NAME

 bare_vqg

FUNCTION

 Compute bare coulomb term in G-space on the FFT mesh i.e. 4pi/(G+q)**2

INPUTS

  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2))
  divgq0= value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq. Used if q = Gamma
  gmet(3,3)=metrix tensor in G space in Bohr**-2.
  izero=if 1, unbalanced components of V(q,g) are set to zero
  hyb_mixing=hybrid mixing coefficient for the Fock contribution
  hyb_mixing_sr=hybrid mixing coefficient for the short-range Fock contribution
  hyb_range_fock=hybrid range for separation
  nfft=Total number of FFT grid points.
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft

OUTPUT

  vqg(nfft)=4pi/(G+q)**2, G=0 component is set to divgq0/pi if q = Gamma.

NOTES

  This routine operates on the full FFT mesh. DO NOT PASS MPI_TYPE
  One can easily implemente MPI-FFT by just calling this routine and then
  extracting the G-vectors treated by the node.

PARENTS

      fock_getghc

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

2258 subroutine bare_vqg(qphon,gsqcut,gmet,izero,hyb_mixing,hyb_mixing_sr,hyb_range_fock,nfft,nkpt_bz,ngfft,ucvol,vqg)
2259 
2260 
2261 !This section has been created automatically by the script Abilint (TD).
2262 !Do not modify the following lines by hand.
2263 #undef ABI_FUNC
2264 #define ABI_FUNC 'bare_vqg'
2265 !End of the abilint section
2266 
2267  implicit none
2268 
2269 !Arguments ------------------------------------
2270 !scalars
2271  integer,intent(in) :: izero,nfft,nkpt_bz
2272  real(dp),intent(in) :: gsqcut,hyb_mixing,hyb_mixing_sr,hyb_range_fock,ucvol
2273 !arrays
2274  integer,intent(in) :: ngfft(18)
2275  real(dp),intent(in) :: gmet(3,3),qphon(3)
2276  real(dp),intent(out) ::  vqg(nfft)
2277 
2278 !Local variables-------------------------------
2279 !scalars
2280  integer,parameter :: cplex1=1
2281  integer :: i1,i2,i23,i3,id1,id2,id3
2282  integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max
2283  integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05
2284  real(dp),parameter :: tolfix=1.000000001e0_dp ! Same value as the one used in hartre
2285  real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,rcut,divgq0
2286  character(len=100) :: msg
2287 !arrays
2288  integer :: id(3)
2289  real(dp),allocatable :: gq(:,:)
2290 
2291 ! *************************************************************************
2292 
2293  if (abs(hyb_mixing_sr)>tol8.and.abs(hyb_range_fock)<tol8) then
2294    msg='SR mixing<>0 while range separation=0!'
2295    MSG_BUG(msg)
2296  end if
2297 
2298 !Treatment of the divergence at q+g=zero
2299 !For the time being, only Spencer-Alavi scheme...
2300  rcut= (three*nkpt_bz*ucvol/four_pi)**(one/three)
2301  divgq0= two_pi*rcut**two
2302 !divgq0=zero
2303 !Initialize a few quantities
2304  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
2305  cutoff=gsqcut*tolfix
2306  vqg=zero
2307 
2308 !Some peculiar values of q
2309  qeq0=0; if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
2310  qeq05=0
2311  if (qeq0==0) then
2312    if (abs(abs(qphon(1))-half)<tol12.or.abs(abs(qphon(2))-half)<tol12.or. &
2313 &   abs(abs(qphon(3))-half)<tol12) qeq05=1
2314  end if
2315 
2316 !In order to speed the routine, precompute the components of g+q
2317 !Also check if the booked space was large enough...
2318  ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
2319  do ii=1,3
2320    id(ii)=ngfft(ii)/2+2
2321    do ing=1,ngfft(ii)
2322      ig=ing-(ing/id(ii))*ngfft(ii)-1
2323      gq(ii,ing)=ig+qphon(ii)
2324    end do
2325  end do
2326  ig1max=-1;ig2max=-1;ig3max=-1
2327  ig1min=n1;ig2min=n2;ig3min=n3
2328 
2329  id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
2330 
2331  ! Triple loop on each dimension
2332  do i3=1,n3
2333    ig3=i3-(i3/id3)*n3-1
2334    ! Precompute some products that do not depend on i2 and i1
2335    gs3=gq(3,i3)*gq(3,i3)*gmet(3,3)
2336    gqgm23=gq(3,i3)*gmet(2,3)*2
2337    gqgm13=gq(3,i3)*gmet(1,3)*2
2338 
2339    do i2=1,n2
2340      ig2=i2-(i2/id2)*n2-1
2341      gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23)
2342      gqgm12=gq(2,i2)*gmet(1,2)*2
2343      gqg2p3=gqgm13+gqgm12
2344 
2345      i23=n1*(i2-1 +(n2)*(i3-1))
2346      ! Do the test that eliminates the Gamma point outside of the inner loop
2347      ii1=1
2348      if (i23==0 .and. qeq0==1  .and. ig2==0 .and. ig3==0) then
2349        ii1=2
2350        ! value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq
2351        vqg(1+i23)=hyb_mixing*divgq0
2352 
2353 !      Note the combination of Spencer-Alavi and Erfc screening
2354        if (abs(hyb_range_fock)>tol8)then
2355          vqg(1+i23)=vqg(1+i23)+hyb_mixing_sr*(pi/hyb_range_fock**2)
2356 !        This would give a combination of Spencer-Alavi and Erfc screening,
2357 !        unfortunately, it modifies also the tests for pure HSE06, so was not retained.
2358 !        vqg(1+i23)=vqg(1+i23)+hyb_mixing_sr*min(divgq0,pi/(hyb_range_fock**2))
2359        endif
2360 
2361      end if
2362 
2363      ! Final inner loop on the first dimension (note the lower limit)
2364      do i1=ii1,n1
2365        gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3)
2366        ii=i1+i23
2367 
2368        if(gs<=cutoff)then
2369          ! Identify min/max indexes (to cancel unbalanced contributions later)
2370          ! Count (q+g)-vectors with similar norm
2371          if ((qeq05==1).and.(izero==1)) then
2372            ig1=i1-(i1/id1)*n1-1
2373            ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1)
2374            ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2)
2375            ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3)
2376          end if
2377 
2378          den=piinv/gs
2379 
2380 !        Spencer-Alavi screening
2381          if (abs(hyb_mixing)>tol8)then
2382            vqg(ii)=vqg(ii)+hyb_mixing*den*(one-cos(rcut*sqrt(four_pi/den)))
2383 !&         vqg(ii)=vqg(ii)+hyb_mixing*den
2384          endif
2385 !        Erfc screening
2386          if (abs(hyb_mixing_sr)>tol8) then
2387            vqg(ii)=vqg(ii)+hyb_mixing_sr*den*(one-exp(-pi/(den*hyb_range_fock**2)))
2388 !          This other possibility combines Erfc and Spencer-Alavi screening in case rcut is too small or hyb_range_fock too large
2389 !          if(divgq0<pi/(hyb_range_fock**2))then
2390 !            vqg(ii)=vqg(ii)+hyb_mixing_sr*den*&
2391 !&             (one-exp(-pi/(den*hyb_range_fock**2)))*(one-cos(rcut*sqrt(four_pi/den)))
2392 !          endif
2393          endif
2394 
2395        end if ! Cut-off
2396      end do ! End loop on i1
2397    end do ! End loop on i2
2398  end do ! End loop on i3
2399 
2400  if (izero==1) then
2401    ! Set contribution of unbalanced components to zero
2402    if (qeq0==1) then !q=0
2403      call zerosym(vqg,cplex1,n1,n2,n3)
2404    else if (qeq05==1) then
2405      !q=1/2; this doesn't work in parallel
2406      ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
2407      ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
2408      ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
2409      if (abs(abs(qphon(1))-half)<tol12) then
2410        if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
2411        if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
2412      end if
2413      if (abs(abs(qphon(2))-half)<tol12) then
2414        if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
2415        if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
2416      end if
2417      if (abs(abs(qphon(3))-half)<tol12) then
2418        if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
2419        if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
2420      end if
2421      call zerosym(vqg,cplex1,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
2422    end if
2423  end if
2424 
2425  ABI_DEALLOCATE(gq)
2426 
2427 end subroutine bare_vqg

ABINIT/m_fock [ Modules ]

[ Top ] [ Modules ]

NAME

  m_fock

FUNCTION

  This module provides the definition of
  the fock_type used to store data for the calculation of Fock exact exchange term
  and the procedures to perform this calculation.

COPYRIGHT

  Copyright (C) 2012-2018 ABINIT group (CMartins,FJ,FA,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

23 #if defined HAVE_CONFIG_H
24 #include "config.h"
25 #endif
26 
27 #include "abi_common.h"
28 
29 module m_fock
30 
31  use defs_basis
32  use defs_datatypes
33  use defs_abitypes
34  use m_abicore
35  use m_errors
36  use m_mpinfo
37  use m_xmpi
38  use libxc_functionals
39  use m_pawang
40  use m_pawtab
41  use m_pawfgr
42  use m_pawfgrtab
43  use m_pawcprj
44  use m_cgtools
45 
46  use m_time,            only : timab
47  use m_fstrings,        only : itoa, ftoa, sjoin
48  use m_symtk,           only : mati3inv, matr3inv
49  use m_fftcore,         only : sphereboundary
50  use m_fft,             only : zerosym, fourwf
51  use m_kg,              only : ph1d3d, getph
52  use m_kpts,            only : listkk
53 
54  implicit none
55 
56  private

ABINIT/strfock [ Functions ]

[ Top ] [ Functions ]

NAME

 strfock

FUNCTION

 Compute Fock energy contribution to stress tensor (Cartesian coordinates).

INPUTS

  gsqcut=cutoff value on $G^2$ for (large) sphere inside fft box.
  $gsqcut=(boxcut^2)*ecut/(2._dp*(\pi^2))$
  gprimd(3,3)=reciprocal space dimensional primitive translations
  hyb_mixing=hybrid mixing coefficient for the Fock contribution
  hyb_mixing_sr=hybrid mixing coefficient for the short-range Fock contribution
  hyb_range_fock=hybrid range for separation
  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
  nkpt_bz= number of k points in the BZ
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  rhog(2,nfft)=Fourier transform of charge density (bohr^-3)
  rhog2(2,nfft)= optional argument: Fourier transform of a second charge density (bohr^-3)
  ucvol=unit cell volume (bohr^3)
  vqg(nfft)=4pi/(G+q)**2

OUTPUT

  fockstr(6)=components of Fock part of stress tensor
   (Cartesian coordinates, symmetric tensor) in hartree/bohr^3
   Definition of symmetric tensor storage: store 6 unique components
   in the order 11, 22, 33, 32, 31, 21 (suggested by Xavier Gonze).

PARENTS

      fock_getghc

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

2469 subroutine strfock(gprimd,gsqcut,fockstr,hyb_mixing,hyb_mixing_sr,hyb_range_fock,mpi_enreg,nfft,ngfft,&
2470 &                  nkpt_bz,rhog,ucvol,qphon,&
2471 &                 rhog2) ! optional argument
2472 
2473 
2474 !This section has been created automatically by the script Abilint (TD).
2475 !Do not modify the following lines by hand.
2476 #undef ABI_FUNC
2477 #define ABI_FUNC 'strfock'
2478 !End of the abilint section
2479 
2480  implicit none
2481 
2482 !Arguments ------------------------------------
2483 !scalars
2484  integer,intent(in) :: nfft,nkpt_bz
2485  real(dp),intent(in) :: gsqcut,hyb_mixing,hyb_mixing_sr,hyb_range_fock,ucvol
2486  type(MPI_type),intent(in) :: mpi_enreg
2487 !arrays
2488  integer,intent(in) :: ngfft(18)
2489  real(dp),intent(in) :: gprimd(3,3),rhog(2,nfft),qphon(3)
2490  real(dp),intent(in),optional :: rhog2(2,nfft)
2491  real(dp),intent(out) :: fockstr(6)
2492 
2493 !Local variables-------------------------------
2494 !scalars
2495  integer,parameter :: im=2,re=1
2496  integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,me_fft,n1,n2,n3,nproc_fft
2497  real(dp) :: arg,cutoff,gsquar,rcut,rhogsq,tolfix=1.000000001_dp,tot,tot1,divgq0
2498  character(len=100) :: msg
2499 !arrays
2500  real(dp) :: gcart(3),tsec(2)
2501  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
2502  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
2503 
2504 ! *************************************************************************
2505 
2506  call timab(568,1,tsec)
2507 
2508  if (abs(hyb_mixing_sr)>tol8.and.abs(hyb_range_fock)<tol8) then
2509    msg='strfock: SR mixing<>0 while range separation=0!'
2510    MSG_BUG(msg)
2511  end if
2512 
2513  fockstr(:)=zero
2514  rcut= (three*nkpt_bz*ucvol/four_pi)**(one/three)
2515  irho2=0;if (present(rhog2)) irho2=1
2516  divgq0=two_pi/three*rcut**2
2517 
2518 !Conduct looping over all fft grid points to find G vecs inside gsqcut
2519 !Include G**2 on surface of cutoff sphere as well as inside:
2520  cutoff=gsqcut*tolfix
2521  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
2522  me_fft=ngfft(11)
2523  nproc_fft=ngfft(10)
2524  id1=n1/2+2
2525  id2=n2/2+2
2526  id3=n3/2+2
2527 
2528 
2529  ii=0
2530  ! Get the distrib associated with this fft_grid
2531  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2532 
2533  do i3=1,n3
2534    ig3=i3-(i3/id3)*n3-1
2535    do i2=1,n2
2536      ig2=i2-(i2/id2)*n2-1
2537      if (fftn2_distrib(i2)==me_fft) then
2538        do i1=1,n1
2539          tot=zero; tot1=zero
2540          ig1=i1-(i1/id1)*n1-1
2541          ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
2542 
2543 !        Compute cartesian components of G
2544          gcart(1)=gprimd(1,1)*(dble(ig1)+qphon(1))+gprimd(1,2)*(dble(ig2)+qphon(2))+gprimd(1,3)*(dble(ig3)+qphon(3))
2545          gcart(2)=gprimd(2,1)*(dble(ig1)+qphon(1))+gprimd(2,2)*(dble(ig2)+qphon(2))+gprimd(2,3)*(dble(ig3)+qphon(3))
2546          gcart(3)=gprimd(3,1)*(dble(ig1)+qphon(1))+gprimd(3,2)*(dble(ig2)+qphon(2))+gprimd(3,3)*(dble(ig3)+qphon(3))
2547 !        Compute |G+q|^2
2548          gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2
2549 !        take |rho(G)|^2 for complex rhog
2550          if (irho2==0) then
2551            rhogsq=rhog(re,ii)**2+rhog(im,ii)**2
2552          else
2553            rhogsq=rhog(re,ii)*rhog2(re,ii)+rhog(im,ii)*rhog2(im,ii)
2554          end if
2555 !        Case G=0:
2556          if(gsquar<tol10) then
2557            if (abs(hyb_mixing_sr)>tol8) cycle
2558            if (abs(hyb_mixing)>tol8) then
2559              fockstr(1)=fockstr(1)+hyb_mixing*divgq0*rhogsq
2560              fockstr(2)=fockstr(2)+hyb_mixing*divgq0*rhogsq
2561              fockstr(3)=fockstr(3)+hyb_mixing*divgq0*rhogsq
2562              cycle
2563            end if
2564          end if
2565 
2566 !        Spencer-Alavi screening
2567          if (abs(hyb_mixing)>tol8) then
2568            arg=two_pi*rcut*sqrt(gsquar)
2569            tot=hyb_mixing*rhogsq*piinv/(gsquar**2)*(1-cos(arg)-arg*sin(arg)/two)
2570            tot1=hyb_mixing*rhogsq/three*rcut*sin(arg)/sqrt(gsquar)
2571          end if
2572 
2573 !        Erfc screening
2574          if (abs(hyb_mixing_sr)>tol8) then
2575            arg=-gsquar*pi**2/(hyb_range_fock**2)
2576            tot=tot+hyb_mixing_sr*rhogsq*piinv/(gsquar**2)*(1.d0-exp(arg)*(1-arg))
2577          end if
2578          fockstr(1)=fockstr(1)+tot*gcart(1)*gcart(1)+tot1
2579          fockstr(2)=fockstr(2)+tot*gcart(2)*gcart(2)+tot1
2580          fockstr(3)=fockstr(3)+tot*gcart(3)*gcart(3)+tot1
2581          fockstr(4)=fockstr(4)+tot*gcart(3)*gcart(2)
2582          fockstr(5)=fockstr(5)+tot*gcart(3)*gcart(1)
2583          fockstr(6)=fockstr(6)+tot*gcart(2)*gcart(1)
2584        end do
2585      end if
2586    end do
2587  end do
2588 
2589 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel
2590 #ifdef FC_IBM
2591  write(std_out,*)' strfock : before mpi_comm, fockstr=',fockstr
2592 #endif
2593 
2594 !Init mpi_comm
2595  if(mpi_enreg%nproc_fft>1)then
2596    call timab(48,1,tsec)
2597    call xmpi_sum(fockstr,mpi_enreg%comm_fft ,ierr)
2598    call timab(48,2,tsec)
2599  end if
2600 
2601 #ifdef FC_IBM
2602 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel
2603  write(std_out,*)' strfock : after mpi_comm, fockstr=',fockstr
2604 #endif
2605 
2606 !Normalize and add term -efock/ucvol on diagonal
2607 !efock has been set to zero because it is not yet known. It will be added later.
2608  fockstr(1)=-fockstr(1)
2609  fockstr(2)=-fockstr(2)
2610  fockstr(3)=-fockstr(3)
2611  fockstr(4)=-fockstr(4)
2612  fockstr(5)=-fockstr(5)
2613  fockstr(6)=-fockstr(6)
2614 
2615  call timab(568,2,tsec)
2616 
2617 end subroutine strfock

m_fock/fock_ACE_destroy [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_ACE_destroy

FUNCTION

  Clean and destroy fock datastructure.

INPUTS

  fockACE <type(fock_ACE_type)>= all the quantities to calculate Fock exact exchange in the ACE context

PARENTS

      scfcv

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1497 subroutine fock_ACE_destroy(fockACE)
1498 
1499 
1500 !This section has been created automatically by the script Abilint (TD).
1501 !Do not modify the following lines by hand.
1502 #undef ABI_FUNC
1503 #define ABI_FUNC 'fock_ACE_destroy'
1504 !End of the abilint section
1505 
1506  implicit none
1507 
1508 !Arguments ------------------------------------
1509  type(fock_ACE_type),pointer :: fockACE(:,:)
1510 !Local variables-------------------------------
1511  integer :: dim1,dim2,ii,jj
1512 ! *************************************************************************
1513  DBG_ENTER("COLL")
1514 
1515  dim1=size(fockACE,1)
1516  dim2=size(fockACE,2)
1517  do jj=1,dim2
1518    do ii=1,dim1
1519      if (allocated(fockACE(ii,jj)%xi)) then
1520        ABI_DEALLOCATE(fockACE(ii,jj)%xi)
1521      end if
1522    end do
1523  end do
1524  DBG_EXIT("COLL")
1525 
1526 end subroutine fock_ACE_destroy

m_fock/fock_calc_ene [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_calc_ene

FUNCTION

  Calculate the Fock contribution to the total energy

INPUTS

  fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
  ikpt= reduced planewave coordinates.

OUTPUT

  none

SIDE EFFECTS

  energies <type(energies_type)>=storage for energies computed here :
   | e_exactX = Fock contribution to the total energy (Hartree)

NOTES

 If the cgocc_bz are not updated at each iteration, be careful to calculate Fock energy at the same frequency.
 TO CHECK == CHANGE IN SOME DEFINTIONS

PARENTS

      vtorho

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1560 subroutine fock_calc_ene(dtset,fock,fock_energy,ikpt,nband,occ)
1561 
1562 
1563 !This section has been created automatically by the script Abilint (TD).
1564 !Do not modify the following lines by hand.
1565 #undef ABI_FUNC
1566 #define ABI_FUNC 'fock_calc_ene'
1567 !End of the abilint section
1568 
1569  implicit none
1570 
1571 !Arguments ------------------------------------
1572 !scalars
1573  integer,intent(in) :: ikpt,nband
1574  real(dp),intent(inout) :: fock_energy
1575  type(dataset_type),intent(in) :: dtset
1576  type(fock_common_type),pointer :: fock
1577 !arrays
1578  real(dp),intent(in) :: occ(nband)
1579 
1580 !Local variables-------------------------------
1581  integer :: iband
1582 
1583 ! *************************************************************************
1584 
1585  do iband=1,nband
1586 
1587    ! Select only the occupied states (such that fock%occ_bz > 10^-8)
1588    if (abs(occ(iband))>tol8) then
1589 !     fock_energy=fock_energy + half*fock%eigen_ikpt(iband)*occ(iband)*dtset%wtk(ikpt)
1590      !* Sum the contribution of each occupied states at point k_i
1591      !* No need to multiply %wtk by ucvol since there is no factor 1/ucvol in the definition of %wtk
1592 
1593 !* accumulate Fock contributions to the forces.
1594 !     if (fock%optfor) then
1595        fock%forces(:,:)=fock%forces(:,:)+occ(iband)*dtset%wtk(ikpt)*fock%forces_ikpt(:,:,iband)
1596 !     endif
1597    end if
1598  end do
1599 
1600 
1601 end subroutine fock_calc_ene

m_fock/fock_destroy [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_destroy

FUNCTION

  Clean and destroy fock datastructure.

INPUTS

  fock <type(fock_type)>= all the quantities to calculate Fock exact exchange

PARENTS

      scfcv

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1287 subroutine fock_destroy(fock)
1288 
1289 
1290 !This section has been created automatically by the script Abilint (TD).
1291 !Do not modify the following lines by hand.
1292 #undef ABI_FUNC
1293 #define ABI_FUNC 'fock_destroy'
1294 !End of the abilint section
1295 
1296  implicit none
1297 
1298 !Arguments ------------------------------------
1299  type(fock_type),pointer :: fock
1300 
1301 ! *************************************************************************
1302 
1303  DBG_ENTER("COLL")
1304  if (fock%fock_common%use_ACE/=0) then
1305    ABI_DATATYPE_DEALLOCATE(fock%fockACE)
1306  end if
1307  ABI_DATATYPE_DEALLOCATE(fock%fock_common)
1308  ABI_DATATYPE_DEALLOCATE(fock%fock_BZ)
1309  ABI_DATATYPE_DEALLOCATE(fock)
1310 
1311  DBG_EXIT("COLL")
1312 end subroutine fock_destroy
1313 
1314 subroutine fock_common_destroy(fock)
1315 
1316 
1317 !This section has been created automatically by the script Abilint (TD).
1318 !Do not modify the following lines by hand.
1319 #undef ABI_FUNC
1320 #define ABI_FUNC 'fock_common_destroy'
1321 !End of the abilint section
1322 
1323  implicit none
1324 
1325 !Arguments ------------------------------------
1326  type(fock_common_type),pointer :: fock
1327 
1328 ! *************************************************************************
1329 
1330  DBG_ENTER("COLL")
1331 
1332   if (allocated(fock%atindx)) then
1333    ABI_DEALLOCATE(fock%atindx)
1334  endif
1335  if (allocated(fock%typat)) then
1336     ABI_DEALLOCATE(fock%typat)
1337  endif
1338  ! real arrays
1339  if (allocated(fock%forces)) then
1340    ABI_DEALLOCATE(fock%forces)
1341  endif
1342  if (allocated(fock%nband)) then
1343    ABI_DEALLOCATE(fock%nband)
1344  endif
1345  if (allocated(fock%forces_ikpt)) then
1346    ABI_DEALLOCATE(fock%forces_ikpt)
1347  endif
1348  if (allocated(fock%stress_ikpt)) then
1349    ABI_DEALLOCATE(fock%stress_ikpt)
1350  endif
1351  if (allocated(fock%eigen_ikpt)) then
1352     ABI_DEALLOCATE(fock%eigen_ikpt)
1353  endif
1354 !*Deallocate datatypes
1355  if (allocated(fock%pawfgrtab)) then
1356     call pawfgrtab_free(fock%pawfgrtab)
1357     ABI_DATATYPE_DEALLOCATE(fock%pawfgrtab)
1358  endif
1359 
1360  ! Put the integer to 0
1361  fock%ieigen=0
1362  fock%ikpt=0
1363  fock%isppol=0
1364 
1365  if (allocated(fock%symrec)) then
1366     ABI_DEALLOCATE(fock%symrec)
1367  endif
1368 
1369 !* [description of divergence in |q+G|=0]
1370 !* Put the real (dp) to 0
1371  fock%gsqcut=zero
1372  fock%hyb_mixing=zero
1373  fock%hyb_mixing_sr=zero
1374  fock%hyb_range_dft=zero
1375  fock%hyb_range_fock=zero
1376 
1377  DBG_EXIT("COLL")
1378 end subroutine fock_common_destroy
1379 
1380 
1381 subroutine fock_BZ_destroy(fock)
1382 
1383 
1384 !This section has been created automatically by the script Abilint (TD).
1385 !Do not modify the following lines by hand.
1386 #undef ABI_FUNC
1387 #define ABI_FUNC 'fock_BZ_destroy'
1388 !End of the abilint section
1389 
1390  implicit none
1391 
1392 !Arguments ------------------------------------
1393  type(fock_BZ_type),pointer :: fock
1394 
1395 ! *************************************************************************
1396 
1397  DBG_ENTER("COLL")
1398 
1399  if (allocated(fock%cwaveocc_bz)) then
1400    ABI_DEALLOCATE(fock%cwaveocc_bz)
1401  endif
1402  if (allocated(fock%cgocc)) then
1403    ABI_DEALLOCATE(fock%cgocc)
1404  endif
1405  if (allocated(fock%npwarr)) then
1406    ABI_DEALLOCATE(fock%npwarr)
1407  endif
1408  if (allocated(fock%occ_bz)) then
1409    ABI_DEALLOCATE(fock%occ_bz)
1410  endif
1411  if (allocated(fock%cwaveocc_prj)) then
1412    call pawcprj_free(fock%cwaveocc_prj)
1413    ABI_DATATYPE_DEALLOCATE(fock%cwaveocc_prj)
1414  endif
1415  ! Deallocate integer arrays
1416 
1417  if (allocated(fock%kg_bz)) then
1418    ABI_DEALLOCATE(fock%kg_bz)
1419  endif
1420  if (allocated(fock%nbandocc_bz)) then
1421    ABI_DEALLOCATE(fock%nbandocc_bz)
1422  endif
1423  if (allocated(fock%istwfk_bz)) then
1424    ABI_DEALLOCATE(fock%istwfk_bz)
1425  endif
1426  if (allocated(fock%calc_phase)) then
1427     ABI_DEALLOCATE(fock%calc_phase)
1428  endif
1429  if (allocated(fock%timerev)) then
1430     ABI_DEALLOCATE(fock%timerev)
1431  endif
1432  if (allocated(fock%tab_ibg)) then
1433     ABI_DEALLOCATE(fock%tab_ibg)
1434  endif
1435  if (allocated(fock%tab_icg)) then
1436     ABI_DEALLOCATE(fock%tab_icg)
1437  endif
1438  if (allocated(fock%tab_icp)) then
1439     ABI_DEALLOCATE(fock%tab_icp)
1440  endif
1441  if (allocated(fock%tab_ikpt)) then
1442     ABI_DEALLOCATE(fock%tab_ikpt)
1443  endif
1444  if (allocated(fock%tab_symkpt)) then
1445     ABI_DEALLOCATE(fock%tab_symkpt)
1446  endif
1447 
1448 !* [description of IBZ and BZ]
1449 !* Deallocate real arrays
1450  if (allocated(fock%wtk_bz)) then
1451    ABI_DEALLOCATE(fock%wtk_bz)
1452  endif
1453  if (allocated(fock%kptns_bz)) then
1454     ABI_DEALLOCATE(fock%kptns_bz)
1455  endif
1456  if (allocated(fock%phase)) then
1457     ABI_DEALLOCATE(fock%phase)
1458  endif
1459 !* Put the integer to 0
1460    fock%nkpt_bz=0
1461 
1462 !* Deallocate real arrays
1463 
1464 !* Deallocate integer arrays
1465  if (allocated(fock%gbound_bz)) then
1466     ABI_DEALLOCATE(fock%gbound_bz)
1467  endif
1468 
1469 !* [description of size of arrays/pointers]
1470 !* Put the integer to 0
1471  fock%mkpt=0
1472  fock%mkptband=0
1473  call destroy_mpi_enreg(fock%mpi_enreg)
1474 
1475  DBG_EXIT("COLL")
1476 
1477 end subroutine fock_BZ_destroy

m_fock/fock_get_getghc_call [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_get_getghc_call

FUNCTION

  Returns the value of fock%getghc_call_

PARENTS

SOURCE

2116 pure integer function fock_get_getghc_call(fock)
2117 
2118 
2119 !This section has been created automatically by the script Abilint (TD).
2120 !Do not modify the following lines by hand.
2121 #undef ABI_FUNC
2122 #define ABI_FUNC 'fock_get_getghc_call'
2123 !End of the abilint section
2124 
2125  implicit none
2126 
2127 !Arguments ------------------------------------
2128 !scalars
2129  type(fock_common_type),intent(in) :: fock
2130 
2131 ! *************************************************************************
2132 
2133  fock_get_getghc_call = fock%getghc_call_
2134 
2135 end function fock_get_getghc_call

m_fock/fock_init [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_init

FUNCTION

  Init all scalars, arrays and pointers in the structure.

INPUTS

  cg(2,mcg)= wavefunctions
  dtset <type(dataset_type)>= all input variables for this dataset
  gsqcut= Fourier cutoff on G^2 used to calculate charge density
  kg(3,mpw*mkmem)= reduced planewave coordinates.
  mcg= size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mpi_enreg=information about MPI parallelization
  npwarr_bz(nkpt)= number of planewaves in basis at this k point
  occ(mband*nkpt*nsppol)= occupation number for each band (often 2) at each k point

OUTPUT

  none

SIDE EFFECTS

  fock <type(fock_type)>= all the quantities to calculate Fock exact exchange are initialized

NOTES

  ############################
  ### Not fully tested yet ###
  ############################

  The current version is restricted to the case nsym=1, nspinor=1 and mkmem/=0.

PARENTS

      scfcv

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

 485 subroutine fock_init(atindx,cplex,dtset,fock,gsqcut,kg,mpi_enreg,nattyp,npwarr,pawang,pawfgr,pawtab,rprimd)
 486 
 487 
 488 !This section has been created automatically by the script Abilint (TD).
 489 !Do not modify the following lines by hand.
 490 #undef ABI_FUNC
 491 #define ABI_FUNC 'fock_init'
 492 !End of the abilint section
 493 
 494  implicit none
 495 
 496 !Arguments ------------------------------------
 497 !scalars
 498  integer, intent(in) :: cplex
 499  real(dp),intent(in) :: gsqcut
 500  type(dataset_type),intent(in) :: dtset
 501  type(MPI_type),intent(in) :: mpi_enreg
 502  type(fock_type),intent(inout),pointer :: fock
 503  type(pawfgr_type),intent(in),target :: pawfgr
 504  type(pawang_type),intent(in),target :: pawang
 505 !arrays
 506  integer, intent(in) :: atindx(dtset%natom),nattyp(dtset%ntypat), npwarr(dtset%nkpt)
 507  integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem)
 508  real(dp), intent(in) :: rprimd(3,3)
 509  type(pawtab_type), intent(in),target :: pawtab(dtset%ntypat*dtset%usepaw)
 510 !Local variables-------------------------------
 511 !scalars
 512  integer :: iatom,ibg,icg,icp,ier,ik,ikg,ikpt,isppol,isym,itypat,jkpt,jpw,jsym,mband,mgfft,mkpt,mkptband
 513  integer :: n1,n2,n3,n4,n5,n6,nband,ncpgr,nkpt_bz,nproc_hf,npwj,timrev,use_ACE,v1,v2,v3
 514  integer :: my_jkpt,jkg_this_proc,my_nsppol,my_nspinor
 515  real(dp) :: dksqmax,arg
 516  character(len=500) :: msg
 517 !arrays
 518  integer :: indx(1),l_size_atm(dtset%natom),shiftg(3),symm(3,3),ident(3,3),symrec(3,3,dtset%nsym)
 519  real(dp) :: gmet(3,3),gprimd(3,3),tau_nons(3),phktnons(2,1),tsec(2),Rtnons(3,dtset%nsym)
 520  integer,allocatable :: dimcprj(:),indkk(:,:),kg_tmp(:),my_ikgtab(:),my_ibgtab(:,:),my_icgtab(:,:),my_icptab(:,:),invsym(:)
 521  real(dp),allocatable :: kptns_hf(:,:), phase1d(:,:)
 522  type(fock_common_type),pointer :: fockcommon
 523  type(fock_BZ_type),pointer :: fockbz
 524 
 525 ! *************************************************************************
 526 
 527  DBG_ENTER("COLL")
 528 
 529  call timab(1500,1,tsec)
 530 
 531  if (dtset%nspinor/=1) then
 532    msg='Hartree-Fock option can be used only with option nspinor=1.'
 533    MSG_ERROR(msg)
 534  end if
 535 
 536 
 537 ! =====================================
 538 ! === Define useful local variables ===
 539 ! =====================================
 540 
 541  nkpt_bz=dtset%nkpthf
 542  nproc_hf=mpi_enreg%nproc_hf
 543  mband=dtset%nbandhf
 544 
 545  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
 546  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
 547 
 548 !* Allocations
 549  ABI_ALLOCATE(kptns_hf,(3,nkpt_bz))
 550  kptns_hf=zero
 551  ABI_ALLOCATE(indkk,(nkpt_bz,6))
 552  indkk=0
 553  ABI_ALLOCATE(phase1d,(2,(2*n1+1)*(2*n2+1)*(2*n3+1)))
 554  phase1d=zero
 555  ABI_ALLOCATE(kg_tmp,(3*dtset%mpw))
 556 
 557 !* Initialize the array my_ikgtab = shifts in arrays kg(ikg) associated to ikpt
 558  ABI_ALLOCATE(my_ikgtab,(dtset%nkpt))
 559  ikg=0
 560  do ikpt=1,dtset%nkpt
 561    nband=dtset%nband(ikpt)
 562    if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband,-1,mpi_enreg%me_kpt))) then
 563 !* The point ikpt is treated on this processor.
 564      my_ikgtab(ikpt)=ikg
 565 !* The array kg is distributed, the shift ikg is incremented only on this proc.
 566      ikg=ikg+npwarr(ikpt)
 567    else
 568      my_ikgtab(ikpt)=-1
 569 !* Default value is -1.
 570    end if
 571  end do
 572 
 573 !* Initialize the array my_ibgtab = shifts in arrays occ(ibg) associated to ikpt
 574 !* Initialize the array my_icgtab = shifts in arrays cg(icg) associated to ikpt
 575  ABI_ALLOCATE(my_ibgtab,(dtset%nkpt,dtset%nsppol))
 576  ABI_ALLOCATE(my_icgtab,(dtset%nkpt,dtset%nsppol))
 577  ABI_ALLOCATE(my_icptab,(dtset%nkpt,dtset%nsppol))
 578  ibg=0; icg=0 ;icp=0
 579  do isppol=1,dtset%nsppol
 580    do ikpt=1,dtset%nkpt
 581      nband=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 582      if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband,isppol,mpi_enreg%me_kpt))) then
 583 !* The states with (ikpt,isppol) are stored on this processor.
 584        my_icgtab(ikpt,isppol)=icg
 585        my_icptab(ikpt,isppol)=icp
 586 !* The array cg is distributed, the shift icg is incremented only on this proc.
 587        icg=icg+npwarr(ikpt)*nband
 588        icp=icp+nband
 589      else
 590        my_icgtab(ikpt,isppol)=-1
 591        my_icptab(ikpt,isppol)=-1
 592 !* Otherwise, the states with (ikpt,isspol) are not stored on this processor and default value is -1.
 593      end if
 594 !* The array occ is shared among the proc, the shift ibg is always incremented.
 595      my_ibgtab(ikpt,isppol)=ibg
 596      ibg=ibg+nband
 597    end do
 598  end do
 599 
 600  if (.not.(associated(fock))) then
 601 
 602 ! =================================
 603 ! === Create the fock structure ===
 604 ! =================================
 605    ABI_DATATYPE_ALLOCATE(fock,)
 606    ABI_DATATYPE_ALLOCATE(fock%fock_common,)
 607    ABI_DATATYPE_ALLOCATE(fock%fock_BZ,)
 608 ! ========================================================
 609 ! === Set all the other state-dependent fields to zero ===
 610 ! ========================================================
 611    fockcommon=>fock%fock_common
 612    fockbz=> fock%fock_BZ
 613    ABI_ALLOCATE(fockcommon%nband,(dtset%nkpt*dtset%nsppol))
 614    do ikpt=1,dtset%nkpt*dtset%nsppol
 615      fockcommon%nband(ikpt)=dtset%nband(ikpt)
 616    end do
 617 
 618    nband=dtset%mband
 619    fockcommon%ikpt= 0
 620 !* Will contain the k-point ikpt of the current state
 621    fockcommon%isppol= 0
 622 !* Will contain the spin isppol of the current state
 623    fockcommon%ieigen=0
 624 !* Will contain the band index of the current state
 625 !* if the value is 0, the Fock contribution to the eigenvalue is not calculated.
 626    ABI_ALLOCATE(fockcommon%eigen_ikpt,(nband))
 627    fockcommon%eigen_ikpt=0.d0
 628 !* Will contain the Fock contributions to the eigenvalue of the current state
 629 
 630 !* Compute the dimension of arrays in "spin" w.r.t parallelism
 631    my_nsppol=dtset%nsppol
 632    if (mpi_enreg%nproc_kpt>1) my_nsppol=1
 633 !* my_nsppol=1 when nsppol=1 or nsppol=2 and only one spin is treated by the processor.
 634 !* my_nsppol=2 when nsppol=2 and no parallelization over kpt (both spins are treated by the processor).
 635 
 636 !* Compute mkpt the size of arrays/pointers for k points w.r.t. parallelism
 637 !* Compute mkptband the size of arrays/pointers for occupied states w.r.t. parallelism
 638    if (nproc_hf<nkpt_bz) then
 639 !* Parallelization over kpts only
 640      mkpt=nkpt_bz/nproc_hf
 641      if (mod(nkpt_bz,nproc_hf) /=0) mkpt=mkpt+1
 642      mkptband=mkpt*mband
 643    else
 644 !* Parallelization over occupied states
 645      if (nproc_hf<nkpt_bz*mband) then
 646        mkptband=(nkpt_bz*mband)/nproc_hf
 647        if (mod((nkpt_bz*mband),nproc_hf) /=0) mkptband=mkptband+1
 648        mkpt=1
 649        if (mod(nproc_hf,nkpt_bz) /=0) mkpt=2
 650      else
 651        mkptband=1
 652        mkpt=1
 653      end if
 654    end if
 655 
 656 ! mpi_enreg settings
 657    call copy_mpi_enreg(mpi_enreg,fockbz%mpi_enreg)
 658    fockbz%mpi_enreg%me_kpt=mpi_enreg%me_hf
 659    fockbz%mpi_enreg%comm_kpt=mpi_enreg%comm_hf
 660    fockbz%mpi_enreg%nproc_kpt=mpi_enreg%nproc_hf
 661    if (allocated(fockbz%mpi_enreg%proc_distrb)) then
 662      ABI_DEALLOCATE(fockbz%mpi_enreg%proc_distrb)
 663    end if
 664    ABI_ALLOCATE(fockbz%mpi_enreg%proc_distrb,(nkpt_bz,mband,1))
 665    do jkpt=1,nkpt_bz
 666      fockbz%mpi_enreg%proc_distrb(jkpt,:,1)=fockbz%mpi_enreg%me_kpt
 667    end do
 668 
 669    mgfft=dtset%mgfft
 670    fockcommon%usepaw=dtset%usepaw
 671    if (fockcommon%usepaw==1)then
 672      mgfft=dtset%mgfftdg
 673      n4=dtset%ngfftdg(4) ; n5=dtset%ngfftdg(5) ; n6=dtset%ngfftdg(6)
 674    end if
 675    fockcommon%optfor=.false.; fockcommon%optstr=.false.
 676    if(dtset%optforces==1) fockcommon%optfor=.true.
 677    if (fockcommon%optfor) then
 678      ABI_ALLOCATE(fockcommon%forces_ikpt,(3,dtset%natom,nband))
 679      ABI_ALLOCATE(fockcommon%forces,(3,dtset%natom))
 680      fockcommon%forces=zero
 681    endif
 682    use_ACE=1 ! Default. Normal users do not have access to this variable, although the next line allows experts to make tests.
 683    if(dtset%userie==1729)use_ACE=0 ! Hidden possibility to disable ACE
 684 
 685    fockcommon%use_ACE=use_ACE
 686    call fockbz_create(fockbz,mgfft,dtset%mpw,mkpt,mkptband,my_nsppol,dtset%natom,n4,n5,n6,use_ACE)
 687 
 688 !* Initialize %mband, %mkpt, %mkptband = size of arrays
 689    fockcommon%mband=mband
 690    fockbz%mkpt=mkpt
 691    fockbz%mkptband=mkptband
 692    fockcommon%my_nsppol = my_nsppol
 693    fockcommon%nsppol = dtset%nsppol
 694    if (fockcommon%use_ACE/=0) then
 695      ABI_DATATYPE_ALLOCATE(fock%fockACE,(dtset%nkpt,dtset%nsppol))
 696      my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 697      do isppol=1,dtset%nsppol
 698        do ikpt=1,dtset%nkpt
 699          nband=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 700          ABI_ALLOCATE(fock%fockACE(ikpt,isppol)%xi,(2,npwarr(ikpt)*my_nspinor,nband))
 701        end do
 702      end do
 703    end if
 704 !========Initialze PAW data========
 705    fockcommon%ntypat=dtset%ntypat
 706    fockcommon%natom=dtset%natom
 707    if (fockcommon%usepaw==1) then
 708      fockbz%mcprj=mkptband*my_nsppol
 709      fockcommon%pawfgr => pawfgr
 710      fockbz%pawang => pawang
 711      fockcommon%pawtab => pawtab
 712      ABI_DATATYPE_ALLOCATE(fockcommon%pawfgrtab,(dtset%natom))
 713      do iatom = 1, dtset%natom
 714        itypat=dtset%typat(iatom)
 715        l_size_atm(iatom) = pawtab(itypat)%lcut_size
 716      end do
 717      call pawfgrtab_init(fockcommon%pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat)
 718      call pawfgrtab_nullify(fockcommon%pawfgrtab)
 719      ABI_DATATYPE_ALLOCATE(fockbz%cwaveocc_prj,(dtset%natom,fockbz%mcprj))
 720      ABI_ALLOCATE(dimcprj,(dtset%natom))
 721      call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
 722      ncpgr = 0
 723      if (dtset%optforces== 1) ncpgr = 3
 724 !     if (dtset%optstress /= 0) ncpgr = 6
 725 !     ncpgr=3*dtset%optforces+6*dtset%optstress
 726      call pawcprj_alloc(fockbz%cwaveocc_prj,ncpgr,dimcprj)
 727      ABI_DEALLOCATE(dimcprj)
 728      ABI_ALLOCATE(fockcommon%atindx,(dtset%natom))
 729      fockcommon%atindx=atindx
 730      ABI_ALLOCATE(fockcommon%typat,(dtset%natom))
 731      fockcommon%typat=dtset%typat
 732    end if
 733 ! ==========================================
 734 ! === Initialize the convergence options ===
 735 ! ==========================================
 736    write(msg,'(2a)') ch10,'Fock_init: initialization of Fock operator parameters:'
 737    call wrtout(std_out,msg,'COLL')
 738 
 739    fockcommon%fock_converged=.false.
 740    fockcommon%scf_converged=.false.
 741 
 742 !* Number of iterations with fixed occupied states when calculating the exact exchange contribution.
 743    if (dtset%nnsclohf<0) then
 744      msg='The parameter nnsclohf must be a non-negative integer.'
 745      MSG_ERROR(msg)
 746    end if
 747    if (dtset%nnsclohf==0) then
 748      fockcommon%nnsclo_hf=1
 749      msg=' - The parameter nnsclohf is set to its default value 1.'
 750      call wrtout(std_out,msg,'COLL')
 751 !* Default value is set to 1 (updating cgocc at each step)
 752 !* May be useful to put default to 3
 753    else
 754      fockcommon%nnsclo_hf=dtset%nnsclohf
 755      write(msg,'(a,i3)') ' - The parameter nnsclohf is set to the value:', dtset%nnsclohf
 756      call wrtout(std_out,msg,'COLL')
 757 !* value chosen by the user
 758    end if
 759 
 760 ! =========================================
 761 ! === Initialize the hybrid coefficient ===
 762 ! =========================================
 763    fockcommon%ixc = dtset%ixc
 764 !  By convention, positive values are the default values for the ixc,
 765 !  while negative values have been set by the user (and stored as negative numbers)
 766    fockcommon%hyb_mixing=abs(dtset%hyb_mixing)
 767    fockcommon%hyb_mixing_sr=abs(dtset%hyb_mixing_sr)
 768    fockcommon%hyb_range_dft=abs(dtset%hyb_range_dft)
 769    fockcommon%hyb_range_fock=abs(dtset%hyb_range_fock)
 770 
 771 !  Set the hybrid parameters if functional from libxc for which parameters can be changed, or if the user asked to do so.
 772 !  Usually, these parameters were obtained from libxc,
 773 !  but the user might have possibly modified them. By the way, must define them here for the usual changeable fonctionals,
 774 !  since otherwise might inherit them from the previous dataset !
 775    if(dtset%ixc<0)then
 776      if (dtset%ixc==-406.or.dtset%ixc==-427.or.dtset%ixc==-428 .or. &
 777 &      min(dtset%hyb_mixing,dtset%hyb_mixing_sr,dtset%hyb_range_dft,dtset%hyb_range_fock)<-tol8)then
 778        call libxc_functionals_set_hybridparams(hyb_mixing=fockcommon%hyb_mixing,&
 779 &                                              hyb_mixing_sr=fockcommon%hyb_mixing_sr,&
 780 &                                              hyb_range=fockcommon%hyb_range_dft)
 781      end if
 782    end if
 783 
 784 
 785 ! ======================================================
 786 ! === Initialize the data relative to Poisson solver ===
 787 ! ======================================================
 788 
 789 !* gsqcut = cutoff value on G^2 for sphere inside the fft box (input for vhartre).
 790    fockcommon%gsqcut= gsqcut
 791 
 792 ! =======================================================
 793 ! === Initialize the properties of the k-points in BZ ===
 794 ! =======================================================
 795 !* Initialize %nkpt_bz = nb of k point in BZ for the calculation of exchange
 796    fockbz%nkpt_bz=nkpt_bz
 797 !* Initialize the array %wtk_bz = weight assigned to each k point.
 798    fockbz%wtk_bz=1.0_dp/dble(nkpt_bz)
 799 
 800 
 801    if (dtset%kptopt>=1 .and. dtset%kptopt<=4) then
 802 ! ============================================
 803 ! === Initialize the set of k-points in BZ ===
 804 ! ============================================
 805      kptns_hf(:,1:nkpt_bz)=dtset%kptns_hf(:,1:nkpt_bz)
 806 !* kptns_hf contains the special k points obtained by the Monkhorst & Pack method, in reduced coordinates. (output)
 807 
 808 ! =======================================================
 809 ! === Compute the transformation to go from IBZ to BZ ===
 810 ! =======================================================
 811 !* Compute the reciprocal space metric.
 812      call matr3inv(rprimd,gprimd)
 813      gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
 814 
 815 !* Calculate the array indkk which describes how to get IBZ from BZ
 816 !* dksqmax=maximal value of the norm**2 of the difference between a kpt2 vector and the closest k-point found from the kptns1 set, using symmetries. (output)
 817 !* sppoldbl=1, no spin-polarisation doubling is required.
 818      timrev=1 ; if (dtset%kptopt==3 .or. dtset%kptopt==4) timrev=0
 819 !* timrev=1 if the use of time-reversal is allowed ; 0 otherwise
 820      if (dtset%kptopt==2 .or. dtset%kptopt==3) then
 821 !* No space symmetry is used, if kptopt==2 time reversal symmetry is used.
 822        symm=0 ; symm(1,1)=1 ; symm(2,2)=1 ; symm(3,3)=1
 823        call listkk(dksqmax,gmet,indkk(1:nkpt_bz,:),dtset%kptns,kptns_hf,dtset%nkpt, &
 824 &          nkpt_bz,1,1,indx,symm,timrev)
 825      else
 826 !* As in getkgrid, no use of antiferromagnetic symmetries thans to the option sppoldbl=1
 827        call listkk(dksqmax,gmet,indkk(1:nkpt_bz,:),dtset%kptns,kptns_hf,dtset%nkpt, &
 828 &          nkpt_bz,dtset%nsym,1,dtset%symafm,dtset%symrel,timrev)
 829      end if
 830 !* indkk(nkpt_bz,6) describes the k point of IBZ that generates each k point of BZ
 831 !*    indkk(:,1)   = k point of IBZ, kpt_ibz
 832 !*    indkk(:,2)   = symmetry operation to apply to kpt_ibz to give the k point of BZ
 833 !*                   (if 0, means no symmetry operation, equivalent to identity )
 834 !*    indkk(:,3:5) = Umklapp vectors to apply to remain in BZ
 835 !*    indkk(:,6)   = 1 if time-reversal was used to generate the k point of BZ, 0 otherwise
 836 !* No use of symafm to generate spin down wfs from spin up wfs for the moment
 837 
 838    else
 839      if (dtset%kptopt==0) then
 840 !* kptopt =0 : read directly nkpt, kpt, kptnrm and wtk in the input file
 841 !*              => this case is not allowed for the moment
 842        msg='Hartree-Fock option can not be used with option kptopt=0.'
 843        MSG_ERROR(msg)
 844      else
 845 !* kptopt <0 : rely on kptbounds, and ndivk to set up a band structure calculation
 846 !*              => a band structure calculation is not yet allowed.
 847        msg='Hartree-Fock option can not be used with option kptopt<0.'
 848        MSG_ERROR(msg)
 849      end if
 850    end if
 851 
 852 !! =======================================================
 853 !! === Initialize the properties of the k-points in BZ ===
 854 !! =======================================================
 855 !       jkg=0
 856 !!* Initialize the arrays %npwarr_bz, %kg_j, %phase_j, %gbound_j
 857 !       do jkpt=1,nkpt_bz
 858 !         ikpt=indkk(jkpt,1)
 859 !!* ikpt = the point of IBZ that jkpt is an image of in BZ
 860 !         npwj=npwarr(ikpt)
 861 !!* npwj = number of planewaves in basis at point jkpt = at point ikpt
 862 !         jsym=indkk(jkpt,2)
 863 !!* jsym = symmetry operation to apply to get jkpt from ikpt
 864 !         shiftg(:)=indkk(jkpt,3:5)
 865 !!* shiftg = Bravais vector G0 to add to remain in BZ
 866 !         if (jsym/=0) then
 867 !           symm(:,:)=dtset%symrel(:,:,jsym)
 868 !           tau_nons(:)=dtset%tnons(:,jsym)
 869 !!* The symmetry operation in k-space (symm) and the non-symorphic translation (tau_nons) are now defined.
 870 !           if(sum(tau_nons(:)**2)>tol8) then
 871 !!* Initialize %calc_phase(jkpt) to 1
 872 !             fock%calc_phase(jkpt)=1
 873 !!* Compute the phase factor exp(i*2*pi*G.tau) for all G.
 874 !             indx(1)=1
 875 !             phase1d=zero
 876 !             call getph(indx,1,n1,n2,n3,phase1d,tau_nons)
 877 !!* Although the routine getph is orignally written for atomic phase factors, it does precisely what we want
 878 !             arg=two_pi*(dtset%kptns(1,ikpt)*tau_nons(1) + dtset%kptns(2,ikpt)*tau_nons(2) &
 879 !&                + dtset%kptns(3,ikpt)*tau_nons(3))
 880 !             phktnons(1,1)=cos(arg)
 881 !             phktnons(2,1)=sin(arg)
 882 !!              phktnons(1,1)=one
 883 !!              phktnons(2,1)=zero
 884 !!* Convert 1D phase factors to 3D phase factors exp(i*2*pi*(k+G).tau) and store it in %phase_j
 885 !             call ph1d3d(1,1,kg(:,1+tab_indikpt(1,ikpt):npwj+tab_indikpt(1,ikpt)),1,1,npwj,n1, &
 886 !&              n2,n3,phktnons,phase1d,fock%phase(:,1+jkg:npwj+jkg))
 887 !           end if
 888 !         else
 889 !           symm=0 ; symm(1,1)=1 ; symm(2,2)=1 ; symm(3,3)=1
 890 !           tau_nons(:)=zero
 891 !           shiftg(:)=0
 892 !         end if
 893 !!* Apply time-reversal symmetry if required
 894 !         if(indkk(jkpt,6)/=0) then
 895 !!* Initialize %timerev(jkpt) to 1
 896 !           fock%timerev(jkpt)=1
 897 !           symm(:,:)=-symm(:,:)
 898 !         end if
 899 
 900 !!* Initialize %istwfk_bz(jkpt) to
 901 !         fock%istwfk_bz(jkpt)=dtset%istwfk(ikpt)
 902 
 903 !!* Initialize %tab_ikpt and %tab_ibgcg
 904 !         fock%tab_ikpt(jkpt)=ikpt
 905 !         fock%tab_ibgcg(1:dtset%nsppol,jkpt)=tab_indikpt(2:1+dtset%nsppol,ikpt)
 906 !         fock%tab_ibgcg(1+dtset%nsppol:2*dtset%nsppol,jkpt)= &
 907 !&          tab_indikpt(2+dtset%nsppol:2*dtset%nsppol+1,ikpt)
 908 
 909 !!* Initialize %npwarr_bz
 910 !         fock%npwarr_bz(jkpt)=npwj
 911 
 912 !!* Initialize %kg_bz
 913 !         do jpw=1,npwj
 914 !           v1=kg(1,jpw+tab_indikpt(1,ikpt)) ; v2=kg(2,jpw+tab_indikpt(1,ikpt)) ; v3=kg(3,jpw+tab_indikpt(1,ikpt))
 915 !           fock%kg_bz(1,jpw+jkg)=-shiftg(1)+symm(1,1)*v1+symm(2,1)*v2+symm(3,1)*v3
 916 !           fock%kg_bz(2,jpw+jkg)=-shiftg(2)+symm(1,2)*v1+symm(2,2)*v2+symm(3,2)*v3
 917 !           fock%kg_bz(3,jpw+jkg)=-shiftg(3)+symm(1,3)*v1+symm(2,3)*v2+symm(3,3)*v3
 918 !!* The symmetry operation symm must be transposed when used. (cf. docs about wfconv)
 919 !         end do
 920 
 921 !!* Initialize %gbound_bz
 922 !         call sphereboundary(fock%gbound_bz(:,:,jkpt),fock%istwfk_bz(jkpt), &
 923 !&          fock%kg_bz(:,1+jkg:npwj+jkg),dtset%mgfft,npwj)
 924 
 925 !!* Update of the shift to be applied
 926 !         jkg=jkg+npwj
 927 !       end do
 928 
 929 ! ==========================================================
 930 ! === Initialize the k-points in BZ and their properties ===
 931 ! ==========================================================
 932 !   jkg=0;
 933 
 934    do isym=1,dtset%nsym
 935          call mati3inv(dtset%symrel(:,:,isym),symrec(:,:,isym))
 936          Rtnons (:,isym)= MATMUL(TRANSPOSE(symrec(:,:,isym)),dtset%tnons(:,isym))
 937    end do
 938    ABI_ALLOCATE(fockcommon%symrec,(3,3,dtset%nsym))
 939    fockcommon%symrec=symrec
 940 
 941    ABI_ALLOCATE(invsym,(dtset%nsym))
 942    invsym=0
 943    ident(1,:3)=(/1,0,0/)
 944    ident(2,:3)=(/0,1,0/)
 945    ident(3,:3)=(/0,0,1/)
 946    do isym=1,dtset%nsym
 947      symm(:,:)=MATMUL(dtset%symrel(:,:,isym),dtset%symrel(:,:,isym))
 948      if (all(symm(:,:)==ident(:,:))) then
 949        invsym(isym)=isym
 950      else
 951        do jsym=1,dtset%nsym
 952          symm(:,:)=MATMUL(dtset%symrel(:,:,isym),dtset%symrel(:,:,jsym))
 953          if (all(symm(:,:)==ident(:,:))) then
 954             invsym(isym)=jsym
 955             cycle
 956          end if
 957        end do
 958      end if
 959      if(invsym(isym)==0) then
 960        msg='No inverse has been found for isym'
 961        MSG_ERROR(msg)
 962      end if
 963    end do
 964 
 965    jkg_this_proc=0;my_jkpt=0
 966 !indkk(1:nkpt_bz,2)=(/1,1,3,1,11,7,9,1/)
 967    do jkpt=1,nkpt_bz
 968 
 969 !* If this processor does not calculate exchange with the k point jkpt, skip the rest of the k-point loop.
 970      if (proc_distrb_cycle(mpi_enreg%distrb_hf,jkpt,1,mband,1,mpi_enreg%me_hf)) cycle
 971 !       if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,jkpt,1,dtset%nbandhf,1,mpi_enreg%me_kpt))) then
 972 !* The processor does own a copy of the array kg of ikpt ; increment the shift.
 973 !         jkg=jkg+npwj
 974 !        end if
 975 ! Skip the rest of the k-point loop
 976 !       cycle
 977 !     end if
 978      my_jkpt=my_jkpt+1
 979 
 980      ikpt=indkk(jkpt,1)
 981 !* ikpt = the point of IBZ that jkpt is an image of in BZ
 982      npwj=npwarr(ikpt)
 983      fockbz%npwarr(my_jkpt)=npwarr(ikpt)
 984 !* npwj = number of planewaves in basis at point jkpt = at point ikpt
 985      jsym=indkk(jkpt,2)
 986 !* jsym = symmetry operation to apply to get jkpt from ikpt
 987      fockbz%tab_symkpt(my_jkpt)=invsym(jsym)
 988      shiftg(:)=indkk(jkpt,3:5)
 989 !* shiftg = Bravais vector G0 to add to remain in BZ
 990 
 991 !* Initialize the array %kptns_bz = the k points in full BZ
 992      fockbz%kptns_bz(:,my_jkpt)=kptns_hf(:,jkpt)
 993 
 994 !* Initialize the array %jstwfk = how is stored the wavefunction at each k point
 995      if (dtset%istwfk(ikpt)/=1) then
 996        fockbz%istwfk_bz(my_jkpt)=set_istwfk(kptns_hf(:,jkpt))
 997      end if
 998 
 999 !* One can take advantage of the time-reversal symmetry in this case.
1000 !* Initialize the array %wtk_bz = weight assigned to each k point.
1001 !     fock%wtk_bz(my_jkpt)=dtset%wtk(jkpt)/ucvol
1002 !* Caution, the definition takes into account "ucvol" !
1003 
1004 !* Initialize the array %npwarr_bz = number of planewaves in basis at each k point
1005 !     fock%npwarr_bz(my_jkpt)=npwj
1006 
1007 !!* Initialize the array %tab_ikpt = indices of k-point in IBZ ikpt for each k point jkpt in BZ (here,ikpt=jkpt)
1008      fockbz%tab_ikpt(my_jkpt)=ikpt
1009 
1010 
1011 !!* Initialize the array %tab_ibgcg = indices of cprj(ikpt)/occ(ikpt) and cg(ikpt) for each k point jkpt
1012 !     if (my_nsppol==2) then
1013 !!* In this case, my_nsppol=dtset%nsppol=2
1014 !       fock%tab_ibgcg(1:2,my_jkpt)=tab_indikpt(2:3,ikpt)
1015 !       fock%tab_ibgcg(3:4,my_jkpt)=tab_indikpt(4:5,ikpt)
1016 !     else
1017 !       if(mpi_enreg%my_isppoltab(1)==1) then
1018 !!* In this case, my_nsppol=1 and the up spin is treated (dtset%nsppol= 1 or 2)
1019 !         fock%tab_ibgcg(1,my_jkpt)=tab_indikpt(2,ikpt)
1020 !         fock%tab_ibgcg(2,my_jkpt)=tab_indikpt(2+dtset%nsppol,ikpt)
1021 !       else
1022 !!* In this case, my_nsppol=1 and the dn spin is treated (so dtset%nsppol=2)
1023 !         fock%tab_ibgcg(1,my_jkpt)=tab_indikpt(3,ikpt)
1024 !         fock%tab_ibgcg(2,my_jkpt)=tab_indikpt(5,ikpt)
1025 !       end if
1026 !     end if
1027 
1028 !* Initialize the array %kg_bz = reduced planewave coordinates at each k point
1029      if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me_kpt))) then
1030 !* We perform the test with isppol=-1 (both spins) and the occupied band (dtset%nbandhf).
1031 !* We assume that paral_kgb==0 (a k-point may not be present on several proc.)
1032 !* The array kg for ikpt is stored on this processor and copied in kg_tmp.
1033        ikg=my_ikgtab(ikpt)
1034 !* ikg = the shift in kg to get the G-vectors associated to ikpt
1035        do ik=1,3
1036 !         kg_tmp(1+(ik-1)*npwj:ik*npwj)=kg(ik,1+tab_indikpt(1,ikpt):npwj+tab_indikpt(1,ikpt))
1037          kg_tmp(1+(ik-1)*npwj:ik*npwj)=kg(ik,1+ikg:npwj+ikg)
1038        end do
1039 !       jkg=jkg+npwj
1040      end if
1041 !* Broadcast the array kg_tmp to all the processors of comm_kpt.
1042 !* Since paral_kgb==0, all the bands of a k-point are treated on the same proc.
1043      call xmpi_bcast(kg_tmp,mpi_enreg%proc_distrb(ikpt,1,1),mpi_enreg%comm_kpt,ier)
1044      do ik=1,3
1045        fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)=kg_tmp(1+(ik-1)*npwj:ik*npwj)
1046      end do
1047 
1048 !* Apply a symmetry operation on kg_bz if necessary
1049      if (jsym/=0) then
1050        symm(:,:)=dtset%symrel(:,:,jsym)
1051 !      tau_nons(:)=dtset%tnons(:,jsym)
1052        tau_nons(:)=-Rtnons(:,invsym(jsym))
1053 !* The symmetry operation in k-space (symm) and the non-symorphic translation (tau_nons) are now defined.
1054        if(sum(tau_nons(:)**2)>tol8) then
1055 !* Initialize %calc_phase(jkpt) to 1
1056          fockbz%calc_phase(my_jkpt)=1
1057 !* Compute the phase factor exp(i*2*pi*G.tau) for all G.
1058          indx(1)=1
1059          phase1d=zero
1060          call getph(indx,1,n1,n2,n3,phase1d,tau_nons)
1061 !* Although the routine getph is orignally written for atomic phase factors, it does precisely what we want
1062          arg=two_pi*(dtset%kptns(1,ikpt)*tau_nons(1) + dtset%kptns(2,ikpt)*tau_nons(2) &
1063 &            + dtset%kptns(3,ikpt)*tau_nons(3))
1064          phktnons(1,1)=cos(arg)
1065          phktnons(2,1)=sin(arg)
1066 !          phktnons(1,1)=one
1067 !          phktnons(2,1)=zero
1068 !* Convert 1D phase factors to 3D phase factors exp(i*2*pi*(k+G).tau) and store it in %phase_j
1069          call ph1d3d(1,1,fockbz%kg_bz(:,1+jkg_this_proc:npwj+jkg_this_proc),1,1,npwj,n1,n2,n3, &
1070 &          phktnons,phase1d,fockbz%phase(:,1+jkg_this_proc:npwj+jkg_this_proc))
1071        end if
1072 !* Apply time-reversal symmetry if required
1073        if(indkk(jkpt,6)/=0) then
1074 !* Initialize %timerev(jkpt) to 1
1075          fockbz%timerev(my_jkpt)=1
1076          symm(:,:)=-symm(:,:)
1077        end if
1078 !* Initialize %kg_bz
1079        do jpw=1,npwj
1080          v1=fockbz%kg_bz(1,jpw+jkg_this_proc) ; v2=fockbz%kg_bz(2,jpw+jkg_this_proc) ; v3=fockbz%kg_bz(3,jpw+jkg_this_proc)
1081          fockbz%kg_bz(1,jpw+jkg_this_proc)=-shiftg(1)+symm(1,1)*v1+symm(2,1)*v2+symm(3,1)*v3
1082          fockbz%kg_bz(2,jpw+jkg_this_proc)=-shiftg(2)+symm(1,2)*v1+symm(2,2)*v2+symm(3,2)*v3
1083          fockbz%kg_bz(3,jpw+jkg_this_proc)=-shiftg(3)+symm(1,3)*v1+symm(2,3)*v2+symm(3,3)*v3
1084 !* The symmetry operation symm must be transposed when used. (cf. docs about wfconv)
1085        end do
1086      else
1087 !* Ths symmetry operation is the identity.
1088 !* Apply time-reversal symmetry if required
1089        if(indkk(jkpt,6)/=0) then
1090 !* Initialize %timerev(jkpt) to 1
1091          fockbz%timerev(my_jkpt)=1
1092          fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)=-fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)
1093        end if
1094      end if
1095 
1096 !* Initialize the array %gbound_bz = boundary of the basis sphere of G vectors at each k point
1097      call sphereboundary(fockbz%gbound_bz(:,:,my_jkpt),fockbz%istwfk_bz(my_jkpt),&
1098 &      fockbz%kg_bz(:,1+jkg_this_proc:npwj+jkg_this_proc),mgfft,npwj)
1099 
1100      jkg_this_proc=jkg_this_proc+npwj
1101 
1102 !* Initialize the arrays %tab_ibg = shifts in arrays cprj and occ (ibg) for each k point jkpt
1103 !* Initialize the arrays %tab_icg = shifts in arrays cg(icg) for each k point jkpt
1104      if (my_nsppol==1) then
1105          fockbz%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,1+mpi_enreg%my_isppoltab(2))
1106          fockbz%tab_icg(my_jkpt,1)=my_icgtab(ikpt,1+mpi_enreg%my_isppoltab(2))
1107          fockbz%tab_icp(my_jkpt,1)=my_icptab(ikpt,1+mpi_enreg%my_isppoltab(2))
1108 !* if mpy_isppoltab(2)=0, the up spin is treated (dtset%nsppol= 1 or 2)
1109 !* if mpy_isppoltab(2)=1, the dn spin is treated (so dtset%nsppol=2)
1110 
1111 !       if(mpi_enreg%my_isppoltab(2)==1) then
1112 !* In this case, my_nsppol=1 and the up spin is treated (dtset%nsppol= 1 or 2)
1113 !         fock%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,1)
1114 !         fock%tab_icg(my_jkpt,1)=my_icgtab(ikpt,1)
1115 !       else
1116 !* In this case, my_nsppol=1 and the dn spin is treated (so dtset%nsppol=2)
1117 !         fock%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,2)
1118 !         fock%tab_icg(my_jkpt,1)=my_icgtab(ikpt,2)
1119 !       end if
1120      else
1121 !* In this case, my_nsppol=dtset%nsppol=2
1122        fockbz%tab_ibg(my_jkpt,:)=my_ibgtab(ikpt,:)
1123        fockbz%tab_icg(my_jkpt,:)=my_icgtab(ikpt,:)
1124        fockbz%tab_icp(my_jkpt,:)=my_icptab(ikpt,:)
1125      end if
1126 
1127    enddo
1128 
1129 !* Deallocation
1130    ABI_DEALLOCATE(invsym)
1131 
1132  end if
1133  ABI_DEALLOCATE(indkk)
1134  ABI_DEALLOCATE(kg_tmp)
1135  ABI_DEALLOCATE(kptns_hf)
1136  ABI_DEALLOCATE(my_ibgtab)
1137  ABI_DEALLOCATE(my_icgtab)
1138  ABI_DEALLOCATE(my_icptab)
1139  ABI_DEALLOCATE(my_ikgtab)
1140  ABI_DEALLOCATE(phase1d)
1141  call fock_print(fockcommon,fockbz,unit=std_out)
1142 
1143  call timab(1500,2,tsec)
1144 
1145 DBG_EXIT("COLL")
1146 
1147 end subroutine fock_init

m_fock/fock_print [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_print

FUNCTION

  Print info on the fock_type data type

INPUTS

  fock<crystal_t>=The object
  [unit]=Unit number for output
  [prtvol]=Verbosity level
  [mode_paral]=Either "COLL" or "PERS"
  [header]=String to be printed as header for additional info.

OUTPUT

  Only printing

PARENTS

      m_fock

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

2165 subroutine fock_print(fockcommon,fockbz,header,unit,mode_paral,prtvol)
2166 
2167 
2168 !This section has been created automatically by the script Abilint (TD).
2169 !Do not modify the following lines by hand.
2170 #undef ABI_FUNC
2171 #define ABI_FUNC 'fock_print'
2172 !End of the abilint section
2173 
2174  implicit none
2175 
2176 !Arguments ------------------------------------
2177 !scalars
2178  integer,optional,intent(in) :: unit,prtvol
2179  character(len=4),optional,intent(in) :: mode_paral
2180  character(len=*),optional,intent(in) :: header
2181  type(fock_common_type),intent(in) :: fockcommon
2182  type(fock_BZ_type),intent(in) :: fockbz
2183 !Local variables-------------------------------
2184  integer :: my_unt,my_prtvol
2185  character(len=4) :: my_mode
2186  character(len=500) :: msg
2187 
2188 ! *********************************************************************
2189 
2190  my_unt=std_out; if (PRESENT(unit)) my_unt=unit
2191  my_prtvol=0 ; if (PRESENT(prtvol)) my_prtvol=prtvol
2192  my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral
2193 
2194  msg=' ==== Info on fock_type ==== '
2195  if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
2196  call wrtout(my_unt,msg,my_mode)
2197 
2198  ! Important dimensions
2199  call wrtout(my_unt,sjoin(" my_nsppol ...",itoa(fockcommon%my_nsppol)),my_mode)
2200  call wrtout(my_unt,sjoin(" nkpt_bz .....",itoa(fockbz%nkpt_bz)),my_mode)
2201 
2202  ! Options
2203  call wrtout(my_unt,sjoin(" nnsclo_hf .......",itoa(fockcommon%nnsclo_hf)),my_mode)
2204  call wrtout(my_unt,sjoin(" ixc .............",itoa(fockcommon%ixc)),my_mode)
2205  call wrtout(my_unt,sjoin(" hybrid mixing....",ftoa(fockcommon%hyb_mixing)),my_mode)
2206  call wrtout(my_unt,sjoin(" hybrid SR mixing ",ftoa(fockcommon%hyb_mixing_sr)),my_mode)
2207  call wrtout(my_unt,sjoin(" hybrid range DFT ",ftoa(fockcommon%hyb_range_dft)),my_mode)
2208  call wrtout(my_unt,sjoin(" hybrid range Fock",ftoa(fockcommon%hyb_range_fock)),my_mode)
2209 
2210 ! write(msg,"(a,f12.1,a)")" Memory required for HF u(r) states: ",product(shape(fockbz%cwaveocc_bz)) * dp * b2Mb, " [Mb]"
2211 ! call wrtout(my_unt,msg,my_mode)
2212 
2213  ! Extra info.
2214  if (my_prtvol > 0) then
2215    call wrtout(my_unt,"Extra info not available",my_mode)
2216  end if
2217 
2218 end subroutine fock_print

m_fock/fock_set_getghc_call [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_set_getghc_call

FUNCTION

  Set the value of fock%getghc_call, Returns the old value

PARENTS

SOURCE

2079 integer function fock_set_getghc_call(fock, new) result(old)
2080 
2081 
2082 !This section has been created automatically by the script Abilint (TD).
2083 !Do not modify the following lines by hand.
2084 #undef ABI_FUNC
2085 #define ABI_FUNC 'fock_set_getghc_call'
2086 !End of the abilint section
2087 
2088  implicit none
2089 
2090 !Arguments ------------------------------------
2091 !scalars
2092  integer,intent(in) :: new
2093  type(fock_common_type),intent(inout) :: fock
2094 
2095 ! *************************************************************************
2096 
2097  old = fock%getghc_call_
2098  fock%getghc_call_ = new
2099 
2100 end function fock_set_getghc_call

m_fock/fock_set_ieigen [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_set_ieigen

FUNCTION

  Set the value of ieigen to the value given in argument.

INPUTS

  fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
  iband= index of the band iband

OUTPUT

  none

PARENTS

      cgwf

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1238 subroutine fock_set_ieigen(fock,iband)
1239 
1240 
1241 !This section has been created automatically by the script Abilint (TD).
1242 !Do not modify the following lines by hand.
1243 #undef ABI_FUNC
1244 #define ABI_FUNC 'fock_set_ieigen'
1245 !End of the abilint section
1246 
1247  implicit none
1248 
1249 !Arguments ------------------------------------
1250  integer, intent(in) :: iband
1251  type(fock_common_type),pointer :: fock
1252 
1253 ! *************************************************************************
1254 
1255 !Nothing to do if fock pointer is not associated...
1256 
1257 ! ======================================================
1258 ! === Update the data relative to the current states ===
1259 ! ======================================================
1260 
1261 !* Copy of the value iband in the field ieigen
1262  if (associated(fock)) then
1263    fock%ieigen=iband
1264    fock%iband=iband
1265  end if
1266 
1267 end subroutine fock_set_ieigen

m_fock/fock_type [ Types ]

[ Top ] [ m_fock ] [ Types ]

NAME

  fock_type

FUNCTION

   This object stores the occupied wavefunctions and other quantities
   needed to calculate Fock exact exchange

SOURCE

 69  type, public :: fock_type
 70   type(fock_common_type), pointer :: fock_common=> null()
 71   type(fock_BZ_type), pointer :: fock_BZ=> null()
 72   type(fock_ACE_type), pointer :: fockACE(:,:)=> null()
 73  end type fock_type
 74 
 75 
 76  type, public :: fock_common_type
 77 
 78 ! Integer scalars
 79   !integer :: mcgocc_bz,mkg_bz,mocc
 80   !integer :: natom,ntypat
 81 
 82   integer :: usepaw
 83     ! 0 if norm-conserving psps, 1 for PAW (not implemented)
 84 
 85   integer :: ikpt,isppol,ieigen,iband
 86     ! data relative to the current states.
 87 
 88   integer :: mband
 89     ! maximum number of bands
 90 
 91   integer :: my_nsppol
 92    ! my_nsppol=1 when nsppol=1 or nsppol=2 and only one spin is treated by the processor.
 93    ! my_nsppol=2 when nsppol=2 and no parallelization over kpt (both spins are treated by the processor).
 94 
 95   integer :: natom
 96    ! Number of atoms, input variable
 97 
 98   integer :: nsppol
 99    ! Number of independent spin polarizations, input variable
100    ! Note that this value does not take into account the MPI distribution of the wavefunctions.
101 
102   integer :: ntypat
103    ! Number of type of atoms
104 
105   integer :: nnsclo_hf
106     ! Number of iterations with fixed occupied states when calculating the exact exchange contribution.
107 
108   integer :: ixc
109     ! XC option (abinit input variable)
110 
111   integer :: use_ACE
112     ! option to use the ACE method of Lin Lin
113     !==0 if the normal Fock operator is to be created and/or used
114     !==1 if the ACE operator is to be created and/or used
115 
116   integer ABI_PRIVATE :: getghc_call_ = 1
117   ! 1 if fock_getghc should be called in getghc, 0 otherwise
118 
119 ! Logical
120   logical :: optfor
121     ! option to calculate forces
122 
123   logical :: optstr
124     ! option to calculate stresses
125 
126   logical :: fock_converged
127     ! .false. if the Fock cycle (with changing Fock/ACE operator) is not converged
128     ! .true. if the Fock cycle (with changing Fock/ACE operator) has converged
129 
130   logical :: scf_converged
131     ! .false. if the SCF cycle (with fixed Fock/ACE operator) is not converged
132     ! .true. if the SCF cycle (with fixed Fock/ACE operator) has converged
133 
134 ! Real(dp) scalars
135 
136   real(dp) :: gsqcut
137     !  cutoff value on G**2 for sphere inside fft box.
138     !   (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)). Used in hartre
139 
140   real(dp) :: hyb_mixing
141     ! hybrid mixing coefficient for the Fock contribution
142 
143   real(dp) :: hyb_mixing_sr
144     ! hybrid mixing coefficient for the short-range Fock contribution
145 
146   real(dp) :: hyb_range_dft
147     ! hybrid range for separation, used in the DFT functional
148     ! (should be equal to hyb_range_fock, but this is not true for HSE03)
149 
150   real(dp) :: hyb_range_fock
151     ! hybrid range for separation, used in the fock contribution
152 
153   integer, allocatable :: atindx(:)
154     !  atindx(natom)=index table for atoms (see gstate.f)
155 
156   integer, allocatable  :: nband(:)
157    ! nband(nkpt)
158    ! Number of bands for each k point
159 
160   integer, allocatable :: symrec(:,:,:)
161 
162   integer,allocatable :: typat(:)
163    ! typat(natom)
164    ! type of each atom
165 
166 ! Real(dp) arrays
167   real(dp) :: stress(6)
168     ! stress(6)
169     ! contribution of the fock term to stresses
170 
171   real(dp), allocatable  :: stress_ikpt(:,:)
172     ! stress(6,nband)
173     ! contribution of the fock term to stresses for the current band
174 
175   real(dp), allocatable :: forces_ikpt(:,:,:)
176     ! forces(3,natom,nband))
177     ! contribution of the fock term to forces for the current band
178 
179   real(dp), allocatable :: forces(:,:)
180     ! forces(3,natom))
181     ! contribution of the fock term to forces
182 
183   real(dp), allocatable :: eigen_ikpt(:)
184     ! eigen_ikpt,(nband))
185     !  Will contain the band index of the current state
186     !  if the value is 0, the Fock contribution to the eigenvalue is not calculated.
187 
188 ! Pointers to PAW-types (associated only if usepaw==1)
189 ! Note that these are references to already existing objects.
190 
191   type(pawtab_type), pointer :: pawtab(:)
192   type(pawfgr_type),pointer :: pawfgr
193   type(pawfgrtab_type),allocatable :: pawfgrtab(:)
194 
195  end type fock_common_type
196 
197  type, public :: fock_BZ_type
198 
199   integer :: mcprj
200     ! dimension of cwaveocc_cprj
201 
202   integer :: mkpt
203     ! maximum number of k-points for Fock treated by this node
204 
205   integer :: mkptband
206     ! size of occupied states stored by this node.
207 
208   integer :: nkpt_bz
209     ! Number of k-points in the BZ for Fock operator
210 
211   integer, allocatable   :: gbound_bz(:,:,:)
212     ! gbound_bz(2*mgfft+8,2,mkpt)
213     ! Tables for zero-padded FFT of wavefunctions.
214 
215   integer, allocatable :: kg_bz(:,:)
216     ! kg_bz(3,mpw*mkpt)
217     ! G-vectors for each k-point in the BZ treate by this node
218 
219   integer, allocatable :: nbandocc_bz(:,:)
220     ! nbandocc_bz,(mkpt,my_nsppol))
221     ! nb of bands at each k point
222 
223   integer, allocatable :: npwarr(:)
224     ! npwarr(mkpt)
225 
226   integer, allocatable :: istwfk_bz(:)
227     ! istwfk_bz,(mkpt))
228     ! storage mode of the wavefunction at each k-point
229 
230   integer, allocatable :: calc_phase(:)
231     ! calc_phase,(mkpt))
232     ! 1 if a phase factor must be considered (0 otherwise) at each k point
233 
234   integer, allocatable :: tab_symkpt(:)
235     ! tab_symkpt,(mkpt))
236     ! indices of symmetry operation to apply to get jkpt in full BZ from ikpt in IBZ
237 
238   integer, allocatable :: timerev(:)
239     ! timerev,(mkpt))
240     ! 1 if time reversal symmetry must be used (0 otherwise) at each k point
241 
242   integer, allocatable :: tab_ibg(:,:)
243     ! tab_ibg,(mkpt,my_nsppol))
244     ! indices of cprj(ikpt)/occ(ikpt) in the arrays cprj/occ for each k-point jkpt
245 
246   integer, allocatable :: tab_icg(:,:)
247     ! tab_icg,(mkpt,my_nsppol))
248     ! indices of cg(ikpt) in the arrays cg for each k-point jkpt
249 
250   integer, allocatable :: tab_icp(:,:)
251     ! tab_icg,(mkpt,my_nsppol))
252     ! indices of cprj(ikpt) in the arrays cprj for each k-point jkpt
253 
254 
255   integer, allocatable :: tab_ikpt(:)
256     ! tab_ikpt,(mkpt))
257     ! indices of k-point ikpt in IBZ which corresponds to each k-point jkpt in full BZ
258 
259   real(dp), allocatable :: cgocc(:,:,:)
260     ! cgocc(2,npw*mkptband,my_nsppol)
261     ! wavefunction in the G-space
262 
263   real(dp), allocatable :: cwaveocc_bz(:,:,:,:,:,:)
264     ! (2,n4,n5,n6,mkptband,my_nsppol))
265     ! occupied states of each bands at each k point (used to construct Fock operator), in the real space
266   real(dp), allocatable :: occ_bz(:,:)
267     ! occ_bz(mkptband,my_nsppol))
268     ! occupancy of each bands at each k point
269 
270   real(dp), allocatable :: wtk_bz(:)
271     ! wtk_bz,(mkpt))
272     ! weights assigned to each k point in the BZ
273     ! Caution, the definition takes into account "ucvol" !
274 
275   real(dp), allocatable :: kptns_bz(:,:)
276     ! kptns_bz(3,mkpt)
277     ! k-points in full BZ
278 
279   real(dp), allocatable :: phase(:,:)
280     ! phase(2,mpw*mkpt))
281     ! phase factor the cg array will be multiplied with at each k point
282 
283   type(MPI_type) :: mpi_enreg
284   type(pawang_type),pointer :: pawang
285   type(pawcprj_type), allocatable :: cwaveocc_prj(:,:)
286  end type fock_BZ_type
287 !----------------------------------------------------------------------
288 
289  type,public :: fock_ACE_type
290 
291 ! ===== Real pointers
292   real(dp), allocatable :: xi(:,:,:)
293 
294  end type fock_ACE_type
295 !----------------------------------------------------------------------
296 
297  public :: fock_init                  ! Initialize the object.
298  public :: fock_set_ieigen            ! Set the value of ieigen to the value given in argument.
299  public :: fock_updateikpt            ! Update the value of energies%e_xc and energies%e_xcdc with Fock contribution.
300  public :: fock_destroy               ! Free memory.
301  public :: fock_ACE_destroy           ! Free memory.
302  public :: fock_common_destroy        ! Free memory.
303  public :: fock_bz_destroy            ! Free memory.
304  public :: fock_calc_ene              ! Calculate the Fock contribution to the total energy.
305  public :: fock_update_exc            ! Update the value of energies%e_xc and energies%e_xcdc with Fock contribution.
306  public :: fock_updatecwaveocc        ! Update in the fock datastructure the fields relative to the occupied states.
307  public :: fock_set_getghc_call       ! Enable/disable the call to fock_getghc in getghc.
308  public :: fock_get_getghc_call       ! Return the value of the flag used to enable/disable the call to fock_getghc in getghc.
309  public :: fock_print                 ! Print info on the object.

m_fock/fock_update_exc [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_update_exc

FUNCTION

  Update the value of energies%e_xc and energies%e_xcdc with Fock contribution

INPUTS

OUTPUT

  none

  energies <type(energies_type)>=storage for energies computed here :
   | e_fock= Fock contribution to the total energy (Hartree)

NOTES

   If the cgocc_bz are not updated at each iteration, be careful to calculate Fock energy at the same frequency.

PARENTS

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1629 subroutine fock_update_exc(fock_energy,xc_energy,xcdc_energy)
1630 
1631 
1632 !This section has been created automatically by the script Abilint (TD).
1633 !Do not modify the following lines by hand.
1634 #undef ABI_FUNC
1635 #define ABI_FUNC 'fock_update_exc'
1636 !End of the abilint section
1637 
1638  implicit none
1639 
1640 !Arguments ------------------------------------
1641  real(dp),intent(in) :: fock_energy
1642  real(dp),intent(inout) :: xc_energy,xcdc_energy
1643 
1644 ! *************************************************************************
1645 
1646 !xc_energy = fock%hyb_mixing*fock_energy
1647 !xcdc_energy = two*fock%hyb_mixing*fock_energy
1648  xc_energy =  fock_energy
1649  xcdc_energy = two*fock_energy
1650 !CMartins : For an atom, ewald should be set to zero (at the beginning of the loop) and
1651 !the contribution in !|q+G|=0 should be an approximation to the missing component of Vloc in G=0
1652 !energies%e_ewald=energies%e_ewald-half*fock%divgq0*fock%wtk_bz(1)*piinv
1653 
1654 end subroutine fock_update_exc

m_fock/fock_updatecwaveocc [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_updatecwaveocc

FUNCTION

  Update in the fock datastructure the fields relative to the occupied states.

INPUTS

  cg(2,mcg)= Input wavefunctions
  cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
  dtset <type(dataset_type)>=all input variables for this dataset
  fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
  indsym(4,nsym,natom) :: 1:3 shift, and 4 final atom, of symmetry isym operating on iatom
                            (S^{-1}(R - t) = r0 + L, see symatm.F90
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mpi_enreg=information about MPI parallelization
  npwarr(nkpt)=number of planewaves in basis at this k point
  occ(mband*nkpt*nsppol)= occupation number for each band (often 2) at each k point
  ucvol= unit cell volume ($\textrm{bohr}^{3}$)

OUTPUT

SIDE EFFECTS

   The field fock%cgocc_bz contains the table cg at the end.
   The fields kg_bz, occ_bz and fock%cwaveocc_prj are simultaneously updated.

NOTES

  ############################
  ### Not fully tested yet ###
  ############################

 May be improved by selecting only the occupied states with the same spin isppol.

PARENTS

      scfcv

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1702 subroutine fock_updatecwaveocc(cg,cprj,dtset,fock,indsym,mcg,mcprj,&
1703 &                              mpi_enreg,nattyp,npwarr,occ,ucvol)
1704 
1705 
1706 !This section has been created automatically by the script Abilint (TD).
1707 !Do not modify the following lines by hand.
1708 #undef ABI_FUNC
1709 #define ABI_FUNC 'fock_updatecwaveocc'
1710 !End of the abilint section
1711 
1712  implicit none
1713 
1714 !scalars
1715  integer, intent(in) :: mcg,mcprj
1716  real(dp), intent(in) :: ucvol
1717  type(dataset_type),intent(in) :: dtset
1718  type(fock_type),intent(inout),pointer :: fock
1719  type(MPI_type),intent(in) :: mpi_enreg
1720 !arrays
1721  integer, intent(in) :: indsym(4,dtset%nsym,dtset%natom),nattyp(dtset%ntypat),npwarr(dtset%nkpt)
1722  real(dp),intent(in) :: cg(2,mcg),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
1723  type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj)
1724 
1725 !Local variables-------------------------------
1726 !scalars
1727  integer,parameter :: tim_fourwf0=0
1728  integer :: iatm,iatom,iband,iband0,iband_cprj,ibg,icg,icp,ier,ikpt,ilmn,isize,ispinor,isppol,itypat,jbg,jcg,jkg,jkpt,jpw,jstwfk!,ii1,ii2
1729  integer :: lmnmax,mband,mband0,mgfft,mkpt,mpw,my_jsppol,my_jband,my_jkpt
1730  integer :: nband,ncpgr,n4,n5,n6,nkpt_bz,npwj,nsppol,nspinor
1731  real(dp),parameter :: weight1=one
1732  real(dp) :: cgre,cgim,invucvol
1733  character(len=500) :: message
1734 ! arrays
1735  integer :: ngfft(18)
1736  integer, ABI_CONTIGUOUS pointer :: gbound_k(:,:),kg_k(:,:)
1737  integer,allocatable :: dimlmn(:),indlmn(:,:,:),indsym_(:,:,:),typat_srt(:)
1738  real(dp) :: tsec(2),tsec2(2),dcp(3)
1739  real(dp),allocatable :: cgocc_tmp(:),cgocc(:,:),dummytab2(:,:),dummytab3(:,:,:),phase_jkpt(:,:)
1740  type(pawcprj_type),allocatable :: cprj_tmp(:,:)
1741  type(fock_common_type),pointer :: fockcommon
1742  type(fock_BZ_type),pointer :: fockbz
1743 ! *************************************************************************
1744 
1745  call timab(1502,1,tsec)
1746 
1747  ABI_CHECK(associated(fock),"fock must be associated")
1748 
1749  if (associated(fock)) then
1750 
1751    fockcommon=>fock%fock_common
1752    fockbz=> fock%fock_BZ
1753 
1754      invucvol=1.d0/sqrt(ucvol)
1755 ! Local variables = useful dimensions
1756      mband=fockcommon%mband
1757      mkpt=fockbz%mkpt
1758      mpw=dtset%mpw
1759      mgfft=dtset%mgfft
1760      ngfft=dtset%ngfft
1761      nkpt_bz=fockbz%nkpt_bz
1762      nsppol=dtset%nsppol
1763      nspinor=1
1764      ncpgr=0
1765 
1766 ! Local variables : useful arrays
1767      ABI_ALLOCATE(cgocc,(2,mpw))
1768      cgocc=zero
1769      ABI_ALLOCATE(cgocc_tmp,(2*mpw+1))
1770      cgocc_tmp=zero
1771      if (fockcommon%usepaw==1) then
1772        mgfft=dtset%mgfftdg
1773        ngfft=dtset%ngfftdg
1774        ABI_DATATYPE_ALLOCATE(cprj_tmp,(dtset%natom,nspinor))
1775        ABI_ALLOCATE(dimlmn,(dtset%natom))
1776        call pawcprj_getdim(dimlmn,dtset%natom,nattyp,dtset%ntypat,dtset%typat,fockcommon%pawtab,"O")
1777        ncpgr = 0
1778        if (dtset%optforces== 1) ncpgr = 3
1779 !       if (dtset%optstress /= 0) ncpgr = 6
1780        call pawcprj_alloc(cprj_tmp,ncpgr,dimlmn)
1781 
1782        lmnmax=maxval(fockcommon%pawtab(:)%lmn_size)
1783        ABI_ALLOCATE(indlmn,(6,lmnmax,dtset%ntypat))
1784        do itypat=1,dtset%ntypat
1785          isize=size(fockcommon%pawtab(itypat)%indlmn,2)
1786          indlmn(:,1:isize,itypat)=fockcommon%pawtab(itypat)%indlmn(:,1:isize)
1787        end do
1788        ABI_ALLOCATE(indsym_,(4,dtset%nsym,dtset%natom))
1789        ABI_ALLOCATE(typat_srt,(dtset%natom))
1790 
1791        if (dtset%nsym==1) then
1792          indsym_=0
1793          do iatom=1,dtset%natom
1794            iatm=fockcommon%atindx(iatom)
1795            typat_srt(iatm)=dtset%typat(iatom)
1796            indsym_(4,:,iatom)=iatom
1797          end do
1798        else
1799          do iatom=1,dtset%natom
1800            iatm=fockcommon%atindx(iatom)
1801            typat_srt(iatm)=dtset%typat(iatom)
1802            indsym_(1:3,:,iatm)=indsym(1:3,:,iatom)
1803            indsym_(4,:,iatm)=fockcommon%atindx(indsym(4,:,iatom))
1804        end do
1805        end if
1806      end if
1807 
1808 ! Local variables to perform FFT
1809      n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
1810      ABI_ALLOCATE(dummytab3,(n4,n5,n6))
1811 
1812 
1813      if(ANY(fockbz%calc_phase(:)/=0)) then
1814        ABI_ALLOCATE(phase_jkpt,(2,mpw))
1815        phase_jkpt=zero
1816      end if
1817 
1818 
1819 ! =======================================================
1820 ! === Update the data relative to the occupied states ===
1821 ! =======================================================
1822 !* The arrays cgocc_bz, kg_bz, occ_bz and npwarr_bz are already allocated with the maximal size.
1823 !     if ((dtset%kptopt>=1).and.(dtset%kptopt<=4)) then
1824 !       if (dtset%kptopt/=3) then
1825 
1826      do isppol=1,nsppol
1827        jbg=0 ; jcg=0 ; jkg=0 ; icp=0
1828        my_jsppol=isppol
1829        if ((isppol==2).and.(mpi_enreg%nproc_kpt/=1)) my_jsppol=1
1830 !* Both spins are treated on the same proc., only in the case where nproc_kpt=1;
1831 !* otherwise each proc. treats only one spin.
1832 
1833        ! MG: This loop is not effient!
1834        ! Ok the number of k-points in the BZ is usually `small` when hybrids are used
1835        ! but what happens if we use a 12x12x12.
1836        ! One should loop over the IBZ, broadcast and reconstruct the star of the k-point.
1837        my_jkpt=0
1838        do jkpt=1,nkpt_bz
1839 
1840          if (proc_distrb_cycle(mpi_enreg%distrb_hf,jkpt,1,mband,1,mpi_enreg%me_hf)) cycle
1841 
1842 !* In this case, the processor does not calculate the exchange with any occupied state on jkpt.
1843 
1844 !               if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,jkpt,1,dtset%nbandhf,jsppol,mpi_enreg%me_kpt))) then
1845 !* The state (jkpt,jband,jsppol) is stored in the array cg of this processor and copied in cgocc_tmp.
1846 !                 icg=icg+dtset%nband(jkpt)*npwj
1847 !               end if
1848 !               ibg=ibg+dtset%nband(jkpt)
1849 !               cycle
1850 !             end if
1851          my_jkpt=my_jkpt+1
1852 
1853          ikpt=fockbz%tab_ikpt(my_jkpt)
1854 
1855 !* ikpt = the point of IBZ that jkpt is an image of in BZ
1856          npwj=npwarr(ikpt)
1857 !* npwj= number of plane wave in basis for the wavefunction
1858          jstwfk=fockbz%istwfk_bz(my_jkpt)
1859 !* jstwfk= how is stored the wavefunction
1860          ibg=fockbz%tab_ibg(my_jkpt,my_jsppol)
1861 !* ibg = shift to be applied on the location of data in the array occ
1862          icg=fockbz%tab_icg(my_jkpt,my_jsppol)
1863 !* icg = shift to be applied on the location of data in the array cg
1864          icp=fockbz%tab_icp(my_jkpt,my_jsppol)
1865 !* icp = shift to be applied on the location of data in the array cprj
1866          gbound_k => fockbz%gbound_bz(:,:,my_jkpt)
1867 !* boundary of the basis sphere of G vectors
1868          kg_k => fockbz%kg_bz(:,1+jkg:npwj+jkg)
1869 !* reduced plean wave coordinates
1870          if (fockbz%calc_phase(my_jkpt)==1) then
1871            phase_jkpt(:,1:npwj)=fockbz%phase(:,1+jkg:npwj+jkg)
1872          end if
1873 !* phase factor at k-point j
1874 
1875 !* Initialize the band counter
1876          my_jband=0
1877          do iband=1,dtset%nband(ikpt)
1878 
1879 
1880            cgocc_tmp=zero
1881            if (fockcommon%usepaw==1) then
1882              call pawcprj_set_zero(cprj_tmp)
1883            end if
1884 
1885            if(ABS(occ(iband+ibg))>tol8) then
1886 !* If the band is occupied
1887 
1888 !* To avoid segmentation fault, my_jband should not be greater than nbandhf
1889              if ((my_jband+1)>mband) then
1890                write(message,*) 'The number of occupied band',my_jband+1,' at k-point',&
1891 &                 ikpt,' is greater than the value of nbandhf ', mband
1892                MSG_ERROR(message)
1893              end if
1894 
1895 !* If the processor does not calculate the exchange with the occupied state (jkpt,my_jband), cycle
1896 !             if (mpi_enreg%distrb_hf(jkpt,(my_jband+1),1)/=mpi_enreg%me_hf) cycle
1897              if (mpi_enreg%distrb_hf(jkpt,iband,1)/=mpi_enreg%me_hf) cycle
1898 !                   if (mpi_enreg%proc_distrb(jkpt,jband,jsppol)==mpi_enreg%me_kpt) then
1899 !* The state (jkpt,jband,jsppol) is stored in the array cg of this processor ; shift are incremented.
1900 !                     icg=icg+npwj
1901 !                   end if
1902 !                   ibg=ibg+1
1903 !* Skip the end of the loop
1904 !                   cycle
1905 !                 end if
1906 
1907 !* increment the number of occupied bands treated on this processor
1908              my_jband = my_jband+1
1909 
1910 !* In this case, the processor calculates the exchange with the occupied state (jkpt,my_jband).
1911              if (mpi_enreg%proc_distrb(ikpt,iband,isppol)==mpi_enreg%me_kpt) then
1912 !* The state (ikpt,iband,isppol) is stored in the array cg of this processor and copied in cgocc_tmp.
1913                if(icg==-1) then
1914                  write(100,*) 'icg=-1',mpi_enreg%me,isppol,my_jsppol,jkpt,my_jkpt,ikpt,iband
1915                end if
1916              ! MG: Why packing re and im part?
1917                cgocc_tmp(1)=occ(iband+ibg)
1918                cgocc_tmp(2:npwj+1)=cg(1,1+(iband-1)*npwj+icg:iband*npwj+icg)
1919                cgocc_tmp(npwj+2:2*npwj+1)=cg(2,1+(iband-1)*npwj+icg:iband*npwj+icg)
1920                if (fockcommon%usepaw==1) then
1921                  call pawcprj_copy(cprj(:,icp+iband:icp+iband+nspinor-1),cprj_tmp)
1922                end if
1923              end if
1924 
1925 !* Broadcast the state (ikpt,iband,isppol) to all the processors of comm_kpt for cgocc
1926              call timab(1503,1,tsec2)
1927              call xmpi_bcast(cgocc_tmp,mpi_enreg%proc_distrb(ikpt,iband,isppol),mpi_enreg%comm_kpt,ier)
1928 
1929 !* Broadcast the state (ikpt,iband,isppol) to all the processors of comm_kpt for cprj
1930              if (fockcommon%usepaw==1) then
1931                call pawcprj_bcast(cprj_tmp,dtset%natom,nspinor,dimlmn,ncpgr,mpi_enreg%proc_distrb(ikpt,iband,isppol),&
1932 &               mpi_enreg%comm_kpt,ier)
1933              end if
1934              call timab(1503,2,tsec2)
1935 !* Keep the processors in %comm_kpt which needs the values in cgocc_tmp to build their own %cwaveocc and %occ_bz.
1936              if ((mpi_enreg%nproc_kpt/=1).and.(nsppol==2)) then
1937                if (fockbz%timerev(my_jkpt)==mpi_enreg%my_isppoltab(isppol)) cycle
1938 !* In the case of a parallel spin-polarized calculation
1939 !* when time reversal symmetry is applied at this k-point (timrev==1), only the processors with the opposite spin (my_isppoltab==0) are kept.
1940 !* when time reversal symmetry is not applied at this k-point (timrev==0), only the processors with the same spin (my_isppoltab==1) are kept.
1941 
1942 !             if (fock%timerev(my_jkpt)==1)) then
1943 !               if (mpi_enreg%my_isppoltab(isppol)==1) cycle
1944 !* In the case of a parallel spin-polarized calculation and when time reversal symmetry is applied at this k-point,
1945 !* only the processors with the opposite spin are kept.
1946 !             else
1947 !               if (mpi_enreg%my_isppoltab(isppol)==0) cycle
1948 !* only the processors with isppol are kept.
1949 !             end if
1950              end if
1951 
1952 !* Copy the values of cgocc_tmp in the arrays cgocc and %occ_bz
1953              fockbz%occ_bz(my_jband+jbg,my_jsppol) = cgocc_tmp(1)
1954              cgocc(1,1:npwj) = cgocc_tmp(2:npwj+1)
1955              cgocc(2,1:npwj) = cgocc_tmp(npwj+2:2*npwj+1)
1956 
1957 !* calculate cg and store it in cgocc_bz
1958              if (fockbz%calc_phase(my_jkpt)==1) then
1959                do jpw=1,npwj
1960                  cgre=cgocc(1,jpw) ; cgim=cgocc(2,jpw)
1961                  cgocc(1,jpw) = phase_jkpt(1,jpw)*cgre - phase_jkpt(2,jpw)*cgim
1962                  cgocc(2,jpw) = phase_jkpt(1,jpw)*cgim + phase_jkpt(2,jpw)*cgre
1963                end do
1964              end if ! phase
1965 
1966 !* apply time reversal symmetry if necessary
1967              if (fockbz%timerev(my_jkpt)==1) then
1968                cgocc(2,:) = - cgocc(2,:)
1969                if((mpi_enreg%nproc_kpt==1).and.(nsppol==2)) my_jsppol=mod(my_jsppol,2)+1
1970 !* exchange spin (1 ->2 ; 2-> 1) in the sequential case.
1971              end if
1972 
1973 !* apply FFT to get cwaveocc in real space
1974 
1975              if (allocated(fockbz%cwaveocc_bz)) then
1976 
1977                ABI_ALLOCATE(dummytab2,(2,npwj))
1978                call fourwf(1,dummytab3,cgocc(:,1:npwj),dummytab2,fockbz%cwaveocc_bz(:,:,:,:,my_jband+jbg,my_jsppol), &
1979 &               gbound_k,gbound_k,jstwfk,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,&
1980 &               npwj,npwj,n4,n5,n6,tim_fourwf0,dtset%paral_kgb,0,weight1,weight1,use_gpu_cuda=dtset%use_gpu_cuda)
1981                ABI_DEALLOCATE(dummytab2)
1982 
1983              else
1984                fockbz%cgocc(:,jcg+1+(my_jband-1)*npwj:jcg+my_jband*npwj,my_jsppol)=cgocc(:,1:npwj)
1985              end if
1986 
1987 !* calculate cprj and store it in cwaveocc_prj
1988              if (fockcommon%usepaw==1) then
1989                iband_cprj=(my_jsppol-1)*fockbz%mkptband+jbg+my_jband
1990                nband=1;mband0=1;iband0=1
1991                call pawcprj_symkn(fockbz%cwaveocc_prj(:,iband_cprj:iband_cprj+nspinor-1),cprj_tmp(:,1:nspinor),&
1992 &               indsym_,dimlmn,iband0,indlmn,&
1993 &               fockbz%tab_symkpt(my_jkpt),fockbz%timerev(my_jkpt),dtset%kptns(:,ikpt),fockbz%pawang%l_max-1,lmnmax,&
1994 &               mband0,dtset%natom,nband,nspinor,dtset%nsym,dtset%ntypat,typat_srt,fockbz%pawang%zarot)
1995 
1996                if(dtset%optforces==1) then
1997                  do iatom=1,dtset%natom
1998                    iatm=fockcommon%atindx(iatom)
1999                    do ispinor=iband_cprj,iband_cprj+nspinor-1
2000                      do ilmn=1,fockcommon%pawtab(dtset%typat(iatom))%lmn_size
2001                        dcp(:)= MATMUL(TRANSPOSE(fockcommon%symrec(:,:,fockbz%tab_symkpt(my_jkpt))),&
2002 &                                               fockbz%cwaveocc_prj(iatm,ispinor)%dcp(1,:,ilmn))
2003                        fockbz%cwaveocc_prj(iatm,ispinor)%dcp(1,:,ilmn)=dcp(:)
2004                        dcp(:)= MATMUL(TRANSPOSE(fockcommon%symrec(:,:,fockbz%tab_symkpt(my_jkpt))),&
2005 &                                               fockbz%cwaveocc_prj(iatm,ispinor)%dcp(2,:,ilmn))
2006                        fockbz%cwaveocc_prj(iatm,ispinor)%dcp(2,:,ilmn)=dcp(:)
2007                      end do
2008                    end do
2009                  end do
2010                end if
2011 
2012              end if
2013 
2014            end if ! band occupied
2015 
2016 !* update the shift to apply to occ in all case because this array is not distributed among the proc.
2017 !               ibg=ibg+1
2018 
2019          end do ! iband
2020 
2021 !* Save the true number of occupied bands in the array %nbandocc_bz
2022          fockbz%nbandocc_bz(my_jkpt,my_jsppol) = my_jband
2023 
2024 !* update the shifts to apply
2025          jbg=jbg+my_jband
2026          jcg=jcg+npwj*my_jband
2027          jkg=jkg+npwj
2028        end do ! ikpt
2029      end do ! isppol
2030      if (allocated(fockbz%cwaveocc_bz)) then
2031        fockbz%cwaveocc_bz=fockbz%cwaveocc_bz*invucvol
2032      end if
2033 
2034      ABI_DEALLOCATE(cgocc_tmp)
2035      ABI_DEALLOCATE(cgocc)
2036      if (fockcommon%usepaw==1) then
2037        ABI_DEALLOCATE(indlmn)
2038        ABI_DEALLOCATE(indsym_)
2039        ABI_DEALLOCATE(typat_srt)
2040        ABI_DEALLOCATE(dimlmn)
2041        call pawcprj_free(cprj_tmp)
2042        ABI_DATATYPE_DEALLOCATE(cprj_tmp)
2043      end if
2044      if(allocated(phase_jkpt)) then
2045        ABI_DEALLOCATE(phase_jkpt)
2046      end if
2047      ABI_DEALLOCATE(dummytab3)
2048 
2049 
2050 ! Restricted or unrestricted HF
2051      if (nsppol==1) then
2052 !* Update the array %occ_bz => May be limited to the occupied states only
2053        fockbz%occ_bz(:,:)=half*fockbz%occ_bz(:,:)
2054 
2055 ! If nsppol=1, this is a restricted Hartree-Fock calculation.
2056 ! If nsppol=2, this is an unrestricted Hartree-Fock calculation.
2057      end if
2058 
2059  end if
2060 
2061  call timab(1502,2,tsec)
2062 
2063 end subroutine fock_updatecwaveocc

m_fock/fock_updateikpt [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fock_updateikpt

FUNCTION

  Update the value of ikpt,isppol for the next exact exchange calculation.

INPUTS

  fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
  gs_ham <type(gs_hamiltonian_type)>= all data for the Hamiltonian to be applied
  ikpt= reduced planewave coordinates.
  isppol= number of planewaves in basis at this k point

OUTPUT

  none

SIDE EFFECTS

   The field fock%eigen_ikpt is also set to 0.d0.

NOTES

  May be improved to calculate the star of ikpt. => I think NO finally

PARENTS

      vtorho

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1180 subroutine fock_updateikpt(fock,ikpt,isppol)
1181 
1182 
1183 !This section has been created automatically by the script Abilint (TD).
1184 !Do not modify the following lines by hand.
1185 #undef ABI_FUNC
1186 #define ABI_FUNC 'fock_updateikpt'
1187 !End of the abilint section
1188 
1189  implicit none
1190 
1191 !Arguments ------------------------------------
1192  integer, intent(in) :: ikpt,isppol
1193  type(fock_common_type),pointer :: fock
1194 
1195 ! *************************************************************************
1196 
1197  !write (std_out,*) ' fock_updateikpt : enter'
1198 
1199 ! ======================================================
1200 ! === Update the data relative to the current states ===
1201 ! ======================================================
1202 !* Copy of the value ikpt in the field ikpt
1203    fock%ikpt=ikpt
1204 !* Copy of the value isppol in the field isppol
1205    fock%isppol=isppol
1206 !* Set all the Fock contributions to the eigenvalues to 0.d0.
1207    fock%eigen_ikpt=zero
1208 !* Set all the Fock contributions to the forces to 0.d0.
1209    if ((fock%optfor).and.(fock%use_ACE==0)) then
1210      fock%forces_ikpt=zero
1211    endif
1212 
1213 end subroutine fock_updateikpt

m_fock/fockbz_create [ Functions ]

[ Top ] [ m_fock ] [ Functions ]

NAME

  fockbz_create

FUNCTION

  Create a fock_type structure.

INPUTS

OUTPUT

  none

SIDE EFFECTS

  fock <type(fock_type)>= all the quantities to calculate Fock exact exchange are allocated

NOTES

  ############################
  ### Not fully tested yet ###
  ############################

  The current version is restricted to the case nsym=1, nspinor=1 and mkmem/=0.

PARENTS

      m_fock

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

350 subroutine fockbz_create(fockbz,mgfft,mpw,mkpt,mkptband,my_nsppol,natom,n4,n5,n6,use_ACE)
351 
352 
353 !This section has been created automatically by the script Abilint (TD).
354 !Do not modify the following lines by hand.
355 #undef ABI_FUNC
356 #define ABI_FUNC 'fockbz_create'
357 !End of the abilint section
358 
359  implicit none
360 
361 !Arguments ------------------------------------
362 !scalars
363  integer, intent(in) :: mgfft,mpw,mkpt,mkptband,my_nsppol,natom,n4,n5,n6,use_ACE
364  type(fock_BZ_type) , intent(inout) :: fockbz
365 
366 !Local variables-------------------------------
367 !scalars
368 !arrays
369 ! character(len=500) :: message                   ! to be uncommented, if needed
370 
371 ! *************************************************************************
372 
373  !write (std_out,*) ' fockbz_create : enter'
374 
375 !* Create the array %kptns_bz = the k points in full BZ
376  ABI_ALLOCATE(fockbz%kptns_bz,(3,mkpt))
377  fockbz%kptns_bz=zero
378 !* Create the array %jstwfk = how is stored the wavefunction at each k-point
379 !* By default, the table is initialized to 1 (do NOT take advantage of the time-reversal symmetry)
380  ABI_ALLOCATE(fockbz%istwfk_bz,(mkpt))
381  fockbz%istwfk_bz=1
382 !* Create the array %wtk_bz = weight assigned to each k point.
383  ABI_ALLOCATE(fockbz%wtk_bz,(mkpt))
384  fockbz%wtk_bz=zero
385 !* Create the array %npwarr_bz = number of planewaves in basis at each k-point
386 !    ABI_ALLOCATE(fockbz%npwarr_bz,(mkpt))
387 !    fockbz%npwarr_bz=0
388 !* Create the array %kg_bz = reduced planewave coordinates at each k-point
389  ABI_ALLOCATE(fockbz%kg_bz,(3,mpw*mkpt))
390  fockbz%kg_bz=0
391 !* Create the array %gbound_bz = boundary of the basis sphere of G vectors at each k-point
392  ABI_ALLOCATE(fockbz%gbound_bz,(2*mgfft+8,2,mkpt))
393  fockbz%gbound_bz=0
394 
395 !* Create the array %tab_ikpt = indices of k-point ikpt in IBZ which corresponds to each k-point jkpt in full BZ
396  ABI_ALLOCATE(fockbz%tab_ikpt,(mkpt))
397  fockbz%tab_ikpt=0
398 !* Create the array %tab_symkpt =indices of symmetry operation to apply to get jkpt in full BZ from ikpt in IBZ
399  ABI_ALLOCATE(fockbz%tab_symkpt,(mkpt))
400  fockbz%tab_symkpt=0
401 !* Create the array %tab_ibg = indices of occ(ikpt) in the arrays cprj/occ for each k-point jkpt
402  ABI_ALLOCATE(fockbz%tab_ibg,(mkpt,my_nsppol))
403  fockbz%tab_ibg=0
404 !* Create the array %tab_icp = indices of cprj(ikpt) in the arrays cprj/occ for each k-point jkpt
405  ABI_ALLOCATE(fockbz%tab_icp,(mkpt,my_nsppol))
406  fockbz%tab_icp=0
407 !* Create the array %tab_icg = indices of cg(ikpt) in the arrays cg for each k-point jkpt
408  ABI_ALLOCATE(fockbz%tab_icg,(mkpt,my_nsppol))
409  fockbz%tab_icg=0
410 
411 !* Create the array %calc_phase = 1 if a phase factor must be considered (0 otherwise) at each k point
412  ABI_ALLOCATE(fockbz%calc_phase,(mkpt))
413  fockbz%calc_phase=0
414 !* Create the array %phase = phase factor the cg array will be multiplied with at each k point
415  ABI_ALLOCATE(fockbz%phase,(2,mpw*mkpt))
416  fockbz%phase=zero
417 
418 !* Create the array %timerev i= 1 if time reversal symmetry must be used (0 otherwise) at each k point
419  ABI_ALLOCATE(fockbz%timerev,(mkpt))
420  fockbz%timerev=0
421 
422 !* Create the array %cwaveocc_bz = wavefunctions of each bands at each k point
423 
424  if (use_ACE==1) then
425    ABI_ALLOCATE(fockbz%cgocc,(2,mpw*mkptband,my_nsppol))
426    fockbz%cgocc=zero
427  else
428    ABI_ALLOCATE(fockbz%cwaveocc_bz,(2,n4,n5,n6,mkptband,my_nsppol))
429    fockbz%cwaveocc_bz=zero
430  end if
431 !* Create the array %occ_bz = occupancy of each bands at each k point => will be limited to only the occupied states
432  ABI_ALLOCATE(fockbz%occ_bz,(mkptband,my_nsppol))
433  fockbz%occ_bz=zero
434 !* Create the array %nbandocc_bz = nb of bands at each k point
435  ABI_ALLOCATE(fockbz%nbandocc_bz,(mkpt,my_nsppol))
436  fockbz%nbandocc_bz=0
437 
438  ABI_ALLOCATE(fockbz%npwarr,(mkpt))
439  fockbz%npwarr=0
440 
441 
442 
443 end subroutine fockbz_create