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