Branch data Line data Source code
1 : : !!****m* ABINIT/m_gbt_fullsoc
2 : : !! NAME
3 : : !! m_gbt_fullsoc
4 : : !!
5 : : !! FUNCTION
6 : : !! Driver-level helper for the gstate post-SCF full-SOC
7 : : !! k-string correction.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2026 ABINIT group
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_gbt_fullsoc
24 : :
25 : : use defs_basis
26 : : use m_errors
27 : : use m_xmpi
28 : : use m_abicore
29 : : use defs_datatypes, only : pseudopotential_type
30 : : use defs_abitypes, only : MPI_type
31 : : use m_dtset, only : dataset_type
32 : : use m_pawtab, only : pawtab_type
33 : : use m_pawcprj, only : pawcprj_type
34 : : use m_geometry, only : metric
35 : : use m_wfk, only : wfk_read_my_kptbands
36 : : use m_hide_lapack, only : xheev
37 : : use m_bz_mesh, only : isamek
38 : : use m_hamiltonian, only : gs_hamiltonian_type
39 : : use m_mkffnl, only : mkffnl
40 : : use m_nonlop, only : nonlop
41 : : use m_cgtools, only : dotprod_g
42 : :
43 : : implicit none
44 : :
45 : : private
46 : : !!***
47 : :
48 : : public :: gbt_fullsoc_kstring
49 : : !!***
50 : :
51 : : contains
52 : : !!***
53 : :
54 : : !!****f* m_gbt_fullsoc/gbt_fullsoc_kstring
55 : : !! NAME
56 : : !! gbt_fullsoc_kstring
57 : : !!
58 : : !! FUNCTION
59 : : !! Post-SCF force-theorem correction for use_gbt=3 and use_gbt=4.
60 : : !! use_gbt=3 starts from GBT+PSOC and adds spin-flip SOC. use_gbt=4
61 : : !! starts from a SOC-free GBT reference and adds all SOC matrix elements.
62 : : !!
63 : : !! SOURCE
64 : :
65 : 4 : subroutine gbt_fullsoc_kstring(cg,dtset,eigen,kg,mcg,mpi_enreg,&
66 [ - + ]: 2 : & npwarr,occ,pawtab,psps,etotal_psoc,rprimd,wfk_path,xred)
67 : :
68 : : !Arguments ------------------------------------
69 : : !scalars
70 : : integer,intent(in) :: mcg
71 : : real(dp),intent(in) :: etotal_psoc
72 : : character(len=*),intent(in) :: wfk_path
73 : : type(dataset_type),intent(in) :: dtset
74 : : type(MPI_type),intent(in) :: mpi_enreg
75 : : type(pseudopotential_type),intent(in) :: psps
76 : : !arrays
77 : : integer,intent(in) :: kg(:,:),npwarr(dtset%nkpt)
78 : : real(dp),intent(in) :: cg(2,mcg),eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
79 : : real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
80 : : real(dp),intent(in) :: rprimd(3,3),xred(3,dtset%natom)
81 : : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
82 : :
83 : : !Local variables-------------------------------
84 : : !scalars
85 : : integer :: bd_index,cg_index,dim_block,iband,ikg_index,ikpt,ikpt0,istr,mcg_work,me_cell
86 : : integer :: nneeded_k
87 : : integer :: nclosed,nclosed_local,nclosed_used,nproc_cell,nskipped,nstr,string_id,xmpierr
88 : : logical :: use_wfk_restart
89 : : real(dp) :: delta_fullsoc,delta_local,ecut_eff_restart,etot_fullsoc
90 : : character(len=500) :: msg
91 : 2 : type(gs_hamiltonian_type) :: gs_hamk_full
92 : : !arrays
93 : 2 : integer,allocatable :: bd_start(:),cg_start(:),kg_start(:),next_k(:),string(:)
94 : 2 : integer,allocatable :: kg_work(:,:)
95 : 2 : logical,allocatable :: distrb_flags(:,:,:),needed_k(:),visited(:)
96 : 2 : real(dp),allocatable :: cg_work(:,:),eigen_work(:),occ_work(:)
97 : : ! *************************************************************************
98 : :
99 : : DBG_ENTER("COLL")
100 : :
101 : 2 : delta_fullsoc=zero
102 : 2 : delta_local=zero
103 : :
104 : 2 : me_cell=xmpi_comm_rank(mpi_enreg%comm_cell)
105 : 2 : nproc_cell=xmpi_comm_size(mpi_enreg%comm_cell)
106 : 2 : use_wfk_restart=(dtset%nstep == 0)
107 : :
108 [ - + ]: 2 : ABI_CHECK(dtset%use_gbt == 3 .or. dtset%use_gbt == 4, "Invalid full-SOC GBT mode")
109 : :
110 : : if (psps%usepaw /= 0 .or. psps%useylm /= 0 .or. dtset%paral_kgb /= 0 .or. &
111 [ + + + - : 18 : & dtset%nspinor /= 2 .or. dtset%nsppol /= 1 .or. any(dtset%istwfk(1:dtset%nkpt) /= 1)) then
+ - + - +
- + - + -
- + ]
112 : 0 : msg=' Full-SOC GBT k-string correction cannot be computed: only NC, useylm=0,'
113 [ # # ]: 0 : msg=trim(msg)//' paral_kgb=0, nspinor=2, nsppol=1, istwfk=1 are implemented.'
114 : 0 : ABI_ERROR(msg)
115 : : end if
116 : :
117 [ + - - - : 2 : if (.not. use_wfk_restart .and. (nproc_cell > 1 .or. dtset%mkmem /= dtset%nkpt)) then
- - ]
118 : 0 : msg=' Full-SOC GBT correction during the reference calculation cannot be computed:'
119 [ # # ]: 0 : msg=trim(msg)//' use one MPI rank with mkmem = nkpt, or restart with nstep = 0 from a PSOC WFK.'
120 : 0 : ABI_ERROR(msg)
121 : : end if
122 : :
123 [ - + ]: 2 : if (use_wfk_restart .and. len_trim(wfk_path) == 0) then
124 : 0 : msg=' Full-SOC GBT correction with nstep=0 requires a reference WFK restart path.'
125 : 0 : ABI_ERROR(msg)
126 : : end if
127 : :
128 [ + - - + ]: 2 : if (use_wfk_restart .and. (xmpi_comm_size(mpi_enreg%comm_spinorfft) /= 1 .or. &
129 : : & xmpi_comm_size(mpi_enreg%comm_atom) /= 1)) then
130 : 0 : msg=' Full-SOC GBT nstep=0 correction currently supports k-string parallelism only;'
131 [ # # ]: 0 : msg=trim(msg)//' band/FFT/spinor/atom parallel subcommunicators must be serial.'
132 : 0 : ABI_ERROR(msg)
133 : : end if
134 : :
135 [ - + + - : 6 : ABI_MALLOC(next_k,(dtset%nkpt))
+ - + - -
+ - + ]
136 [ + + ]: 18 : do ikpt=1,dtset%nkpt
137 [ + + ]: 66 : next_k(ikpt)=gbt_find_kpoint(dtset%kptns(:,ikpt)+dtset%qgbt,dtset%kptns,dtset%nkpt)
138 : : end do
139 : :
140 [ + + + - : 18 : if (any(next_k == 0)) then
- + ]
141 : 0 : write(msg,'(a,3es16.8)')' Full-SOC GBT k-string correction cannot be computed: qgbt=',&
142 : 0 : & dtset%qgbt
143 [ # # ]: 0 : msg=trim(msg)//' does not map the present k-point mesh onto closed k strings.'
144 [ # # ]: 0 : msg=trim(msg)//' Use a commensurate k mesh so every k+qgbt point is present.'
145 [ # # ]: 0 : ABI_FREE(next_k)
146 : 0 : ABI_ERROR(msg)
147 : : end if
148 : :
149 [ - + + - : 4 : ABI_MALLOC(kg_start,(dtset%nkpt))
- + - + ]
150 [ - + + - : 4 : ABI_MALLOC(cg_start,(dtset%nkpt))
- + - + ]
151 [ - + + - : 4 : ABI_MALLOC(bd_start,(dtset%nkpt))
- + - + ]
152 [ + + ]: 18 : kg_start(:)=-1
153 [ + + ]: 18 : cg_start(:)=-1
154 : : bd_index=0
155 [ + + ]: 18 : do ikpt=1,dtset%nkpt
156 : 16 : bd_start(ikpt)=bd_index
157 : 18 : bd_index=bd_index+dtset%nband(ikpt)
158 : : end do
159 : :
160 [ - + + - : 4 : ABI_MALLOC(visited,(dtset%nkpt))
- + - + ]
161 [ - + + - : 4 : ABI_MALLOC(string,(dtset%nkpt))
- + - + ]
162 [ - + + - : 4 : ABI_MALLOC(needed_k,(dtset%nkpt))
- + - + ]
163 [ + + ]: 18 : visited(:)=.false.
164 [ + + ]: 18 : needed_k(:)=.false.
165 : 2 : nclosed=0; nskipped=0
166 : :
167 [ + - ]: 2 : if (use_wfk_restart) then
168 [ - + + - : 10 : ABI_MALLOC(distrb_flags,(dtset%nkpt,dtset%mband,dtset%nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
169 [ + + + + : 436 : distrb_flags(:,:,:)=.false.
+ + ]
170 : : end if
171 : :
172 [ + + ]: 18 : do ikpt0=1,dtset%nkpt
173 [ + + ]: 16 : if (visited(ikpt0)) cycle
174 : 8 : nstr=0; ikpt=ikpt0
175 : : do
176 [ - + ]: 16 : if (visited(ikpt)) exit
177 : 16 : nstr=nstr+1
178 : 16 : string(nstr)=ikpt
179 : 16 : visited(ikpt)=.true.
180 : 16 : ikpt=next_k(ikpt)
181 [ + + ]: 16 : if (ikpt == ikpt0) exit
182 [ + - ]: 16 : if (nstr >= dtset%nkpt) exit
183 : : end do
184 : :
185 [ - + ]: 8 : if (ikpt /= ikpt0) then
186 : 0 : nskipped=nskipped+1
187 : 0 : cycle
188 : : end if
189 : :
190 : 8 : dim_block=0
191 [ + + ]: 24 : do istr=1,nstr
192 : 24 : dim_block=dim_block+dtset%nband(string(istr))
193 : : end do
194 [ - + ]: 8 : if (dim_block <= 0) cycle
195 : :
196 : 18 : nclosed=nclosed+1
197 : : end do
198 : :
199 [ - + ]: 2 : if (.not. use_wfk_restart) then
200 : : ikg_index=0; cg_index=0
201 [ # # ]: 0 : do ikpt=1,dtset%nkpt
202 : 0 : kg_start(ikpt)=ikg_index
203 : 0 : cg_start(ikpt)=cg_index
204 : 0 : ikg_index=ikg_index+npwarr(ikpt)
205 : 0 : cg_index=cg_index+npwarr(ikpt)*dtset%nspinor*dtset%nband(ikpt)
206 : : end do
207 : : end if
208 : :
209 : : call gs_hamk_full%init(psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,&
210 : : & dtset%natom,dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
211 : : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
212 : : & mpi_spintab=mpi_enreg%my_isppoltab,gpu_option=dtset%gpu_option,&
213 : 2 : & spinaxis=dtset%spinaxis,use_gbt=0)
214 : :
215 [ + + ]: 18 : visited(:)=.false.
216 : : string_id=0
217 : : nclosed_local=0
218 [ + + ]: 18 : do ikpt0=1,dtset%nkpt
219 [ + + ]: 16 : if (visited(ikpt0)) cycle
220 : 8 : nstr=0; ikpt=ikpt0
221 : : do
222 [ - + ]: 16 : if (visited(ikpt)) exit
223 : 16 : nstr=nstr+1
224 : 16 : string(nstr)=ikpt
225 : 16 : visited(ikpt)=.true.
226 : 16 : ikpt=next_k(ikpt)
227 [ + + ]: 16 : if (ikpt == ikpt0) exit
228 [ + - ]: 16 : if (nstr >= dtset%nkpt) exit
229 : : end do
230 [ - + ]: 8 : if (ikpt /= ikpt0) cycle
231 : :
232 : 8 : dim_block=0
233 [ + + ]: 24 : do istr=1,nstr
234 : 24 : dim_block=dim_block+dtset%nband(string(istr))
235 : : end do
236 [ - + ]: 8 : if (dim_block <= 0) cycle
237 : :
238 : 8 : string_id=string_id+1
239 [ + - - + ]: 8 : if (use_wfk_restart .and. mod(string_id-1,nproc_cell) /= me_cell) cycle
240 : :
241 : : if (use_wfk_restart) then
242 : : ! wfk_read_my_kptbands packs selected k-points in increasing global-k order.
243 : : ! Build the same offsets here, but retain coefficients for this string only.
244 [ + + ]: 72 : needed_k(:)=.false.
245 [ + + + + : 1744 : distrb_flags(:,:,:)=.false.
+ + ]
246 [ + + ]: 72 : kg_start(:)=-1
247 [ + + ]: 72 : cg_start(:)=-1
248 [ + + ]: 24 : do istr=1,nstr
249 : 16 : ikpt=string(istr)
250 : 16 : needed_k(ikpt)=.true.
251 [ + + ]: 408 : do iband=1,dtset%nband(ikpt)
252 : 400 : distrb_flags(ikpt,iband,1)=.true.
253 : : end do
254 : : end do
255 : :
256 : : ikg_index=0; cg_index=0
257 [ + + ]: 72 : do ikpt=1,dtset%nkpt
258 [ + + ]: 64 : if (.not. needed_k(ikpt)) cycle
259 : 16 : kg_start(ikpt)=ikg_index
260 : 16 : cg_start(ikpt)=cg_index
261 : 16 : ikg_index=ikg_index+npwarr(ikpt)
262 : 72 : cg_index=cg_index+npwarr(ikpt)*dtset%nspinor*dtset%nband(ikpt)
263 : : end do
264 : :
265 [ + + + + ]: 72 : nneeded_k=count(needed_k)
266 [ - + ]: 8 : ABI_CHECK(nneeded_k == nstr, "Wrong number of k-points in full-SOC string")
267 : 8 : mcg_work=dtset%mpw*dtset%nspinor*dtset%mband*nneeded_k*dtset%nsppol
268 [ - + + - : 24 : ABI_MALLOC(kg_work,(3,dtset%mpw*nneeded_k))
+ - + - -
+ - + -
+ ]
269 [ - + + - : 24 : ABI_MALLOC(cg_work,(2,mcg_work))
+ - + - -
+ - + ]
270 [ - + + - : 24 : ABI_MALLOC(eigen_work,(dtset%mband*dtset%nkpt*dtset%nsppol))
+ - + - -
+ - + ]
271 [ - + + - : 16 : ABI_MALLOC(occ_work,(dtset%mband*dtset%nkpt*dtset%nsppol))
- + - + ]
272 : 8 : ecut_eff_restart=dtset%ecut*dtset%dilatmx**2
273 : : call wfk_read_my_kptbands(wfk_path,distrb_flags,xmpi_comm_self,ecut_eff_restart,&
274 : : & 0,dtset%istwfk,dtset%kptns,mcg_work,dtset%mband,dtset%mband,nneeded_k,&
275 : : & dtset%mpw,dtset%natom,dtset%nkpt,npwarr,dtset%nspinor,dtset%nsppol,psps%usepaw,&
276 : 8 : & cg_work,kg=kg_work,eigen=eigen_work,occ=occ_work,ask_accurate_=1)
277 : :
278 : : call gbt_fullsoc_one_string(cg_work,bd_start,cg_start,delta_local,dtset,eigen_work,&
279 : : & gs_hamk_full,kg_work,kg_start,mcg_work,mpi_enreg,next_k,npwarr,occ_work,psps,rprimd,&
280 : 8 : & string,nstr)
281 : :
282 [ - + ]: 8 : ABI_FREE(cg_work)
283 [ - + ]: 8 : ABI_FREE(kg_work)
284 [ - + ]: 8 : ABI_FREE(eigen_work)
285 [ - + ]: 8 : ABI_FREE(occ_work)
286 : : else
287 : : call gbt_fullsoc_one_string(cg,bd_start,cg_start,delta_local,dtset,eigen,&
288 : : & gs_hamk_full,kg,kg_start,mcg,mpi_enreg,next_k,npwarr,occ,psps,rprimd,&
289 : 0 : & string,nstr)
290 : : end if
291 : 18 : nclosed_local=nclosed_local+1
292 : : end do
293 : :
294 : 2 : delta_fullsoc=delta_local
295 : 2 : nclosed_used=nclosed_local
296 [ - + ]: 2 : if (nproc_cell > 1) then
297 : 0 : call xmpi_sum(delta_fullsoc,mpi_enreg%comm_cell,xmpierr)
298 : 0 : call xmpi_sum(nclosed_used,mpi_enreg%comm_cell,xmpierr)
299 : : end if
300 : :
301 : 2 : etot_fullsoc=etotal_psoc+delta_fullsoc
302 [ + - ]: 2 : if (me_cell == 0) then
303 : 2 : write(msg,'(a)')' '
304 [ + + ]: 6 : call wrtout([std_out, ab_out], msg)
305 : 2 : write(msg,'(a,i0,a)')'use_gbt=',dtset%use_gbt,' full-SOC k-string correction:'
306 [ + + ]: 6 : call wrtout([std_out, ab_out], msg)
307 : 2 : write(msg,'(a,3es16.8)')' qgbt: ',dtset%qgbt
308 [ + + ]: 6 : call wrtout([std_out, ab_out], msg)
309 : 2 : write(msg,'(a,i0)')' closed k strings found: ',nclosed
310 [ + + ]: 6 : call wrtout([std_out, ab_out], msg)
311 : 2 : write(msg,'(a,i0)')' closed k strings used: ',nclosed_used
312 [ + + ]: 6 : call wrtout([std_out, ab_out], msg)
313 [ - + ]: 2 : if (nskipped > 0) then
314 : 0 : write(msg,'(a,i0,a)')' WARNING: skipped ',nskipped,' incomplete k strings.'
315 [ # # ]: 0 : call wrtout([std_out, ab_out], msg)
316 : : end if
317 : 2 : write(msg,'(a,t32,es22.14,1x,a)')' Delta E_fullSOC:',delta_fullsoc,'Ha'
318 [ + + ]: 6 : call wrtout([std_out, ab_out], msg)
319 [ + - ]: 2 : if (.not. use_wfk_restart) then
320 [ # # ]: 0 : if (dtset%use_gbt == 3) then
321 : 0 : write(msg,'(a,t32,es22.14,1x,a)')' E_total PSOC reference:',etotal_psoc,'Ha'
322 : : else
323 : 0 : write(msg,'(a,t32,es22.14,1x,a)')' E_total SOC-free reference:',etotal_psoc,'Ha'
324 : : end if
325 [ # # ]: 0 : call wrtout([std_out, ab_out], msg)
326 : 0 : write(msg,'(a,t32,es22.14,1x,a)')' E_total reference+Delta:',etot_fullsoc,'Ha'
327 [ # # ]: 0 : call wrtout([std_out, ab_out], msg)
328 : : end if
329 : : if (use_wfk_restart) then
330 [ + + ]: 2 : if (dtset%use_gbt == 3) then
331 : 1 : write(msg,'(a)')' Add Delta E_fullSOC to the converged use_gbt=2 PSOC total energy.'
332 : : else
333 : 1 : write(msg,'(a)')' Add Delta E_fullSOC to the converged SOC-free GBT total energy.'
334 : : end if
335 : : else
336 : 0 : write(msg,'(a)')' Reference+Delta is reported here; the SCF total_energy is not replaced.'
337 : : end if
338 [ + + ]: 6 : call wrtout([std_out, ab_out], msg)
339 : : end if
340 : :
341 : 2 : call gs_hamk_full%free()
342 [ + - ]: 2 : if (allocated(distrb_flags)) then
343 : 2 : ABI_FREE(distrb_flags)
344 : : end if
345 [ - + ]: 2 : ABI_FREE(needed_k)
346 [ - + ]: 2 : ABI_FREE(visited)
347 [ - + ]: 2 : ABI_FREE(string)
348 [ - + ]: 2 : ABI_FREE(kg_start)
349 [ - + ]: 2 : ABI_FREE(cg_start)
350 [ - + ]: 2 : ABI_FREE(bd_start)
351 [ - + ]: 2 : ABI_FREE(next_k)
352 : :
353 : : DBG_EXIT("COLL")
354 : :
355 [ - + - + : 2 : end subroutine gbt_fullsoc_kstring
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ ]
356 : : !!***
357 : :
358 : : !!****f* m_gbt_fullsoc/gbt_fullsoc_one_string
359 : : !! NAME
360 : : !! gbt_fullsoc_one_string
361 : : !!
362 : : !! FUNCTION
363 : : !! Build and diagonalize one closed q-string block.
364 : : !!
365 : : !! SOURCE
366 : :
367 : 8 : subroutine gbt_fullsoc_one_string(cg,bd_start,cg_start,delta_fullsoc,dtset,eigen,&
368 [ - + ]: 8 : & gs_hamk_full,kg,kg_start,mcg,mpi_enreg,next_k,npwarr,occ,psps,rprimd,string,nstr)
369 : :
370 : : !Arguments ------------------------------------
371 : : !scalars
372 : : integer,intent(in) :: mcg,nstr
373 : : real(dp),intent(inout) :: delta_fullsoc
374 : : type(dataset_type),intent(in) :: dtset
375 : : type(MPI_type),intent(in) :: mpi_enreg
376 : : type(gs_hamiltonian_type),intent(inout) :: gs_hamk_full
377 : : type(pseudopotential_type),intent(in) :: psps
378 : : !arrays
379 : : integer,intent(in) :: bd_start(dtset%nkpt),cg_start(dtset%nkpt),kg(:,:)
380 : : integer,intent(in) :: kg_start(dtset%nkpt),next_k(dtset%nkpt),npwarr(dtset%nkpt),string(nstr)
381 : : real(dp),intent(in) :: cg(2,mcg),eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
382 : : real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),rprimd(3,3)
383 : :
384 : : !Local variables-------------------------------
385 : : !scalars
386 : : integer :: dim_block,iband,ikpt,istr,irow,offset
387 : : real(dp) :: delta_string,weight
388 : : character(len=500) :: msg
389 : : !arrays
390 : 8 : integer,allocatable :: state_offset(:)
391 : 8 : real(dp),allocatable :: eig_full(:),occ_ref(:),eps_ref(:)
392 : 8 : complex(dp),allocatable :: hblock(:,:)
393 : : ! *************************************************************************
394 : :
395 : 8 : dim_block=0
396 [ + + ]: 24 : do istr=1,nstr
397 : 24 : dim_block=dim_block+dtset%nband(string(istr))
398 : : end do
399 : :
400 [ - + + - : 24 : ABI_MALLOC(state_offset,(dtset%nkpt))
+ - + - -
+ - + ]
401 [ - + + - : 32 : ABI_MALLOC(hblock,(dim_block,dim_block))
- + + - +
- + - - +
- + ]
402 [ - + + - : 24 : ABI_MALLOC(eig_full,(dim_block))
+ - + - -
+ - + ]
403 [ - + + - : 16 : ABI_MALLOC(occ_ref,(dim_block))
- + - + ]
404 [ - + + - : 16 : ABI_MALLOC(eps_ref,(dim_block))
- + - + ]
405 [ + + + + ]: 18824 : hblock(:,:)=czero
406 [ + + ]: 72 : state_offset(:)=-1
407 : :
408 : : offset=0
409 [ + + ]: 24 : do istr=1,nstr
410 : 16 : ikpt=string(istr)
411 : 16 : state_offset(ikpt)=offset
412 [ + + ]: 400 : do iband=1,dtset%nband(ikpt)
413 : 384 : irow=offset+iband
414 : 384 : eps_ref(irow)=eigen(bd_start(ikpt)+iband)
415 : 384 : occ_ref(irow)=occ(bd_start(ikpt)+iband)
416 : 400 : hblock(irow,irow)=cmplx(eps_ref(irow),zero,kind=dp)
417 : : end do
418 : 24 : offset=offset+dtset%nband(ikpt)
419 : : end do
420 : :
421 [ + + ]: 24 : do istr=1,nstr
422 : 16 : ikpt=string(istr)
423 [ + + ]: 16 : if (dtset%use_gbt == 4) then
424 : : call gbt_add_fullsoc_samek(cg,cg_start,dtset,gs_hamk_full,hblock,ikpt,&
425 : 8 : & kg,kg_start,mcg,mpi_enreg,npwarr,psps,rprimd,state_offset(ikpt))
426 : : end if
427 : : call gbt_add_fullsoc_link(cg,cg_start,dtset,gs_hamk_full,hblock,ikpt,&
428 : : & kg,kg_start,mcg,mpi_enreg,next_k(ikpt),npwarr,psps,rprimd,&
429 : 24 : & state_offset(ikpt),state_offset(next_k(ikpt)))
430 : : end do
431 : :
432 : 8 : call gbt_sort_eig_occ(eps_ref,occ_ref,dim_block)
433 : 8 : call xheev('N','U',dim_block,hblock,eig_full)
434 : :
435 : 8 : weight=dtset%wtk(string(1))
436 [ + + ]: 16 : do istr=2,nstr
437 [ - + ]: 16 : if (abs(dtset%wtk(string(istr))-weight) > tol10) then
438 : 0 : write(msg,'(a)')' Full-SOC GBT correction: nonuniform weights in one k string;'
439 [ # # ]: 0 : msg=trim(msg)//' using the first k-point weight for this string.'
440 : 0 : ABI_WARNING(msg)
441 : 0 : exit
442 : : end if
443 : : end do
444 : :
445 : 8 : delta_string=zero
446 [ + + ]: 392 : do irow=1,dim_block
447 : 392 : delta_string=delta_string+occ_ref(irow)*(eig_full(irow)-eps_ref(irow))
448 : : end do
449 : 8 : delta_fullsoc=delta_fullsoc+weight*delta_string
450 : :
451 [ - + ]: 8 : ABI_FREE(state_offset)
452 [ - + ]: 8 : ABI_FREE(hblock)
453 [ - + ]: 8 : ABI_FREE(eig_full)
454 [ - + ]: 8 : ABI_FREE(occ_ref)
455 [ - + ]: 8 : ABI_FREE(eps_ref)
456 : :
457 : 8 : end subroutine gbt_fullsoc_one_string
458 : : !!***
459 : :
460 : 24 : subroutine gbt_add_fullsoc_samek(cg,cg_start,dtset,gs_hamk_full,hblock,ikpt,&
461 [ - + - + ]: 8 : & kg,kg_start,mcg,mpi_enreg,npwarr,psps,rprimd,offset)
462 : :
463 : : ! Add the spin-conserving full-SOC blocks at logical k. The two GBT spin
464 : : ! components live at the physical wavevectors k-q/2 and k+q/2, respectively.
465 : :
466 : : integer,intent(in) :: ikpt,mcg,offset
467 : : type(dataset_type),intent(in) :: dtset
468 : : type(MPI_type),intent(in) :: mpi_enreg
469 : : type(gs_hamiltonian_type),intent(inout) :: gs_hamk_full
470 : : type(pseudopotential_type),intent(in) :: psps
471 : : integer,intent(in) :: cg_start(dtset%nkpt),kg(:,:),kg_start(dtset%nkpt),npwarr(dtset%nkpt)
472 : : real(dp),intent(in) :: cg(2,mcg),rprimd(3,3)
473 : : complex(dp),intent(inout) :: hblock(:,:)
474 : :
475 : : integer,parameter :: choice=1,cpopt=-1,ider=0,idir=0,nkpg=0,nnlout=1
476 : : integer,parameter :: only_so=1,paw_opt=0,signs=2,tim_nonlop=3
477 : : integer :: ib_bra,ib_ket,idx,npw,ispin,irow,jcol
478 : : real(dp) :: dot_i,dot_r,lambda(1),rmet(3,3),ucvol
479 : : real(dp) :: gmet(3,3),gprimd(3,3),kphys(3),kpg_dum(0,0)
480 : 8 : type(pawcprj_type),allocatable :: cprj_dum(:,:)
481 : 8 : integer,allocatable,target :: kg_k(:,:)
482 : 8 : real(dp),allocatable :: enlout(:),nonlop_dum(:,:),vectin(:,:),vectout(:,:)
483 : 8 : real(dp),allocatable,target :: ffnl(:,:,:,:),ph3d(:,:,:)
484 : 8 : real(dp),allocatable :: ylm(:,:),ylmgr_dum(:,:,:)
485 : :
486 : 8 : npw=npwarr(ikpt)
487 [ + - ]: 8 : if (npw <= 0) return
488 : 8 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
489 [ - + + - : 24 : ABI_MALLOC(kg_k,(3,npw))
+ - - + -
+ - + ]
490 [ + + + + ]: 20268 : kg_k(:,:)=kg(:,kg_start(ikpt)+1:kg_start(ikpt)+npw)
491 [ - + + - : 32 : ABI_MALLOC(ylm,(npw,psps%mpsang*psps%mpsang*psps%useylm))
+ - - - +
- - + - +
- + ]
492 [ - + ]: 8 : ABI_MALLOC(ylmgr_dum,(0,0,0))
493 [ - + + - : 40 : ABI_MALLOC(ffnl,(npw,1,psps%lmnmax,psps%ntypat))
- + + - -
+ + - + -
+ - - + -
+ - + ]
494 [ - + + - : 32 : ABI_MALLOC(ph3d,(2,npw,gs_hamk_full%matblk))
- + + - +
- + - - +
- + - + ]
495 [ - + + - : 24 : ABI_MALLOC(vectin,(2,2*npw))
+ - - + -
+ - + ]
496 [ - + + - : 16 : ABI_MALLOC(vectout,(2,2*npw))
- + - + -
+ ]
497 [ - + - + ]: 8 : ABI_MALLOC(enlout,(1))
498 [ - + - + ]: 8 : ABI_MALLOC(nonlop_dum,(1,1))
499 [ - + - + ]: 8 : ABI_MALLOC(cprj_dum,(0,0))
500 : 8 : lambda(1)=zero
501 : :
502 [ + + ]: 24 : do ispin=1,2
503 [ + + ]: 64 : kphys(:)=dtset%kptns(:,ikpt)+(real(ispin,dp)-1.5_dp)*dtset%qgbt(:)
504 : : call mkffnl(psps%dimekb,1,psps%ekb,ffnl,psps%ffspl,gmet,gprimd,&
505 : : & ider,idir,psps%indlmn,kg_k,kpg_dum,kphys,psps%lmnmax,psps%lnmax,&
506 : : & psps%mpsang,psps%mqgrid_ff,nkpg,npw,psps%ntypat,psps%pspso,&
507 : 16 : & psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm,ylmgr_dum)
508 : : call gs_hamk_full%load_k(kpt_k=kphys,istwf_k=1,npw_k=npw,kg_k=kg_k,ffnl_k=ffnl,&
509 : 16 : & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
510 : :
511 [ + + ]: 408 : do ib_ket=1,dtset%nband(ikpt)
512 [ + + + + ]: 1459104 : vectin(:,:)=zero
513 : 384 : idx=cg_start(ikpt)+(ib_ket-1)*dtset%nspinor*npw+(ispin-1)*npw
514 [ + + + + ]: 729744 : vectin(:,(ispin-1)*npw+1:ispin*npw)=cg(:,idx+1:idx+npw)
515 [ + + + + : 1459488 : vectout(:,:)=zero; enlout(:)=zero
+ + ]
516 : : call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk_full,idir,lambda,mpi_enreg,&
517 : 384 : & 1,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop,vectin,vectout,only_SO=only_so)
518 [ + + ]: 9616 : do ib_bra=1,dtset%nband(ikpt)
519 : 9216 : idx=cg_start(ikpt)+(ib_bra-1)*dtset%nspinor*npw+(ispin-1)*npw
520 : : call dotprod_g(dot_r,dot_i,1,npw,2,cg(:,idx+1:idx+npw),&
521 : 9216 : & vectout(:,(ispin-1)*npw+1:ispin*npw),mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
522 : 9216 : irow=offset+ib_bra; jcol=offset+ib_ket
523 : 9600 : hblock(irow,jcol)=hblock(irow,jcol)+cmplx(dot_r,dot_i,kind=dp)
524 : : end do
525 : : end do
526 : : end do
527 : :
528 [ + - - + : 8 : ABI_FREE(cprj_dum)
- - - - -
+ ]
529 [ - + ]: 8 : ABI_FREE(nonlop_dum)
530 [ - + ]: 8 : ABI_FREE(enlout)
531 [ - + ]: 8 : ABI_FREE(vectout)
532 [ - + ]: 8 : ABI_FREE(vectin)
533 [ - + ]: 8 : ABI_FREE(ph3d)
534 [ - + ]: 8 : ABI_FREE(ffnl)
535 [ - + ]: 8 : ABI_FREE(ylmgr_dum)
536 [ - + ]: 8 : ABI_FREE(ylm)
537 [ - + ]: 8 : ABI_FREE(kg_k)
538 : :
539 [ - + - + : 8 : end subroutine gbt_add_fullsoc_samek
- + - + -
+ - + - +
- + - + -
+ - - - -
- - ]
540 : : !!***
541 : :
542 : : !!****f* m_gbt_fullsoc/gbt_add_fullsoc_link
543 : : !! NAME
544 : : !! gbt_add_fullsoc_link
545 : : !!
546 : : !! FUNCTION
547 : : !! Add <k|H_SO^flip|k+q> matrix elements to one k-string block.
548 : : !!
549 : : !! SOURCE
550 : :
551 : 48 : subroutine gbt_add_fullsoc_link(cg,cg_start,dtset,gs_hamk_full,hblock,ik_bra,&
552 [ - + - + ]: 16 : & kg,kg_start,mcg,mpi_enreg,ik_ket,npwarr,psps,rprimd,offset_bra,offset_ket)
553 : :
554 : : !Arguments ------------------------------------
555 : : !scalars
556 : : integer,intent(in) :: ik_bra,ik_ket,mcg,offset_bra,offset_ket
557 : : type(dataset_type),intent(in) :: dtset
558 : : type(MPI_type),intent(in) :: mpi_enreg
559 : : type(gs_hamiltonian_type),intent(inout) :: gs_hamk_full
560 : : type(pseudopotential_type),intent(in) :: psps
561 : : !arrays
562 : : integer,intent(in) :: cg_start(dtset%nkpt),kg(:,:)
563 : : integer,intent(in) :: kg_start(dtset%nkpt),npwarr(dtset%nkpt)
564 : : real(dp),intent(in) :: cg(2,mcg),rprimd(3,3)
565 : : complex(dp),intent(inout) :: hblock(:,:)
566 : :
567 : : !Local variables-------------------------------
568 : : !scalars
569 : : integer,parameter :: choice=1,cpopt=-1,ider=0,idir=0,nkpg=0,nnlout=1
570 : : integer,parameter :: only_so=1,paw_opt=0,signs=2,tim_nonlop=3
571 : : integer :: ib_bra,ib_ket,idx_bra,idx_ket,irow,jcol,npw_bra,npw_ket
572 : : real(dp) :: dot_i,dot_r,lambda(1),rmet(3,3),ucvol
573 : : real(dp) :: gmet(3,3),gprimd(3,3),k_bra(3),k_ket(3),kpg_dum(0,0)
574 : 16 : type(pawcprj_type),allocatable :: cprj_dum(:,:)
575 : : !arrays
576 : 16 : integer,allocatable,target :: kg_bra(:,:),kg_ket(:,:)
577 : 16 : real(dp),allocatable :: enlout(:),nonlop_dum(:,:),vectin(:,:),vectout(:,:)
578 : 16 : real(dp),allocatable,target :: ffnl_bra(:,:,:,:),ffnl_ket(:,:,:,:),ph3d_bra(:,:,:)
579 : 16 : real(dp),allocatable,target :: ph3d_ket(:,:,:)
580 : 16 : real(dp),allocatable :: ylm_bra(:,:),ylm_ket(:,:),ylmgr_dum(:,:,:)
581 : : ! *************************************************************************
582 : :
583 : 16 : npw_bra=npwarr(ik_bra)
584 : 16 : npw_ket=npwarr(ik_ket)
585 [ + - ]: 16 : if (npw_bra <= 0 .or. npw_ket <= 0) return
586 : :
587 [ + + ]: 64 : k_ket(:)=dtset%kptns(:,ik_ket)-half*dtset%qgbt(:)
588 [ + + ]: 64 : k_bra(:)=dtset%kptns(:,ik_bra)+half*dtset%qgbt(:)
589 : :
590 : 16 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
591 : :
592 [ - + + - : 48 : ABI_MALLOC(kg_ket,(3,npw_ket))
+ - - + -
+ - + ]
593 [ - + + - : 48 : ABI_MALLOC(kg_bra,(3,npw_bra))
+ - - + -
+ - + ]
594 [ + + + + ]: 40536 : kg_ket(:,:)=kg(:,kg_start(ik_ket)+1:kg_start(ik_ket)+npw_ket)
595 [ + + + + ]: 40536 : kg_bra(:,:)=kg(:,kg_start(ik_bra)+1:kg_start(ik_bra)+npw_bra)
596 : :
597 [ - + + - : 64 : ABI_MALLOC(ylm_ket,(npw_ket,psps%mpsang*psps%mpsang*psps%useylm))
+ - - - +
- - + - +
- + ]
598 [ - + + - : 64 : ABI_MALLOC(ylm_bra,(npw_bra,psps%mpsang*psps%mpsang*psps%useylm))
+ - - - +
- - + - +
- + ]
599 [ - + ]: 16 : ABI_MALLOC(ylmgr_dum,(0,0,0))
600 [ - + + - : 80 : ABI_MALLOC(ffnl_ket,(npw_ket,1,psps%lmnmax,psps%ntypat))
- + + - -
+ + - + -
+ - - + -
+ - + ]
601 [ - + + - : 80 : ABI_MALLOC(ffnl_bra,(npw_bra,1,psps%lmnmax,psps%ntypat))
- + + - -
+ + - + -
+ - - + -
+ - + ]
602 : :
603 : : call mkffnl(psps%dimekb,1,psps%ekb,ffnl_ket,psps%ffspl,gmet,gprimd,&
604 : : & ider,idir,psps%indlmn,kg_ket,kpg_dum,k_ket,psps%lmnmax,psps%lnmax,&
605 : : & psps%mpsang,psps%mqgrid_ff,nkpg,npw_ket,psps%ntypat,psps%pspso,&
606 : 16 : & psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_ket,ylmgr_dum)
607 : : call mkffnl(psps%dimekb,1,psps%ekb,ffnl_bra,psps%ffspl,gmet,gprimd,&
608 : : & ider,idir,psps%indlmn,kg_bra,kpg_dum,k_bra,psps%lmnmax,psps%lnmax,&
609 : : & psps%mpsang,psps%mqgrid_ff,nkpg,npw_bra,psps%ntypat,psps%pspso,&
610 : 16 : & psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_bra,ylmgr_dum)
611 : :
612 [ - + + - : 64 : ABI_MALLOC(ph3d_ket,(2,npw_ket,gs_hamk_full%matblk))
- + + - +
- + - - +
- + - + ]
613 [ - + + - : 64 : ABI_MALLOC(ph3d_bra,(2,npw_bra,gs_hamk_full%matblk))
- + + - +
- + - - +
- + - + ]
614 : : call gs_hamk_full%load_k(kpt_k=k_ket,istwf_k=1,npw_k=npw_ket,&
615 : : & kg_k=kg_ket,ffnl_k=ffnl_ket,ph3d_k=ph3d_ket,&
616 : 16 : & compute_ph3d=.true.,compute_gbound=.true.)
617 : : call gs_hamk_full%load_kprime(kpt_kp=k_bra,istwf_kp=1,npw_kp=npw_bra,&
618 : : & kg_kp=kg_bra,ffnl_kp=ffnl_bra,ph3d_kp=ph3d_bra,&
619 : 16 : & compute_ph3d=.true.,compute_gbound=.true.)
620 : :
621 [ - + + - : 48 : ABI_MALLOC(vectin,(2,2*npw_ket))
+ - - + -
+ - + ]
622 [ - + + - : 48 : ABI_MALLOC(vectout,(2,2*npw_bra))
+ - - + -
+ - + ]
623 [ - + - + ]: 16 : ABI_MALLOC(enlout,(1))
624 [ - + - + ]: 16 : ABI_MALLOC(nonlop_dum,(1,1))
625 [ - + - + ]: 16 : ABI_MALLOC(cprj_dum,(0,0))
626 : 16 : lambda(1)=zero
627 : :
628 [ + + ]: 400 : do ib_ket=1,dtset%nband(ik_ket)
629 [ + + + + ]: 1459104 : vectin(:,:)=zero
630 : 384 : idx_ket=cg_start(ik_ket)+(ib_ket-1)*dtset%nspinor*npw_ket
631 [ + + + + ]: 729744 : vectin(:,1:npw_ket)=cg(:,idx_ket+1:idx_ket+npw_ket)
632 [ + + + + ]: 1459104 : vectout(:,:)=zero
633 [ + + ]: 768 : enlout(:)=zero
634 : :
635 : : call nonlop(choice,cpopt,cprj_dum,enlout,gs_hamk_full,idir,lambda,mpi_enreg,&
636 : 384 : & 1,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop,vectin,vectout,only_SO=only_so)
637 : :
638 [ + + ]: 9616 : do ib_bra=1,dtset%nband(ik_bra)
639 : 9216 : idx_bra=cg_start(ik_bra)+(ib_bra-1)*dtset%nspinor*npw_bra+npw_bra
640 : : call dotprod_g(dot_r,dot_i,1,npw_bra,2,cg(:,idx_bra+1:idx_bra+npw_bra),&
641 : 9216 : & vectout(:,npw_bra+1:2*npw_bra),mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
642 : 9216 : irow=offset_bra+ib_bra
643 : 9216 : jcol=offset_ket+ib_ket
644 : 9216 : hblock(irow,jcol)=hblock(irow,jcol)+cmplx(dot_r,dot_i,kind=dp)
645 : 9600 : hblock(jcol,irow)=hblock(jcol,irow)+cmplx(dot_r,-dot_i,kind=dp)
646 : : end do
647 : : end do
648 : :
649 [ - + ]: 16 : ABI_FREE(vectin)
650 [ - + ]: 16 : ABI_FREE(vectout)
651 [ - + ]: 16 : ABI_FREE(enlout)
652 [ - + ]: 16 : ABI_FREE(nonlop_dum)
653 [ + - - + : 16 : ABI_FREE(cprj_dum)
- - - - -
+ ]
654 [ - + ]: 16 : ABI_FREE(ph3d_ket)
655 [ - + ]: 16 : ABI_FREE(ph3d_bra)
656 [ - + ]: 16 : ABI_FREE(ffnl_ket)
657 [ - + ]: 16 : ABI_FREE(ffnl_bra)
658 [ - + ]: 16 : ABI_FREE(ylm_ket)
659 [ - + ]: 16 : ABI_FREE(ylm_bra)
660 [ - + ]: 16 : ABI_FREE(ylmgr_dum)
661 [ - + ]: 16 : ABI_FREE(kg_ket)
662 [ - + ]: 16 : ABI_FREE(kg_bra)
663 : :
664 [ - + - + : 16 : end subroutine gbt_add_fullsoc_link
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
- - - -
- ]
665 : : !!***
666 : :
667 : 16 : integer function gbt_find_kpoint(target,kptns,nkpt) result(ik_found)
668 : :
669 : : !Arguments ------------------------------------
670 : : integer,intent(in) :: nkpt
671 : : real(dp),intent(in) :: target(3),kptns(3,nkpt)
672 : :
673 : : !Local variables-------------------------------
674 : : integer :: g0(3),ikpt
675 : : ! *************************************************************************
676 : :
677 : 16 : ik_found=0
678 [ + - ]: 72 : do ikpt=1,nkpt
679 [ + + ]: 72 : if (isamek(target,kptns(:,ikpt),g0)) then
680 : 16 : ik_found=ikpt
681 : : return
682 : : end if
683 : : end do
684 : :
685 : : end function gbt_find_kpoint
686 : : !!***
687 : :
688 : 8 : subroutine gbt_sort_eig_occ(eps,occs,nstate)
689 : :
690 : : !Arguments ------------------------------------
691 : : integer,intent(in) :: nstate
692 : : real(dp),intent(inout) :: eps(nstate),occs(nstate)
693 : :
694 : : !Local variables-------------------------------
695 : : integer :: ii,jj
696 : : real(dp) :: eps_tmp,occ_tmp
697 : : ! *************************************************************************
698 : :
699 [ + + ]: 384 : do ii=2,nstate
700 : 376 : eps_tmp=eps(ii)
701 : 376 : occ_tmp=occs(ii)
702 : 376 : jj=ii-1
703 [ + + ]: 2684 : do while (jj >= 1)
704 [ + + ]: 2680 : if (eps(jj) <= eps_tmp) exit
705 : 2308 : eps(jj+1)=eps(jj)
706 : 2308 : occs(jj+1)=occs(jj)
707 : 2680 : jj=jj-1
708 : : end do
709 : 376 : eps(jj+1)=eps_tmp
710 : 384 : occs(jj+1)=occ_tmp
711 : : end do
712 : :
713 : 8 : end subroutine gbt_sort_eig_occ
714 : : !!***
715 : :
716 : : end module m_gbt_fullsoc
717 : : !!***
|
