Line data Source code
1 : !!****m* ABINIT/m_fock
2 : !! NAME
3 : !! m_fock
4 : !!
5 : !! FUNCTION
6 : !! This module provides the definition of
7 : !! the fock_type used to store data for the calculation of Fock exact exchange term
8 : !! and the procedures to perform this calculation.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2012-2026 ABINIT group (CMartins,FJ,FA,MT)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_fock
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_errors
29 : use m_mpinfo
30 : use m_xmpi
31 : use libxc_functionals
32 : use m_pawang
33 : use m_pawtab
34 : use m_pawfgr
35 : use m_pawfgrtab
36 : use m_pawcprj
37 : use m_cgtools
38 : use m_nctk
39 : use m_dtset
40 :
41 : use defs_abitypes, only : MPI_type
42 : use m_time, only : timab
43 : use m_fstrings, only : itoa, ftoa, sjoin
44 : use m_matrix, only : mati3inv, matr3inv
45 : use m_fftcore, only : sphereboundary
46 : use m_fft, only : zerosym, fourwf
47 : use m_kg, only : ph1d3d, getph
48 : use m_kpts, only : listkk
49 : use m_barevcoul, only : barevcoul
50 :
51 : implicit none
52 :
53 : private
54 : !!***
55 :
56 : !!****t* m_fock/fock_type
57 : !! NAME
58 : !! fock_type
59 : !!
60 : !! FUNCTION
61 : !! This object stores the occupied wavefunctions and other quantities
62 : !! needed to calculate Fock exact exchange
63 : !!
64 : !! SOURCE
65 :
66 : type, public :: fock_type
67 : type(fock_common_type), pointer :: fock_common=> null()
68 : type(fock_BZ_type), pointer :: fock_BZ=> null()
69 : type(fock_ACE_type), pointer :: fockACE(:,:)=> null()
70 : end type fock_type
71 :
72 : type, public :: fock_common_type
73 :
74 : ! Integer scalars
75 : !integer :: mcgocc_bz,mkg_bz,mocc
76 : !integer :: natom,ntypat
77 :
78 : integer :: usepaw
79 : ! 0 if norm-conserving psps, 1 for PAW (not implemented)
80 :
81 : integer :: ikpt,isppol,ieigen,iband
82 : ! data relative to the current states.
83 :
84 : integer :: mband
85 : ! maximum number of bands
86 :
87 : integer :: my_nsppol
88 : ! my_nsppol=1 when nsppol=1 or nsppol=2 and only one spin is treated by the processor.
89 : ! my_nsppol=2 when nsppol=2 and no parallelization over kpt (both spins are treated by the processor).
90 :
91 : integer :: natom
92 : ! Number of atoms, input variable
93 :
94 : integer :: nsppol
95 : ! Number of independent spin polarizations, input variable
96 : ! Note that this value does not take into account the MPI distribution of the wavefunctions.
97 :
98 : integer :: ntypat
99 : ! Number of type of atoms
100 :
101 : integer :: nnsclo_hf
102 : ! Number of iterations with fixed occupied states when calculating the exact exchange contribution.
103 :
104 : integer :: ixc
105 : ! XC option (abinit input variable)
106 :
107 : integer :: use_ACE
108 : ! option to use the ACE method of Lin Lin
109 : !==0 if the normal Fock operator is to be created and/or used
110 : !==1 if the ACE operator is to be created and/or used
111 :
112 : integer :: fock_icutcoul
113 : ! contains input variable fock_icutcoul that specifies the treatment of the Coulomb interaction (cutoff in real-space,
114 : ! divergence treatment etc.)
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 : real(dp) :: rcut
154 : ! contains input variable rcut that specifies the cutoff radius for spherical cutoff
155 :
156 : real(dp) :: e_fock0
157 : ! contribution of the Fock term to energy (computed and stored here in case of ACE)
158 :
159 : integer, allocatable :: atindx(:)
160 : ! atindx(natom)=index table for atoms (see gstate.f)
161 :
162 : integer, allocatable :: nband(:)
163 : ! nband(nkpt)
164 : ! Number of bands for each k point
165 :
166 : integer, allocatable :: symrec(:,:,:)
167 :
168 : integer,allocatable :: typat(:)
169 : ! typat(natom)
170 : ! type of each atom
171 :
172 : ! Real(dp) arrays
173 : real(dp) :: stress(6)
174 : ! stress(6)
175 : ! contribution of the fock term to stresses
176 :
177 : real(dp), allocatable :: stress_ikpt(:,:)
178 : ! stress(6,nband)
179 : ! contribution of the fock term to stresses for the current band
180 :
181 : real(dp), allocatable :: forces_ikpt(:,:,:)
182 : ! forces(3,natom,nband))
183 : ! contribution of the fock term to forces for the current band
184 :
185 : real(dp), allocatable :: forces(:,:)
186 : ! forces(3,natom))
187 : ! contribution of the fock term to forces
188 :
189 : real(dp), allocatable :: eigen_ikpt(:)
190 : ! eigen_ikpt,(nband))
191 : ! Will contain the band index of the current state
192 : ! if the value is 0, the Fock contribution to the eigenvalue is not calculated.
193 :
194 : ! Pointers to PAW-types (associated only if usepaw==1)
195 : ! Note that these are references to already existing objects.
196 :
197 : type(pawtab_type), pointer :: pawtab(:) => null()
198 : type(pawfgr_type),pointer :: pawfgr => null()
199 : type(pawfgrtab_type),allocatable :: pawfgrtab(:)
200 :
201 : end type fock_common_type
202 :
203 : type, public :: fock_BZ_type
204 :
205 : integer :: mcprj
206 : ! dimension of cwaveocc_cprj
207 :
208 : integer :: mkpt
209 : ! maximum number of k-points for Fock treated by this node
210 :
211 : integer :: mkptband
212 : ! size of occupied states stored by this node.
213 :
214 : integer :: nkpt_bz
215 : ! Number of k-points in the BZ for Fock operator
216 :
217 : integer, allocatable :: gbound_bz(:,:,:)
218 : ! gbound_bz(2*mgfft+8,2,mkpt)
219 : ! Tables for zero-padded FFT of wavefunctions.
220 :
221 : integer, allocatable :: kg_bz(:,:)
222 : ! kg_bz(3,mpw*mkpt)
223 : ! G-vectors for each k-point in the BZ treate by this node
224 :
225 : integer, allocatable :: nbandocc_bz(:,:)
226 : ! nbandocc_bz,(mkpt,my_nsppol))
227 : ! nb of bands at each k point
228 :
229 : integer, allocatable :: npwarr(:)
230 : ! npwarr(mkpt)
231 :
232 : integer, allocatable :: istwfk_bz(:)
233 : ! istwfk_bz,(mkpt))
234 : ! storage mode of the wavefunction at each k-point
235 :
236 : integer, allocatable :: calc_phase(:)
237 : ! calc_phase,(mkpt))
238 : ! 1 if a phase factor must be considered (0 otherwise) at each k point
239 :
240 : integer, allocatable :: tab_symkpt(:)
241 : ! tab_symkpt,(mkpt))
242 : ! indices of symmetry operation to apply to get jkpt in full BZ from ikpt in IBZ
243 :
244 : integer, allocatable :: timerev(:)
245 : ! timerev,(mkpt))
246 : ! 1 if time reversal symmetry must be used (0 otherwise) at each k point
247 :
248 : integer, allocatable :: tab_ibg(:,:)
249 : ! tab_ibg,(mkpt,my_nsppol))
250 : ! indices of cprj(ikpt)/occ(ikpt) in the arrays cprj/occ for each k-point jkpt
251 :
252 : integer, allocatable :: tab_icg(:,:)
253 : ! tab_icg,(mkpt,my_nsppol))
254 : ! indices of cg(ikpt) in the arrays cg for each k-point jkpt
255 :
256 : integer, allocatable :: tab_icp(:,:)
257 : ! tab_icg,(mkpt,my_nsppol))
258 : ! indices of cprj(ikpt) in the arrays cprj for each k-point jkpt
259 :
260 : integer, allocatable :: tab_ikpt(:)
261 : ! tab_ikpt,(mkpt))
262 : ! indices of k-point ikpt in IBZ which corresponds to each k-point jkpt in full BZ
263 :
264 : real(dp), allocatable :: cgocc(:,:,:)
265 : ! cgocc(2,npw*mkptband,my_nsppol)
266 : ! wavefunction in the G-space
267 :
268 : real(dp), allocatable :: cwaveocc_bz(:,:,:,:,:,:)
269 : ! (2,n4,n5,n6,mkptband,my_nsppol))
270 : ! occupied states of each bands at each k point (used to construct Fock operator), in the real space
271 : real(dp), allocatable :: occ_bz(:,:)
272 : ! occ_bz(mkptband,my_nsppol))
273 : ! occupancy of each bands at each k point
274 :
275 : real(dp), allocatable :: wtk_bz(:)
276 : ! wtk_bz,(mkpt))
277 : ! weights assigned to each k point in the BZ
278 : ! Caution, the definition takes into account "ucvol" !
279 :
280 : real(dp), allocatable :: kptns_bz(:,:)
281 : ! kptns_bz(3,mkpt)
282 : ! k-points in full BZ
283 :
284 : real(dp), allocatable :: phase(:,:)
285 : ! phase(2,mpw*mkpt))
286 : ! phase factor the cg array will be multiplied with at each k point
287 :
288 : type(MPI_type) :: mpi_enreg
289 : type(pawang_type),pointer :: pawang
290 : type(pawcprj_type), allocatable :: cwaveocc_prj(:,:)
291 :
292 : end type fock_BZ_type
293 : !----------------------------------------------------------------------
294 :
295 : type,public :: fock_ACE_type
296 :
297 : real(dp), allocatable :: xi(:,:,:)
298 :
299 : end type fock_ACE_type
300 : !----------------------------------------------------------------------
301 :
302 : public :: fock_init ! Initialize the object.
303 : !public :: fock_from_wfk ! Initialize the object from external WFK file.
304 : public :: fock_set_ieigen ! Set the value of ieigen to the value given in argument.
305 : public :: fock_updateikpt ! Update the value of energies%e_xc and energies%e_xcdc with Fock contribution.
306 : public :: fock_destroy ! Free memory.
307 : public :: fock_ACE_destroy ! Free memory.
308 : public :: fock_common_destroy ! Free memory.
309 : public :: fock_bz_destroy ! Free memory.
310 : public :: fock_calc_ene ! Calculate the Fock contribution to the total energy.
311 : public :: fock_update_exc ! Update the value of energies%e_xc and energies%e_xcdc with Fock contribution.
312 : public :: fock_updatecwaveocc ! Update in the fock datastructure the fields relative to the occupied states.
313 : public :: fock_set_getghc_call ! Enable/disable the call to fock_getghc in getghc.
314 : public :: fock_get_getghc_call ! Return the value of the flag used to enable/disable the call to fock_getghc in getghc.
315 : public :: fock_print ! Print info on the object.
316 : !!***
317 :
318 : ! Help functions
319 : public :: bare_vqg
320 : public :: strfock
321 :
322 : contains
323 : !!***
324 :
325 : !!****f* m_fock/fockbz_create
326 : !! NAME
327 : !! fockbz_create
328 : !!
329 : !! FUNCTION
330 : !! Create a fock__BZ_type structure.
331 : !!
332 : !! INPUTS
333 : !!
334 : !! NOTES
335 : !!
336 : !! ############################
337 : !! ### Not fully tested yet ###
338 : !! ############################
339 : !!
340 : !! The current version is restricted to the case nsym=1, nspinor=1 and mkmem/=0.
341 : !!
342 : !! SOURCE
343 :
344 123 : subroutine fockbz_create(fockbz,mgfft,mpw,mkpt,mkptband,my_nsppol,n4,n5,n6,use_ACE)
345 :
346 : !Arguments ------------------------------------
347 : type(fock_BZ_type), intent(inout) :: fockbz
348 : integer, intent(in) :: mgfft,mpw,mkpt,mkptband,my_nsppol,n4,n5,n6,use_ACE
349 : ! *************************************************************************
350 :
351 : !write (std_out,*) ' fockbz_create : enter'
352 :
353 : !* Create the array %kptns_bz = the k points in full BZ
354 369 : ABI_MALLOC(fockbz%kptns_bz,(3,mkpt))
355 3047 : fockbz%kptns_bz=zero
356 : !* Create the array %jstwfk = how is stored the wavefunction at each k-point
357 : !* By default, the table is initialized to 1 (do NOT take advantage of the time-reversal symmetry)
358 369 : ABI_MALLOC(fockbz%istwfk_bz,(mkpt))
359 854 : fockbz%istwfk_bz=1
360 : !* Create the array %wtk_bz = weight assigned to each k point.
361 369 : ABI_MALLOC(fockbz%wtk_bz,(mkpt))
362 854 : fockbz%wtk_bz=zero
363 : !* Create the array %npwarr_bz = number of planewaves in basis at each k-point
364 : ! ABI_MALLOC(fockbz%npwarr_bz,(mkpt))
365 : ! fockbz%npwarr_bz=0
366 : !* Create the array %kg_bz = reduced planewave coordinates at each k-point
367 369 : ABI_MALLOC(fockbz%kg_bz,(3,mpw*mkpt))
368 576543 : fockbz%kg_bz=0
369 : !* Create the array %gbound_bz = boundary of the basis sphere of G vectors at each k-point
370 615 : ABI_MALLOC(fockbz%gbound_bz,(2*mgfft+8,2,mkpt))
371 61116 : fockbz%gbound_bz=0
372 :
373 : !* Create the array %tab_ikpt = indices of k-point ikpt in IBZ which corresponds to each k-point jkpt in full BZ
374 246 : ABI_MALLOC(fockbz%tab_ikpt,(mkpt))
375 854 : fockbz%tab_ikpt=0
376 : !* Create the array %tab_symkpt =indices of symmetry operation to apply to get jkpt in full BZ from ikpt in IBZ
377 246 : ABI_MALLOC(fockbz%tab_symkpt,(mkpt))
378 854 : fockbz%tab_symkpt=0
379 : !* Create the array %tab_ibg = indices of occ(ikpt) in the arrays cprj/occ for each k-point jkpt
380 492 : ABI_MALLOC(fockbz%tab_ibg,(mkpt,my_nsppol))
381 1031 : fockbz%tab_ibg=0
382 : !* Create the array %tab_icp = indices of cprj(ikpt) in the arrays cprj/occ for each k-point jkpt
383 369 : ABI_MALLOC(fockbz%tab_icp,(mkpt,my_nsppol))
384 1031 : fockbz%tab_icp=0
385 : !* Create the array %tab_icg = indices of cg(ikpt) in the arrays cg for each k-point jkpt
386 369 : ABI_MALLOC(fockbz%tab_icg,(mkpt,my_nsppol))
387 1031 : fockbz%tab_icg=0
388 :
389 : !* Create the array %calc_phase = 1 if a phase factor must be considered (0 otherwise) at each k point
390 246 : ABI_MALLOC(fockbz%calc_phase,(mkpt))
391 854 : fockbz%calc_phase=0
392 : !* Create the array %phase = phase factor the cg array will be multiplied with at each k point
393 369 : ABI_MALLOC(fockbz%phase,(2,mpw*mkpt))
394 432438 : fockbz%phase=zero
395 :
396 : !* Create the array %timerev i= 1 if time reversal symmetry must be used (0 otherwise) at each k point
397 246 : ABI_MALLOC(fockbz%timerev,(mkpt))
398 854 : fockbz%timerev=0
399 :
400 : !* Create the array %cwaveocc_bz = wavefunctions of each bands at each k point
401 :
402 123 : if (use_ACE==1) then
403 492 : ABI_MALLOC(fockbz%cgocc,(2,mpw*mkptband,my_nsppol))
404 2199326 : fockbz%cgocc=zero
405 : else
406 0 : ABI_MALLOC(fockbz%cwaveocc_bz,(2,n4,n5,n6,mkptband,my_nsppol))
407 0 : fockbz%cwaveocc_bz=zero
408 : end if
409 : !* Create the array %occ_bz = occupancy of each bands at each k point => will be limited to only the occupied states
410 492 : ABI_MALLOC(fockbz%occ_bz,(mkptband,my_nsppol))
411 3513 : fockbz%occ_bz=zero
412 : !* Create the array %nbandocc_bz = nb of bands at each k point
413 369 : ABI_MALLOC(fockbz%nbandocc_bz,(mkpt,my_nsppol))
414 1031 : fockbz%nbandocc_bz=0
415 :
416 246 : ABI_MALLOC(fockbz%npwarr,(mkpt))
417 854 : fockbz%npwarr=0
418 :
419 123 : end subroutine fockbz_create
420 : !!***
421 :
422 : !!****f* m_fock/fock_init
423 : !! NAME
424 : !! fock_init
425 : !!
426 : !! FUNCTION
427 : !! Init fock_t object
428 : !!
429 : !! INPUTS
430 : !! cg(2,mcg)= wavefunctions
431 : !! dtset <type(dataset_type)>= all input variables for this dataset
432 : !! gsqcut= Fourier cutoff on G^2 used to calculate charge density
433 : !! kg(3,mpw*mkmem)= reduced planewave coordinates.
434 : !! mcg= size of wave-functions array (cg) = mpw*nspinor*mband*mkmem*nsppol
435 : !! mpi_enreg=information about MPI parallelization
436 : !! npwarr_bz(nkpt)= number of planewaves in basis at this k point
437 : !! occ(mband*nkpt*nsppol)= occupation number for each band (often 2) at each k point
438 : !!
439 : !! SIDE EFFECTS
440 : !! fock <type(fock_type)>= all the quantities to calculate Fock exact exchange are initialized
441 : !!
442 : !! NOTES
443 : !!
444 : !! ############################
445 : !! ### Not fully tested yet ###
446 : !! ############################
447 : !!
448 : !! The current version is restricted to the case nsym=1, nspinor=1 and mkmem/=0.
449 : !!
450 : !! SOURCE
451 :
452 123 : subroutine fock_init(atindx,cplex,dtset,fock,gsqcut,kg,mpi_enreg,nattyp,npwarr,pawang,pawfgr,pawtab,rprimd)
453 :
454 : !Arguments ------------------------------------
455 : !scalars
456 : integer, intent(in) :: cplex
457 : real(dp),intent(in) :: gsqcut
458 : type(dataset_type),intent(in) :: dtset
459 : type(MPI_type),intent(in) :: mpi_enreg
460 : type(fock_type),intent(inout),pointer :: fock
461 : type(pawfgr_type),intent(in),target :: pawfgr
462 : type(pawang_type),intent(in),target :: pawang
463 : !arrays
464 : integer, intent(in) :: atindx(dtset%natom),nattyp(dtset%ntypat), npwarr(dtset%nkpt)
465 : integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem)
466 : real(dp), intent(in) :: rprimd(3,3)
467 : type(pawtab_type), intent(in),target :: pawtab(dtset%ntypat*dtset%usepaw)
468 : !Local variables-------------------------------
469 : !scalars
470 : integer :: iatom,ibg,icg,icp,ier,ik,ikg,ikpt,isppol,isym,itypat,jkpt,jpw,jsym,mband,mgfft,mkpt,mkptband
471 : integer :: n1,n2,n3,n4,n5,n6,nband,ncpgr,nkpt_bz,nproc_hf,npwj,timrev,use_ACE,v1,v2,v3
472 : integer :: my_jkpt,jkg_this_proc,my_nsppol,my_nspinor
473 : real(dp) :: dksqmax,arg
474 : character(len=500) :: msg
475 : !arrays
476 246 : integer :: indx(1),l_size_atm(dtset%natom),shiftg(3),symm(3,3),ident(3,3),symrec(3,3,dtset%nsym)
477 246 : real(dp) :: gmet(3,3),gprimd(3,3),tau_nons(3),phktnons(2,1),tsec(2),Rtnons(3,dtset%nsym)
478 123 : integer,allocatable :: dimcprj(:),indkk(:,:),kg_tmp(:),my_ikgtab(:),my_ibgtab(:,:),my_icgtab(:,:),my_icptab(:,:),invsym(:)
479 123 : real(dp),allocatable :: kptns_hf(:,:), phase1d(:,:)
480 : type(fock_common_type),pointer :: fockcommon
481 : type(fock_BZ_type),pointer :: fockbz
482 : ! *************************************************************************
483 :
484 : DBG_ENTER("COLL")
485 :
486 123 : call timab(1501,1,tsec)
487 123 : ABI_CHECK_IEQ(dtset%nspinor, 1, 'Hartree-Fock option can be used only with option nspinor = 1')
488 :
489 : ! =====================================
490 : ! === Define useful local variables ===
491 : ! =====================================
492 :
493 123 : nkpt_bz=dtset%nkpthf
494 123 : nproc_hf=mpi_enreg%nproc_hf
495 123 : mband=dtset%nbandhf
496 :
497 123 : n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
498 123 : n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
499 :
500 : !* Allocations
501 369 : ABI_MALLOC(kptns_hf,(3,nkpt_bz))
502 3047 : kptns_hf=zero
503 369 : ABI_MALLOC(indkk,(nkpt_bz,6))
504 5247 : indkk=0
505 369 : ABI_MALLOC(phase1d,(2,(2*n1+1)*(2*n2+1)*(2*n3+1)))
506 27690318 : phase1d=zero
507 369 : ABI_MALLOC(kg_tmp,(3*dtset%mpw))
508 :
509 : !* Initialize the array my_ikgtab = shifts in arrays kg(ikg) associated to ikpt
510 369 : ABI_MALLOC(my_ikgtab,(dtset%nkpt))
511 123 : ikg=0
512 803 : do ikpt=1,dtset%nkpt
513 680 : nband=dtset%nband(ikpt)
514 803 : if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband,-1,mpi_enreg%me_kpt))) then
515 : !* The point ikpt is treated on this processor.
516 520 : my_ikgtab(ikpt)=ikg
517 : !* The array kg is distributed, the shift ikg is incremented only on this proc.
518 520 : ikg=ikg+npwarr(ikpt)
519 : else
520 160 : my_ikgtab(ikpt)=-1
521 : !* Default value is -1.
522 : end if
523 : end do
524 :
525 : !* Initialize the array my_ibgtab = shifts in arrays occ(ibg) associated to ikpt
526 : !* Initialize the array my_icgtab = shifts in arrays cg(icg) associated to ikpt
527 492 : ABI_MALLOC(my_ibgtab,(dtset%nkpt,dtset%nsppol))
528 369 : ABI_MALLOC(my_icgtab,(dtset%nkpt,dtset%nsppol))
529 369 : ABI_MALLOC(my_icptab,(dtset%nkpt,dtset%nsppol))
530 123 : ibg=0; icg=0 ;icp=0
531 260 : do isppol=1,dtset%nsppol
532 957 : do ikpt=1,dtset%nkpt
533 697 : nband=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
534 697 : if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband,isppol,mpi_enreg%me_kpt))) then
535 : !* The states with (ikpt,isppol) are stored on this processor.
536 537 : my_icgtab(ikpt,isppol)=icg
537 537 : my_icptab(ikpt,isppol)=icp
538 : !* The array cg is distributed, the shift icg is incremented only on this proc.
539 537 : icg=icg+npwarr(ikpt)*nband
540 537 : icp=icp+nband
541 : else
542 160 : my_icgtab(ikpt,isppol)=-1
543 160 : my_icptab(ikpt,isppol)=-1
544 : !* Otherwise, the states with (ikpt,isspol) are not stored on this processor and default value is -1.
545 : end if
546 : !* The array occ is shared among the proc, the shift ibg is always incremented.
547 697 : my_ibgtab(ikpt,isppol)=ibg
548 834 : ibg=ibg+nband
549 : end do
550 : end do
551 :
552 123 : if (.not.(associated(fock))) then
553 :
554 : ! =================================
555 : ! === Create the fock structure ===
556 : ! =================================
557 123 : ABI_MALLOC(fock,)
558 123 : ABI_MALLOC(fock%fock_common,)
559 123 : ABI_MALLOC(fock%fock_BZ,)
560 : ! ========================================================
561 : ! === Set all the other state-dependent fields to zero ===
562 : ! ========================================================
563 123 : fockcommon=>fock%fock_common
564 123 : fockbz=> fock%fock_BZ
565 369 : ABI_MALLOC(fockcommon%nband,(dtset%nkpt*dtset%nsppol))
566 820 : do ikpt=1,dtset%nkpt*dtset%nsppol
567 820 : fockcommon%nband(ikpt)=dtset%nband(ikpt)
568 : end do
569 :
570 123 : nband=dtset%mband
571 123 : fockcommon%ikpt= 0
572 : !* Will contain the k-point ikpt of the current state
573 123 : fockcommon%isppol= 0
574 : !* Will contain the spin isppol of the current state
575 123 : fockcommon%ieigen=0
576 : !* Will contain the band index of the current state
577 : !* if the value is 0, the Fock contribution to the eigenvalue is not calculated.
578 369 : ABI_MALLOC(fockcommon%eigen_ikpt,(nband))
579 748 : fockcommon%eigen_ikpt=0.d0
580 : !* Will contain the Fock contributions to the eigenvalue of the current state
581 :
582 : !* Compute the dimension of arrays in "spin" w.r.t parallelism
583 123 : my_nsppol=dtset%nsppol
584 123 : if (mpi_enreg%nproc_spkpt>1) my_nsppol=1
585 : !* my_nsppol=1 when nsppol=1 or nsppol=2 and only one spin is treated by the processor.
586 : !* my_nsppol=2 when nsppol=2 and no parallelization over kpt (both spins are treated by the processor).
587 :
588 : !* Compute mkpt the size of arrays/pointers for k points w.r.t. parallelism
589 : !* Compute mkptband the size of arrays/pointers for occupied states w.r.t. parallelism
590 123 : if (nproc_hf<nkpt_bz) then
591 : !* Parallelization over kpts only
592 87 : mkpt=nkpt_bz/nproc_hf
593 87 : if (mod(nkpt_bz,nproc_hf) /=0) mkpt=mkpt+1
594 87 : mkptband=mkpt*mband
595 : else
596 : !* Parallelization over occupied states
597 36 : if (nproc_hf<nkpt_bz*mband) then
598 35 : mkptband=(nkpt_bz*mband)/nproc_hf
599 35 : if (mod((nkpt_bz*mband),nproc_hf) /=0) mkptband=mkptband+1
600 35 : mkpt=1
601 35 : if (mod(nproc_hf,nkpt_bz) /=0) mkpt=2
602 : else
603 1 : mkptband=1
604 1 : mkpt=1
605 : end if
606 : end if
607 :
608 : ! mpi_enreg settings
609 123 : call copy_mpi_enreg(mpi_enreg,fockbz%mpi_enreg)
610 123 : fockbz%mpi_enreg%me_kpt=mpi_enreg%me_hf
611 123 : fockbz%mpi_enreg%comm_kpt=mpi_enreg%comm_hf
612 123 : fockbz%mpi_enreg%nproc_spkpt=mpi_enreg%nproc_hf
613 123 : ABI_SFREE(fockbz%mpi_enreg%proc_distrb)
614 492 : ABI_MALLOC(fockbz%mpi_enreg%proc_distrb,(nkpt_bz,mband,1))
615 854 : do jkpt=1,nkpt_bz
616 3960 : fockbz%mpi_enreg%proc_distrb(jkpt,:,1)=fockbz%mpi_enreg%me_kpt
617 : end do
618 :
619 123 : mgfft=dtset%mgfft
620 123 : fockcommon%usepaw=dtset%usepaw
621 123 : if (fockcommon%usepaw==1)then
622 26 : mgfft=dtset%mgfftdg
623 26 : n4=dtset%ngfftdg(4) ; n5=dtset%ngfftdg(5) ; n6=dtset%ngfftdg(6)
624 : end if
625 123 : fockcommon%optfor=.false.; fockcommon%optstr=.false.
626 123 : if(dtset%optforces==1) fockcommon%optfor=.true.
627 123 : if (fockcommon%optfor) then
628 92 : ABI_MALLOC(fockcommon%forces_ikpt,(3,dtset%natom,nband))
629 69 : ABI_MALLOC(fockcommon%forces,(3,dtset%natom))
630 207 : fockcommon%forces=zero
631 : end if
632 123 : use_ACE=1 ! Default. Normal users do not have access to this variable, although the next line allows experts to make tests.
633 123 : if(dtset%userie==1729)use_ACE=0 ! Hidden possibility to disable ACE
634 :
635 123 : fockcommon%use_ACE=use_ACE
636 123 : call fockbz_create(fockbz,mgfft,dtset%mpw,mkpt,mkptband,my_nsppol,n4,n5,n6,use_ACE)
637 :
638 123 : fockcommon%fock_icutcoul = dtset%fock_icutcoul
639 123 : fockcommon%rcut = dtset%fock_rcut
640 :
641 : !* Initialize %mband, %mkpt, %mkptband = size of arrays
642 123 : fockcommon%mband=mband
643 123 : fockbz%mkpt=mkpt
644 123 : fockbz%mkptband=mkptband
645 123 : fockcommon%my_nsppol = my_nsppol
646 123 : fockcommon%nsppol = dtset%nsppol
647 123 : if (fockcommon%use_ACE/=0) then
648 1326 : ABI_MALLOC(fock%fockACE,(dtset%nkpt,dtset%nsppol))
649 123 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
650 260 : do isppol=1,dtset%nsppol
651 957 : do ikpt=1,dtset%nkpt
652 697 : nband=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
653 2925 : ABI_MALLOC(fock%fockACE(ikpt,isppol)%xi,(2,npwarr(ikpt)*my_nspinor,nband))
654 : end do
655 : end do
656 : end if
657 : !========Initialze PAW data========
658 123 : fockcommon%ntypat=dtset%ntypat
659 123 : fockcommon%natom=dtset%natom
660 123 : if (fockcommon%usepaw==1) then
661 26 : fockbz%mcprj=mkptband*my_nsppol
662 26 : fockcommon%pawfgr => pawfgr
663 26 : fockbz%pawang => pawang
664 26 : fockcommon%pawtab => pawtab
665 130 : ABI_MALLOC(fockcommon%pawfgrtab,(dtset%natom))
666 78 : do iatom = 1, dtset%natom
667 52 : itypat=dtset%typat(iatom)
668 78 : l_size_atm(iatom) = pawtab(itypat)%lcut_size
669 : end do
670 26 : call pawfgrtab_init(fockcommon%pawfgrtab,cplex,l_size_atm,dtset%nspden,dtset%typat)
671 26 : call pawfgrtab_nullify(fockcommon%pawfgrtab)
672 1922 : ABI_MALLOC(fockbz%cwaveocc_prj,(dtset%natom,fockbz%mcprj))
673 78 : ABI_MALLOC(dimcprj,(dtset%natom))
674 26 : call pawcprj_getdim(dimcprj,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
675 26 : ncpgr = 0
676 26 : if (dtset%optforces== 1) ncpgr = 3
677 : ! if (dtset%optstress /= 0) ncpgr = 6
678 : ! ncpgr=3*dtset%optforces+6*dtset%optstress
679 26 : call pawcprj_alloc(fockbz%cwaveocc_prj,ncpgr,dimcprj)
680 26 : ABI_FREE(dimcprj)
681 52 : ABI_MALLOC(fockcommon%atindx,(dtset%natom))
682 104 : fockcommon%atindx=atindx
683 52 : ABI_MALLOC(fockcommon%typat,(dtset%natom))
684 104 : fockcommon%typat=dtset%typat
685 : end if
686 : ! ==========================================
687 : ! === Initialize the convergence options ===
688 : ! ==========================================
689 123 : write(msg,'(2a)') ch10,'Fock_init: initialization of Fock operator parameters:'
690 123 : call wrtout(std_out,msg)
691 :
692 123 : fockcommon%fock_converged=.false.
693 123 : fockcommon%scf_converged=.false.
694 :
695 : !* Number of iterations with fixed occupied states when calculating the exact exchange contribution.
696 123 : if (dtset%nnsclohf<0) then
697 0 : ABI_ERROR('The parameter nnsclohf must be a non-negative integer.')
698 : end if
699 123 : if (dtset%nnsclohf==0) then
700 116 : fockcommon%nnsclo_hf=1
701 116 : msg=' - The parameter nnsclohf is set to its default value 1.'
702 116 : call wrtout(std_out,msg)
703 : !* Default value is set to 1 (updating cgocc at each step)
704 : !* May be useful to put default to 3
705 : else
706 7 : fockcommon%nnsclo_hf=dtset%nnsclohf
707 7 : write(msg,'(a,i3)') ' - The parameter nnsclohf is set to the value:', dtset%nnsclohf
708 7 : call wrtout(std_out,msg)
709 : !* value chosen by the user
710 : end if
711 :
712 : ! =========================================
713 : ! === Initialize the hybrid coefficient ===
714 : ! =========================================
715 123 : fockcommon%ixc = dtset%ixc
716 : ! By convention, positive values are the default values for the ixc,
717 : ! while negative values have been set by the user (and stored as negative numbers)
718 123 : fockcommon%hyb_mixing=abs(dtset%hyb_mixing)
719 123 : fockcommon%hyb_mixing_sr=abs(dtset%hyb_mixing_sr)
720 123 : fockcommon%hyb_range_dft=abs(dtset%hyb_range_dft)
721 123 : fockcommon%hyb_range_fock=abs(dtset%hyb_range_fock)
722 :
723 : ! Set the hybrid parameters if functional from libxc for which parameters can be changed, or if the user asked to do so.
724 : ! Usually, these parameters were obtained from libxc,
725 : ! but the user might have possibly modified them. By the way, must define them here for the usual changeable fonctionals,
726 : ! since otherwise might inherit them from the previous dataset !
727 123 : if(dtset%ixc<0)then
728 52 : if (dtset%ixc==-406.or.dtset%ixc==-427.or.dtset%ixc==-428 .or. &
729 : & min(dtset%hyb_mixing,dtset%hyb_mixing_sr,dtset%hyb_range_dft,dtset%hyb_range_fock)<-tol8)then
730 : call libxc_functionals_set_hybridparams(hyb_mixing=fockcommon%hyb_mixing,&
731 : & hyb_mixing_sr=fockcommon%hyb_mixing_sr,&
732 48 : & hyb_range=fockcommon%hyb_range_dft)
733 : end if
734 : end if
735 :
736 :
737 : ! ======================================================
738 : ! === Initialize the data relative to Poisson solver ===
739 : ! ======================================================
740 :
741 : !* gsqcut = cutoff value on G^2 for sphere inside the fft box (input for vhartre).
742 123 : fockcommon%gsqcut= gsqcut
743 :
744 : ! =======================================================
745 : ! === Initialize the properties of the k-points in BZ ===
746 : ! =======================================================
747 : !* Initialize %nkpt_bz = nb of k point in BZ for the calculation of exchange
748 123 : fockbz%nkpt_bz=nkpt_bz
749 : !* Initialize the array %wtk_bz = weight assigned to each k point.
750 854 : fockbz%wtk_bz=1.0_dp/dble(nkpt_bz)
751 :
752 :
753 123 : if (dtset%kptopt>=1 .and. dtset%kptopt<=4) then
754 : ! ============================================
755 : ! === Initialize the set of k-points in BZ ===
756 : ! ============================================
757 3047 : kptns_hf(:,1:nkpt_bz)=dtset%kptns_hf(:,1:nkpt_bz)
758 : !* kptns_hf contains the special k points obtained by the Monkhorst & Pack method, in reduced coordinates. (output)
759 :
760 : ! =======================================================
761 : ! === Compute the transformation to go from IBZ to BZ ===
762 : ! =======================================================
763 : !* Compute the reciprocal space metric.
764 123 : call matr3inv(rprimd,gprimd)
765 4920 : gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
766 :
767 : !* Calculate the array indkk which describes how to get IBZ from BZ
768 : !* 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)
769 : !* sppoldbl=1, no spin-polarisation doubling is required.
770 123 : timrev=1 ; if (dtset%kptopt==3 .or. dtset%kptopt==4) timrev=0
771 : !* timrev=1 if the use of time-reversal is allowed ; 0 otherwise
772 123 : if (dtset%kptopt==2 .or. dtset%kptopt==3) then
773 : !* No space symmetry is used, if kptopt==2 time reversal symmetry is used.
774 0 : symm=0 ; symm(1,1)=1 ; symm(2,2)=1 ; symm(3,3)=1
775 : call listkk(dksqmax,gmet,indkk(1:nkpt_bz,:),dtset%kptns,kptns_hf,dtset%nkpt, &
776 0 : & nkpt_bz,1,1,indx,symm,timrev,xmpi_comm_self)
777 : else
778 : !* As in getkgrid, no use of antiferromagnetic symmetries thans to the option sppoldbl=1
779 : call listkk(dksqmax,gmet,indkk(1:nkpt_bz,:),dtset%kptns,kptns_hf,dtset%nkpt, &
780 123 : & nkpt_bz,dtset%nsym,1,dtset%symafm,dtset%symrel,timrev, xmpi_comm_self)
781 : end if
782 : !* indkk(nkpt_bz,6) describes the k point of IBZ that generates each k point of BZ
783 : !* indkk(:,1) = k point of IBZ, kpt_ibz
784 : !* indkk(:,2) = symmetry operation to apply to kpt_ibz to give the k point of BZ
785 : !* (if 0, means no symmetry operation, equivalent to identity )
786 : !* indkk(:,3:5) = Umklapp vectors to apply to remain in BZ
787 : !* indkk(:,6) = 1 if time-reversal was used to generate the k point of BZ, 0 otherwise
788 : !* No use of symafm to generate spin down wfs from spin up wfs for the moment
789 :
790 : else
791 0 : if (dtset%kptopt==0) then
792 : !* kptopt =0 : read directly nkpt, kpt, kptnrm and wtk in the input file
793 : !* => this case is not allowed for the moment
794 0 : ABI_ERROR('Hartree-Fock option can not be used with option kptopt=0.')
795 : else
796 : !* kptopt <0 : rely on kptbounds, and ndivk to set up a band structure calculation
797 : !* => a band structure calculation is not yet allowed.
798 0 : ABI_ERROR('Hartree-Fock option can not be used with option kptopt<0.')
799 : end if
800 : end if
801 :
802 : !! =======================================================
803 : !! === Initialize the properties of the k-points in BZ ===
804 : !! =======================================================
805 : ! jkg=0
806 : !!* Initialize the arrays %npwarr_bz, %kg_j, %phase_j, %gbound_j
807 : ! do jkpt=1,nkpt_bz
808 : ! ikpt=indkk(jkpt,1)
809 : !!* ikpt = the point of IBZ that jkpt is an image of in BZ
810 : ! npwj=npwarr(ikpt)
811 : !!* npwj = number of planewaves in basis at point jkpt = at point ikpt
812 : ! jsym=indkk(jkpt,2)
813 : !!* jsym = symmetry operation to apply to get jkpt from ikpt
814 : ! shiftg(:)=indkk(jkpt,3:5)
815 : !!* shiftg = Bravais vector G0 to add to remain in BZ
816 : ! if (jsym/=0) then
817 : ! symm(:,:)=dtset%symrel(:,:,jsym)
818 : ! tau_nons(:)=dtset%tnons(:,jsym)
819 : !!* The symmetry operation in k-space (symm) and the non-symorphic translation (tau_nons) are now defined.
820 : ! if(sum(tau_nons(:)**2)>tol8) then
821 : !!* Initialize %calc_phase(jkpt) to 1
822 : ! fock%calc_phase(jkpt)=1
823 : !!* Compute the phase factor exp(i*2*pi*G.tau) for all G.
824 : ! indx(1)=1
825 : ! phase1d=zero
826 : ! call getph(indx,1,n1,n2,n3,phase1d,tau_nons)
827 : !!* Although the routine getph is orignally written for atomic phase factors, it does precisely what we want
828 : ! arg=two_pi*(dtset%kptns(1,ikpt)*tau_nons(1) + dtset%kptns(2,ikpt)*tau_nons(2) &
829 : !& + dtset%kptns(3,ikpt)*tau_nons(3))
830 : ! phktnons(1,1)=cos(arg)
831 : ! phktnons(2,1)=sin(arg)
832 : !! phktnons(1,1)=one
833 : !! phktnons(2,1)=zero
834 : !!* Convert 1D phase factors to 3D phase factors exp(i*2*pi*(k+G).tau) and store it in %phase_j
835 : ! call ph1d3d(1,1,kg(:,1+tab_indikpt(1,ikpt):npwj+tab_indikpt(1,ikpt)),1,1,npwj,n1, &
836 : !& n2,n3,phktnons,phase1d,fock%phase(:,1+jkg:npwj+jkg))
837 : ! end if
838 : ! else
839 : ! symm=0 ; symm(1,1)=1 ; symm(2,2)=1 ; symm(3,3)=1
840 : ! tau_nons(:)=zero
841 : ! shiftg(:)=0
842 : ! end if
843 : !!* Apply time-reversal symmetry if required
844 : ! if(indkk(jkpt,6)/=0) then
845 : !!* Initialize %timerev(jkpt) to 1
846 : ! fock%timerev(jkpt)=1
847 : ! symm(:,:)=-symm(:,:)
848 : ! end if
849 :
850 : !!* Initialize %istwfk_bz(jkpt) to
851 : ! fock%istwfk_bz(jkpt)=dtset%istwfk(ikpt)
852 :
853 : !!* Initialize %tab_ikpt and %tab_ibgcg
854 : ! fock%tab_ikpt(jkpt)=ikpt
855 : ! fock%tab_ibgcg(1:dtset%nsppol,jkpt)=tab_indikpt(2:1+dtset%nsppol,ikpt)
856 : ! fock%tab_ibgcg(1+dtset%nsppol:2*dtset%nsppol,jkpt)= &
857 : !& tab_indikpt(2+dtset%nsppol:2*dtset%nsppol+1,ikpt)
858 :
859 : !!* Initialize %npwarr_bz
860 : ! fock%npwarr_bz(jkpt)=npwj
861 :
862 : !!* Initialize %kg_bz
863 : ! do jpw=1,npwj
864 : ! v1=kg(1,jpw+tab_indikpt(1,ikpt)) ; v2=kg(2,jpw+tab_indikpt(1,ikpt)) ; v3=kg(3,jpw+tab_indikpt(1,ikpt))
865 : ! fock%kg_bz(1,jpw+jkg)=-shiftg(1)+symm(1,1)*v1+symm(2,1)*v2+symm(3,1)*v3
866 : ! fock%kg_bz(2,jpw+jkg)=-shiftg(2)+symm(1,2)*v1+symm(2,2)*v2+symm(3,2)*v3
867 : ! fock%kg_bz(3,jpw+jkg)=-shiftg(3)+symm(1,3)*v1+symm(2,3)*v2+symm(3,3)*v3
868 : !!* The symmetry operation symm must be transposed when used. (cf. docs about wfconv)
869 : ! end do
870 :
871 : !!* Initialize %gbound_bz
872 : ! call sphereboundary(fock%gbound_bz(:,:,jkpt),fock%istwfk_bz(jkpt), &
873 : !& fock%kg_bz(:,1+jkg:npwj+jkg),dtset%mgfft,npwj)
874 :
875 : !!* Update of the shift to be applied
876 : ! jkg=jkg+npwj
877 : ! end do
878 :
879 : ! ==========================================================
880 : ! === Initialize the k-points in BZ and their properties ===
881 : ! ==========================================================
882 : ! jkg=0;
883 :
884 9627 : do isym=1,dtset%nsym
885 9504 : call mati3inv(dtset%symrel(:,:,isym),symrec(:,:,isym))
886 266235 : Rtnons (:,isym)= MATMUL(TRANSPOSE(symrec(:,:,isym)),dtset%tnons(:,isym))
887 : end do
888 369 : ABI_MALLOC(fockcommon%symrec,(3,3,dtset%nsym))
889 123798 : fockcommon%symrec=symrec
890 :
891 369 : ABI_MALLOC(invsym,(dtset%nsym))
892 9627 : invsym=0
893 492 : ident(1,:3)=(/1,0,0/)
894 492 : ident(2,:3)=(/0,1,0/)
895 492 : ident(3,:3)=(/0,0,1/)
896 9627 : do isym=1,dtset%nsym
897 636768 : symm(:,:)=MATMUL(dtset%symrel(:,:,isym),dtset%symrel(:,:,isym))
898 61206 : if (all(symm(:,:)==ident(:,:))) then
899 4042 : invsym(isym)=isym
900 : else
901 765542 : do jsym=1,dtset%nsym
902 50925360 : symm(:,:)=MATMUL(dtset%symrel(:,:,isym),dtset%symrel(:,:,jsym))
903 1562638 : if (all(symm(:,:)==ident(:,:))) then
904 15878 : invsym(isym)=jsym
905 15878 : cycle
906 : end if
907 : end do
908 : end if
909 9627 : if(invsym(isym)==0) then
910 0 : ABI_ERROR('No inverse has been found for isym')
911 : end if
912 : end do
913 :
914 123 : jkg_this_proc=0;my_jkpt=0
915 : !indkk(1:nkpt_bz,2)=(/1,1,3,1,11,7,9,1/)
916 854 : do jkpt=1,nkpt_bz
917 :
918 : !* If this processor does not calculate exchange with the k point jkpt, skip the rest of the k-point loop.
919 731 : if (proc_distrb_cycle(mpi_enreg%distrb_hf,jkpt,1,mband,1,mpi_enreg%me_hf)) cycle
920 : ! if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,jkpt,1,dtset%nbandhf,1,mpi_enreg%me_kpt))) then
921 : !* The processor does own a copy of the array kg of ikpt ; increment the shift.
922 : ! jkg=jkg+npwj
923 : ! end if
924 : ! Skip the rest of the k-point loop
925 : ! cycle
926 : ! end if
927 731 : my_jkpt=my_jkpt+1
928 :
929 731 : ikpt=indkk(jkpt,1)
930 : !* ikpt = the point of IBZ that jkpt is an image of in BZ
931 731 : npwj=npwarr(ikpt)
932 731 : fockbz%npwarr(my_jkpt)=npwarr(ikpt)
933 : !* npwj = number of planewaves in basis at point jkpt = at point ikpt
934 731 : jsym=indkk(jkpt,2)
935 : !* jsym = symmetry operation to apply to get jkpt from ikpt
936 731 : fockbz%tab_symkpt(my_jkpt)=invsym(jsym)
937 2924 : shiftg(:)=indkk(jkpt,3:5)
938 : !* shiftg = Bravais vector G0 to add to remain in BZ
939 :
940 : !* Initialize the array %kptns_bz = the k points in full BZ
941 2924 : fockbz%kptns_bz(:,my_jkpt)=kptns_hf(:,jkpt)
942 :
943 : !* Initialize the array %jstwfk = how is stored the wavefunction at each k point
944 731 : if (dtset%istwfk(ikpt)/=1) then
945 1 : fockbz%istwfk_bz(my_jkpt)=set_istwfk(kptns_hf(:,jkpt))
946 : end if
947 :
948 : !* One can take advantage of the time-reversal symmetry in this case.
949 : !* Initialize the array %wtk_bz = weight assigned to each k point.
950 : ! fock%wtk_bz(my_jkpt)=dtset%wtk(jkpt)/ucvol
951 : !* Caution, the definition takes into account "ucvol" !
952 :
953 : !* Initialize the array %npwarr_bz = number of planewaves in basis at each k point
954 : ! fock%npwarr_bz(my_jkpt)=npwj
955 :
956 : !!* Initialize the array %tab_ikpt = indices of k-point in IBZ ikpt for each k point jkpt in BZ (here,ikpt=jkpt)
957 731 : fockbz%tab_ikpt(my_jkpt)=ikpt
958 :
959 :
960 : !!* Initialize the array %tab_ibgcg = indices of cprj(ikpt)/occ(ikpt) and cg(ikpt) for each k point jkpt
961 : ! if (my_nsppol==2) then
962 : !!* In this case, my_nsppol=dtset%nsppol=2
963 : ! fock%tab_ibgcg(1:2,my_jkpt)=tab_indikpt(2:3,ikpt)
964 : ! fock%tab_ibgcg(3:4,my_jkpt)=tab_indikpt(4:5,ikpt)
965 : ! else
966 : ! if(mpi_enreg%my_isppoltab(1)==1) then
967 : !!* In this case, my_nsppol=1 and the up spin is treated (dtset%nsppol= 1 or 2)
968 : ! fock%tab_ibgcg(1,my_jkpt)=tab_indikpt(2,ikpt)
969 : ! fock%tab_ibgcg(2,my_jkpt)=tab_indikpt(2+dtset%nsppol,ikpt)
970 : ! else
971 : !!* In this case, my_nsppol=1 and the dn spin is treated (so dtset%nsppol=2)
972 : ! fock%tab_ibgcg(1,my_jkpt)=tab_indikpt(3,ikpt)
973 : ! fock%tab_ibgcg(2,my_jkpt)=tab_indikpt(5,ikpt)
974 : ! end if
975 : ! end if
976 :
977 : !* Initialize the array %kg_bz = reduced planewave coordinates at each k point
978 731 : if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,-1,mpi_enreg%me_kpt))) then
979 : !* We perform the test with isppol=-1 (both spins) and the occupied band (dtset%nbandhf).
980 : !* We assume that paral_kgb==0 (a k-point may not be present on several proc.)
981 : !* The array kg for ikpt is stored on this processor and copied in kg_tmp.
982 687 : ikg=my_ikgtab(ikpt)
983 : !* ikg = the shift in kg to get the G-vectors associated to ikpt
984 2748 : do ik=1,3
985 : ! kg_tmp(1+(ik-1)*npwj:ik*npwj)=kg(ik,1+tab_indikpt(1,ikpt):npwj+tab_indikpt(1,ikpt))
986 378546 : kg_tmp(1+(ik-1)*npwj:ik*npwj)=kg(ik,1+ikg:npwj+ikg)
987 : end do
988 : ! jkg=jkg+npwj
989 : end if
990 : !* Broadcast the array kg_tmp to all the processors of comm_kpt.
991 : !* Since paral_kgb==0, all the bands of a k-point are treated on the same proc.
992 731 : call xmpi_bcast(kg_tmp,mpi_enreg%proc_distrb(ikpt,1,1),mpi_enreg%comm_kpt,ier)
993 2924 : do ik=1,3
994 419666 : fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)=kg_tmp(1+(ik-1)*npwj:ik*npwj)
995 : end do
996 :
997 : !* Apply a symmetry operation on kg_bz if necessary
998 731 : if (jsym/=0) then
999 9503 : symm(:,:)=dtset%symrel(:,:,jsym)
1000 : ! tau_nons(:)=dtset%tnons(:,jsym)
1001 2924 : tau_nons(:)=-Rtnons(:,invsym(jsym))
1002 : !* The symmetry operation in k-space (symm) and the non-symorphic translation (tau_nons) are now defined.
1003 2924 : if(sum(tau_nons(:)**2)>tol8) then
1004 : !* Initialize %calc_phase(jkpt) to 1
1005 193 : fockbz%calc_phase(my_jkpt)=1
1006 : !* Compute the phase factor exp(i*2*pi*G.tau) for all G.
1007 193 : indx(1)=1
1008 20042938 : phase1d=zero
1009 193 : call getph(indx,1,n1,n2,n3,phase1d,tau_nons)
1010 : !* Although the routine getph is orignally written for atomic phase factors, it does precisely what we want
1011 : arg=two_pi*(dtset%kptns(1,ikpt)*tau_nons(1) + dtset%kptns(2,ikpt)*tau_nons(2) &
1012 193 : & + dtset%kptns(3,ikpt)*tau_nons(3))
1013 193 : phktnons(1,1)=cos(arg)
1014 193 : phktnons(2,1)=sin(arg)
1015 : ! phktnons(1,1)=one
1016 : ! phktnons(2,1)=zero
1017 : !* Convert 1D phase factors to 3D phase factors exp(i*2*pi*(k+G).tau) and store it in %phase_j
1018 : call ph1d3d(1,1,fockbz%kg_bz(:,1+jkg_this_proc:npwj+jkg_this_proc),1,1,npwj,n1,n2,n3, &
1019 193 : & phktnons,phase1d,fockbz%phase(:,1+jkg_this_proc:npwj+jkg_this_proc))
1020 : end if
1021 : !* Apply time-reversal symmetry if required
1022 731 : if(indkk(jkpt,6)/=0) then
1023 : !* Initialize %timerev(jkpt) to 1
1024 0 : fockbz%timerev(my_jkpt)=1
1025 0 : symm(:,:)=-symm(:,:)
1026 : end if
1027 : !* Initialize %kg_bz
1028 139645 : do jpw=1,npwj
1029 138914 : 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)
1030 138914 : fockbz%kg_bz(1,jpw+jkg_this_proc)=-shiftg(1)+symm(1,1)*v1+symm(2,1)*v2+symm(3,1)*v3
1031 138914 : fockbz%kg_bz(2,jpw+jkg_this_proc)=-shiftg(2)+symm(1,2)*v1+symm(2,2)*v2+symm(3,2)*v3
1032 139645 : fockbz%kg_bz(3,jpw+jkg_this_proc)=-shiftg(3)+symm(1,3)*v1+symm(2,3)*v2+symm(3,3)*v3
1033 : !* The symmetry operation symm must be transposed when used. (cf. docs about wfconv)
1034 : end do
1035 : else
1036 : !* Ths symmetry operation is the identity.
1037 : !* Apply time-reversal symmetry if required
1038 0 : if(indkk(jkpt,6)/=0) then
1039 : !* Initialize %timerev(jkpt) to 1
1040 0 : fockbz%timerev(my_jkpt)=1
1041 0 : fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)=-fockbz%kg_bz(ik,1+jkg_this_proc:npwj+jkg_this_proc)
1042 : end if
1043 : end if
1044 :
1045 : !* Initialize the array %gbound_bz = boundary of the basis sphere of G vectors at each k point
1046 : call sphereboundary(fockbz%gbound_bz(:,:,my_jkpt),fockbz%istwfk_bz(my_jkpt),&
1047 731 : & fockbz%kg_bz(:,1+jkg_this_proc:npwj+jkg_this_proc),mgfft,npwj)
1048 :
1049 731 : jkg_this_proc=jkg_this_proc+npwj
1050 :
1051 : !* Initialize the arrays %tab_ibg = shifts in arrays cprj and occ (ibg) for each k point jkpt
1052 : !* Initialize the arrays %tab_icg = shifts in arrays cg(icg) for each k point jkpt
1053 1585 : if (my_nsppol==1) then
1054 691 : fockbz%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,1+mpi_enreg%my_isppoltab(2))
1055 691 : fockbz%tab_icg(my_jkpt,1)=my_icgtab(ikpt,1+mpi_enreg%my_isppoltab(2))
1056 691 : fockbz%tab_icp(my_jkpt,1)=my_icptab(ikpt,1+mpi_enreg%my_isppoltab(2))
1057 : !* if mpy_isppoltab(2)=0, the up spin is treated (dtset%nsppol= 1 or 2)
1058 : !* if mpy_isppoltab(2)=1, the dn spin is treated (so dtset%nsppol=2)
1059 :
1060 : ! if(mpi_enreg%my_isppoltab(2)==1) then
1061 : !* In this case, my_nsppol=1 and the up spin is treated (dtset%nsppol= 1 or 2)
1062 : ! fock%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,1)
1063 : ! fock%tab_icg(my_jkpt,1)=my_icgtab(ikpt,1)
1064 : ! else
1065 : !* In this case, my_nsppol=1 and the dn spin is treated (so dtset%nsppol=2)
1066 : ! fock%tab_ibg(my_jkpt,1)=my_ibgtab(ikpt,2)
1067 : ! fock%tab_icg(my_jkpt,1)=my_icgtab(ikpt,2)
1068 : ! end if
1069 : else
1070 : !* In this case, my_nsppol=dtset%nsppol=2
1071 120 : fockbz%tab_ibg(my_jkpt,:)=my_ibgtab(ikpt,:)
1072 120 : fockbz%tab_icg(my_jkpt,:)=my_icgtab(ikpt,:)
1073 120 : fockbz%tab_icp(my_jkpt,:)=my_icptab(ikpt,:)
1074 : end if
1075 :
1076 : enddo
1077 :
1078 : !* Deallocation
1079 123 : ABI_FREE(invsym)
1080 :
1081 : end if
1082 123 : ABI_FREE(indkk)
1083 123 : ABI_FREE(kg_tmp)
1084 123 : ABI_FREE(kptns_hf)
1085 123 : ABI_FREE(my_ibgtab)
1086 123 : ABI_FREE(my_icgtab)
1087 123 : ABI_FREE(my_icptab)
1088 123 : ABI_FREE(my_ikgtab)
1089 123 : ABI_FREE(phase1d)
1090 123 : call fock_print(fockcommon,fockbz,unit=std_out)
1091 :
1092 123 : call timab(1501,2,tsec)
1093 :
1094 : DBG_EXIT("COLL")
1095 :
1096 123 : end subroutine fock_init
1097 : !!***
1098 :
1099 : !!****f* m_fock/fock_updateikpt
1100 : !! NAME
1101 : !! fock_updateikpt
1102 : !!
1103 : !! FUNCTION
1104 : !! Update the value of ikpt,isppol for the next exact exchange calculation.
1105 : !!
1106 : !! INPUTS
1107 : !! fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
1108 : !! ikpt= k-point index
1109 : !! isppol= Spin index
1110 : !!
1111 : !! SIDE EFFECTS
1112 : !! The field fock%eigen_ikpt is also set to 0.d0.
1113 : !!
1114 : !! NOTES
1115 : !! May be improved to calculate the star of ikpt. => I think NO finally
1116 : !!
1117 : !! SOURCE
1118 :
1119 6063 : subroutine fock_updateikpt(fock,ikpt,isppol)
1120 :
1121 : !Arguments ------------------------------------
1122 : type(fock_common_type),pointer :: fock
1123 : integer, intent(in) :: ikpt,isppol
1124 : ! *************************************************************************
1125 :
1126 : !write (std_out,*) ' fock_updateikpt : enter'
1127 :
1128 : ! ======================================================
1129 : ! === Update the data relative to the current states ===
1130 : ! ======================================================
1131 : !* Copy of the value ikpt in the field ikpt
1132 6063 : fock%ikpt=ikpt
1133 : !* Copy of the value isppol in the field isppol
1134 6063 : fock%isppol=isppol
1135 : !* Set all the Fock contributions to the eigenvalues to 0.d0.
1136 35071 : fock%eigen_ikpt=zero
1137 : !* Set all the Fock contributions to the forces to 0.d0.
1138 6063 : if ((fock%optfor).and.(fock%use_ACE==0)) then
1139 0 : fock%forces_ikpt=zero
1140 : end if
1141 :
1142 6063 : end subroutine fock_updateikpt
1143 : !!***
1144 :
1145 : !!****f* m_fock/fock_set_ieigen
1146 : !! NAME
1147 : !! fock_set_ieigen
1148 : !!
1149 : !! FUNCTION
1150 : !! Set the value of ieigen to the value given in argument.
1151 : !!
1152 : !! INPUTS
1153 : !! fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
1154 : !! iband= index of the band iband
1155 : !!
1156 : !! OUTPUT
1157 : !! none
1158 : !!
1159 : !! SOURCE
1160 :
1161 3149627 : subroutine fock_set_ieigen(fock,iband)
1162 :
1163 : !Arguments ------------------------------------
1164 : type(fock_common_type),pointer :: fock
1165 : integer, intent(in) :: iband
1166 : ! *************************************************************************
1167 :
1168 : !Nothing to do if fock pointer is not associated...
1169 :
1170 : ! ======================================================
1171 : ! === Update the data relative to the current states ===
1172 : ! ======================================================
1173 :
1174 : !* Copy of the value iband in the field ieigen
1175 3149627 : if (associated(fock)) then
1176 34259 : fock%ieigen=iband
1177 34259 : fock%iband=iband
1178 : end if
1179 :
1180 3149627 : end subroutine fock_set_ieigen
1181 : !!***
1182 :
1183 : !!****f* m_fock/fock_destroy
1184 : !! NAME
1185 : !! fock_destroy
1186 : !!
1187 : !! FUNCTION
1188 : !! Clean and destroy fock datastructure.
1189 : !!
1190 : !! INPUTS
1191 : !! fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
1192 : !!
1193 : !! SOURCE
1194 123 : subroutine fock_destroy(fock)
1195 :
1196 : !Arguments ------------------------------------
1197 : type(fock_type),pointer :: fock
1198 : ! *************************************************************************
1199 :
1200 123 : if (fock%fock_common%use_ACE/=0) then
1201 123 : ABI_FREE(fock%fockACE)
1202 : end if
1203 123 : ABI_FREE(fock%fock_common)
1204 123 : ABI_FREE(fock%fock_BZ)
1205 123 : ABI_FREE(fock)
1206 :
1207 123 : end subroutine fock_destroy
1208 : !!***
1209 :
1210 123 : subroutine fock_common_destroy(fock)
1211 :
1212 : !Arguments ------------------------------------
1213 : type(fock_common_type),pointer :: fock
1214 : ! *************************************************************************
1215 :
1216 : DBG_ENTER("COLL")
1217 :
1218 123 : ABI_SFREE(fock%atindx)
1219 123 : ABI_SFREE(fock%typat)
1220 :
1221 : ! real arrays
1222 123 : ABI_SFREE(fock%forces)
1223 123 : ABI_SFREE(fock%nband)
1224 123 : ABI_SFREE(fock%forces_ikpt)
1225 123 : ABI_SFREE(fock%stress_ikpt)
1226 123 : ABI_SFREE(fock%eigen_ikpt)
1227 :
1228 : ! Deallocate datatypes
1229 123 : if (allocated(fock%pawfgrtab)) then
1230 26 : call pawfgrtab_free(fock%pawfgrtab)
1231 78 : ABI_FREE(fock%pawfgrtab)
1232 : end if
1233 :
1234 : ! Put the integer to 0
1235 123 : fock%ieigen=0
1236 123 : fock%ikpt=0
1237 123 : fock%isppol=0
1238 :
1239 123 : ABI_SFREE(fock%symrec)
1240 :
1241 : !* [description of divergence in |q+G|=0]
1242 : !* Put the real (dp) to 0
1243 123 : fock%gsqcut=zero
1244 123 : fock%hyb_mixing=zero
1245 123 : fock%hyb_mixing_sr=zero
1246 123 : fock%hyb_range_dft=zero
1247 123 : fock%hyb_range_fock=zero
1248 :
1249 : DBG_EXIT("COLL")
1250 :
1251 123 : end subroutine fock_common_destroy
1252 : !!***
1253 :
1254 123 : subroutine fock_BZ_destroy(fock)
1255 :
1256 : !Arguments ------------------------------------
1257 : type(fock_BZ_type),pointer :: fock
1258 : ! *************************************************************************
1259 :
1260 : DBG_ENTER("COLL")
1261 :
1262 123 : ABI_SFREE(fock%cwaveocc_bz)
1263 123 : ABI_SFREE(fock%cgocc)
1264 123 : ABI_SFREE(fock%npwarr)
1265 123 : ABI_SFREE(fock%occ_bz)
1266 123 : if (allocated(fock%cwaveocc_prj)) then
1267 26 : call pawcprj_free(fock%cwaveocc_prj)
1268 1238 : ABI_FREE(fock%cwaveocc_prj)
1269 : end if
1270 : ! Deallocate integer arrays
1271 :
1272 123 : ABI_SFREE(fock%kg_bz)
1273 123 : ABI_SFREE(fock%nbandocc_bz)
1274 123 : ABI_SFREE(fock%istwfk_bz)
1275 123 : ABI_SFREE(fock%calc_phase)
1276 123 : ABI_SFREE(fock%timerev)
1277 123 : ABI_SFREE(fock%tab_ibg)
1278 123 : ABI_SFREE(fock%tab_icg)
1279 123 : ABI_SFREE(fock%tab_icp)
1280 123 : ABI_SFREE(fock%tab_ikpt)
1281 123 : ABI_SFREE(fock%tab_symkpt)
1282 :
1283 : !* [description of IBZ and BZ]
1284 : !* Deallocate real arrays
1285 123 : ABI_SFREE(fock%wtk_bz)
1286 123 : ABI_SFREE(fock%kptns_bz)
1287 123 : ABI_SFREE(fock%phase)
1288 : !* Put the integer to 0
1289 123 : fock%nkpt_bz=0
1290 :
1291 : !* Deallocate real arrays
1292 :
1293 : !* Deallocate integer arrays
1294 123 : ABI_SFREE(fock%gbound_bz)
1295 :
1296 : !* [description of size of arrays/pointers]
1297 : !* Put the integer to 0
1298 123 : fock%mkpt=0
1299 123 : fock%mkptband=0
1300 123 : call destroy_mpi_enreg(fock%mpi_enreg)
1301 :
1302 : DBG_EXIT("COLL")
1303 :
1304 123 : end subroutine fock_BZ_destroy
1305 : !!***
1306 :
1307 : !!****f* m_fock/fock_ACE_destroy
1308 : !! NAME
1309 : !! fock_ACE_destroy
1310 : !!
1311 : !! FUNCTION
1312 : !! Clean and destroy fock datastructure.
1313 : !!
1314 : !! INPUTS
1315 : !! fockACE <type(fock_ACE_type)>= all the quantities to calculate Fock exact exchange in the ACE context
1316 : !!
1317 : !! SOURCE
1318 :
1319 123 : subroutine fock_ACE_destroy(fockACE)
1320 :
1321 : !Arguments ------------------------------------
1322 : type(fock_ACE_type),pointer :: fockACE(:,:)
1323 :
1324 : !Local variables-------------------------------
1325 : integer :: dim1,dim2,ii,jj
1326 : ! *************************************************************************
1327 :
1328 : DBG_ENTER("COLL")
1329 :
1330 123 : dim1=size(fockACE,1)
1331 123 : dim2=size(fockACE,2)
1332 260 : do jj=1,dim2
1333 957 : do ii=1,dim1
1334 834 : if (allocated(fockACE(ii,jj)%xi)) then
1335 697 : ABI_FREE(fockACE(ii,jj)%xi)
1336 : end if
1337 : end do
1338 : end do
1339 : DBG_EXIT("COLL")
1340 :
1341 123 : end subroutine fock_ACE_destroy
1342 : !!***
1343 :
1344 :
1345 : !!****f* m_fock/fock_calc_ene
1346 : !! NAME
1347 : !! fock_calc_ene
1348 : !!
1349 : !! FUNCTION
1350 : !! Calculate the Fock contribution to the total energy
1351 : !!
1352 : !! INPUTS
1353 : !! fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
1354 : !! ikpt= reduced planewave coordinates.
1355 : !!
1356 : !! OUTPUT
1357 : !! none
1358 : !!
1359 : !! SIDE EFFECTS
1360 : !! energies <type(energies_type)>=storage for energies computed here :
1361 : !! | e_exactX = Fock contribution to the total energy (Hartree)
1362 : !!
1363 : !! NOTES
1364 : !! If the cgocc_bz are not updated at each iteration, be careful to calculate Fock energy at the same frequency.
1365 : !! TO CHECK == CHANGE IN SOME DEFINTIONS
1366 : !!
1367 : !! SOURCE
1368 :
1369 0 : subroutine fock_calc_ene(dtset,fock,fock_energy,ikpt,nband,occ)
1370 :
1371 : !Arguments ------------------------------------
1372 : !scalars
1373 : integer,intent(in) :: ikpt,nband
1374 : real(dp),intent(inout) :: fock_energy
1375 : type(dataset_type),intent(in) :: dtset
1376 : type(fock_common_type),pointer :: fock
1377 : !arrays
1378 : real(dp),intent(in) :: occ(nband)
1379 :
1380 : !Local variables-------------------------------
1381 : integer :: iband
1382 : ! *************************************************************************
1383 :
1384 : ABI_UNUSED(fock_energy)
1385 :
1386 0 : do iband=1,nband
1387 :
1388 : ! Select only the occupied states (such that fock%occ_bz > 10^-8)
1389 0 : if (abs(occ(iband))>tol8) then
1390 : ! fock_energy=fock_energy + half*fock%eigen_ikpt(iband)*occ(iband)*dtset%wtk(ikpt)
1391 : !* Sum the contribution of each occupied states at point k_i
1392 : !* No need to multiply %wtk by ucvol since there is no factor 1/ucvol in the definition of %wtk
1393 :
1394 : !* accumulate Fock contributions to the forces.
1395 : ! if (fock%optfor) then
1396 0 : fock%forces(:,:)=fock%forces(:,:)+occ(iband)*dtset%wtk(ikpt)*fock%forces_ikpt(:,:,iband)
1397 : ! end if
1398 : end if
1399 : end do
1400 :
1401 0 : end subroutine fock_calc_ene
1402 : !!***
1403 :
1404 : !!****f* m_fock/fock_update_exc
1405 : !! NAME
1406 : !! fock_update_exc
1407 : !!
1408 : !! FUNCTION
1409 : !! Update the value of energies%e_xc and energies%e_xcdc with Fock contribution
1410 : !!
1411 : !! INPUTS
1412 : !!
1413 : !! OUTPUT
1414 : !! none
1415 : !!
1416 : !! energies <type(energies_type)>=storage for energies computed here :
1417 : !! | e_fock= Fock contribution to the total energy (Hartree)
1418 : !!
1419 : !! NOTES
1420 : !! If the cgocc_bz are not updated at each iteration, be careful to calculate Fock energy at the same frequency.
1421 : !!
1422 : !! SOURCE
1423 :
1424 0 : subroutine fock_update_exc(fock_energy,xc_energy,xcdc_energy)
1425 :
1426 : !Arguments ------------------------------------
1427 : real(dp),intent(in) :: fock_energy
1428 : real(dp),intent(inout) :: xc_energy,xcdc_energy
1429 : ! *************************************************************************
1430 :
1431 : !xc_energy = fock%hyb_mixing*fock_energy
1432 : !xcdc_energy = two*fock%hyb_mixing*fock_energy
1433 0 : xc_energy = fock_energy
1434 0 : xcdc_energy = two*fock_energy
1435 : !CMartins : For an atom, ewald should be set to zero (at the beginning of the loop) and
1436 : !the contribution in !|q+G|=0 should be an approximation to the missing component of Vloc in G=0
1437 : !energies%e_ewald=energies%e_ewald-half*fock%divgq0*fock%wtk_bz(1)*piinv
1438 :
1439 0 : end subroutine fock_update_exc
1440 : !!***
1441 :
1442 : !----------------------------------------------------------------------
1443 :
1444 : !!****f* m_fock/fock_updatecwaveocc
1445 : !! NAME
1446 : !! fock_updatecwaveocc
1447 : !!
1448 : !! FUNCTION
1449 : !! Update in the fock datastructure the fields relative to the occupied states.
1450 : !!
1451 : !! INPUTS
1452 : !! cg(2,mcg)= Input wavefunctions
1453 : !! cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
1454 : !! dtset <type(dataset_type)>=all input variables for this dataset
1455 : !! fock <type(fock_type)>= all the quantities to calculate Fock exact exchange
1456 : !! indsym(4,nsym,natom) :: 1:3 shift, and 4 final atom, of symmetry isym operating on iatom
1457 : !! (S^{-1}(R - t) = r0 + L, see symatm.F90
1458 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
1459 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
1460 : !! mpi_enreg=information about MPI parallelization
1461 : !! npwarr(nkpt)=number of planewaves in basis at this k point
1462 : !! occ(mband*nkpt*nsppol)= occupation number for each band (often 2) at each k point
1463 : !! ucvol= unit cell volume ($\textrm{bohr}^{3}$)
1464 : !!
1465 : !! OUTPUT
1466 : !!
1467 : !! SIDE EFFECTS
1468 : !! The field fock%cgocc_bz contains the table cg at the end.
1469 : !! The fields kg_bz, occ_bz and fock%cwaveocc_prj are simultaneously updated.
1470 : !!
1471 : !! NOTES
1472 : !!
1473 : !! ############################
1474 : !! ### Not fully tested yet ###
1475 : !! ############################
1476 : !!
1477 : !! May be improved by selecting only the occupied states with the same spin isppol.
1478 : !!
1479 : !! SOURCE
1480 :
1481 1238 : subroutine fock_updatecwaveocc(cg,cprj,dtset,fock,indsym,mcg,mcprj,&
1482 1238 : mpi_enreg,nattyp,npwarr,occ,ucvol)
1483 :
1484 : !scalars
1485 : integer, intent(in) :: mcg,mcprj
1486 : real(dp), intent(in) :: ucvol
1487 : type(dataset_type),intent(in) :: dtset
1488 : type(fock_type),intent(inout),pointer :: fock
1489 : type(MPI_type),intent(in) :: mpi_enreg
1490 : !arrays
1491 : integer, intent(in) :: indsym(4,dtset%nsym,dtset%natom),nattyp(dtset%ntypat),npwarr(dtset%nkpt)
1492 : real(dp),intent(in) :: cg(2,mcg),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
1493 : type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj)
1494 :
1495 : !Local variables-------------------------------
1496 : !scalars
1497 : integer,parameter :: tim_fourwf0=0
1498 : 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
1499 : integer :: lmnmax,mband,mband0,mgfft,mkpt,mpw,my_jsppol,my_jband,my_jkpt
1500 : integer :: nband,ncpgr,n4,n5,n6,nkpt_bz,npwj,nsppol,nspinor
1501 : real(dp),parameter :: weight1=one
1502 : real(dp) :: cgre,cgim,invucvol
1503 : character(len=500) :: message
1504 : ! arrays
1505 : integer :: ngfft(18)
1506 1238 : integer, contiguous, pointer :: gbound_k(:,:),kg_k(:,:)
1507 1238 : integer,allocatable :: dimlmn(:),indlmn(:,:,:),indsym_(:,:,:),typat_srt(:)
1508 : real(dp) :: tsec(2),tsec2(2),dcp(3)
1509 1238 : real(dp),allocatable :: cgocc_tmp(:),cgocc(:,:),dummytab2(:,:),dummytab3(:,:,:),phase_jkpt(:,:)
1510 1238 : type(pawcprj_type),allocatable :: cprj_tmp(:,:)
1511 : type(fock_common_type),pointer :: fockcommon
1512 : type(fock_BZ_type),pointer :: fockbz
1513 : ! *************************************************************************
1514 :
1515 1238 : call timab(1502,1,tsec)
1516 :
1517 1238 : ABI_CHECK(associated(fock),"fock must be associated")
1518 :
1519 1238 : if (associated(fock)) then
1520 :
1521 1238 : fockcommon=>fock%fock_common
1522 1238 : fockbz=> fock%fock_BZ
1523 :
1524 1238 : invucvol=1.d0/sqrt(ucvol)
1525 : ! Local variables = useful dimensions
1526 1238 : mband=fockcommon%mband
1527 1238 : mkpt=fockbz%mkpt
1528 1238 : mpw=dtset%mpw
1529 1238 : mgfft=dtset%mgfft
1530 23522 : ngfft=dtset%ngfft
1531 1238 : nkpt_bz=fockbz%nkpt_bz
1532 1238 : nsppol=dtset%nsppol
1533 1238 : nspinor=1
1534 1238 : ncpgr=0
1535 :
1536 : ! Local variables : useful arrays
1537 3714 : ABI_MALLOC(cgocc,(2,mpw))
1538 1466249 : cgocc=zero
1539 3714 : ABI_MALLOC(cgocc_tmp,(2*mpw+1))
1540 979150 : cgocc_tmp=zero
1541 1238 : if (fockcommon%usepaw==1) then
1542 204 : mgfft=dtset%mgfftdg
1543 3876 : ngfft=dtset%ngfftdg
1544 1224 : ABI_MALLOC(cprj_tmp,(dtset%natom,nspinor))
1545 612 : ABI_MALLOC(dimlmn,(dtset%natom))
1546 204 : call pawcprj_getdim(dimlmn,dtset%natom,nattyp,dtset%ntypat,dtset%typat,fockcommon%pawtab,"O")
1547 : ncpgr = 0
1548 204 : if (dtset%optforces== 1) ncpgr = 3
1549 : ! if (dtset%optstress /= 0) ncpgr = 6
1550 204 : call pawcprj_alloc(cprj_tmp,ncpgr,dimlmn)
1551 :
1552 504 : lmnmax=maxval(fockcommon%pawtab(:)%lmn_size)
1553 816 : ABI_MALLOC(indlmn,(6,lmnmax,dtset%ntypat))
1554 504 : do itypat=1,dtset%ntypat
1555 300 : isize=size(fockcommon%pawtab(itypat)%indlmn,2)
1556 13272 : indlmn(:,1:isize,itypat)=fockcommon%pawtab(itypat)%indlmn(:,1:isize)
1557 : end do
1558 816 : ABI_MALLOC(indsym_,(4,dtset%nsym,dtset%natom))
1559 408 : ABI_MALLOC(typat_srt,(dtset%natom))
1560 :
1561 204 : if (dtset%nsym==1) then
1562 0 : indsym_=0
1563 0 : do iatom=1,dtset%natom
1564 0 : iatm=fockcommon%atindx(iatom)
1565 0 : typat_srt(iatm)=dtset%typat(iatom)
1566 0 : indsym_(4,:,iatom)=iatom
1567 : end do
1568 : else
1569 612 : do iatom=1,dtset%natom
1570 408 : iatm=fockcommon%atindx(iatom)
1571 408 : typat_srt(iatm)=dtset%typat(iatom)
1572 35000 : indsym_(1:3,:,iatm)=indsym(1:3,:,iatom)
1573 9260 : indsym_(4,:,iatm)=fockcommon%atindx(indsym(4,:,iatom))
1574 : end do
1575 : end if
1576 : end if
1577 :
1578 : ! Local variables to perform FFT
1579 1238 : n4=ngfft(4) ; n5=ngfft(5) ; n6=ngfft(6)
1580 6190 : ABI_MALLOC(dummytab3,(n4,n5,n6))
1581 :
1582 :
1583 4772 : if(ANY(fockbz%calc_phase(:)/=0)) then
1584 518 : ABI_MALLOC(phase_jkpt,(2,mpw))
1585 205240 : phase_jkpt=zero
1586 : end if
1587 :
1588 :
1589 : ! =======================================================
1590 : ! === Update the data relative to the occupied states ===
1591 : ! =======================================================
1592 : !* The arrays cgocc_bz, kg_bz, occ_bz and npwarr_bz are already allocated with the maximal size.
1593 : ! if ((dtset%kptopt>=1).and.(dtset%kptopt<=4)) then
1594 : ! if (dtset%kptopt/=3) then
1595 :
1596 2611 : do isppol=1,nsppol
1597 1373 : jbg=0 ; jcg=0 ; jkg=0 ; icp=0
1598 1373 : my_jsppol=isppol
1599 1373 : if ((isppol==2).and.(mpi_enreg%nproc_spkpt/=1)) my_jsppol=1
1600 : !* Both spins are treated on the same proc., only in the case where nproc_spkpt=1;
1601 : !* otherwise each proc. treats only one spin.
1602 :
1603 : ! MG: This loop is not effient!
1604 : ! Ok the number of k-points in the BZ is usually `small` when hybrids are used
1605 : ! but what happens if we use a 12x12x12.
1606 : ! One should loop over the IBZ, broadcast and reconstruct the star of the k-point.
1607 1373 : my_jkpt=0
1608 9068 : do jkpt=1,nkpt_bz
1609 :
1610 6457 : if (proc_distrb_cycle(mpi_enreg%distrb_hf,jkpt,1,mband,1,mpi_enreg%me_hf)) cycle
1611 :
1612 : !* In this case, the processor does not calculate the exchange with any occupied state on jkpt.
1613 :
1614 : ! if (.NOT.(proc_distrb_cycle(mpi_enreg%proc_distrb,jkpt,1,dtset%nbandhf,jsppol,mpi_enreg%me_kpt))) then
1615 : !* The state (jkpt,jband,jsppol) is stored in the array cg of this processor and copied in cgocc_tmp.
1616 : ! icg=icg+dtset%nband(jkpt)*npwj
1617 : ! end if
1618 : ! ibg=ibg+dtset%nband(jkpt)
1619 : ! cycle
1620 : ! end if
1621 6457 : my_jkpt=my_jkpt+1
1622 :
1623 6457 : ikpt=fockbz%tab_ikpt(my_jkpt)
1624 :
1625 : !* ikpt = the point of IBZ that jkpt is an image of in BZ
1626 6457 : npwj=npwarr(ikpt)
1627 : !* npwj= number of plane wave in basis for the wavefunction
1628 6457 : jstwfk=fockbz%istwfk_bz(my_jkpt)
1629 : !* jstwfk= how is stored the wavefunction
1630 6457 : ibg=fockbz%tab_ibg(my_jkpt,my_jsppol)
1631 : !* ibg = shift to be applied on the location of data in the array occ
1632 6457 : icg=fockbz%tab_icg(my_jkpt,my_jsppol)
1633 : !* icg = shift to be applied on the location of data in the array cg
1634 6457 : icp=fockbz%tab_icp(my_jkpt,my_jsppol)
1635 : !* icp = shift to be applied on the location of data in the array cprj
1636 6457 : gbound_k => fockbz%gbound_bz(:,:,my_jkpt)
1637 : !* boundary of the basis sphere of G vectors
1638 6457 : kg_k => fockbz%kg_bz(:,1+jkg:npwj+jkg)
1639 : !* reduced plean wave coordinates
1640 6457 : if (fockbz%calc_phase(my_jkpt)==1) then
1641 737484 : phase_jkpt(:,1:npwj)=fockbz%phase(:,1+jkg:npwj+jkg)
1642 : end if
1643 : !* phase factor at k-point j
1644 :
1645 : !* Initialize the band counter
1646 6457 : my_jband=0
1647 : !FBru: Here run over all the nbandhf bands instead of just the truly occupied states
1648 : ! this is a very tiny waste, but solves parallelization issues
1649 33837 : do iband=1,dtset%nbandhf
1650 :
1651 :
1652 12867116 : cgocc_tmp=zero
1653 27380 : if (fockcommon%usepaw==1) then
1654 2470 : call pawcprj_set_zero(cprj_tmp)
1655 : end if
1656 :
1657 : !* To avoid segmentation fault, my_jband should not be greater than nbandhf
1658 : !FBru: this error should never happen in practice since we now limit the loop to nbandhf
1659 27380 : if ((my_jband+1)>mband) then
1660 0 : write(message,*) 'The number of occupied band',my_jband+1,' at k-point',&
1661 0 : & ikpt,' is greater than the value of nbandhf ', mband
1662 0 : ABI_ERROR(message)
1663 : end if
1664 :
1665 : !* If the processor does not calculate the exchange with the occupied state (jkpt,my_jband), cycle
1666 : ! if (mpi_enreg%distrb_hf(jkpt,(my_jband+1),1)/=mpi_enreg%me_hf) cycle
1667 27380 : if (mpi_enreg%distrb_hf(jkpt,iband,1)/=mpi_enreg%me_hf) cycle
1668 : ! if (mpi_enreg%proc_distrb(jkpt,jband,jsppol)==mpi_enreg%me_kpt) then
1669 : !* The state (jkpt,jband,jsppol) is stored in the array cg of this processor ; shift are incremented.
1670 : ! icg=icg+npwj
1671 : ! end if
1672 : ! ibg=ibg+1
1673 : !* Skip the end of the loop
1674 : ! cycle
1675 : ! end if
1676 :
1677 : !* increment the number of occupied bands treated on this processor
1678 27380 : my_jband = my_jband+1
1679 :
1680 : !* In this case, the processor calculates the exchange with the occupied state (jkpt,my_jband).
1681 27380 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol)==mpi_enreg%me_kpt) then
1682 : !* The state (ikpt,iband,isppol) is stored in the array cg of this processor and copied in cgocc_tmp.
1683 23444 : if(icg==-1) then
1684 0 : write(100,*) 'icg=-1',mpi_enreg%me,isppol,my_jsppol,jkpt,my_jkpt,ikpt,iband
1685 : end if
1686 : ! MG: Why packing re and im part?
1687 23444 : cgocc_tmp(1)=occ(iband+ibg)
1688 4988621 : cgocc_tmp(2:npwj+1)=cg(1,1+(iband-1)*npwj+icg:iband*npwj+icg)
1689 4988621 : cgocc_tmp(npwj+2:2*npwj+1)=cg(2,1+(iband-1)*npwj+icg:iband*npwj+icg)
1690 23444 : if (fockcommon%usepaw==1) then
1691 2470 : call pawcprj_copy(cprj(:,icp+iband:icp+iband+nspinor-1),cprj_tmp)
1692 : end if
1693 : end if
1694 :
1695 : !* Broadcast the state (ikpt,iband,isppol) to all the processors of comm_kpt for cgocc
1696 27380 : call timab(1503,1,tsec2)
1697 27380 : call xmpi_bcast(cgocc_tmp,mpi_enreg%proc_distrb(ikpt,iband,isppol),mpi_enreg%comm_kpt,ier)
1698 :
1699 : !* Broadcast the state (ikpt,iband,isppol) to all the processors of comm_kpt for cprj
1700 27380 : if (fockcommon%usepaw==1) then
1701 : call pawcprj_bcast(cprj_tmp,dtset%natom,nspinor,dimlmn,ncpgr,mpi_enreg%proc_distrb(ikpt,iband,isppol),&
1702 2470 : & mpi_enreg%comm_kpt,ier)
1703 : end if
1704 27380 : call timab(1503,2,tsec2)
1705 : !* Keep the processors in %comm_kpt which needs the values in cgocc_tmp to build their own %cwaveocc and %occ_bz.
1706 27380 : if ((mpi_enreg%nproc_spkpt/=1).and.(nsppol==2)) then
1707 0 : if (fockbz%timerev(my_jkpt)==mpi_enreg%my_isppoltab(isppol)) cycle
1708 : !* In the case of a parallel spin-polarized calculation
1709 : !* when time reversal symmetry is applied at this k-point (timrev==1), only the processors with the opposite spin (my_isppoltab==0) are kept.
1710 : !* 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.
1711 :
1712 : ! if (fock%timerev(my_jkpt)==1)) then
1713 : ! if (mpi_enreg%my_isppoltab(isppol)==1) cycle
1714 : !* In the case of a parallel spin-polarized calculation and when time reversal symmetry is applied at this k-point,
1715 : !* only the processors with the opposite spin are kept.
1716 : ! else
1717 : ! if (mpi_enreg%my_isppoltab(isppol)==0) cycle
1718 : !* only the processors with isppol are kept.
1719 : ! end if
1720 : end if
1721 :
1722 : !* Copy the values of cgocc_tmp in the arrays cgocc and %occ_bz
1723 27380 : fockbz%occ_bz(my_jband+jbg,my_jsppol) = cgocc_tmp(1)
1724 6214285 : cgocc(1,1:npwj) = cgocc_tmp(2:npwj+1)
1725 6214285 : cgocc(2,1:npwj) = cgocc_tmp(npwj+2:2*npwj+1)
1726 :
1727 : !* calculate cg and store it in cgocc_bz
1728 27380 : if (fockbz%calc_phase(my_jkpt)==1) then
1729 1229388 : do jpw=1,npwj
1730 1222956 : cgre=cgocc(1,jpw) ; cgim=cgocc(2,jpw)
1731 1222956 : cgocc(1,jpw) = phase_jkpt(1,jpw)*cgre - phase_jkpt(2,jpw)*cgim
1732 1229388 : cgocc(2,jpw) = phase_jkpt(1,jpw)*cgim + phase_jkpt(2,jpw)*cgre
1733 : end do
1734 : end if ! phase
1735 :
1736 : !* apply time reversal symmetry if necessary
1737 27380 : if (fockbz%timerev(my_jkpt)==1) then
1738 0 : cgocc(2,:) = - cgocc(2,:)
1739 0 : if((mpi_enreg%nproc_spkpt==1).and.(nsppol==2)) my_jsppol=mod(my_jsppol,2)+1
1740 : !* exchange spin (1 ->2 ; 2-> 1) in the sequential case.
1741 : end if
1742 :
1743 : !* apply FFT to get cwaveocc in real space
1744 :
1745 27380 : if (allocated(fockbz%cwaveocc_bz)) then
1746 :
1747 0 : ABI_MALLOC(dummytab2,(2,npwj))
1748 : call fourwf(1,dummytab3,cgocc(:,1:npwj),dummytab2,fockbz%cwaveocc_bz(:,:,:,:,my_jband+jbg,my_jsppol), &
1749 : & gbound_k,gbound_k,jstwfk,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,&
1750 0 : & npwj,npwj,n4,n5,n6,tim_fourwf0,0,weight1,weight1,gpu_option=dtset%gpu_option)
1751 0 : ABI_FREE(dummytab2)
1752 :
1753 : else
1754 18588095 : fockbz%cgocc(:,jcg+1+(my_jband-1)*npwj:jcg+my_jband*npwj,my_jsppol)=cgocc(:,1:npwj)
1755 : end if
1756 :
1757 : !* calculate cprj and store it in cwaveocc_prj
1758 61217 : if (fockcommon%usepaw==1) then
1759 2470 : iband_cprj=(my_jsppol-1)*fockbz%mkptband+jbg+my_jband
1760 2470 : nband=1;mband0=1;iband0=1
1761 : call pawcprj_symkn(fockbz%cwaveocc_prj(:,iband_cprj:iband_cprj+nspinor-1),cprj_tmp(:,1:nspinor),&
1762 : & indsym_,dimlmn,iband0,indlmn,&
1763 : & fockbz%tab_symkpt(my_jkpt),fockbz%timerev(my_jkpt),dtset%kptns(:,ikpt),fockbz%pawang%l_max-1,lmnmax,&
1764 2470 : & mband0,dtset%natom,nband,nspinor,dtset%nsym,dtset%ntypat,typat_srt,fockbz%pawang%zarot)
1765 :
1766 2470 : if(dtset%optforces==1) then
1767 2310 : do iatom=1,dtset%natom
1768 1540 : iatm=fockcommon%atindx(iatom)
1769 3850 : do ispinor=iband_cprj,iband_cprj+nspinor-1
1770 12220 : do ilmn=1,fockcommon%pawtab(dtset%typat(iatom))%lmn_size
1771 : dcp(:)= MATMUL(TRANSPOSE(fockcommon%symrec(:,:,fockbz%tab_symkpt(my_jkpt))),&
1772 228500 : & fockbz%cwaveocc_prj(iatm,ispinor)%dcp(1,:,ilmn))
1773 36560 : fockbz%cwaveocc_prj(iatm,ispinor)%dcp(1,:,ilmn)=dcp(:)
1774 : dcp(:)= MATMUL(TRANSPOSE(fockcommon%symrec(:,:,fockbz%tab_symkpt(my_jkpt))),&
1775 228500 : & fockbz%cwaveocc_prj(iatm,ispinor)%dcp(2,:,ilmn))
1776 38100 : fockbz%cwaveocc_prj(iatm,ispinor)%dcp(2,:,ilmn)=dcp(:)
1777 : end do
1778 : end do
1779 : end do
1780 : end if
1781 :
1782 : end if
1783 :
1784 : !* update the shift to apply to occ in all case because this array is not distributed among the proc.
1785 : ! ibg=ibg+1
1786 :
1787 : end do ! iband
1788 :
1789 : !* Save the true number of occupied bands in the array %nbandocc_bz
1790 6457 : fockbz%nbandocc_bz(my_jkpt,my_jsppol) = my_jband
1791 :
1792 : !* update the shifts to apply
1793 6457 : jbg=jbg+my_jband
1794 6457 : jcg=jcg+npwj*my_jband
1795 7830 : jkg=jkg+npwj
1796 : end do ! ikpt
1797 : end do ! isppol
1798 1238 : if (allocated(fockbz%cwaveocc_bz)) then
1799 0 : fockbz%cwaveocc_bz=fockbz%cwaveocc_bz*invucvol
1800 : end if
1801 :
1802 1238 : ABI_FREE(cgocc_tmp)
1803 1238 : ABI_FREE(cgocc)
1804 1238 : if (fockcommon%usepaw==1) then
1805 204 : ABI_FREE(indlmn)
1806 204 : ABI_FREE(indsym_)
1807 204 : ABI_FREE(typat_srt)
1808 204 : ABI_FREE(dimlmn)
1809 204 : call pawcprj_free(cprj_tmp)
1810 612 : ABI_FREE(cprj_tmp)
1811 : end if
1812 1238 : if(allocated(phase_jkpt)) then
1813 259 : ABI_FREE(phase_jkpt)
1814 : end if
1815 1238 : ABI_FREE(dummytab3)
1816 :
1817 :
1818 : ! Restricted or unrestricted HF
1819 1238 : if (nsppol==1) then
1820 : !* Update the array %occ_bz => May be limited to the occupied states only
1821 27516 : fockbz%occ_bz(:,:)=half*fockbz%occ_bz(:,:)
1822 :
1823 : ! If nsppol=1, this is a restricted Hartree-Fock calculation.
1824 : ! If nsppol=2, this is an unrestricted Hartree-Fock calculation.
1825 : end if
1826 :
1827 : end if
1828 :
1829 1238 : call timab(1502,2,tsec)
1830 :
1831 2476 : end subroutine fock_updatecwaveocc
1832 : !!***
1833 :
1834 : !----------------------------------------------------------------------
1835 :
1836 : !!****f* m_fock/fock_set_getghc_call
1837 : !! NAME
1838 : !! fock_set_getghc_call
1839 : !!
1840 : !! FUNCTION
1841 : !! Set the value of fock%getghc_call, Returns the old value
1842 : !!
1843 : !! SOURCE
1844 :
1845 0 : integer function fock_set_getghc_call(fock, new) result(old)
1846 :
1847 : !Arguments ------------------------------------
1848 : !scalars
1849 : type(fock_common_type),intent(inout) :: fock
1850 : integer,intent(in) :: new
1851 : ! *************************************************************************
1852 :
1853 0 : old = fock%getghc_call_
1854 0 : fock%getghc_call_ = new
1855 :
1856 0 : end function fock_set_getghc_call
1857 : !!***
1858 :
1859 : !----------------------------------------------------------------------
1860 :
1861 : !!****f* m_fock/fock_get_getghc_call
1862 : !! NAME
1863 : !! fock_get_getghc_call
1864 : !!
1865 : !! FUNCTION
1866 : !! Returns the value of fock%getghc_call_
1867 : !!
1868 : !! SOURCE
1869 :
1870 155370 : pure integer function fock_get_getghc_call(fock)
1871 :
1872 : !Arguments ------------------------------------
1873 : type(fock_common_type),intent(in) :: fock
1874 : ! *************************************************************************
1875 :
1876 155370 : fock_get_getghc_call = fock%getghc_call_
1877 :
1878 155370 : end function fock_get_getghc_call
1879 : !!***
1880 :
1881 : !----------------------------------------------------------------------
1882 :
1883 : !!****f* m_fock/fock_print
1884 : !! NAME
1885 : !! fock_print
1886 : !!
1887 : !! FUNCTION
1888 : !! Print info on the fock_type data type
1889 : !!
1890 : !! INPUTS
1891 : !! fock<crystal_t>=The object
1892 : !! [unit]=Unit number for output
1893 : !! [prtvol]=Verbosity level
1894 : !! [mode_paral]=Either "COLL" or "PERS"
1895 : !! [header]=String to be printed as header for additional info.
1896 : !!
1897 : !! OUTPUT
1898 : !! Only printing
1899 : !!
1900 : !! SOURCE
1901 :
1902 123 : subroutine fock_print(fockcommon,fockbz,header,unit,mode_paral,prtvol)
1903 :
1904 : !Arguments ------------------------------------
1905 : !scalars
1906 : integer,optional,intent(in) :: unit,prtvol
1907 : character(len=4),optional,intent(in) :: mode_paral
1908 : character(len=*),optional,intent(in) :: header
1909 : type(fock_common_type),intent(in) :: fockcommon
1910 : type(fock_BZ_type),intent(in) :: fockbz
1911 :
1912 : !Local variables-------------------------------
1913 : integer :: my_unt,my_prtvol
1914 : character(len=4) :: my_mode
1915 : character(len=500) :: msg
1916 : ! *********************************************************************
1917 :
1918 123 : my_unt=std_out; if (PRESENT(unit)) my_unt=unit
1919 123 : my_prtvol=0 ; if (PRESENT(prtvol)) my_prtvol=prtvol
1920 123 : my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode=mode_paral
1921 :
1922 123 : msg=' ==== Info on fock_type ==== '
1923 123 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
1924 123 : call wrtout(my_unt,msg,my_mode)
1925 :
1926 : ! Important dimensions
1927 123 : call wrtout(my_unt,sjoin(" my_nsppol ...",itoa(fockcommon%my_nsppol)),my_mode)
1928 123 : call wrtout(my_unt,sjoin(" nkpt_bz .....",itoa(fockbz%nkpt_bz)),my_mode)
1929 :
1930 : ! Options
1931 123 : call wrtout(my_unt,sjoin(" nnsclo_hf .......",itoa(fockcommon%nnsclo_hf)),my_mode)
1932 123 : call wrtout(my_unt,sjoin(" ixc .............",itoa(fockcommon%ixc)),my_mode)
1933 123 : call wrtout(my_unt,sjoin(" hybrid mixing....",ftoa(fockcommon%hyb_mixing)),my_mode)
1934 123 : call wrtout(my_unt,sjoin(" hybrid SR mixing ",ftoa(fockcommon%hyb_mixing_sr)),my_mode)
1935 123 : call wrtout(my_unt,sjoin(" hybrid range DFT ",ftoa(fockcommon%hyb_range_dft)),my_mode)
1936 123 : call wrtout(my_unt,sjoin(" hybrid range Fock",ftoa(fockcommon%hyb_range_fock)),my_mode)
1937 :
1938 : ! write(msg,"(a,f12.1,a)")" Memory required for HF u(r) states: ",product(shape(fockbz%cwaveocc_bz)) * dp * b2Mb, " [Mb]"
1939 : ! call wrtout(my_unt,msg,my_mode)
1940 :
1941 : ! Extra info.
1942 : !if (my_prtvol > 0) then
1943 : ! call wrtout(my_unt,"Extra info not available",my_mode)
1944 : !end if
1945 :
1946 123 : end subroutine fock_print
1947 : !!***
1948 :
1949 : !----------------------------------------------------------------------
1950 :
1951 : !!****f* ABINIT/bare_vqg
1952 : !! NAME
1953 : !! bare_vqg
1954 : !!
1955 : !! FUNCTION
1956 : !! Compute bare coulomb term in G-space on the FFT mesh i.e. 4pi/(G+q)**2 for a specified qpoint
1957 : !!
1958 : !! INPUTS
1959 : !! qpoint(3)=reduced coordinates for the phonon wavelength
1960 : !! fockcommon=all the technical details of the Fock operator
1961 : !! gmet(3,3)=metrix tensor in G space in Bohr**-2.
1962 : !! nfft=Total number of FFT grid points.
1963 : !! nkpt_bz=total number of kpoints in the full BZ
1964 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1965 : !! ucvol=unitcell volume
1966 : !!
1967 : !! OUTPUT
1968 : !! vqg(nfft)=4pi/(G+q)**2, G=0 component is set to an analytic value if q = Gamma.
1969 : !!
1970 : !! NOTES
1971 : !! This routine operates on the full FFT mesh. DO NOT PASS MPI_TYPE
1972 : !! One can easily implemente MPI-FFT by just calling this routine and then
1973 : !! extracting the G-vectors treated by the node.
1974 : !!
1975 : !! SOURCE
1976 :
1977 215486 : subroutine bare_vqg(qpoint,fockcommon,gmet,nfft,nkpt_bz,ngfft,ucvol,vqg)
1978 :
1979 : !Arguments ------------------------------------
1980 : !scalars
1981 : type(fock_common_type),intent(in) :: fockcommon
1982 : integer,intent(in) :: nfft,nkpt_bz
1983 : real(dp),intent(in) :: ucvol
1984 : !arrays
1985 : integer,intent(in) :: ngfft(18)
1986 : real(dp),intent(in) :: qpoint(3)
1987 : real(dp),intent(inout) :: gmet(3,3)
1988 : real(dp),intent(out) :: vqg(nfft)
1989 :
1990 : !Local variables-------------------------------
1991 : !scalars
1992 : integer :: izero
1993 : real(dp) :: rcut !,gs2,gs3 ! gqgm12,gqgm13,gqgm23,
1994 430972 : real(dp) :: vqg_sr(nfft)
1995 : ! *************************************************************************
1996 :
1997 215486 : if (abs(fockcommon%hyb_mixing_sr)>tol8.and.abs(fockcommon%hyb_range_fock)<tol8) then
1998 0 : ABI_BUG('SR mixing<>0 while range separation=0!')
1999 : end if
2000 :
2001 215486 : izero = fockcommon%usepaw
2002 :
2003 910108658 : vqg = zero
2004 :
2005 215486 : if (abs(fockcommon%hyb_mixing)>tol8) then
2006 132282 : call barevcoul(fockcommon%rcut,fockcommon%fock_icutcoul,qpoint,fockcommon%gsqcut,gmet,nfft,nkpt_bz,ngfft,ucvol,izero,vqg)
2007 : ! Rescale the interaction with the factor hyb_mixing
2008 579933498 : vqg = vqg * fockcommon%hyb_mixing
2009 : end if
2010 :
2011 215486 : if (abs(fockcommon%hyb_mixing_sr)>tol8) then
2012 81956 : rcut= one / fockcommon%hyb_range_fock
2013 328017368 : vqg_sr = zero
2014 81956 : call barevcoul(rcut,fockcommon%fock_icutcoul,qpoint,fockcommon%gsqcut,gmet,nfft,nkpt_bz,ngfft,ucvol,izero,vqg_sr,shortrange=.true.)
2015 : ! Rescale the interaction with the factor hyb_mixing_sr and add it to the full range one
2016 328017368 : vqg = vqg + vqg_sr * fockcommon%hyb_mixing_sr
2017 : end if
2018 :
2019 215486 : end subroutine bare_vqg
2020 : !!***
2021 :
2022 : !!****f* ABINIT/strfock
2023 : !! NAME
2024 : !! strfock
2025 : !!
2026 : !! FUNCTION
2027 : !! Compute Fock energy contribution to stress tensor (Cartesian coordinates).
2028 : !!
2029 : !! INPUTS
2030 : !! fockcommon= basic information for fock calculations
2031 : !! gprimd(3,3)=reciprocal space dimensional primitive translations
2032 : !! mpi_enreg=information about MPI parallelization
2033 : !! nfft=(effective) number of FFT grid points (for this processor)
2034 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
2035 : !! nkpt_bz= number of k points in the BZ
2036 : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
2037 : !! rhog(2,nfft)=Fourier transform of charge density (bohr^-3)
2038 : !! rhog2(2,nfft)= optional argument: Fourier transform of a second charge density (bohr^-3)
2039 : !! ucvol=unit cell volume (bohr^3)
2040 : !!
2041 : !! OUTPUT
2042 : !! fockstr(6)=components of Fock part of stress tensor
2043 : !! (Cartesian coordinates, symmetric tensor) in hartree/bohr^3
2044 : !! Definition of symmetric tensor storage: store 6 unique components
2045 : !! in the order 11, 22, 33, 32, 31, 21 (suggested by Xavier Gonze).
2046 : !!
2047 : !! SOURCE
2048 :
2049 12960 : subroutine strfock(fockcommon,gprimd,fockstr,mpi_enreg,nfft,ngfft,&
2050 12960 : nkpt_bz,ndat,rhog,ucvol,qphon,&
2051 0 : rhog2,gpu_option) ! optional arguments
2052 :
2053 : !Arguments ------------------------------------
2054 : !scalars
2055 : type(fock_common_type),intent(in) :: fockcommon
2056 : integer,intent(in) :: nfft,nkpt_bz,ndat
2057 : integer,intent(in),optional :: gpu_option
2058 : real(dp),intent(in) :: ucvol
2059 : type(MPI_type),intent(in) :: mpi_enreg
2060 : !arrays
2061 : integer,intent(in) :: ngfft(18)
2062 : real(dp),intent(in) :: gprimd(3,3),rhog(2,nfft,ndat),qphon(3)
2063 : real(dp),intent(in),optional :: rhog2(2,nfft,ndat)
2064 : real(dp),intent(out) :: fockstr(6,ndat)
2065 :
2066 : !Local variables-------------------------------
2067 : !scalars
2068 : integer,parameter :: im=2,re=1
2069 : integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,idat,me_fft,n1,n2,n3,nproc_fft
2070 : real(dp) :: arg,gsquar,rcut,rhogsq,tot,tot1 !tolfix=1.000000001_dp,
2071 : logical :: rcut_spencer_alavi
2072 : #ifdef HAVE_OPENMP_OFFLOAD
2073 : real(dp) :: hyb_mixing,hyb_mixing_sr,hyb_range_fock
2074 : ! Cray has trouble with reduction on array, so we use 6 scalars instead
2075 : real(dp) :: fockstr1,fockstr2,fockstr3,fockstr4,fockstr5,fockstr6
2076 : #endif
2077 : !character(len=500) :: msg
2078 : !arrays
2079 25920 : real(dp) :: gcart(3),tsec(2),gmet(3,3),vqg(nfft)
2080 12960 : real(dp), allocatable :: v_gcart(:,:,:,:)
2081 12960 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
2082 12960 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
2083 : ! *************************************************************************
2084 :
2085 12960 : call timab(568,1,tsec)
2086 :
2087 12960 : if (abs(fockcommon%hyb_mixing_sr)>tol8.and.abs(fockcommon%hyb_range_fock)<tol8) then
2088 0 : ABI_BUG('strfock: SR mixing<>0 while range separation=0!')
2089 : end if
2090 :
2091 518400 : gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
2092 12960 : call bare_vqg(qphon,fockcommon,gmet,nfft,nkpt_bz,ngfft,ucvol,vqg)
2093 :
2094 : !if(gpu_option==ABI_GPU_DISABLED) then
2095 441080 : fockstr(:,:)=zero
2096 : !else(gpu_option==ABI_GPU_OPENMP) then
2097 : ! gpu_set_to_zero(fockstr, 6*ndat)
2098 : !end if
2099 :
2100 : ! fockcommon%rcut is zero, rcut is a function of the cell volume (Spencer-Alavi scheme)
2101 : ! Therefore gives a contribution to the stress
2102 12960 : rcut_spencer_alavi = fockcommon%rcut<tol8
2103 12960 : if(rcut_spencer_alavi) then
2104 12864 : rcut = (three*nkpt_bz*ucvol/four_pi)**(one/three)
2105 : else
2106 : rcut = fockcommon%rcut
2107 : endif
2108 12960 : irho2=0;if (present(rhog2)) irho2=1
2109 :
2110 : !Conduct looping over all fft grid points to find G vecs inside fockcommon%gsqcut
2111 :
2112 12960 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
2113 12960 : me_fft=ngfft(11)
2114 12960 : nproc_fft=ngfft(10)
2115 12960 : id1=n1/2+2
2116 12960 : id2=n2/2+2
2117 12960 : id3=n3/2+2
2118 :
2119 :
2120 12960 : ii=0
2121 : ! Get the distrib associated with this fft_grid
2122 12960 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
2123 :
2124 12960 : if(gpu_option==ABI_GPU_DISABLED) then
2125 : !$OMP PARALLEL DO &
2126 : !$OMP& PRIVATE(idat,i3,i2,i1,ig3,ig2,ig1,tot,tot1,ii,gcart,gsquar,rhogsq,arg)
2127 74120 : do idat=1,ndat
2128 1070035 : do i3=1,n3
2129 995915 : ig3=i3-(i3/id3)*n3-1
2130 17745894 : do i2=1,n2
2131 16688819 : ig2=i2-(i2/id2)*n2-1
2132 17684734 : if (fftn2_distrib(i2)==me_fft) then
2133 317020057 : do i1=1,n1
2134 300331238 : tot=zero; tot1=zero
2135 300331238 : ig1=i1-(i1/id1)*n1-1
2136 300331238 : ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
2137 :
2138 : ! Compute cartesian components of G
2139 300331238 : gcart(1)=gprimd(1,1)*(dble(ig1)+qphon(1))+gprimd(1,2)*(dble(ig2)+qphon(2))+gprimd(1,3)*(dble(ig3)+qphon(3))
2140 300331238 : gcart(2)=gprimd(2,1)*(dble(ig1)+qphon(1))+gprimd(2,2)*(dble(ig2)+qphon(2))+gprimd(2,3)*(dble(ig3)+qphon(3))
2141 300331238 : gcart(3)=gprimd(3,1)*(dble(ig1)+qphon(1))+gprimd(3,2)*(dble(ig2)+qphon(2))+gprimd(3,3)*(dble(ig3)+qphon(3))
2142 : ! Compute |G+q|^2
2143 300331238 : gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2
2144 :
2145 :
2146 :
2147 : ! take |rho(G)|^2 for complex rhog
2148 300331238 : if (irho2==0) then
2149 300331238 : rhogsq=rhog(re,ii,idat)**2+rhog(im,ii,idat)**2
2150 : else
2151 0 : rhogsq=rhog(re,ii,idat)*rhog2(re,ii,idat)+rhog(im,ii,idat)*rhog2(im,ii,idat)
2152 : end if
2153 : ! Case G=0:
2154 300331238 : if(gsquar<tol10) then
2155 3790 : if (abs(fockcommon%hyb_mixing_sr)>tol8) cycle
2156 2835 : if (abs(fockcommon%hyb_mixing)>tol8) then
2157 2799 : if (rcut_spencer_alavi) then
2158 : ! vqg(1) already contains the factor fockcommon%hyb_mixing
2159 2751 : fockstr(1,idat)=fockstr(1,idat)+vqg(1)/3.0_dp*rhogsq
2160 2751 : fockstr(2,idat)=fockstr(2,idat)+vqg(1)/3.0_dp*rhogsq
2161 2751 : fockstr(3,idat)=fockstr(3,idat)+vqg(1)/3.0_dp*rhogsq
2162 : endif
2163 : cycle
2164 : end if
2165 : end if
2166 :
2167 : ! Spherical cutoff screening
2168 300327484 : if (abs(fockcommon%hyb_mixing)>tol8) then
2169 222776987 : arg=two_pi*rcut*sqrt(gsquar)
2170 222776987 : tot=fockcommon%hyb_mixing*rhogsq*piinv/(gsquar**2)*(1-cos(arg)-arg*sin(arg)/two)
2171 222776987 : if (rcut_spencer_alavi) then
2172 221204171 : tot1=fockcommon%hyb_mixing*rhogsq/three*rcut*sin(arg)/sqrt(gsquar)
2173 : end if
2174 : end if
2175 :
2176 : ! Erfc screening
2177 300327484 : if (abs(fockcommon%hyb_mixing_sr)>tol8) then
2178 76555169 : arg=-gsquar*pi**2/(fockcommon%hyb_range_fock**2)
2179 76555169 : tot=tot+fockcommon%hyb_mixing_sr*rhogsq*piinv/(gsquar**2)*(1.d0-exp(arg)*(1-arg))
2180 : end if
2181 300327484 : fockstr(1,idat)=fockstr(1,idat)+tot*gcart(1)*gcart(1)+tot1
2182 300327484 : fockstr(2,idat)=fockstr(2,idat)+tot*gcart(2)*gcart(2)+tot1
2183 300327484 : fockstr(3,idat)=fockstr(3,idat)+tot*gcart(3)*gcart(3)+tot1
2184 300327484 : fockstr(4,idat)=fockstr(4,idat)+tot*gcart(3)*gcart(2)
2185 300327484 : fockstr(5,idat)=fockstr(5,idat)+tot*gcart(3)*gcart(1)
2186 317017258 : fockstr(6,idat)=fockstr(6,idat)+tot*gcart(2)*gcart(1)
2187 : end do
2188 : end if
2189 : end do
2190 : end do
2191 : end do !ndat
2192 0 : else if(gpu_option==ABI_GPU_OPENMP) then
2193 0 : ABI_MALLOC(v_gcart, (3,n1,n2,n3))
2194 : #ifdef HAVE_OPENMP_OFFLOAD
2195 : hyb_mixing = fockcommon%hyb_mixing
2196 : hyb_mixing_sr = fockcommon%hyb_mixing_sr
2197 : hyb_range_fock = fockcommon%hyb_range_fock
2198 :
2199 : !$OMP TARGET ENTER DATA MAP(alloc:v_gcart)
2200 : !$OMP TARGET ENTER DATA MAP(to:gprimd,qphon,vqg,ffti2_local,fftn2_distrib)
2201 :
2202 : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO &
2203 : !$OMP& PRIVATE(idat,i3,i2,i1,ig3,ig2,ig1,tot,tot1,ii,gcart,gsquar,rhogsq,arg) &
2204 : !$OMP& MAP(to:v_gcart,gprimd,qphon)
2205 : do i3=1,n3
2206 : do i2=1,n2
2207 : do i1=1,n1
2208 : ig3=i3-(i3/id3)*n3-1
2209 : ig2=i2-(i2/id2)*n2-1
2210 : ig1=i1-(i1/id1)*n1-1
2211 :
2212 : ! Compute cartesian components of G
2213 : v_gcart(1,i1,i2,i3)=gprimd(1,1)*(dble(ig1)+qphon(1))+gprimd(1,2)*(dble(ig2)+qphon(2))+gprimd(1,3)*(dble(ig3)+qphon(3))
2214 : v_gcart(2,i1,i2,i3)=gprimd(2,1)*(dble(ig1)+qphon(1))+gprimd(2,2)*(dble(ig2)+qphon(2))+gprimd(2,3)*(dble(ig3)+qphon(3))
2215 : v_gcart(3,i1,i2,i3)=gprimd(3,1)*(dble(ig1)+qphon(1))+gprimd(3,2)*(dble(ig2)+qphon(2))+gprimd(3,3)*(dble(ig3)+qphon(3))
2216 : end do
2217 : end do
2218 : end do
2219 :
2220 : !$OMP TARGET TEAMS DISTRIBUTE &
2221 : !$OMP& PRIVATE(idat,fockstr1,fockstr2,fockstr3,fockstr4,fockstr5,fockstr6) &
2222 : !$OMP& MAP(to:rhog,v_gcart,vqg,ffti2_local,fftn2_distrib) MAP(tofrom:fockstr)
2223 : do idat=1,ndat
2224 : fockstr1=zero
2225 : fockstr2=zero
2226 : fockstr3=zero
2227 : fockstr4=zero
2228 : fockstr5=zero
2229 : fockstr6=zero
2230 : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i3,i2,i1,tot,tot1,ii,gsquar,rhogsq,arg) &
2231 : !$OMP& REDUCTION(+:fockstr1,fockstr2,fockstr3,fockstr4,fockstr5,fockstr6)
2232 : do i3=1,n3
2233 : do i2=1,n2
2234 : do i1=1,n1
2235 : tot=zero; tot1=zero
2236 : ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1))
2237 : if (fftn2_distrib(i2)==me_fft) then
2238 : ! Compute |G+q|^2
2239 : gsquar=v_gcart(1,i1,i2,i3)**2+v_gcart(2,i1,i2,i3)**2+v_gcart(3,i1,i2,i3)**2
2240 : ! take |rho(G)|^2 for complex rhog
2241 : rhogsq=rhog(re,ii,idat)**2+rhog(im,ii,idat)**2
2242 : ! Case G=0:
2243 : if(gsquar<tol10) then
2244 : if (abs(hyb_mixing)>tol8 .and. abs(hyb_mixing_sr)<tol8) then
2245 : ! vqg(1) already contains the factor fockcommon%hyb_mixing
2246 : if (rcut_spencer_alavi) then
2247 : fockstr1=fockstr1+vqg(1)/3.0_dp*rhogsq
2248 : fockstr2=fockstr2+vqg(1)/3.0_dp*rhogsq
2249 : fockstr3=fockstr3+vqg(1)/3.0_dp*rhogsq
2250 : end if
2251 : end if
2252 :
2253 : else
2254 :
2255 : ! Spherical cutoff screening
2256 : if (abs(hyb_mixing)>tol8) then
2257 : arg=two_pi*rcut*sqrt(gsquar)
2258 : tot=hyb_mixing*rhogsq*piinv/(gsquar**2)*(1-cos(arg)-arg*sin(arg)/two)
2259 : if (rcut_spencer_alavi) then
2260 : tot1=hyb_mixing*rhogsq/three*rcut*sin(arg)/sqrt(gsquar)
2261 : end if
2262 : end if
2263 :
2264 : ! Erfc screening
2265 : if (abs(hyb_mixing_sr)>tol8) then
2266 : arg=-gsquar*pi**2/(hyb_range_fock**2)
2267 : tot=tot+hyb_mixing_sr*rhogsq*piinv/(gsquar**2)*(1.d0-exp(arg)*(1-arg))
2268 : end if
2269 : fockstr1=fockstr1+tot*v_gcart(1,i1,i2,i3)*v_gcart(1,i1,i2,i3)+tot1
2270 : fockstr2=fockstr2+tot*v_gcart(2,i1,i2,i3)*v_gcart(2,i1,i2,i3)+tot1
2271 : fockstr3=fockstr3+tot*v_gcart(3,i1,i2,i3)*v_gcart(3,i1,i2,i3)+tot1
2272 : fockstr4=fockstr4+tot*v_gcart(3,i1,i2,i3)*v_gcart(2,i1,i2,i3)
2273 : fockstr5=fockstr5+tot*v_gcart(3,i1,i2,i3)*v_gcart(1,i1,i2,i3)
2274 : fockstr6=fockstr6+tot*v_gcart(2,i1,i2,i3)*v_gcart(1,i1,i2,i3)
2275 : end if
2276 : end if
2277 : end do
2278 : end do
2279 : end do
2280 : fockstr(1,idat)=fockstr1
2281 : fockstr(2,idat)=fockstr2
2282 : fockstr(3,idat)=fockstr3
2283 : fockstr(4,idat)=fockstr4
2284 : fockstr(5,idat)=fockstr5
2285 : fockstr(6,idat)=fockstr6
2286 : end do !ndat
2287 :
2288 : !$OMP TARGET EXIT DATA MAP(delete:v_gcart,gprimd,qphon,vqg,ffti2_local,fftn2_distrib)
2289 : #endif
2290 0 : ABI_FREE(v_gcart)
2291 : end if
2292 :
2293 : !Init mpi_comm
2294 12960 : if(mpi_enreg%nproc_fft>1)then
2295 0 : call timab(48,1,tsec)
2296 0 : call xmpi_sum(fockstr,mpi_enreg%comm_fft ,ierr)
2297 0 : call timab(48,2,tsec)
2298 : end if
2299 :
2300 :
2301 : !Normalize and add term -efock/ucvol on diagonal
2302 : !efock has been set to zero because it is not yet known. It will be added later.
2303 74120 : fockstr(1,:)=-fockstr(1,:)
2304 74120 : fockstr(2,:)=-fockstr(2,:)
2305 74120 : fockstr(3,:)=-fockstr(3,:)
2306 74120 : fockstr(4,:)=-fockstr(4,:)
2307 74120 : fockstr(5,:)=-fockstr(5,:)
2308 74120 : fockstr(6,:)=-fockstr(6,:)
2309 :
2310 12960 : call timab(568,2,tsec)
2311 :
2312 12960 : end subroutine strfock
2313 : !!***
2314 :
2315 18280 : end module m_fock
2316 : !!***
|
