Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_calc_hybrid
8 : USE m_judft
9 : use m_store_load_hybrid
10 : private
11 : public calc_hybrid
12 : CONTAINS
13 :
14 186 : SUBROUTINE calc_hybrid(fi,mpdata,hybdat,fmpi,nococonv,stars,enpara,&
15 : results,xcpot,v,iter, iterHF)
16 : use m_work_package
17 : use m_set_coul_participation
18 : USE m_types_hybdat
19 : USE m_types
20 : USE m_mixedbasis
21 : USE m_coulombmatrix
22 : USE m_hf_init
23 : USE m_hf_setup
24 : USE m_hsfock
25 : USE m_io_hybrid
26 : USE m_eig66_io
27 : use m_eig66_mpi
28 : use m_distribute_mpi
29 : use m_create_coul_comms
30 : use m_eigvec_setup
31 : use m_distrib_vx
32 : use, intrinsic :: iso_c_binding
33 : #ifdef CPP_MPI
34 : use mpi
35 : #endif
36 : #ifdef CPP_PROG_THREAD
37 : use m_thread_lib
38 : #endif
39 :
40 : IMPLICIT NONE
41 :
42 : type(t_fleurinput), intent(in) :: fi
43 : type(t_mpdata), intent(inout) :: mpdata
44 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
45 : TYPE(t_mpi), INTENT(IN) :: fmpi
46 : TYPE(t_nococonv), INTENT(IN) :: nococonv
47 : type(t_stars), intent(in) :: stars
48 : TYPE(t_enpara), INTENT(IN) :: enpara
49 : TYPE(t_results), INTENT(INOUT) :: results
50 : TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
51 : TYPE(t_potden), INTENT(IN) :: v
52 : INTEGER, INTENT(INOUT) :: iter, iterHF
53 :
54 : ! local variables
55 : type(t_hybmpi) :: glob_mpi, wp_mpi
56 1332 : type(t_work_package) :: work_pack(fi%input%jspins)
57 : INTEGER :: jsp, nk, err, i, wp_rank, ierr, ik
58 : INTEGER :: j_wp, n_wps, root_comm
59 : INTEGER :: max_band_pack, jq
60 186 : type(t_lapw) :: lapw
61 : LOGICAL :: init_vex = .TRUE. !In first call we have to init v_nonlocal
62 : character(len=999) :: msg
63 186 : REAL, ALLOCATABLE :: eig_irr(:, :)
64 186 : integer, allocatable :: vx_loc(:,:), weights(:)
65 : type(c_ptr) :: threadId
66 186 : type(t_mat), allocatable :: vx_tmp(:,:)
67 :
68 186 : CALL timestart("hybrid code")
69 :
70 : #ifdef CPP_MPI
71 : #ifdef CPP_PROG_THREAD
72 186 : if(fmpi%l_mpi_multithreaded) call start_prog_thread(threadId)
73 : #endif
74 : #endif
75 :
76 186 : IF (fi%kpts%nkptf == 0) THEN
77 : CALL judft_error("kpoint-set of full BZ not available", &
78 0 : hint="to generate fi%kpts in the full BZ you should specify a k-mesh in inp.xml")
79 : END IF
80 :
81 : !Check if new non-local potential shall be generated
82 186 : hybdat%l_subvxc = fi%hybinp%l_hybrid .AND. (.NOT. xcpot%is_name("exx"))
83 : !If this is the first iteration loop we can not calculate a new non-local potential
84 : !hybdat%l_calhf = (results%last_distance >= 0.0) .AND. (results%last_distance < fi%input%minDistance)
85 : !make sure we do at least one PBE first
86 186 : if(iter == 1 .and. iterHF == 0) hybdat%l_calhf = .False.
87 :
88 186 : IF (.NOT. hybdat%l_calhf) THEN
89 262 : hybdat%l_subvxc = hybdat%l_subvxc .AND. hybdat%l_addhf
90 : else
91 12 : call glob_mpi%init(fmpi%mpi_comm)
92 12 : results%te_hfex%core = 0
93 :
94 : !Check if we are converged well enough to calculate a new potential
95 12 : hybdat%l_addhf = .TRUE.
96 :
97 : !In first iteration allocate some memory
98 12 : IF (init_vex) THEN
99 6 : call first_iteration_alloc(fi, hybdat)
100 6 : init_vex = .FALSE.
101 : END IF
102 12 : hybdat%l_subvxc = (hybdat%l_subvxc .AND. hybdat%l_addhf)
103 12 : IF (.NOT. ALLOCATED(results%w_iks)) allocate(results%w_iks(fi%input%neig, fi%kpts%nkpt, fi%input%jspins))
104 :
105 12 : iterHF = iterHF + 1
106 :
107 : !Delete broyd files
108 12 : CALL system("rm -f broyd*")
109 :
110 : !check if z-reflection trick can be used
111 :
112 :
113 12 : CALL timestart("Preparation for hybrid functionals")
114 : !construct the mixed-basis
115 12 : CALL timestart("generation of mixed basis")
116 12 : if(glob_mpi%rank == 0) write (*,*) "iterHF = " // int2str(iterHF)
117 : CALL mixedbasis(fi%atoms, fi%kpts, fi%input, fi%cell, xcpot, fi%mpinp, mpdata, fi%hybinp, hybdat,&
118 12 : enpara, fmpi, v, iterHF)
119 12 : CALL timestop("generation of mixed basis")
120 :
121 : ! setup parallelization
122 12 : n_wps = min(glob_mpi%size, fi%kpts%nkpt)
123 60 : allocate(weights(n_wps), source=0)
124 36 : do j_wp = 1, n_wps
125 36 : do ik = j_wp, fi%kpts%nkpt, n_wps
126 36 : weights(j_wp) = weights(j_wp) + fi%kpts%eibz(ik)%nkpt
127 : enddo
128 : enddo
129 12 : call distribute_mpi(weights, glob_mpi, wp_mpi, wp_rank)
130 12 : call hybdat%set_nobd(fi, results)
131 12 : call hybdat%set_nbands(fi, fmpi, results)
132 28 : do jsp = 1,fi%input%jspins
133 28 : call work_pack(jsp)%init(fi, hybdat, mpdata, wp_mpi, jsp, wp_rank, n_wps)
134 : enddo
135 :
136 102 : if(.not. allocated(hybdat%zmat)) allocate(hybdat%zmat(fi%kpts%nkptf, fi%input%jspins))
137 :
138 28 : DO jsp = 1, fi%input%jspins
139 156 : DO nk = 1,fi%kpts%nkptf
140 : ! IF (hybdat%zmat(nk, jsp)%mat%matsize2 .NE. hybdat%nbands(nk, jsp)) THEN ! This IF caused deadlocks.
141 128 : CALL lapw%init(fi%input, fi%noco, nococonv,fi%kpts, fi%atoms, fi%sym, nk, fi%cell)
142 : call eigvec_setup(hybdat%zmat(nk, jsp), fi, lapw, work_pack, fmpi, &
143 144 : hybdat%nbands(nk, jsp), nk, jsp, hybdat%eig_id)
144 : ! END IF
145 : enddo
146 : enddo
147 12 : call bcast_eigvecs(hybdat, fi, nococonv, fmpi)
148 :
149 42 : if(.not. allocated(hybdat%coul)) allocate(hybdat%coul(fi%kpts%nkpt))
150 12 : call set_coul_participation(hybdat, fi, fmpi, work_pack)
151 12 : call create_coul_comms(hybdat, fi, fmpi)
152 :
153 48 : do i =1,fi%kpts%nkpt
154 48 : if(hybdat%coul(i)%l_participate) then
155 36 : call hybdat%coul(i)%alloc(fi, mpdata%num_radbasfn, mpdata%n_g, i, .false.)
156 : else
157 0 : call hybdat%coul(i)%mini_alloc(fi)
158 : endif
159 : enddo
160 :
161 : ! use jsp=1 for coulomb work-planning
162 12 : CALL coulombmatrix(fmpi, fi, mpdata, hybdat, xcpot)
163 :
164 48 : do i =1,fi%kpts%nkpt
165 48 : if(hybdat%coul(i)%l_participate) then
166 36 : call hybdat%coul(i)%mpi_bcast(fi, hybdat%coul(i)%comm, 0)
167 : endif
168 : enddo
169 :
170 12 : CALL hf_init(mpdata, fi, hybdat)
171 12 : CALL timestop("Preparation for hybrid functionals")
172 :
173 12 : call judft_comm_split(glob_mpi%comm, wp_mpi%rank, 0, root_comm)
174 :
175 12 : CALL timestart("Calculation of non-local HF potential")
176 112 : allocate(vx_loc(fi%kpts%nkpt,fi%input%jspins), source=-1)
177 112 : allocate(vx_tmp(fi%kpts%nkpt, fi%input%jspins))
178 28 : DO jsp = 1, fi%input%jspins
179 : CALL HF_setup(mpdata,fi, fmpi, nococonv, results, jsp, enpara, &
180 16 : hybdat, v%mt(:, 0, :, :), eig_irr)
181 :
182 16 : call timestart("get max_q")
183 16 : hybdat%max_q = 0
184 40 : DO i = 1,work_pack(jsp)%k_packs(1)%size
185 24 : max_band_pack = 0
186 112 : DO jq = 1, size(work_pack(jsp)%k_packs(i)%q_packs)
187 112 : max_band_pack = max(max_band_pack, work_pack(jsp)%k_packs(i)%q_packs(jq)%submpi%size)
188 : enddo
189 40 : hybdat%max_q = hybdat%max_q + size(work_pack(jsp)%k_packs(i)%q_packs) * fi%kpts%nkptf * max_band_pack
190 : END DO
191 : #ifdef CPP_MPI
192 16 : call MPI_Allreduce(MPI_IN_PLACE, hybdat%max_q, 1, MPI_INTEGER, MPI_MAX, MPI_COMM_WORLD, ierr)
193 : #endif
194 16 : hybdat%max_q = hybdat%max_q + 20 ! The 20 is kind of dirty. It is meant as a safety because of the MPI_BARRIER in hsfock, related to the ex_to_vx call.
195 16 : call timestop("get max_q")
196 :
197 40 : DO i = 1,work_pack(jsp)%k_packs(1)%size
198 24 : nk = work_pack(jsp)%k_packs(i)%nk
199 24 : PRINT*, 'kpoint= ', nk
200 24 : CALL lapw%init(fi%input, fi%noco, nococonv,fi%kpts, fi%atoms, fi%sym, nk, fi%cell)
201 : CALL hsfock(fi, work_pack(jsp)%k_packs(i), mpdata, lapw, jsp, hybdat, eig_irr, &
202 24 : nococonv, stars, results, xcpot, fmpi, vx_tmp(nk, jsp))
203 40 : if(work_pack(jsp)%k_packs(i)%submpi%root()) vx_loc(nk, jsp) = fmpi%irank
204 : END DO
205 :
206 : #ifdef CPP_MPI
207 16 : CALL timestart("balancing MPI_Barriers")
208 1056 : DO WHILE (hybdat%max_q > 0)
209 1040 : call MPI_Barrier(MPI_COMM_WORLD, ierr)
210 1040 : hybdat%max_q = hybdat%max_q - 1
211 : END DO
212 16 : CALL timestop("balancing MPI_Barriers")
213 : #endif
214 :
215 28 : call work_pack(jsp)%free()
216 : END DO
217 : #ifdef CPP_MPI
218 12 : call timestart("MPI_Allred te_hfex%core")
219 12 : if(wp_mpi%rank == 0) call MPI_Allreduce(MPI_IN_PLACE, results%te_hfex%core, 1, MPI_DOUBLE_PRECISION, MPI_SUM, root_comm, ierr)
220 12 : call timestop("MPI_Allred te_hfex%core")
221 : #endif
222 12 : CALL timestop("Calculation of non-local HF potential")
223 : #ifdef CPP_MPI
224 36 : call MPI_Allreduce(MPI_IN_PLACE, vx_loc, size(vx_loc), MPI_INTEGER, MPI_MAX, fmpi%mpi_comm, ierr)
225 : #endif
226 12 : call distrib_vx(fi, fmpi, nococonv, vx_loc, vx_tmp, hybdat)
227 12 : call store_hybrid_data(fi, fmpi, hybdat)
228 :
229 : #ifdef CPP_MPI
230 12 : call timestart("Hybrid imbalance")
231 : call MPI_Barrier(fmpi%mpi_comm, err)
232 12 : call timestop("Hybrid imbalance")
233 : #endif
234 :
235 : ENDIF
236 :
237 : #ifdef CPP_MPI
238 : #ifdef CPP_PROG_THREAD
239 186 : if(fmpi%l_mpi_multithreaded) call stop_prog_thread(threadId)
240 : #endif
241 : #endif
242 420 : CALL timestop("hybrid code")
243 : CONTAINS
244 6 : subroutine first_iteration_alloc(fi, hybdat)
245 : implicit none
246 : type(t_fleurinput), intent(in) :: fi
247 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
248 :
249 6 : if(allocated(hybdat%nbands)) then
250 0 : deallocate(hybdat%nbands, stat=err, errmsg=msg)
251 : if(err /= 0) THEN
252 : write (*,*) "errorcode", err
253 : write (*,*) "errormessage", msg
254 : endif
255 : endif
256 :
257 96 : allocate(hybdat%nbands(fi%kpts%nkptf, fi%input%jspins), source=0)
258 :
259 6 : if(allocated(hybdat%nbasm)) deallocate(hybdat%nbasm)
260 66 : allocate(hybdat%nbasm(fi%kpts%nkptf), source=0)
261 :
262 6 : if(allocated(hybdat%div_vv)) deallocate(hybdat%div_vv)
263 1682 : allocate(hybdat%div_vv(fi%input%neig, fi%kpts%nkpt, fi%input%jspins), source=0.0)
264 6 : end subroutine first_iteration_alloc
265 : END SUBROUTINE calc_hybrid
266 : END MODULE m_calc_hybrid
|
