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