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