Line data Source code
1 : !!*** ABINIT/m_orbmag
2 : !! NAME
3 : !! m_orbmag
4 : !!
5 : !! FUNCTION
6 : !! This module contains the declaration of data types and methods
7 : !! used to handle orbital magnetization
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2011-2026 ABINIT group (JWZ)
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 : !! NOTES
16 : !! These routines implement the theory developed in Zwanziger, Torrent, Gonze
17 : !! Phys Rev B 107, 165157 (2023). This paper will be referred to in the comments as ZTG23.
18 : !!
19 : !! SOURCE
20 :
21 : #if defined HAVE_CONFIG_H
22 : #include "config.h"
23 : #endif
24 :
25 : #include "abi_common.h"
26 :
27 : ! note: in a typical index over lmn2_size, think of it as row ilmn, column jlmn for element d_ij.
28 : ! In pawinit line 356, klmn is constructed such that ilmn <= jlmn. Thus we have the upper triangular
29 : ! part of the dij matrix. When looping over both ilmn and jlmn, element dij with i>j is constructed
30 : ! by symmetry from element dji.
31 : #define MATPACK(row,col) (MAX(row,col)*(MAX(row,col)-1)/2 + MIN(row,col))
32 :
33 : #define EPSABG(i,j,k) (((i - j) * (j - k) * (k - i)) / 2)
34 :
35 : module m_orbmag
36 :
37 : use defs_basis
38 : use m_errors
39 : use m_abicore
40 : use m_xmpi
41 : use m_dtset
42 :
43 : use defs_datatypes, only : pseudopotential_type
44 : use defs_abitypes, only : MPI_type
45 : use m_crystal, only : crystal_t
46 : use m_cgprj, only : getcprj
47 : use m_cgtools, only : cg_zdotc,cg_zdotu,projbd
48 : use m_dtfil
49 : use m_ebands
50 : use m_fft, only : fourwf
51 : use m_getghc, only : getghc
52 : use m_getgh1c
53 : use m_hamiltonian
54 : use m_hdr
55 : use m_kg, only : getph,mkkin,mkkpg,ph1d3d
56 : use m_mkffnl, only : mkffnl
57 : use m_mpinfo, only : proc_distrb_cycle,proc_distrb_nband
58 : use m_nctk
59 : use netcdf
60 : use m_nonlop, only : nonlop
61 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free,pawcprj_getdim, pawcprj_get, pawcprj_put
62 : use m_pawdij, only : pawv1
63 : use m_pawfgr, only : pawfgr_type
64 : use m_pawfgrtab, only : pawfgrtab_type
65 : use m_paw_ij, only : paw_ij_type
66 : use m_pawrad, only : nderiv_gen,pawrad_type,pawrad_deducer0,simp_gen,poisson
67 : use m_paw_sphharm, only : setsym_ylm,slxyzs,realgaunt,make_dyadic
68 : use m_pawtab, only : pawtab_type
69 : use m_spacepar, only : make_vectornd
70 : use m_time, only : cwtime, timab
71 :
72 : implicit none
73 :
74 : ! these parameters name the various output terms
75 : integer,parameter :: chern_nterms=2
76 : integer,parameter :: ibcc=1,ibvv1=2
77 : integer,parameter :: orbmag_nterms=6
78 : integer,parameter :: incc=1,invv1=2,invv2=3
79 : integer,parameter :: innl=4,inlr=5,inbm=6
80 :
81 : ! these parameters are constants used repeatedly
82 :
83 : ! accounts for exp(i k.r) in abinit derivatives rather than exp( 2pi i k.r)
84 : real(dp),parameter :: c2=one/(two_pi*two_pi)
85 : complex(dp),parameter :: com=-half*j_dpc ! Orbital magnetism pre-factor
86 : complex(dp),parameter :: cbc=-com ! Berry curvature pre-factor
87 :
88 : ! local datatype for orbmag data on kpt mesh, for eventual output to netcdf
89 : type,private :: orbmag_mesh_type
90 : ! scalars
91 :
92 : integer :: mband, nkpt, nsppol
93 : ! number of bands, kpts, spin polarizations
94 :
95 : integer :: natom, ntypat
96 : ! atoms and types of atoms
97 :
98 : integer :: chern_nterms
99 : ! number of chern terms to store on the kpt mesh
100 : ! CC, VV1, VV2
101 :
102 : integer :: orbmag_nterms
103 : ! number of orbmag terms to store on the kpt mesh
104 : ! CC, VV1, VV2, NL, L_R, B.M
105 :
106 : integer :: n4,n5,n6
107 : ! real space grid dimenions for rmesh
108 :
109 : real(dp),allocatable :: lambsig(:)
110 : ! lambsig(ntypat)
111 :
112 : real(dp),allocatable :: nucdipmom(:,:)
113 : ! nucdipmom(3,natom)
114 :
115 : real(dp),allocatable :: cmesh(:,:,:,:,:)
116 : ! 3 for the 3 directions
117 : ! cmesh(mband,nkpt,nsppol,3,chern_terms)
118 :
119 : real(dp),allocatable :: chern_terms(:,:,:,:)
120 : ! 3 for the 3 directions
121 : ! chern_terms(dtset%mband,dtset%nsppol,3,chern_nterms)
122 :
123 : real(dp),allocatable :: chern_trace(:,:)
124 : ! 3 for the 3 directions
125 : ! chern_trace(3,chern_nterms)
126 :
127 : real(dp),allocatable :: omesh(:,:,:,:,:)
128 : ! 3 for the 3 directions
129 : ! omesh(mband,nkpt,nsppol,3,orbmag_terms)
130 :
131 : real(dp),allocatable :: orbmag_terms(:,:,:,:)
132 : ! 3 for the 3 directions
133 : ! orbmag_terms(dtset%mband,dtset%nsppol,3,orbmag_nterms)
134 :
135 : real(dp),allocatable :: orbmag_trace(:,:)
136 : ! 3 for the 3 directions
137 : ! orbmag_trace(3,orbmag_nterms)
138 :
139 : real(dp),allocatable :: rmesh(:,:,:,:,:)
140 : ! total orbmag on real mesh
141 : ! 3 for the 3 directions
142 : ! rmesh(n4,n5,n6,3,orbmag_nterms)
143 :
144 : real(dp),allocatable :: rmesh_int(:,:)
145 : ! integrated orbmag on real mesh, useful for checking
146 : ! 3 for the 3 directions
147 : ! rmesh_int(3,orbmag_nterms)
148 :
149 : contains
150 :
151 : procedure :: init => orbmag_init
152 : procedure :: free => orbmag_free
153 : procedure :: mpisum => orbmag_mpisum
154 : procedure :: term_scale => orbmag_term_scale
155 : procedure :: output => orbmag_output
156 :
157 : end type orbmag_mesh_type
158 :
159 : ! local datatype for various onsite terms. Probably overkill, but convenient.
160 : type,private :: dterm_type
161 : ! scalars
162 : integer :: lmnmax
163 : integer :: lmn2max
164 : integer :: natom
165 : integer :: ndij
166 : integer :: has_aij=0
167 : integer :: has_qij=0
168 : integer :: has_LR=0
169 : integer :: has_BM=0
170 :
171 : ! sum of \Delta A_ij
172 : ! typically will be just paw_ij
173 : ! aij(natom,lmn2max,ndij)
174 : complex(dp),allocatable :: aij(:,:,:)
175 :
176 : ! <phi|phi> - <tphi|tphi>
177 : ! qij(natom,lmn2max,ndij)
178 : complex(dp),allocatable :: qij(:,:,:)
179 :
180 : ! onsite L_R/2
181 : ! <phi|L_R/2|phi> - <tphi|L_R/2|tphi>
182 : ! ebk_LR(2*lmn2max,natom,ndij,3)
183 : real(dp),allocatable :: ekb_LR(:,:,:,:)
184 :
185 : ! ebk_BM(2*lmn2max,natom,ndij,3)
186 : real(dp),allocatable :: ekb_BM(:,:,:,:)
187 :
188 : contains
189 :
190 : procedure :: init => dterm_init
191 : procedure :: free => dterm_free
192 :
193 : end type dterm_type
194 :
195 : ! Bound methods:
196 :
197 : public :: orbmag
198 : public :: orbmag_ncpp
199 :
200 : private :: orbmag_cc_k
201 : private :: orbmag_cc_k_ncpp
202 : private :: orbmag_vv_k
203 : private :: orbmag_nl_k
204 : private :: nonlop_orbmag_nl
205 : private :: orbmag_nl1_k
206 : private :: me_proj_mesh
207 : ! private :: cprj_test
208 : private :: make_d
209 : private :: dterm_aij
210 : private :: dterm_qij
211 : private :: dterm_LR
212 : private :: dterm_BM
213 : private :: local_fermie
214 :
215 : private :: lamb_core
216 : private :: make_pcg1
217 : private :: gauge_treatment
218 : private :: para_to_diag
219 : private :: orbmag_init
220 : private :: orbmag_free
221 : private :: orbmag_mpisum
222 : private :: orbmag_term_scale
223 : private :: orbmag_output
224 : private :: orbmag_ncwrite ! Write orbmag_mesh contributions to netcdf file.
225 : private :: dterm_init
226 : private :: dterm_free
227 :
228 : CONTAINS !========================================================================================
229 : !!***
230 :
231 : !!****f* ABINIT/orbmag_ncpp
232 : !! NAME
233 : !! orbmag_ncpp
234 : !!
235 : !! FUNCTION
236 : !! This routine computes the orbital magnetization and Berry curvature based on input
237 : !! wavefunctions and DDK wavefuntions, specialized for NCPP
238 : !!
239 : !! INPUTS
240 : !! cg(2,mcg)=all ground state wavefunctions
241 : !! cg1(2,mcg1,3)=all DDK wavefunctions in all 3 directions
242 : !! cprj(dtset%natom,mcprj)<type(pawcprj_type)>=all ground state cprj
243 : !! crystal(crystal_t)=structured datatype holding details about unit cell
244 : !! dtfil <type(datafiles_type)>=variables related to files
245 : !! dtset <type(dataset_type)>=all input variables for this dataset
246 : !! ebands_k(ebands_t)=structured datatype holding GS eigenvalues
247 : !! gsqcut=large sphere cut-off
248 : !! hdr(hdr_type)=structured dataype with header info for eventual output
249 : !! kg(3,mpw*mkmem_rbz)=basis sphere of planewaves at k
250 : !! mcg=dimension of cg
251 : !! mcg1=dimension of cg1
252 : !! mcprj=dimension of cprj
253 : !! mkmem_rbz=kpts in memory
254 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
255 : !! mpw=max number of planewaves at k
256 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
257 : !! ngfftf(18)=FFT grid size information (from pawfgr%ngfft)
258 : !! paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
259 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
260 : !! pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data
261 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
262 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
263 : !! rprimd(3,3)=real space translation vectors
264 : !! usevxctau=1 if kinetic energy density contribution has to be included (mGGA)
265 : !! vtrial(nfftf,dtset%nspden)=GS potential (Hartree)
266 : !! vxctau(nfftf,nspden,4*usevxctau)=derivative of e_xc with respect to kinetic energy density, for mGGA
267 : !! ylm(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)=all ylm's
268 : !! ylmgr(mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm)=gradients of ylm's
269 : !!
270 : !! OUTPUT
271 : !! only printing in call to orbmag_output
272 : !!
273 : !! SIDE EFFECTS
274 : !!
275 : !! TODO
276 : !!
277 : !! NOTES
278 : !! See Zwanziger, Torrent, and Gonze Phys Rev B 107, 165157 (2023), "ZTG23"
279 : !! DDK wavefunctions are used for the derivatives.
280 : !!
281 : !! SOURCE
282 :
283 1 : subroutine orbmag_ncpp(cg,cg1,dtfil,dtset,crystal,ebands_k,kg,gsqcut,hdr,&
284 1 : & mcg,mcg1,mkmem_rbz,mpi_enreg,mpw,nfftf,ngfftf,pawfgr,pawtab,psps,&
285 1 : & usevxctau,vtrial,vxctau)
286 :
287 : !Arguments ------------------------------------
288 : !scalars
289 : integer,intent(in) :: mcg,mcg1,mkmem_rbz,mpw,nfftf,usevxctau
290 : real(dp),intent(in) :: gsqcut
291 : type(crystal_t),intent(in) :: crystal
292 : type(datafiles_type), intent(in) :: dtfil
293 : type(dataset_type),intent(in) :: dtset
294 : type(ebands_t),intent(in) :: ebands_k
295 : type(hdr_type),intent(in) :: hdr
296 : type(MPI_type), intent(inout) :: mpi_enreg
297 : type(pawfgr_type),intent(in) :: pawfgr
298 : type(pseudopotential_type), intent(in) :: psps
299 :
300 : !arrays
301 : integer,intent(in) :: kg(3,mpw*mkmem_rbz),ngfftf(18)
302 : real(dp),intent(in) :: cg(2,mcg),cg1(2,mcg1,3)
303 : real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden)
304 : real(dp),intent(inout) :: vxctau(nfftf,dtset%nspden,4*usevxctau)
305 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
306 :
307 : !Local
308 : !scalars
309 : character(len=fnlen) :: fname
310 : integer :: bdtot_index,dimffnl,exchn2n3d,iatom,icg,ider,idir,ikg,ikg1,ikpt
311 : integer :: indx,isppol,istwf_k,itypat
312 : integer :: mcgk,me,my_nspinor,nband_k,nband_me,ncid
313 : integer :: ngfft1,ngfft2,ngfft3,ngfft4,ngfft5,ngfft6
314 : integer :: nkpg,nn,nproc,npw_k,npwsp,nucdip_dirs,spaceComm,usecprj_local=0
315 : integer,parameter :: master=0
316 : real(dp) :: ecut_eff,fermie
317 : logical :: has_nucdip
318 1 : type(gs_hamiltonian_type) :: gs_hamk
319 1 : type(orbmag_mesh_type) :: orbmag_mesh
320 :
321 : !arrays
322 1 : integer,allocatable :: atindx(:),atindx1(:),kg_k(:,:),nattyp(:)
323 : real(dp) :: kpoint(3),omlamb(3)
324 1 : real(dp),allocatable :: cg1_k(:,:,:),eig_k(:)
325 1 : real(dp),allocatable,target :: cg_k(:,:)
326 1 : real(dp),allocatable :: ffnl_k(:,:,:,:),kinpw(:),kpg_k(:,:),occ_k(:),ph1d(:,:),ph3d(:,:,:)
327 1 : real(dp),allocatable :: trnrm(:),vlocal(:,:,:,:),vectornd(:,:,:),vectornd_pac(:,:,:,:,:)
328 1 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:),ylm_k(:,:),ylmgr_k(:,:,:)
329 1 : type(paw_ij_type),allocatable :: ncpp_paw_ij(:)
330 :
331 : !----------------------------------------------
332 :
333 : ! set up basic FFT parameters
334 1 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
335 1 : istwf_k = 1
336 1 : spaceComm=mpi_enreg%comm_cell
337 1 : nproc=xmpi_comm_size(spaceComm)
338 1 : me = mpi_enreg%me_kpt
339 1 : ngfft1=dtset%ngfft(1) ; ngfft2=dtset%ngfft(2) ; ngfft3=dtset%ngfft(3)
340 1 : ngfft4=dtset%ngfft(4) ; ngfft5=dtset%ngfft(5) ; ngfft6=dtset%ngfft(6)
341 1 : ecut_eff = dtset%ecut*(dtset%dilatmx)**2
342 1 : exchn2n3d = 0; ikg1 = 0
343 :
344 : ! Fermi energy
345 1 : call local_fermie(dtset,ebands_k,fermie,mpi_enreg)
346 :
347 : ! initialize orbmag_mesh datatype
348 1 : call orbmag_mesh%init(dtset)
349 10 : orbmag_mesh%nucdipmom=dtset%nucdipmom
350 : ! if user input lambsig specifically in the input file, use it
351 3 : if ( any ( abs(dtset%lambsig).GT.tol8 ) ) then
352 0 : orbmag_mesh%lambsig=dtset%lambsig
353 : ! else use the value read in to pawtab structure (which might well be zero)
354 1 : else if (psps%usepaw .EQ. 1) then
355 0 : orbmag_mesh%lambsig=pawtab(1:dtset%ntypat)%lamb_shielding
356 : else
357 3 : orbmag_mesh%lambsig=zero
358 : end if
359 :
360 :
361 : !Definition of atindx array
362 : !Generate an index table of atoms, in order for them to be used type after type.
363 3 : ABI_MALLOC(atindx,(dtset%natom))
364 2 : ABI_MALLOC(atindx1,(dtset%natom))
365 3 : ABI_MALLOC(nattyp,(psps%ntypat))
366 3 : indx=1
367 3 : do itypat=1,psps%ntypat
368 2 : nattyp(itypat)=0
369 7 : do iatom=1,dtset%natom
370 6 : if(dtset%typat(iatom)==itypat)then
371 2 : atindx(iatom)=indx
372 2 : atindx1(indx)=iatom
373 2 : indx=indx+1
374 2 : nattyp(itypat)=nattyp(itypat)+1
375 : end if
376 : end do
377 : end do
378 :
379 1 : ABI_MALLOC(ncpp_paw_ij,(0))
380 3 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
381 1 : call getph(atindx,dtset%natom,ngfft1,ngfft2,ngfft3,ph1d,crystal%xred)
382 :
383 : !==== Initialize most of the Hamiltonian ====
384 : !Allocate all arrays and initialize quantities that do not depend on k and spin.
385 : !gs_hamk is the normal hamiltonian at k
386 : call gs_hamk%init(psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,&
387 : & dtset%typat,crystal%xred,dtset%nfft,dtset%mgfft,dtset%ngfft,crystal%rprimd,&
388 : & dtset%nloalg,paw_ij=ncpp_paw_ij,ph1d=ph1d,usecprj=usecprj_local,&
389 1 : & nucdipmom=dtset%nucdipmom)
390 :
391 : ! iterate over spin channels
392 1 : bdtot_index=0
393 1 : icg = 0
394 2 : do isppol = 1, dtset%nsppol
395 :
396 : !========= construct local potential ==================
397 6 : ABI_MALLOC(vlocal,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc))
398 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
399 1 : & dtset%nspden, gs_hamk%nvloc, 1, pawfgr, mpi_enreg, vtrial, vlocal)
400 1 : call gs_hamk%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.)
401 :
402 : !======== compute nuclear dipole vector potential (may be zero) ==========
403 5 : has_nucdip = ANY( ABS(dtset%nucdipmom) .GT. tol8 )
404 1 : if(has_nucdip) then
405 1 : nucdip_dirs=3
406 5 : ABI_MALLOC(vectornd,(nfftf,dtset%nspden,nucdip_dirs))
407 12295 : vectornd = zero
408 : call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,&
409 1 : & dtset%nspden,dtset%nucdipmom,crystal%rprimd,vectornd,crystal%xred)
410 7 : ABI_MALLOC(vectornd_pac,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,nucdip_dirs))
411 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
412 1 : & dtset%nspden, gs_hamk%nvloc, nucdip_dirs, pawfgr, mpi_enreg, vectornd,vectornd_pac)
413 1 : ABI_FREE(vectornd)
414 1 : call gs_hamk%load_spin(isppol,vectornd=vectornd_pac)
415 : else
416 0 : nucdip_dirs=0
417 : end if
418 :
419 : ! !======== compute vxctaulocal if vxctau present =====================
420 :
421 1 : if (usevxctau==1) then
422 0 : ABI_MALLOC(vxctaulocal,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,4))
423 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
424 0 : & dtset%nspden, gs_hamk%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
425 0 : call gs_hamk%load_spin(isppol, vxctaulocal=vxctaulocal)
426 : end if
427 :
428 1 : ikg = 0
429 : !============= BIG FAT KPT LOOP :) ===========================
430 33 : do ikpt = 1, dtset%nkpt
431 :
432 32 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
433 32 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
434 :
435 : ! if the current kpt is not on the current processor, cycle
436 32 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
437 32 : bdtot_index=bdtot_index+nband_k
438 : cycle
439 : end if
440 :
441 128 : kpoint(:)=ebands_k%kptns(:,ikpt)
442 32 : npw_k = ebands_k%npwarr(ikpt)
443 32 : npwsp = npw_k*dtset%nspinor
444 :
445 : ! retrieve kg_k at this k point
446 96 : ABI_MALLOC(kg_k,(3,npw_k))
447 18752 : kg_k(1:3,1:npw_k) = kg(1:3,ikg+1:ikg+npw_k)
448 :
449 : ! retrieve ylm at this k point
450 128 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
451 128 : ABI_MALLOC(ylmgr_k,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm))
452 :
453 : ! retrieve occupation numbers at this k point
454 96 : ABI_MALLOC(occ_k,(nband_k))
455 160 : occ_k(:)=ebands_k%occ(1:nband_k,ikpt,isppol)
456 :
457 : ! Compute kinetic energy at kpt
458 96 : ABI_MALLOC(kinpw,(npw_k))
459 4712 : kinpw(:) = zero
460 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,crystal%gmet,&
461 32 : & kg_k,kinpw,kpoint,npw_k,0,0)
462 :
463 : ! Compute k+G at this k point (relevant only for useylm=1)
464 32 : nkpg=3*dtset%optforces*dtset%nloalg(3)
465 128 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
466 32 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
467 :
468 128 : ABI_MALLOC(ph3d,(2,npw_k,dtset%natom))
469 :
470 : ! Compute nonlocal form factors ffnl at all (k+G):
471 32 : ider=0;idir=0;dimffnl=1
472 160 : ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,dtset%ntypat))
473 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_k,psps%ffspl,&
474 : & crystal%gmet,crystal%gprimd,ider,idir,psps%indlmn,&
475 : & kg_k,kpg_k,kpoint,psps%lmnmax,&
476 : & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
477 : & npw_k,dtset%ntypat,psps%pspso,psps%qgrid_ff,crystal%rmet,&
478 32 : & psps%usepaw,psps%useylm,ylm_k,ylmgr_k,kinpw=kinpw)
479 : ! - Load k-dependent quantities in the Hamiltonian
480 : call gs_hamk%load_k(kpt_k=kpoint(:),istwf_k=istwf_k,npw_k=npw_k,&
481 : & kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl_k,ph3d_k=ph3d,&
482 32 : & compute_ph3d=.TRUE.,compute_gbound=(mpi_enreg%paral_kgb.NE.1))
483 :
484 32 : ABI_SFREE(ylm_k)
485 32 : ABI_SFREE(ylmgr_k)
486 :
487 : ! retrieve ground state wavefunctions at this k point and isppol
488 32 : mcgk = npw_k*nband_k*dtset%nspinor
489 96 : ABI_MALLOC(cg_k,(2,mcgk))
490 56224 : cg_k = cg(1:2,icg+1:icg+mcgk)
491 :
492 : ! retrieve first order wavefunctions at this k point and isppol
493 128 : ABI_MALLOC(cg1_k,(2,mcgk,3))
494 168640 : cg1_k = cg1(1:2,icg+1:icg+mcgk,1:3)
495 :
496 : ! retrieve zeroth order eigenvalues at this k point and isppol
497 64 : ABI_MALLOC(eig_k,(nband_k))
498 160 : eig_k(:)=ebands_k%eig(1:nband_k,ikpt,isppol)
499 :
500 : ! set up normalization factors at this k point
501 64 : ABI_MALLOC(trnrm,(nband_k))
502 160 : trnrm(1:nband_k) = ebands_k%occ(1:nband_k,ikpt,isppol)*dtset%wtk(ikpt)/crystal%ucvol
503 :
504 : !--------------------------------------------------------------------------------
505 : ! Finally ready to compute contributions to orbital magnetism and Berry curvature
506 : !--------------------------------------------------------------------------------
507 :
508 : ! ZTG23 Eq. 36 term 2 and Eq. 46 term 1
509 : call orbmag_cc_k_ncpp(dtset,eig_k,fermie,cg1_k,gs_hamk,ikpt,isppol,&
510 32 : & mcgk,mpi_enreg,nband_k,npw_k,orbmag_mesh)
511 :
512 : ! accumulate terms
513 160 : do nn = 1, nband_k
514 128 : if(abs(trnrm(nn)).LT.tol8) cycle
515 : orbmag_mesh%chern_terms(nn,isppol,1:3,ibcc:ibvv1) = orbmag_mesh%chern_terms(nn,isppol,1:3,ibcc:ibvv1) + &
516 2176 : & trnrm(nn)*orbmag_mesh%cmesh(nn,ikpt,isppol,1:3,ibcc:ibvv1)
517 : orbmag_mesh%orbmag_terms(nn,isppol,1:3,incc:inbm) = orbmag_mesh%orbmag_terms(nn,isppol,1:3,incc:inbm) + &
518 6304 : & trnrm(nn)*orbmag_mesh%omesh(nn,ikpt,isppol,1:3,incc:inbm)
519 : end do ! loop on bands
520 :
521 32 : icg = icg + mcgk
522 32 : ikg = ikg + npw_k
523 32 : bdtot_index=bdtot_index+nband_k
524 :
525 32 : ABI_SFREE(ffnl_k)
526 32 : ABI_SFREE(ph3d)
527 32 : ABI_SFREE(kpg_k)
528 32 : ABI_SFREE(kg_k)
529 32 : ABI_SFREE(kinpw)
530 32 : ABI_SFREE(cg_k)
531 32 : ABI_SFREE(cg1_k)
532 32 : ABI_SFREE(eig_k)
533 32 : ABI_SFREE(occ_k)
534 33 : ABI_SFREE(trnrm)
535 :
536 : end do ! end loop over kpts
537 :
538 1 : ABI_SFREE(vlocal)
539 1 : ABI_SFREE(vectornd_pac)
540 2 : ABI_SFREE(vxctaulocal)
541 :
542 : end do ! end loop over isppol
543 :
544 : ! accumulate data over processors
545 1 : call orbmag_mesh%mpisum(nproc,spaceComm)
546 :
547 : ! prepare terms for output to abo file
548 1 : call orbmag_mesh%term_scale(crystal,dtset)
549 :
550 : ! get the Lamb term
551 1 : call lamb_core(atindx,dtset,omlamb,pawtab,psps%usepaw)
552 :
553 : ! output raw data to netcdf file for more detailed postprocessing
554 1 : if (me == master) then
555 1 : fname = trim(dtfil%filnam_ds(4))//'_ORBMAG.nc'
556 1 : NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
557 1 : call orbmag_ncwrite(crystal,dtset,ebands_k,hdr,ncid,orbmag_mesh)
558 1 : NCF_CHECK(nf90_close(ncid))
559 : end if
560 :
561 : ! output summary to abo file
562 1 : call orbmag_mesh%output(dtset,omlamb)
563 :
564 : !!---------------------------------------------------
565 : !! deallocate memory
566 : !!---------------------------------------------------
567 :
568 1 : call gs_hamk%free()
569 1 : call orbmag_mesh%free()
570 1 : ABI_SFREE(ncpp_paw_ij)
571 1 : ABI_SFREE(ph1d)
572 1 : ABI_SFREE(atindx)
573 1 : ABI_SFREE(atindx1)
574 1 : ABI_SFREE(nattyp)
575 :
576 :
577 2 : end subroutine orbmag_ncpp
578 : !!***
579 :
580 :
581 : !!****f* ABINIT/orbmag
582 : !! NAME
583 : !! orbmag
584 : !!
585 : !! FUNCTION
586 : !! This routine computes the orbital magnetization and Berry curvature based on input
587 : !! wavefunctions and DDK wavefuntions.
588 : !!
589 : !! INPUTS
590 : !! cg(2,mcg)=all ground state wavefunctions
591 : !! cg1(2,mcg1,3)=all DDK wavefunctions in all 3 directions
592 : !! cprj(dtset%natom,mcprj)<type(pawcprj_type)>=all ground state cprj
593 : !! crystal(crystal_t)=structured datatype holding details about unit cell
594 : !! dtfil <type(datafiles_type)>=variables related to files
595 : !! dtset <type(dataset_type)>=all input variables for this dataset
596 : !! ebands_k(ebands_t)=structured datatype holding GS eigenvalues
597 : !! gsqcut=large sphere cut-off
598 : !! hdr(hdr_type)=structured dataype with header info for eventual output
599 : !! kg(3,mpw*mkmem_rbz)=basis sphere of planewaves at k
600 : !! mcg=dimension of cg
601 : !! mcg1=dimension of cg1
602 : !! mcprj=dimension of cprj
603 : !! mkmem_rbz=kpts in memory
604 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
605 : !! mpw=max number of planewaves at k
606 : !! nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90)
607 : !! ngfftf(18)=FFT grid size information (from pawfgr%ngfft)
608 : !! paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
609 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
610 : !! pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data
611 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
612 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
613 : !! rprimd(3,3)=real space translation vectors
614 : !! usevxctau=1 if kinetic energy density contribution has to be included (mGGA)
615 : !! vtrial(nfftf,dtset%nspden)=GS potential (Hartree)
616 : !! vxctau(nfftf,nspden,4*usevxctau)=derivative of e_xc with respect to kinetic energy density, for mGGA
617 : !! ylm(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)=all ylm's
618 : !! ylmgr(mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm)=gradients of ylm's
619 : !!
620 : !! OUTPUT
621 : !! only printing in call to orbmag_output
622 : !!
623 : !! SIDE EFFECTS
624 : !!
625 : !! TODO
626 : !!
627 : !! NOTES
628 : !! See Zwanziger, Torrent, and Gonze Phys Rev B 107, 165157 (2023), "ZTG23"
629 : !! DDK wavefunctions are used for the derivatives.
630 : !!
631 : !! SOURCE
632 :
633 16 : subroutine orbmag(cg,cg1,cprj,crystal,dtfil,dtset,ebands_k,gsqcut,hdr,kg,mcg,mcg1,&
634 16 : & mcprj,mkmem_rbz,mpi_enreg,mpw,nfftf,ngfftf,paw_ij,pawfgr,pawrad,&
635 16 : & pawtab,psps,usevxctau,vtrial,vxctau,ylm,ylmgr)
636 :
637 : !Arguments ------------------------------------
638 : !scalars
639 : integer,intent(in) :: mcprj,mcg,mcg1,mkmem_rbz,mpw,nfftf,usevxctau
640 : real(dp),intent(in) :: gsqcut
641 : type(crystal_t),intent(in) :: crystal
642 : type(datafiles_type), intent(in) :: dtfil
643 : type(dataset_type),intent(in) :: dtset
644 : type(ebands_t),intent(in) :: ebands_k
645 : type(hdr_type),intent(in) :: hdr
646 : type(MPI_type), intent(inout) :: mpi_enreg
647 : type(pawfgr_type),intent(in) :: pawfgr
648 : type(pseudopotential_type), intent(in) :: psps
649 :
650 : !arrays
651 : integer,intent(in) :: kg(3,mpw*mkmem_rbz),ngfftf(18)
652 : real(dp),intent(in) :: cg(2,mcg),cg1(2,mcg1,3)
653 : real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden)
654 : real(dp),intent(inout) :: vxctau(nfftf,dtset%nspden,4*usevxctau)
655 : real(dp),intent(in) :: ylm(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)
656 : real(dp),intent(in) :: ylmgr(mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm)
657 : type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj)
658 : type(paw_ij_type),intent(inout) :: paw_ij(dtset%natom*psps%usepaw)
659 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw)
660 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
661 :
662 : !Local
663 : !scalars
664 : character(len=fnlen) :: fname
665 : integer :: adir,bdtot_index,choice,cpopt,dimffnl,exchn2n3d
666 : integer :: iat,iatom,icg,icprj,ider,idir
667 : integer :: ikg,ikg1,ikpt,ilm,indx,isppol,istwf_k,itypat,lmn2max
668 : integer :: me,mcgk,mcprjk,my_nspinor,nband_k,nband_me,ncid,ngfft1,ngfft2,ngfft3,ngfft4
669 : integer :: ngfft5,ngfft6,nn,nkpg,npw_k,npwsp,nproc,nucdip_dirs,spaceComm
670 : integer,parameter :: master=0
671 : real(dp) :: arg,ecut_eff,fermie
672 : logical :: has_nucdip
673 16 : type(dterm_type) :: dterm
674 16 : type(gs_hamiltonian_type) :: gs_hamk
675 16 : type(orbmag_mesh_type) :: orbmag_mesh
676 :
677 : !arrays
678 16 : integer,allocatable :: atindx(:),atindx1(:),dimlmn(:),kg_k(:,:),nattyp(:)
679 : real(dp) :: kpoint(3),omlamb(3)
680 16 : real(dp),allocatable :: cg1_k(:,:,:),cwavef(:,:),dkinpw(:,:),eig_k(:)
681 16 : real(dp),allocatable :: ffnl_k(:,:,:,:),kinpw(:),kpg_k(:,:),occ_k(:)
682 16 : real(dp),allocatable,target :: cg_k(:,:),gcg1_k(:,:,:)
683 16 : real(dp),allocatable,target :: ph1d(:,:)
684 16 : real(dp),allocatable :: ph3d(:,:,:),phkxred(:,:),trnrm(:)
685 16 : real(dp),allocatable :: vectornd(:,:,:),vectornd_pac(:,:,:,:,:),vlocal(:,:,:,:)
686 16 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:)
687 16 : real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:)
688 16 : type(pawcprj_type),allocatable :: cprj_k(:,:),cprj1_k(:,:,:),cwaveprj(:,:)
689 :
690 : !----------------------------------------------
691 :
692 : ! set up basic FFT parameters
693 16 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
694 16 : istwf_k = 1
695 16 : spaceComm=mpi_enreg%comm_cell
696 16 : nproc=xmpi_comm_size(spaceComm)
697 16 : me = mpi_enreg%me_kpt
698 16 : ngfft1=dtset%ngfft(1) ; ngfft2=dtset%ngfft(2) ; ngfft3=dtset%ngfft(3)
699 16 : ngfft4=dtset%ngfft(4) ; ngfft5=dtset%ngfft(5) ; ngfft6=dtset%ngfft(6)
700 16 : ecut_eff = dtset%ecut*(dtset%dilatmx)**2
701 16 : exchn2n3d = 0; ikg1 = 0
702 :
703 : ! Fermi energy
704 16 : call local_fermie(dtset,ebands_k,fermie,mpi_enreg)
705 :
706 : !Definition of atindx array
707 : !Generate an index table of atoms, in order for them to be used type after type.
708 48 : ABI_MALLOC(atindx,(dtset%natom))
709 32 : ABI_MALLOC(atindx1,(dtset%natom))
710 48 : ABI_MALLOC(nattyp,(psps%ntypat))
711 16 : indx=1
712 44 : do itypat=1,psps%ntypat
713 28 : nattyp(itypat)=0
714 96 : do iatom=1,dtset%natom
715 80 : if(dtset%typat(iatom)==itypat)then
716 28 : atindx(iatom)=indx
717 28 : atindx1(indx)=iatom
718 28 : indx=indx+1
719 28 : nattyp(itypat)=nattyp(itypat)+1
720 : end if
721 : end do
722 : end do
723 :
724 48 : ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom))
725 16 : call getph(atindx,dtset%natom,ngfft1,ngfft2,ngfft3,ph1d,crystal%xred)
726 :
727 48 : ABI_MALLOC(kg_k,(3,mpw))
728 48 : ABI_MALLOC(kinpw,(mpw))
729 16 : if (abs(dtset%orbmag) .EQ. 3) then
730 0 : ABI_MALLOC(dkinpw,(mpw,3))
731 : end if
732 :
733 32 : ABI_MALLOC(dimlmn,(dtset%natom))
734 16 : call pawcprj_getdim(dimlmn,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
735 :
736 120 : ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor))
737 16 : call pawcprj_alloc(cwaveprj,0,dimlmn)
738 :
739 16 : lmn2max = psps%lmnmax*(psps%lmnmax+1)/2
740 : ! note: in make_d, terms will be filled as iatom using atindx
741 16 : call dterm%init(psps%lmnmax,lmn2max,dtset%natom,paw_ij(1)%ndij)
742 16 : call make_d(atindx,dterm,dtset,crystal%gprimd,paw_ij,pawrad,pawtab,psps)
743 :
744 : ! initialize orbmag_mesh datatype
745 16 : call orbmag_mesh%init(dtset)
746 144 : orbmag_mesh%nucdipmom=dtset%nucdipmom
747 : ! if user input lambsig specifically in the input file, use it
748 32 : if ( any ( abs(dtset%lambsig).GT.tol8 ) ) then
749 46 : orbmag_mesh%lambsig=dtset%lambsig
750 : ! else use the value read in to pawtab structure (which might well be zero)
751 : else
752 14 : orbmag_mesh%lambsig=pawtab(1:dtset%ntypat)%lamb_shielding
753 : end if
754 :
755 : !==== Initialize most of the Hamiltonian ====
756 : !Allocate all arrays and initialize quantities that do not depend on k and spin.
757 : !gs_hamk is the normal hamiltonian at k
758 : call gs_hamk%init(psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,&
759 : & dtset%typat,crystal%xred,dtset%nfft,dtset%mgfft,dtset%ngfft,crystal%rprimd,&
760 16 : & dtset%nloalg,nucdipmom=dtset%nucdipmom,paw_ij=paw_ij)
761 :
762 : ! iterate over spin channels
763 16 : bdtot_index=0
764 16 : icg = 0
765 16 : icprj = 0
766 33 : do isppol = 1, dtset%nsppol
767 :
768 : !========= construct local potential ==================
769 102 : ABI_MALLOC(vlocal,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc))
770 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
771 17 : & dtset%nspden, gs_hamk%nvloc, 1, pawfgr, mpi_enreg, vtrial, vlocal)
772 17 : call gs_hamk%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.)
773 :
774 : !======== compute nuclear dipole vector potential (may be zero) ==========
775 69 : has_nucdip = ANY( ABS(dtset%nucdipmom) .GT. tol8 )
776 17 : if(has_nucdip) then
777 17 : nucdip_dirs=3
778 85 : ABI_MALLOC(vectornd,(nfftf,dtset%nspden,nucdip_dirs))
779 4309757 : vectornd = zero
780 : call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,&
781 17 : & dtset%nspden,dtset%nucdipmom,crystal%rprimd,vectornd,crystal%xred)
782 119 : ABI_MALLOC(vectornd_pac,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,nucdip_dirs))
783 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
784 17 : & dtset%nspden, gs_hamk%nvloc, nucdip_dirs, pawfgr, mpi_enreg, vectornd,vectornd_pac)
785 17 : ABI_FREE(vectornd)
786 17 : call gs_hamk%load_spin(isppol,vectornd=vectornd_pac)
787 : else
788 0 : nucdip_dirs=0
789 : end if
790 :
791 : !======== compute vxctaulocal if vxctau present =====================
792 :
793 17 : if (usevxctau==1) then
794 14 : ABI_MALLOC(vxctaulocal,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,4))
795 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
796 2 : & dtset%nspden, gs_hamk%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
797 2 : call gs_hamk%load_spin(isppol, vxctaulocal=vxctaulocal)
798 : end if
799 :
800 17 : ikg = 0
801 : !============= BIG FAT KPT LOOP :) ===========================
802 692 : do ikpt = 1, dtset%nkpt
803 :
804 675 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
805 675 : nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
806 :
807 : ! if the current kpt is not on the current processor, cycle
808 675 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
809 675 : bdtot_index=bdtot_index+nband_k
810 : cycle
811 : end if
812 :
813 2700 : kpoint(:)=ebands_k%kptns(:,ikpt)
814 675 : npw_k = ebands_k%npwarr(ikpt)
815 675 : npwsp = npw_k*dtset%nspinor
816 :
817 : ! retrieve kg_k at this k point
818 516967 : kg_k(1:3,1:npw_k) = kg(1:3,ikg+1:ikg+npw_k)
819 :
820 : ! retrieve ylm at this k point
821 2700 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang))
822 2700 : ABI_MALLOC(ylmgr_k,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm))
823 3375 : do ilm=1,psps%mpsang*psps%mpsang
824 518992 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
825 1560351 : ylmgr_k(1:npw_k,1:3,ilm)=ylmgr(1+ikg:npw_k+ikg,1:3,ilm)
826 : end do
827 :
828 : ! retrieve occupation numbers at this k point
829 2025 : ABI_MALLOC(occ_k,(nband_k))
830 : !occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
831 4911 : occ_k(:)=ebands_k%occ(1:nband_k,ikpt,isppol)
832 :
833 : ! Compute kinetic energy at kpt
834 132460 : kinpw(:) = zero
835 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,crystal%gmet,&
836 675 : & kg_k,kinpw,kpoint,npw_k,0,0)
837 675 : if (abs(dtset%orbmag).EQ.3) then
838 0 : do adir=1,3
839 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,crystal%gmet,&
840 0 : & kg_k,dkinpw(:,adir),kpoint,npw_k,adir,0)
841 : end do
842 : end if
843 :
844 : ! Compute k+G at this k point
845 675 : nkpg = 3
846 2700 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
847 675 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
848 :
849 : ! Make 3d phase factors
850 2025 : ABI_MALLOC(phkxred,(2,dtset%natom))
851 1766 : do iat = 1, dtset%natom
852 1091 : iatom = atindx(iat)
853 4364 : arg=two_pi*DOT_PRODUCT(kpoint,crystal%xred(:,iat))
854 1766 : phkxred(1,iatom)=DCOS(arg);phkxred(2,iatom)=DSIN(arg)
855 : end do
856 2700 : ABI_MALLOC(ph3d,(2,npw_k,dtset%natom))
857 : call ph1d3d(1,dtset%natom,kg_k,dtset%natom,dtset%natom,&
858 675 : & npw_k,ngfft1,ngfft2,ngfft3,phkxred,ph1d,ph3d)
859 :
860 : ! Compute nonlocal form factors ffnl at all (k+G):
861 675 : ider=1 ! ffnl and 1st derivatives
862 675 : idir=4 ! ignored when ider = 0; idir=0 means d ffnl/ dk in reduced units referenced
863 : ! to reciprocal translations
864 : ! idir=4 meand d ffnl / dk in reduced units referenced to real space
865 : ! translations. rfddk = 1 wavefunctions are computed using this convention.
866 675 : dimffnl=4 ! 1 + number of derivatives
867 4050 : ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,dtset%ntypat))
868 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_k,psps%ffspl,&
869 : & crystal%gmet,crystal%gprimd,ider,idir,psps%indlmn,&
870 : & kg_k,kpg_k,kpoint,psps%lmnmax,&
871 : & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
872 : & npw_k,dtset%ntypat,psps%pspso,psps%qgrid_ff,crystal%rmet,&
873 675 : & psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
874 : ! - Load k-dependent quantities in the Hamiltonian
875 : call gs_hamk%load_k(kpt_k=kpoint(:),istwf_k=istwf_k,npw_k=npw_k,&
876 : & kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl_k,ph3d_k=ph3d,&
877 675 : & compute_gbound=.TRUE.)
878 :
879 675 : ABI_SFREE(ylm_k)
880 675 : ABI_SFREE(ylmgr_k)
881 :
882 : ! retrieve ground state wavefunctions at this k point and isppol
883 675 : mcgk = npw_k*nband_k*dtset%nspinor
884 2025 : ABI_MALLOC(cg_k,(2,mcgk))
885 2664978 : cg_k = cg(1:2,icg+1:icg+mcgk)
886 :
887 : ! retrieve first order wavefunctions at this k point and isppol
888 2700 : ABI_MALLOC(cg1_k,(2,mcgk,3))
889 7994259 : cg1_k = cg1(1:2,icg+1:icg+mcgk,1:3)
890 :
891 : ! retrieve zeroth order eigenvalues at this k point and isppol
892 1350 : ABI_MALLOC(eig_k,(nband_k))
893 : !eig_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
894 4911 : eig_k(:)=ebands_k%eig(1:nband_k,ikpt,isppol)
895 :
896 : ! retrieve cprj_k at this k point and isppol
897 675 : mcprjk = nband_k*dtset%nspinor
898 16420 : ABI_MALLOC(cprj_k,(dtset%natom,mcprjk))
899 675 : call pawcprj_alloc(cprj_k,cprj(1,1)%ncpgr,dimlmn)
900 : call pawcprj_get(atindx,cprj_k,cprj,dtset%natom,1,icprj,ikpt,0,isppol,dtset%mband,&
901 675 : & mkmem_rbz,dtset%natom,nband_k,nband_k,dtset%nspinor,dtset%nsppol,0)
902 :
903 : ! gauge treatment of cg1_k
904 1350 : ABI_MALLOC(gcg1_k,(2,mcgk,3))
905 : call gauge_treatment(atindx,cg_k,cg1_k,cprj_k,dimlmn,dkinpw,dtset,eig_k,gcg1_k,gs_hamk,&
906 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,mpw,nband_k,ngfft4,ngfft5,ngfft6,npw_k,&
907 675 : & nucdip_dirs,occ_k,vectornd_pac)
908 :
909 : ! compute <p|gcg1> cprjs
910 46560 : ABI_MALLOC(cprj1_k,(dtset%natom,mcprjk,3))
911 2700 : do adir = 1, 3
912 2700 : call pawcprj_alloc(cprj1_k(:,:,adir),0,dimlmn)
913 : end do
914 675 : choice = 1
915 675 : cpopt = 0
916 675 : idir = 0
917 2025 : ABI_MALLOC(cwavef,(2,npwsp))
918 4911 : do nn = 1, nband_k
919 17619 : do adir = 1, 3
920 8003592 : cwavef(1:2,1:npwsp) = gcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,adir)
921 : call getcprj(choice,cpopt,cwavef,cwaveprj,gs_hamk%ffnl_k,idir,&
922 : & psps%indlmn,istwf_k,kg_k,gs_hamk%kpg_k,kpoint,psps%lmnmax,dtset%mgfft,&
923 : & mpi_enreg,1,dtset%natom,nattyp,dtset%ngfft,dtset%nloalg,npw_k,&
924 : & dtset%nspinor,dtset%ntypat,phkxred,ph1d,gs_hamk%ph3d_k,&
925 12708 : & crystal%ucvol,psps%useylm)
926 : call pawcprj_put(atindx,cwaveprj,cprj1_k(:,:,adir),dtset%natom,&
927 : & nn,0,ikpt,0,isppol,dtset%mband,mkmem_rbz,dtset%natom,1,nband_k,&
928 16944 : & dimlmn,dtset%nspinor,dtset%nsppol,0)
929 : end do
930 : end do
931 675 : ABI_SFREE(cwavef)
932 675 : ABI_SFREE(phkxred)
933 :
934 : ! set up normalization factors at this k point
935 1350 : ABI_MALLOC(trnrm,(nband_k))
936 4911 : trnrm(1:nband_k) = ebands_k%occ(1:nband_k,ikpt,isppol)*dtset%wtk(ikpt)/crystal%ucvol
937 :
938 : !--------------------------------------------------------------------------------
939 : ! Finally ready to compute contributions to orbital magnetism and Berry curvature
940 : !--------------------------------------------------------------------------------
941 :
942 : ! ZTG23 Eq. 36 term 2 and Eq. 46 term 1
943 : call orbmag_cc_k(atindx,cprj1_k,dimlmn,dtset,eig_k,fermie,gcg1_k,gs_hamk,ikpt,isppol,&
944 675 : & mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,trnrm)
945 :
946 : ! ZTG23 Eq. 36 terms 3 and 4 and Eq. 46 term 2
947 : call orbmag_vv_k(atindx,cg_k,cprj_k,dimlmn,dtset,eig_k,fermie,gcg1_k,gs_hamk,ikpt,isppol,&
948 675 : & mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,trnrm)
949 :
950 : ! ZTG23 Eq. 36 term 1
951 : call orbmag_nl_k(atindx,cg_k,cprj_k,dimlmn,dterm,dtset,eig_k,gs_hamk,ikpt,isppol,&
952 675 : & mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,pawtab,trnrm)
953 :
954 : ! ZTG23 text after Eq. 42
955 : ! <L_R> contribution
956 : call orbmag_nl1_k(atindx,cg_k,cprj_k,dimlmn,dterm,dtset,gs_hamk,ikpt,inlr,isppol,&
957 675 : & mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,trnrm)
958 :
959 : ! ZTG23 Eq. 43
960 : ! A0.An contribution
961 : call orbmag_nl1_k(atindx,cg_k,cprj_k,dimlmn,dterm,dtset,gs_hamk,ikpt,inbm,isppol,&
962 675 : & mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,trnrm)
963 :
964 : ! accumulate terms
965 4911 : do nn = 1, nband_k
966 4236 : if(abs(trnrm(nn)).LT.tol8) cycle
967 : orbmag_mesh%chern_terms(nn,isppol,1:3,ibcc:ibvv1) = orbmag_mesh%chern_terms(nn,isppol,1:3,ibcc:ibvv1) + &
968 57290 : & trnrm(nn)*orbmag_mesh%cmesh(nn,ikpt,isppol,1:3,ibcc:ibvv1)
969 : orbmag_mesh%orbmag_terms(nn,isppol,1:3,incc:inbm) = orbmag_mesh%orbmag_terms(nn,isppol,1:3,incc:inbm) + &
970 166671 : & trnrm(nn)*orbmag_mesh%omesh(nn,ikpt,isppol,1:3,incc:inbm)
971 : end do ! loop on bands
972 :
973 675 : icg = icg + mcgk
974 675 : icprj = icprj + mcprjk
975 675 : ikg = ikg + npw_k
976 675 : bdtot_index=bdtot_index+nband_k
977 :
978 675 : ABI_SFREE(ffnl_k)
979 675 : ABI_SFREE(ph3d)
980 675 : ABI_SFREE(kpg_k)
981 675 : ABI_SFREE(cg_k)
982 675 : ABI_SFREE(cg1_k)
983 675 : ABI_SFREE(gcg1_k)
984 675 : ABI_SFREE(eig_k)
985 675 : ABI_SFREE(occ_k)
986 675 : call pawcprj_free(cprj_k)
987 9135 : ABI_SFREE(cprj_k)
988 2700 : do adir = 1, 3
989 2700 : call pawcprj_free(cprj1_k(:,:,adir))
990 : end do
991 26055 : ABI_SFREE(cprj1_k)
992 692 : ABI_SFREE(trnrm)
993 :
994 : end do ! end loop over kpts
995 :
996 17 : ABI_SFREE(vlocal)
997 17 : ABI_SFREE(vectornd_pac)
998 33 : ABI_SFREE(vxctaulocal)
999 :
1000 : end do ! end loop over isppol
1001 :
1002 : ! accumulate data over processors
1003 16 : call orbmag_mesh%mpisum(nproc,spaceComm)
1004 :
1005 : ! prepare terms for output to abo file
1006 16 : call orbmag_mesh%term_scale(crystal,dtset)
1007 :
1008 : ! get the Lamb term
1009 16 : call lamb_core(atindx,dtset,omlamb,pawtab,psps%usepaw)
1010 :
1011 : ! output raw data to netcdf file for more detailed postprocessing
1012 16 : if (me == master) then
1013 16 : fname = trim(dtfil%filnam_ds(4))//'_ORBMAG.nc'
1014 16 : NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
1015 16 : call orbmag_ncwrite(crystal,dtset,ebands_k,hdr,ncid,orbmag_mesh)
1016 16 : NCF_CHECK(nf90_close(ncid))
1017 : end if
1018 :
1019 : ! output summary to abo file
1020 16 : call orbmag_mesh%output(dtset,omlamb)
1021 :
1022 : !---------------------------------------------------
1023 : ! deallocate memory
1024 : !---------------------------------------------------
1025 :
1026 16 : call gs_hamk%free()
1027 :
1028 16 : ABI_SFREE(kg_k)
1029 16 : ABI_SFREE(kinpw)
1030 16 : ABI_SFREE(dkinpw)
1031 16 : ABI_SFREE(ph1d)
1032 :
1033 16 : ABI_SFREE(atindx)
1034 16 : ABI_SFREE(atindx1)
1035 16 : ABI_SFREE(nattyp)
1036 :
1037 16 : ABI_FREE(dimlmn)
1038 16 : call pawcprj_free(cwaveprj)
1039 52 : ABI_FREE(cwaveprj)
1040 :
1041 16 : call dterm%free()
1042 16 : call orbmag_mesh%free()
1043 :
1044 32 : end subroutine orbmag
1045 : !!***
1046 :
1047 : !!****f*m_orbmag/orbmag_mpisum
1048 : !! NAME
1049 : !! orbmag_mpisum
1050 : !!
1051 : !! FUNCTION
1052 : !! accumulate data in orbmag_mesh_type over processes
1053 : !!
1054 : !! INPUTS
1055 : !!
1056 : !! OUTPUT
1057 : !!
1058 : !! SOURCE
1059 :
1060 17 : subroutine orbmag_mpisum(self,nproc,spaceComm)
1061 : !Arguments
1062 : class(orbmag_mesh_type),intent(inout),target :: self
1063 : integer,intent(in) :: nproc,spaceComm
1064 :
1065 : !Local variables
1066 : integer :: buff_size,ierr
1067 17 : real(dp),allocatable :: buffer1(:),buffer2(:)
1068 :
1069 17 : if (nproc > 1) then
1070 0 : if (allocated(self%cmesh)) then
1071 0 : buff_size=size(self%cmesh)
1072 0 : ABI_MALLOC(buffer1,(buff_size))
1073 0 : ABI_MALLOC(buffer2,(buff_size))
1074 0 : buffer1=zero;buffer2=zero
1075 : buffer1(1:buff_size) = &
1076 0 : & reshape(self%cmesh,(/self%mband*self%nkpt*self%nsppol*3*chern_nterms/))
1077 0 : call xmpi_sum(buffer1,buffer2,buff_size,spaceComm,ierr)
1078 : self%cmesh(1:self%mband,1:self%nkpt,1:self%nsppol,1:3,1:chern_nterms)=&
1079 0 : & reshape(buffer2,(/self%mband,self%nkpt,self%nsppol,3,chern_nterms/))
1080 0 : ABI_FREE(buffer1)
1081 0 : ABI_FREE(buffer2)
1082 : end if
1083 0 : if (allocated(self%chern_terms)) then
1084 0 : buff_size=size(self%chern_terms)
1085 0 : ABI_MALLOC(buffer1,(buff_size))
1086 0 : ABI_MALLOC(buffer2,(buff_size))
1087 0 : buffer1=zero;buffer2=zero
1088 : buffer1(1:buff_size) = &
1089 0 : & reshape(self%chern_terms,(/self%mband*self%nsppol*3*chern_nterms/))
1090 0 : call xmpi_sum(buffer1,buffer2,buff_size,spaceComm,ierr)
1091 : self%chern_terms(1:self%mband,1:self%nsppol,1:3,1:chern_nterms)=&
1092 0 : & reshape(buffer2,(/self%mband,self%nsppol,3,chern_nterms/))
1093 0 : ABI_FREE(buffer1)
1094 0 : ABI_FREE(buffer2)
1095 : end if
1096 0 : if (allocated(self%omesh)) then
1097 0 : buff_size=size(self%omesh)
1098 0 : ABI_MALLOC(buffer1,(buff_size))
1099 0 : ABI_MALLOC(buffer2,(buff_size))
1100 0 : buffer1=zero;buffer2=zero
1101 : buffer1(1:buff_size) = &
1102 0 : & reshape(self%omesh,(/self%mband*self%nkpt*self%nsppol*3*orbmag_nterms/))
1103 0 : call xmpi_sum(buffer1,buffer2,buff_size,spaceComm,ierr)
1104 : self%omesh(1:self%mband,1:self%nkpt,1:self%nsppol,1:3,1:orbmag_nterms)=&
1105 0 : & reshape(buffer2,(/self%mband,self%nkpt,self%nsppol,3,orbmag_nterms/))
1106 0 : ABI_FREE(buffer1)
1107 0 : ABI_FREE(buffer2)
1108 : end if
1109 0 : if (allocated(self%orbmag_terms)) then
1110 0 : buff_size=size(self%orbmag_terms)
1111 0 : ABI_MALLOC(buffer1,(buff_size))
1112 0 : ABI_MALLOC(buffer2,(buff_size))
1113 0 : buffer1=zero;buffer2=zero
1114 : buffer1(1:buff_size) = &
1115 0 : & reshape(self%orbmag_terms,(/self%mband*self%nsppol*3*orbmag_nterms/))
1116 0 : call xmpi_sum(buffer1,buffer2,buff_size,spaceComm,ierr)
1117 : self%orbmag_terms(1:self%mband,1:self%nsppol,1:3,1:orbmag_nterms)=&
1118 0 : & reshape(buffer2,(/self%mband,self%nsppol,3,orbmag_nterms/))
1119 0 : ABI_FREE(buffer1)
1120 0 : ABI_FREE(buffer2)
1121 : end if
1122 0 : if (allocated(self%rmesh)) then
1123 0 : buff_size=size(self%rmesh)
1124 0 : ABI_MALLOC(buffer1,(buff_size))
1125 0 : ABI_MALLOC(buffer2,(buff_size))
1126 0 : buffer1=zero;buffer2=zero
1127 : buffer1(1:buff_size) = &
1128 0 : & reshape(self%rmesh,(/self%n4*self%n5*self%n6*3*orbmag_nterms/))
1129 0 : call xmpi_sum(buffer1,buffer2,buff_size,spaceComm,ierr)
1130 : self%rmesh(1:self%n4,1:self%n5,1:self%n6,1:3,1:orbmag_nterms)=&
1131 0 : & reshape(buffer2,(/self%n4,self%n5,self%n6,3,orbmag_nterms/))
1132 0 : ABI_FREE(buffer1)
1133 0 : ABI_FREE(buffer2)
1134 : end if
1135 :
1136 : end if
1137 :
1138 17 : end subroutine orbmag_mpisum
1139 : !!***
1140 :
1141 :
1142 : !!****f*m_orbmag/orbmag_term_scale
1143 : !! NAME
1144 : !! orbmag_term_scale
1145 : !!
1146 : !! FUNCTION
1147 : !! change frames and scale terms as needed
1148 : !!
1149 : !! INPUTS
1150 : !!
1151 : !! OUTPUT
1152 : !!
1153 : !! SOURCE
1154 :
1155 17 : subroutine orbmag_term_scale(self,crystal,dtset)
1156 :
1157 : !Arguments ------------------------------------
1158 : !scalars
1159 : class(orbmag_mesh_type),intent(inout),target :: self
1160 : type(crystal_t),intent(in) :: crystal
1161 : type(dataset_type),intent(in) :: dtset
1162 :
1163 : !arrays
1164 :
1165 : !Local variables -------------------------
1166 : !scalars
1167 : integer :: i4,i5,i6,isppol,iterm,nn
1168 : !arrays
1169 :
1170 : !--------------------------------------------------------------------
1171 :
1172 119 : do iterm = 1, orbmag_nterms
1173 227 : do isppol = 1, dtset%nsppol
1174 762 : do nn = 1, self%mband
1175 660 : if((iterm.EQ.inlr).OR.(iterm.EQ.inbm)) then
1176 : self%orbmag_terms(nn,isppol,1:3,iterm) = &
1177 3864 : & MATMUL(crystal%rprimd,self%orbmag_terms(nn,isppol,1:3,iterm))
1178 : else
1179 : self%orbmag_terms(nn,isppol,1:3,iterm) = &
1180 5888 : & crystal%ucvol*MATMUL(crystal%gprimd,self%orbmag_terms(nn,isppol,1:3,iterm))
1181 : end if
1182 : end do ! nn
1183 : end do !isppol
1184 : end do
1185 :
1186 17 : if (dtset%orbmag .EQ. 4) then
1187 7 : do iterm = 1, orbmag_nterms
1188 277 : do i4=1,self%n4
1189 12426 : do i5=1,self%n5
1190 559170 : do i6=1,self%n6
1191 558900 : if((iterm.EQ.inlr).OR.(iterm.EQ.inbm)) then
1192 : self%rmesh(i4,i5,i6,1:3,iterm) = &
1193 3827250 : & MATMUL(crystal%rprimd,self%rmesh(i4,i5,i6,1:3,iterm))
1194 : else
1195 : self%rmesh(i4,i5,i6,1:3,iterm) = &
1196 5832000 : & crystal%ucvol*MATMUL(crystal%gprimd,self%rmesh(i4,i5,i6,1:3,iterm))
1197 : end if
1198 : end do
1199 : end do
1200 : end do
1201 : end do
1202 : end if
1203 :
1204 17 : if (dtset%orbmag .EQ. 4) then
1205 7 : do iterm = 1, orbmag_nterms
1206 24 : self%rmesh_int(1:3,iterm) = zero
1207 277 : do i4=1,self%n4
1208 12426 : do i5=1,self%n5
1209 559170 : do i6=1,self%n6
1210 : self%rmesh_int(1:3,iterm) = self%rmesh_int(1:3,iterm) + &
1211 2199150 : & self%rmesh(i4,i5,i6,1:3,iterm)
1212 : end do
1213 : end do
1214 : end do
1215 : end do
1216 : end if
1217 :
1218 51 : do iterm = 1, chern_nterms
1219 87 : do isppol = 1, dtset%nsppol
1220 254 : do nn = 1, self%mband
1221 : self%chern_terms(nn,isppol,1:3,iterm) = &
1222 2980 : & crystal%ucvol*MATMUL(crystal%gprimd,self%chern_terms(nn,isppol,1:3,iterm))
1223 : end do ! nn
1224 : end do !isppol
1225 : end do
1226 :
1227 : !! convert orbmag magnetization to orbital moment
1228 : !! Berry curvature terms are ignored
1229 2405 : self%orbmag_terms(:,:,:,incc:inbm)=crystal%ucvol*self%orbmag_terms(:,:,:,incc:inbm)
1230 :
1231 : !! accumulate trace of terms
1232 35 : do isppol = 1, dtset%nsppol
1233 127 : do nn = 1, self%mband
1234 : self%orbmag_trace(1:3,1:orbmag_nterms) = self%orbmag_trace(1:3,1:orbmag_nterms) + &
1235 2300 : & self%orbmag_terms(nn,isppol,1:3,1:orbmag_nterms)
1236 : self%chern_trace(1:3,1:chern_nterms) = self%chern_trace(1:3,1:chern_nterms) + &
1237 846 : & self%chern_terms(nn,isppol,1:3,1:chern_nterms)
1238 : end do ! nn
1239 : end do ! isppol
1240 :
1241 17 : end subroutine orbmag_term_scale
1242 : !!***
1243 :
1244 : !!****f* ABINIT/orbmag_nl1_k
1245 : !! NAME
1246 : !! orbmag_nl1_k
1247 : !!
1248 : !! FUNCTION
1249 : !! make NL(1) term at k
1250 : !!
1251 : !! INPUTS
1252 : !! atindx(natom)=index table for atoms (see gstate.f)
1253 : !! cg_k(2,mcgk) ground state wavefunctions at this k point
1254 : !! cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k
1255 : !! dimlmn(dtset%natom)=cprj lmn dimensions
1256 : !! dterm <type(dterm_type)> data related to onsite interactions
1257 : !! dtset <type(dataset_type)>=all input variables for this dataset
1258 : !! gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k
1259 : !! ikpt=current k pt
1260 : !! isppol=current spin polarization
1261 : !! mcgk=2nd dimension of cg_k
1262 : !! mcprjk=dimension of cprj_k
1263 : !! mkmem_rbz=kpts in memory
1264 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
1265 : !! nband_k=bands at this kpt
1266 : !! oterm=chooses which onsite term to apply
1267 : !! npw_k=planewaves at this k point
1268 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
1269 : !!
1270 : !! OUTPUT
1271 : !!
1272 : !! SIDE EFFECTS
1273 : !! orbmag_mesh%omesh
1274 : !! if nl1_option = 1, orbmag contribution of <L_R> is returned
1275 : !! if nl1_option = 2, orbmag contribution of <A0.An> is returned
1276 : !!
1277 : !! TODO
1278 : !!
1279 : !! NOTES
1280 : !! returns \sum_{Rij}<u|p_i>a_ij<p_j|u> for various a_ij inputs
1281 : !!
1282 : !! SOURCE
1283 :
1284 1350 : subroutine orbmag_nl1_k(atindx,cg_k,cprj_k,dimlmn,dterm,dtset,gs_hamk,ikpt,&
1285 1350 : & oterm,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,trnrm)
1286 :
1287 : !Arguments ------------------------------------
1288 : !scalars
1289 : integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k,oterm
1290 : type(dterm_type),intent(in) :: dterm
1291 : type(dataset_type),intent(in) :: dtset
1292 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1293 : type(MPI_type), intent(inout) :: mpi_enreg
1294 : type(orbmag_mesh_type),intent(inout) :: orbmag_mesh
1295 :
1296 : !arrays
1297 : integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom)
1298 : real(dp),intent(in) :: trnrm(nband_k)
1299 : real(dp),intent(in),target :: cg_k(2,mcgk)
1300 : type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk)
1301 :
1302 : !Local variables -------------------------
1303 : !scalars
1304 : integer :: adir,choice,cpopt,dimekb1,dimekb2,dimekb3
1305 : integer :: n4,n5,n6,ndat,nn,nnlout,npwsp
1306 : integer :: paw_opt,signs,tim_nonlop
1307 : complex(dp) :: prefac_m,ormesh_fac
1308 : logical :: need_ormesh
1309 1350 : type(gs_hamiltonian_type),target :: gs_hamk_local
1310 : !arrays
1311 : real(dp) :: enlout(1),lambda(1),udotu(2)
1312 1350 : real(dp),allocatable :: fofr(:,:,:,:),svectout(:,:)
1313 1350 : real(dp),allocatable,target :: vectout(:,:)
1314 1350 : real(dp),pointer :: cwavef(:,:)
1315 1350 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
1316 : !--------------------------------------------------------------------
1317 :
1318 1350 : npwsp = npw_k*dtset%nspinor
1319 1350 : prefac_m = cone
1320 4050 : ABI_MALLOC(vectout,(2,npwsp))
1321 2700 : ABI_MALLOC(svectout,(2,npwsp))
1322 :
1323 : ! make a deep copy of gs_hamk so we can sub in the LR and BM
1324 : ! Dij values without destroying gs_hamk
1325 1350 : call gs_hamk%copy(gs_hamk_local)
1326 1350 : dimekb1=size(gs_hamk_local%ekb_spin,1)
1327 1350 : dimekb2=size(gs_hamk_local%ekb_spin,2)
1328 1350 : dimekb3=size(gs_hamk_local%ekb_spin,3)
1329 :
1330 1350 : need_ormesh = (dtset%orbmag .EQ. 4)
1331 1350 : if (need_ormesh) then
1332 2 : n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6); ndat=1
1333 10 : ABI_MALLOC(fofr,(2,n4,n5,n6*ndat))
1334 : end if
1335 :
1336 9700 : ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor))
1337 1350 : call pawcprj_alloc(cwaveprj,cprj_k(1,1)%ncpgr,dimlmn)
1338 :
1339 9822 : do nn = 1, nband_k
1340 8472 : cwavef => cg_k(1:2,(nn-1)*npwsp+1:nn*npwsp)
1341 : call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,&
1342 8472 : & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0)
1343 :
1344 35238 : do adir = 1, 3
1345 :
1346 5848488 : gs_hamk_local%ekb_spin = zero
1347 38124 : select case (oterm)
1348 : case ( inlr )
1349 : gs_hamk_local%ekb_spin(1:dimekb1,1:dimekb2,1:dimekb3,1,1) = &
1350 2784396 : & dterm%ekb_LR(1:dimekb1,1:dimekb2,1:dimekb3,adir)
1351 : case ( inbm )
1352 : gs_hamk_local%ekb_spin(1:dimekb1,1:dimekb2,1:dimekb3,1,1) = &
1353 2784396 : & dterm%ekb_BM(1:dimekb1,1:dimekb2,1:dimekb3,adir)
1354 : case default
1355 25416 : gs_hamk_local%ekb_spin(1:dimekb1,1:dimekb2,1:dimekb3,1,1) = zero
1356 : end select
1357 :
1358 : ! use nonlop to construct vectout = \sum_ij |p_i>D_ij<p_j|unk>
1359 25416 : choice = 1; cpopt = 4; paw_opt = 1; signs = 2; nnlout = 1; ndat = 1
1360 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk_local,adir,lambda,mpi_enreg,ndat,nnlout,&
1361 25416 : & paw_opt,signs,svectout,tim_nonlop,cwavef,vectout)
1362 25416 : udotu=cg_zdotc(npwsp,cwavef,vectout)
1363 25416 : orbmag_mesh%omesh(nn,ikpt,isppol,adir,oterm) = udotu(1)
1364 :
1365 33888 : if (need_ormesh) then
1366 24 : ormesh_fac = prefac_m*trnrm(nn)
1367 24 : call me_proj_mesh(cwavef,fofr,gs_hamk,vectout,mpi_enreg,n4,n5,n6,ndat,npw_k,ormesh_fac)
1368 2236704 : orbmag_mesh%rmesh(:,:,:,adir,oterm)=orbmag_mesh%rmesh(:,:,:,adir,oterm)+fofr(1,:,:,:)
1369 : end if
1370 :
1371 : end do !adir
1372 :
1373 : end do !nn
1374 :
1375 1350 : call pawcprj_free(cwaveprj)
1376 4044 : ABI_SFREE(cwaveprj)
1377 1350 : IF(ASSOCIATED(cwavef)) NULLIFY(cwavef)
1378 1350 : ABI_SFREE(fofr)
1379 1350 : ABI_SFREE(vectout)
1380 1350 : ABI_SFREE(svectout)
1381 :
1382 1350 : call gs_hamk_local%free()
1383 :
1384 2700 : end subroutine orbmag_nl1_k
1385 : !!***
1386 :
1387 : !!****f* ABINIT/orbmag_nl_k
1388 : !! NAME
1389 : !! orbmag_nl_k
1390 : !!
1391 : !! FUNCTION
1392 : !! make NL term at k
1393 : !!
1394 : !! INPUTS
1395 : !! atindx(natom)=index table for atoms (see gstate.f)
1396 : !! cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k
1397 : !! dimlmn(dtset%natom)=cprj lmn dimensions
1398 : !! dterm <type(dterm_type)> data related to onsite interactions
1399 : !! dtset <type(dataset_type)>=all input variables for this dataset
1400 : !! eig_k(nband_k)=gs eigenvalues at this kpt
1401 : !! ikpt=current k pt
1402 : !! isppol=current spin polarization
1403 : !! mcprjk=dimension of cprj_k
1404 : !! mkmem_rbz=kpts in memory
1405 : !! nband_k=bands at this kpt
1406 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
1407 : !!
1408 : !! OUTPUT
1409 : !!
1410 : !! SIDE EFFECTS
1411 : !! orbmag_mesh%omesh
1412 : !!
1413 : !! TODO
1414 : !!
1415 : !! NOTES
1416 : !! computes -\frac{i}{2}\sum_{Rij}<u|d_b p_i>D^0_{ij} - E^0s^0_{ij}<d_g p_j|u>
1417 : !! This is ZTG23 Eq. 36 term 1
1418 : !!
1419 : !! SOURCE
1420 :
1421 675 : subroutine orbmag_nl_k(atindx,cg_k,cprj_k,dimlmn,dterm,dtset,eig_k,gs_hamk,ikpt,isppol,&
1422 675 : & mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,pawtab,trnrm)
1423 :
1424 : !Arguments ------------------------------------
1425 : !scalars
1426 : integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k
1427 : type(dterm_type),intent(in) :: dterm
1428 : type(dataset_type),intent(in) :: dtset
1429 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1430 : type(MPI_type), intent(inout) :: mpi_enreg
1431 : type(orbmag_mesh_type),intent(inout) :: orbmag_mesh
1432 :
1433 : !arrays
1434 : integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom)
1435 : real(dp),intent(in),target :: cg_k(2,mcgk)
1436 : real(dp),intent(in) :: eig_k(nband_k),trnrm(nband_k)
1437 : type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk)
1438 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
1439 :
1440 : !Local variables -------------------------
1441 : !scalars
1442 : integer :: adir,bdir,gdir,n4,n5,n6,ndat,nn,npwsp
1443 : complex(dp) :: prefac_m,ormesh_fac,txt
1444 : logical :: need_ormesh
1445 : !arrays
1446 : real(dp) :: nonlop_udotu(2)
1447 675 : real(dp),allocatable :: fofr(:,:,:,:),svectout(:,:)
1448 : real(dp),allocatable,target :: vectout(:,:)
1449 675 : real(dp),pointer :: unk(:,:)
1450 675 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
1451 :
1452 : !--------------------------------------------------------------------
1453 :
1454 675 : npwsp = npw_k*dtset%nspinor
1455 2025 : ABI_MALLOC(vectout,(2,npwsp))
1456 1350 : ABI_MALLOC(svectout,(2,npwsp))
1457 :
1458 675 : need_ormesh = (dtset%orbmag .EQ. 4)
1459 675 : if (need_ormesh) then
1460 1 : n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6); ndat=1
1461 5 : ABI_MALLOC(fofr,(2,n4,n5,n6*ndat))
1462 : end if
1463 :
1464 4850 : ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor))
1465 675 : call pawcprj_alloc(cwaveprj,cprj_k(1,1)%ncpgr,dimlmn)
1466 :
1467 4911 : do nn = 1, nband_k
1468 :
1469 4236 : unk => cg_k(1:2,(nn-1)*npwsp+1:nn*npwsp)
1470 : call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,&
1471 4236 : & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0)
1472 :
1473 4236 : prefac_m = -com*c2
1474 17619 : do adir = 1, 3
1475 : ! even permutation of adir,bdir,gdir
1476 : ! cross product term adir,gdir,bdir leads to (-i/2)(Z-Z*) and
1477 : ! double the original even term
1478 12708 : bdir=modulo(adir,3)+1
1479 12708 : gdir=modulo(bdir,3)+1
1480 :
1481 : call nonlop_orbmag_nl(atindx,cwaveprj,bdir,gdir,dterm,dtset,eig_k(nn),gs_hamk,npw_k,&
1482 12708 : & pawtab,vectout)
1483 :
1484 12708 : nonlop_udotu(1:2)=cg_zdotc(npwsp,unk,vectout)
1485 12708 : txt=prefac_m*CMPLX(nonlop_udotu(1),nonlop_udotu(2))
1486 :
1487 : orbmag_mesh%omesh(nn,ikpt,isppol,adir,innl) = &
1488 12708 : &orbmag_mesh%omesh(nn,ikpt,isppol,adir,innl) + two*REAL(txt)
1489 :
1490 16944 : if (need_ormesh) then
1491 12 : ormesh_fac = trnrm(nn)*prefac_m
1492 12 : call me_proj_mesh(unk,fofr,gs_hamk,vectout,mpi_enreg,n4,n5,n6,ndat,npw_k,ormesh_fac)
1493 : orbmag_mesh%rmesh(:,:,:,adir,innl)=&
1494 1118352 : &orbmag_mesh%rmesh(:,:,:,adir,innl)+two*fofr(1,:,:,:)
1495 : end if
1496 :
1497 : end do ! adir
1498 :
1499 : end do !nn
1500 :
1501 675 : call pawcprj_free(cwaveprj)
1502 2022 : ABI_SFREE(cwaveprj)
1503 675 : IF(ASSOCIATED(unk)) NULLIFY(unk)
1504 675 : ABI_SFREE(fofr)
1505 675 : ABI_SFREE(vectout)
1506 675 : ABI_SFREE(svectout)
1507 :
1508 1350 : end subroutine orbmag_nl_k
1509 : !!***
1510 :
1511 : !!****f* ABINIT/orbmag_cc_k_ncpp
1512 : !! NAME
1513 : !! orbmag_cc_k_ncpp
1514 : !!
1515 : !! FUNCTION
1516 : !! computes <P_c du/dk|H + E*S|P_c du/dk> term in orbital magnetism, for NCPP
1517 : !!
1518 : !! INPUTS
1519 : !! atindx(natom)=index table for atoms (see gstate.f)
1520 : !! cprj1_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for Pc cg1_k
1521 : !! dimlmn(dtset%natom)=cprj lmn dimensions
1522 : !! dtset <type(dataset_type)>=all input variables for this dataset
1523 : !! eig_k(nband_k)=gs eigenvalues at this kpt
1524 : !! fermie=offset energy to use
1525 : !! gcg1_k(2,mcgk,3)=gauge adjusted cg1_k
1526 : !! gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k
1527 : !! ikpt=current k pt
1528 : !! isppol=current spin polarization
1529 : !! mcgk=dimension of cg_k
1530 : !! mcprjk=dimension of cprj_k
1531 : !! mkmem_rbz=kpts in memory
1532 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
1533 : !! nband_k=bands at this kpt
1534 : !! npw_k=number of planewaves at this kpt
1535 : !! occ_k=band occupations at this kpt
1536 : !! ucvol=unit cell volume
1537 : !!
1538 : !! OUTPUT
1539 : !!
1540 : !! SIDE EFFECTS
1541 : !! orbmag_mesh%omesh, orbmag_mesh%cmesh
1542 : !!
1543 : !! TODO
1544 : !!
1545 : !! NOTES
1546 : !! ZTG23 Eq. 36 term 2 and Eq. 46 term 1
1547 : !!
1548 : !! SOURCE
1549 :
1550 32 : subroutine orbmag_cc_k_ncpp(dtset,eig_k,fermie,gcg1_k,gs_hamk,ikpt,isppol,&
1551 : & mcgk,mpi_enreg,nband_k,npw_k,orbmag_mesh)
1552 :
1553 : !Arguments ------------------------------------
1554 : !scalars
1555 : integer,intent(in) :: ikpt,isppol,mcgk,nband_k,npw_k
1556 : real(dp),intent(in) :: fermie
1557 : type(dataset_type),intent(in) :: dtset
1558 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1559 : type(MPI_type), intent(inout) :: mpi_enreg
1560 : type(orbmag_mesh_type),intent(inout) :: orbmag_mesh
1561 :
1562 : !arrays
1563 : real(dp),intent(in) :: eig_k(nband_k)
1564 : real(dp),intent(in),target :: gcg1_k(2,mcgk,3)
1565 :
1566 : !Local variables -------------------------
1567 : !scalars
1568 : integer :: adir,bdir,cpopt,gdir
1569 : integer :: ndat,nn,npwsp,sij_opt,tim_getghc,type_calc
1570 : real(dp) :: lams
1571 : complex(dp) :: prefac_b,prefac_m
1572 : !arrays
1573 : real(dp) bdot(2),mdot(2)
1574 32 : real(dp),allocatable :: gec(:,:),gvnlxc(:,:)
1575 32 : real(dp),allocatable,target :: ghc(:,:)
1576 32 : real(dp),pointer :: du_dbeta(:,:),du_dgamma(:,:)
1577 : complex(dp) :: m1(3),b1(3)
1578 32 : type(pawcprj_type),allocatable :: cwaveprj1(:,:)
1579 : !--------------------------------------------------------------------
1580 :
1581 32 : npwsp = npw_k*dtset%nspinor
1582 :
1583 96 : ABI_MALLOC(ghc,(2,npwsp))
1584 64 : ABI_MALLOC(gec,(2,npwsp))
1585 64 : ABI_MALLOC(gvnlxc,(2,npwsp))
1586 64 : ABI_MALLOC(cwaveprj1,(dtset%natom,0))
1587 :
1588 32 : tim_getghc = 0
1589 32 : lams = zero
1590 32 : ndat = 1
1591 :
1592 32 : prefac_b = cbc*c2
1593 32 : prefac_m = com*c2
1594 160 : do nn = 1, nband_k
1595 128 : cpopt = -1; sij_opt = 0
1596 128 : type_calc = 0 ! apply local and non-local Hamiltonian
1597 :
1598 : m1 = czero
1599 : b1 = czero
1600 :
1601 544 : do adir = 1, 3
1602 :
1603 : ! even permutation of adir,bdir,gdir
1604 : ! cross product term adir,gdir,bdir leads to (-i/2)(Z-Z*) and
1605 : ! double the original even term
1606 384 : bdir=modulo(adir,3)+1
1607 384 : gdir=modulo(bdir,3)+1
1608 :
1609 384 : du_dgamma => gcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,gdir)
1610 :
1611 : ! compute H|Pc d_\gamma u>
1612 : ! [H+(E - 2\mu)]|ket> is needed for orbmag
1613 : ! -2|ket> needed for chern
1614 : call getghc(cpopt,du_dgamma,cwaveprj1,ghc,gec,gs_hamk,gvnlxc,lams,mpi_enreg,&
1615 384 : & ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc)
1616 :
1617 168864 : gec(1:2,1:npwsp) = (eig_k(nn) - two*fermie)*du_dgamma(1:2,1:npwsp)
1618 168864 : ghc(1:2,1:npwsp) = ghc(1:2,1:npwsp) + gec(1:2,1:npwsp)
1619 :
1620 384 : du_dbeta => gcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,bdir)
1621 1536 : mdot = cg_zdotc(npwsp,du_dbeta,ghc); bdot = -two*cg_zdotc(npwsp,du_dbeta,du_dgamma)
1622 :
1623 : ! cross product term adir,gdir,bdir leads to (-i/2)(Z-Z*) and
1624 : ! double the original even term: this is the origin of the factor of two
1625 : orbmag_mesh%omesh(nn,ikpt,isppol,adir,incc) = &
1626 : & orbmag_mesh%omesh(nn,ikpt,isppol,adir,incc) + &
1627 384 : & two*real(prefac_m*CMPLX(mdot(1),mdot(2)))
1628 : orbmag_mesh%cmesh(nn,ikpt,isppol,adir,ibcc) = &
1629 : & orbmag_mesh%cmesh(nn,ikpt,isppol,adir,ibcc) - &
1630 512 : & four*real(prefac_b*CMPLX(bdot(1),bdot(2)))
1631 :
1632 : end do ! adir
1633 :
1634 : end do !nn
1635 :
1636 32 : if(ASSOCIATED(du_dbeta)) NULLIFY(du_dbeta)
1637 32 : if(ASSOCIATED(du_dgamma)) NULLIFY(du_dgamma)
1638 :
1639 32 : ABI_SFREE(ghc)
1640 32 : ABI_SFREE(gec)
1641 32 : ABI_SFREE(gvnlxc)
1642 32 : ABI_SFREE(cwaveprj1)
1643 :
1644 64 : end subroutine orbmag_cc_k_ncpp
1645 : !!***
1646 :
1647 :
1648 : !!****f* ABINIT/orbmag_cc_k
1649 : !! NAME
1650 : !! orbmag_cc_k
1651 : !!
1652 : !! FUNCTION
1653 : !! computes <P_c du/dk|H + E*S|P_c du/dk> term in orbital magnetism
1654 : !!
1655 : !! INPUTS
1656 : !! atindx(natom)=index table for atoms (see gstate.f)
1657 : !! cprj1_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for Pc cg1_k
1658 : !! dimlmn(dtset%natom)=cprj lmn dimensions
1659 : !! dtset <type(dataset_type)>=all input variables for this dataset
1660 : !! eig_k(nband_k)=gs eigenvalues at this kpt
1661 : !! fermie=offset energy to use
1662 : !! gcg1_k(2,mcgk,3)=gauge adjusted cg1_k
1663 : !! gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k
1664 : !! ikpt=current k pt
1665 : !! isppol=current spin polarization
1666 : !! mcgk=dimension of cg_k
1667 : !! mcprjk=dimension of cprj_k
1668 : !! mkmem_rbz=kpts in memory
1669 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
1670 : !! nband_k=bands at this kpt
1671 : !! npw_k=number of planewaves at this kpt
1672 : !! occ_k=band occupations at this kpt
1673 : !! ucvol=unit cell volume
1674 : !!
1675 : !! OUTPUT
1676 : !!
1677 : !! SIDE EFFECTS
1678 : !! orbmag_mesh%omesh, orbmag_mesh%cmesh
1679 : !!
1680 : !! TODO
1681 : !!
1682 : !! NOTES
1683 : !! ZTG23 Eq. 36 term 2 and Eq. 46 term 1
1684 : !!
1685 : !! SOURCE
1686 :
1687 675 : subroutine orbmag_cc_k(atindx,cprj1_k,dimlmn,dtset,eig_k,fermie,&
1688 675 : & gcg1_k,gs_hamk,ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,&
1689 675 : & npw_k,orbmag_mesh,trnrm)
1690 :
1691 : !Arguments ------------------------------------
1692 : !scalars
1693 : integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k
1694 : real(dp),intent(in) :: fermie
1695 : type(dataset_type),intent(in) :: dtset
1696 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1697 : type(MPI_type), intent(inout) :: mpi_enreg
1698 : type(orbmag_mesh_type),intent(inout) :: orbmag_mesh
1699 :
1700 : !arrays
1701 : integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom)
1702 : real(dp),intent(in) :: eig_k(nband_k),trnrm(nband_k)
1703 : real(dp),intent(in),target :: gcg1_k(2,mcgk,3)
1704 : type(pawcprj_type),intent(in) :: cprj1_k(dtset%natom,mcprjk,3)
1705 :
1706 : !Local variables -------------------------
1707 : !scalars
1708 : integer :: adir,bdir,cpopt,dum_dnlbra,dum_dnlket,gdir
1709 : integer :: ndat,n4,n5,n6,nn,npwsp,sij_opt,tim_getghc,type_calc
1710 : real(dp) :: lams
1711 : complex(dp) :: ormesh_fac,prefac_b,prefac_m
1712 : logical :: need_ormesh
1713 : !arrays
1714 : real(dp) bdot(2),mdot(2)
1715 675 : real(dp),allocatable :: fofr(:,:,:,:),gsc(:,:),gvnlxc(:,:)
1716 675 : real(dp),allocatable,target :: ghc(:,:)
1717 675 : real(dp),pointer :: du_dbeta(:,:),du_dgamma(:,:)
1718 : complex(dp) :: m1(3),b1(3)
1719 675 : type(pawcprj_type),allocatable :: cwaveprj1(:,:)
1720 : !--------------------------------------------------------------------
1721 :
1722 675 : npwsp = npw_k*dtset%nspinor
1723 675 : need_ormesh = (dtset%orbmag .EQ. 4)
1724 :
1725 2025 : ABI_MALLOC(ghc,(2,npwsp))
1726 1350 : ABI_MALLOC(gsc,(2,npwsp))
1727 1350 : ABI_MALLOC(gvnlxc,(2,npwsp))
1728 4850 : ABI_MALLOC(cwaveprj1,(dtset%natom,dtset%nspinor))
1729 675 : call pawcprj_alloc(cwaveprj1,0,dimlmn)
1730 :
1731 675 : tim_getghc = 0
1732 675 : lams = zero
1733 675 : ndat = 1
1734 :
1735 675 : if (need_ormesh) then
1736 1 : dum_dnlbra = 0; dum_dnlket = 0
1737 1 : n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6); ndat=1
1738 5 : ABI_MALLOC(fofr,(2,n4,n5,n6*ndat))
1739 : end if
1740 :
1741 675 : prefac_b = cbc*c2
1742 675 : prefac_m = com*c2
1743 4911 : do nn = 1, nband_k
1744 :
1745 : m1 = czero
1746 : b1 = czero
1747 :
1748 17619 : do adir = 1, 3
1749 :
1750 : ! even permutation of adir,bdir,gdir
1751 : ! cross product term adir,gdir,bdir leads to (-i/2)(Z-Z*) and
1752 : ! double the original even term
1753 12708 : bdir=modulo(adir,3)+1
1754 12708 : gdir=modulo(bdir,3)+1
1755 :
1756 12708 : cpopt = 2
1757 12708 : du_dgamma => gcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,gdir)
1758 :
1759 : call pawcprj_get(atindx,cwaveprj1,cprj1_k(:,:,gdir),dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,&
1760 12708 : & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0)
1761 :
1762 : ! compute H|Pc d_\gamma u> and S|Pc d_\gamma u>
1763 : ! [H+E*S - 2\mu*S]|ket> is needed for orbmag
1764 : ! -2*S|ket> needed for Chern
1765 12708 : type_calc = 0 ! apply local and non-local Hamiltonian
1766 12708 : sij_opt = 1 ! compute gsc in addition to ghc
1767 : call getghc(cpopt,du_dgamma,cwaveprj1,ghc,gsc,gs_hamk,gvnlxc,lams,mpi_enreg,&
1768 12708 : & ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc)
1769 :
1770 8003592 : ghc(1:2,1:npwsp) = ghc(1:2,1:npwsp) + gsc(1:2,1:npwsp)*(eig_k(nn) - two*fermie)
1771 :
1772 12708 : du_dbeta => gcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,bdir)
1773 12708 : mdot = cg_zdotc(npwsp,du_dbeta,ghc); bdot = cg_zdotc(npwsp,du_dbeta,gsc)
1774 :
1775 : ! cross product term adir,gdir,bdir leads to (-i/2)(Z-Z*) and
1776 : ! double the original even term: this is the origin of the factor of two
1777 : orbmag_mesh%omesh(nn,ikpt,isppol,adir,incc) = &
1778 : & orbmag_mesh%omesh(nn,ikpt,isppol,adir,incc) + &
1779 12708 : & two*real(prefac_m*CMPLX(mdot(1),mdot(2)))
1780 : orbmag_mesh%cmesh(nn,ikpt,isppol,adir,ibcc) = &
1781 : & orbmag_mesh%cmesh(nn,ikpt,isppol,adir,ibcc) - &
1782 12708 : & four*real(prefac_b*CMPLX(bdot(1),bdot(2)))
1783 :
1784 16944 : if (need_ormesh) then
1785 12 : ormesh_fac = trnrm(nn)*prefac_m
1786 12 : call me_proj_mesh(du_dbeta,fofr,gs_hamk,ghc,mpi_enreg,n4,n5,n6,ndat,npw_k,ormesh_fac)
1787 : orbmag_mesh%rmesh(:,:,:,adir,incc)=&
1788 1118352 : &orbmag_mesh%rmesh(:,:,:,adir,incc)+two*fofr(1,:,:,:)
1789 : end if
1790 :
1791 : end do ! adir
1792 :
1793 : end do !nn
1794 :
1795 675 : if(ASSOCIATED(du_dbeta)) NULLIFY(du_dbeta)
1796 675 : if(ASSOCIATED(du_dgamma)) NULLIFY(du_dgamma)
1797 :
1798 675 : ABI_SFREE(ghc)
1799 675 : ABI_SFREE(gsc)
1800 675 : ABI_SFREE(gvnlxc)
1801 675 : call pawcprj_free(cwaveprj1)
1802 2022 : ABI_SFREE(cwaveprj1)
1803 675 : ABI_SFREE(fofr)
1804 :
1805 1350 : end subroutine orbmag_cc_k
1806 : !!***
1807 :
1808 : !!****f* ABINIT/orbmag_vv_k
1809 : !! NAME
1810 : !! orbmag_vv_k
1811 : !!
1812 : !! FUNCTION
1813 : !! orbmag_vv_k
1814 : !!
1815 : !! INPUTS
1816 : !! atindx(natom)=index table for atoms (see gstate.f)
1817 : !! cg_k(2,mcgk)=ground state wavefunctions at this k point
1818 : !! cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k
1819 : !! dimlmn(dtset%natom)=cprj lmn dimensions
1820 : !! dtset <type(dataset_type)>=all input variables for this dataset
1821 : !! eig_k(nband_k)=gs eigenvalues at this kpt
1822 : !! fermie=offset energy to use
1823 : !! gcg1_k(2,mcgk,3)=gauge treated cg1_k
1824 : !! gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k
1825 : !! ikpt=current k pt
1826 : !! isppol=current spin polarization
1827 : !! mcgk=dimension of cg_k
1828 : !! mcprjk=dimension of cprj_k
1829 : !! mkmem_rbz=kpts in memory
1830 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
1831 : !! nband_k=bands at this kpt
1832 : !! npw_k=number of planewaves at this kpt
1833 : !!
1834 : !! OUTPUT
1835 : !!
1836 : !! SIDE EFFECTS
1837 : !! orbmag_mesh%omesh, orbmag_mesh%cmesh
1838 : !!
1839 : !! TODO
1840 : !!
1841 : !! NOTES
1842 : !! contributions (1) <Pc d_b u|E d_gS|u> + <u|E d_bS|Pc d_g u> and
1843 : !! (2) \sum_n' <u |d_b ES|u_n'><u_n'|d_g ES|u> to orbital magnetization
1844 : !! these are ZTG23 Eq 36 terms 3 and 4, and Eq. 46 term 2
1845 : !!
1846 : !! SOURCE
1847 :
1848 675 : subroutine orbmag_vv_k(atindx,cg_k,cprj_k,dimlmn,dtset,eig_k,fermie,gcg1_k,gs_hamk,&
1849 675 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,orbmag_mesh,trnrm)
1850 :
1851 : !Arguments ------------------------------------
1852 : !scalars
1853 : integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k
1854 : real(dp),intent(in) :: fermie
1855 : type(dataset_type),intent(in) :: dtset
1856 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
1857 : type(MPI_type), intent(inout) :: mpi_enreg
1858 : type(orbmag_mesh_type),intent(inout) :: orbmag_mesh
1859 :
1860 : !arrays
1861 : integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom)
1862 : real(dp),intent(in) :: eig_k(nband_k),trnrm(nband_k)
1863 : real(dp),intent(in),target :: cg_k(2,mcgk),gcg1_k(2,mcgk,3)
1864 : type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk)
1865 :
1866 : !Local variables -------------------------
1867 : !scalars
1868 : integer :: adir,bdir,choice,cpopt,dnl_dum,fourwf_cplex,fourwf_option,gdir
1869 : integer :: n4,n5,n6,ndat,nn,nnlout,np,npwsp,paw_opt,signs,tim_fourwf,tim_getghc
1870 : complex(dp) :: b1,bdotc,bpdotc,gdotc,gpdotc,m1,mv2b,ormesh_fac,prefac_b,prefac_m
1871 : logical :: need_ormesh
1872 : !arrays
1873 : real(dp) :: bdot(2),bpdot(2),gdot(2),gpdot(2),enlout(1),lamv(1)
1874 675 : real(dp),allocatable :: fofr(:,:,:,:),proj_un(:,:),vectout(:,:)
1875 675 : real(dp),allocatable,target :: svectoutb(:,:),svectoutbp(:,:),svectoutg(:,:)
1876 675 : real(dp),pointer :: bra(:,:),du_dbeta(:,:),du_dgamma(:,:),unk(:,:)
1877 675 : type(pawcprj_type),allocatable :: cwaveprj(:,:),vv2_cwaveprj(:,:)
1878 : !--------------------------------------------------------------------
1879 :
1880 675 : npwsp = npw_k*dtset%nspinor
1881 675 : fourwf_cplex = 1
1882 675 : fourwf_option = 0
1883 675 : tim_fourwf = 1
1884 675 : need_ormesh = (dtset%orbmag .EQ. 4)
1885 :
1886 2025 : ABI_MALLOC(vectout,(2,npwsp))
1887 1350 : ABI_MALLOC(svectoutb,(2,npwsp))
1888 1350 : ABI_MALLOC(svectoutg,(2,npwsp))
1889 4850 : ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor))
1890 4175 : ABI_MALLOC(vv2_cwaveprj,(dtset%natom,dtset%nspinor))
1891 675 : call pawcprj_alloc(cwaveprj,cprj_k(1,1)%ncpgr,dimlmn)
1892 675 : call pawcprj_alloc(vv2_cwaveprj,cprj_k(1,1)%ncpgr,dimlmn)
1893 :
1894 675 : tim_getghc = 0
1895 675 : lamv = zero
1896 675 : ndat = 1
1897 675 : nnlout = 1
1898 675 : cpopt = 4 ! cprj and derivs in memory
1899 675 : choice = 5 ! apply dS/dk
1900 675 : paw_opt = 3 ! retain dS/dk|u>
1901 675 : signs = 2
1902 :
1903 675 : if (need_ormesh) then
1904 1 : dnl_dum=0
1905 1 : n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6); ndat=1
1906 2 : ABI_MALLOC(proj_un,(2,npwsp))
1907 2 : ABI_MALLOC(svectoutbp,(2,npwsp))
1908 5 : ABI_MALLOC(fofr,(2,n4,n5,n6*ndat))
1909 : end if
1910 :
1911 675 : prefac_b = cbc*c2
1912 675 : prefac_m = com*c2
1913 4911 : do nn = 1, nband_k
1914 :
1915 : ! extract |u_nk>
1916 4236 : unk => cg_k(1:2,(nn-1)*npwsp+1:nn*npwsp)
1917 : call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,&
1918 4236 : & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0)
1919 :
1920 17619 : do adir=1,3
1921 : ! even permutation of adir,bdir,gdir
1922 : ! cross product term adir,gdir,bdir leads to (-i/2)(Z-Z*) and
1923 : ! double the original even term
1924 12708 : bdir=modulo(adir,3)+1
1925 12708 : gdir=modulo(bdir,3)+1
1926 :
1927 : ! compute dS/dk_b|u_nk>
1928 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,bdir,lamv,mpi_enreg,ndat,nnlout,&
1929 12708 : & paw_opt,signs,svectoutb,tim_getghc,unk,vectout)
1930 :
1931 : ! compute dS/dk_g |u_nk>
1932 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,gdir,lamv,mpi_enreg,ndat,nnlout,&
1933 12708 : & paw_opt,signs,svectoutg,tim_getghc,unk,vectout)
1934 :
1935 : ! extract |Pc du/dk_b>
1936 12708 : du_dbeta => gcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,bdir)
1937 12708 : gdot=cg_zdotc(npwsp,du_dbeta,svectoutg); gdotc=CMPLX(gdot(1),gdot(2))
1938 :
1939 : ! extract |Pc du/dk_g>
1940 12708 : du_dgamma => gcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,gdir)
1941 12708 : bdot=cg_zdotc(npwsp,du_dgamma,svectoutb); bdotc=CMPLX(bdot(1),bdot(2))
1942 :
1943 : ! here is <Pc du/dk_b|dS/dk_g|u_nk>*E_nk
1944 12708 : b1 = -prefac_b*gdotc
1945 12708 : m1 = prefac_m*gdotc*(eig_k(nn) - fermie)
1946 :
1947 : ! add CONJG(<Pc du/dk_b|dS/dk_g|u_nk>)*E_nk
1948 12708 : b1 = b1 - prefac_b*CONJG(bdotc)
1949 12708 : m1 = m1 + prefac_m*CONJG(bdotc)*(eig_k(nn) - fermie)
1950 :
1951 : orbmag_mesh%omesh(nn,ikpt,isppol,adir,invv1) = &
1952 12708 : & orbmag_mesh%omesh(nn,ikpt,isppol,adir,invv1) + two*real(m1)
1953 :
1954 : orbmag_mesh%cmesh(nn,ikpt,isppol,adir,ibvv1) = &
1955 12708 : & orbmag_mesh%cmesh(nn,ikpt,isppol,adir,ibvv1) + two*real(b1)
1956 :
1957 12708 : if (need_ormesh) then
1958 :
1959 12 : ormesh_fac = trnrm(nn)*prefac_m*(eig_k(nn) - fermie)
1960 :
1961 12 : call me_proj_mesh(du_dbeta,fofr,gs_hamk,svectoutg,mpi_enreg,n4,n5,n6,ndat,npw_k,ormesh_fac)
1962 : orbmag_mesh%rmesh(:,:,:,adir,invv1)=&
1963 1118352 : &orbmag_mesh%rmesh(:,:,:,adir,invv1)+two*fofr(1,:,:,:)
1964 :
1965 12 : call me_proj_mesh(svectoutb,fofr,gs_hamk,du_dgamma,mpi_enreg,n4,n5,n6,ndat,npw_k,ormesh_fac)
1966 : orbmag_mesh%rmesh(:,:,:,adir,invv1)=&
1967 1118352 : &orbmag_mesh%rmesh(:,:,:,adir,invv1)+two*fofr(1,:,:,:)
1968 :
1969 : end if
1970 :
1971 12708 : mv2b = czero
1972 50202 : do np = nn+1, nband_k
1973 37494 : bra => cg_k(1:2,(np-1)*npwsp+1:np*npwsp)
1974 37494 : gpdot=cg_zdotc(npwsp,bra,svectoutg); gpdotc=CMPLX(gpdot(1),gpdot(2))
1975 37494 : bpdot=cg_zdotc(npwsp,bra,svectoutb); bpdotc=CMPLX(bpdot(1),bpdot(2))
1976 50202 : mv2b = mv2b + prefac_m*CONJG(bpdotc)*gpdotc*(eig_k(nn) - eig_k(np))
1977 : end do ! np
1978 :
1979 : ! Note that term VV2 has a minus sign in contrast to VV1 (see ZTG Eq. 36)
1980 : orbmag_mesh%omesh(nn,ikpt,isppol,adir,invv2) = &
1981 12708 : & orbmag_mesh%omesh(nn,ikpt,isppol,adir,invv2) - two*real(mv2b)
1982 :
1983 16944 : if (need_ormesh) then
1984 :
1985 : ! project dS/dk_g|u> onto valence space
1986 181488 : proj_un = zero
1987 60 : do np=1,nband_k
1988 48 : bra => cg_k(1:2,(np-1)*npwsp+1:np*npwsp)
1989 48 : gpdot=cg_zdotc(npwsp,bra,svectoutg)
1990 : proj_un(1,1:npwsp) = proj_un(1,1:npwsp) + &
1991 242016 : & bra(1,1:npwsp)*gpdot(1) - bra(2,1:npwsp)*gpdot(2)
1992 : proj_un(2,1:npwsp) = proj_un(2,1:npwsp) + &
1993 242028 : & bra(1,1:npwsp)*gpdot(2) + bra(2,1:npwsp)*gpdot(1)
1994 : end do
1995 :
1996 : ! compute dS/dk_b \sum_' |u'><u'|dS/dk_g|u>
1997 12 : cpopt = -1 ! cprj and derivs computed and not saved
1998 : call nonlop(choice,cpopt,vv2_cwaveprj,enlout,gs_hamk,bdir,lamv,mpi_enreg,ndat,nnlout,&
1999 12 : & paw_opt,signs,svectoutbp,tim_getghc,proj_un,vectout)
2000 12 : cpopt = 4 ! change cpopt back to its usual value in this routine
2001 12 : ormesh_fac = trnrm(nn)*prefac_m*(eig_k(nn) - fermie)
2002 12 : call me_proj_mesh(unk,fofr,gs_hamk,svectoutbp,mpi_enreg,n4,n5,n6,ndat,npw_k,ormesh_fac)
2003 : orbmag_mesh%rmesh(:,:,:,adir,invv2)=&
2004 1118352 : &orbmag_mesh%rmesh(:,:,:,adir,invv2)-two*fofr(1,:,:,:)
2005 :
2006 : end if ! computation of vv2 for rmesh
2007 :
2008 : end do ! adir
2009 :
2010 : end do !nn
2011 :
2012 675 : IF(ASSOCIATED(du_dbeta)) NULLIFY(du_dbeta)
2013 675 : IF(ASSOCIATED(du_dgamma)) NULLIFY(du_dgamma)
2014 675 : IF(ASSOCIATED(unk)) NULLIFY(unk)
2015 675 : IF(ASSOCIATED(bra)) NULLIFY(bra)
2016 :
2017 675 : ABI_SFREE(vectout)
2018 675 : ABI_SFREE(svectoutb)
2019 675 : ABI_SFREE(svectoutg)
2020 675 : call pawcprj_free(cwaveprj)
2021 675 : call pawcprj_free(vv2_cwaveprj)
2022 2022 : ABI_SFREE(cwaveprj)
2023 2022 : ABI_SFREE(vv2_cwaveprj)
2024 675 : ABI_SFREE(fofr)
2025 675 : ABI_SFREE(proj_un)
2026 675 : ABI_SFREE(svectoutbp)
2027 :
2028 1350 : end subroutine orbmag_vv_k
2029 : !!***
2030 :
2031 : !!****f* ABINIT/para_to_diag
2032 : !! NAME
2033 : !! para_to_diag
2034 : !!
2035 : !! FUNCTION
2036 : !! convert cg1_k wavefunction from parallel to diagonal gauge
2037 : !!
2038 : !! INPUTS
2039 : !! atindx(natom)=index table for atoms (see gstate.f)
2040 : !! cg_k(2,mcgk)=ground state wavefunctions at this k point
2041 : !! cg1_k(2,mcgk,3)=DDK wavefunctions at this k point, all 3 directions
2042 : !! cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k
2043 : !! dimlmn(dtset%natom)=cprj lmn dimensions
2044 : !! dtset <type(dataset_type)>=all input variables for this dataset
2045 : !! gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k
2046 : !! ikpt=current k pt
2047 : !! isppol=current spin polarization
2048 : !! mcgk=dimension of cg_k
2049 : !! mcprjk=dimension of cprj_k
2050 : !! mkmem_rbz=kpts in memory
2051 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
2052 : !! nband_k=bands at this kpt
2053 : !! npw_k=number of planewaves at this kpt
2054 : !! occ_k=band occupations at this kpt
2055 : !!
2056 : !! OUTPUT
2057 : !! gcg1_k(2,mcgk,3)=cg1_k converted to requested gauge and/or projection
2058 : !!
2059 : !! NOTES
2060 : !!
2061 : !! SOURCE
2062 :
2063 0 : subroutine para_to_diag(atindx,cg_k,cg1_k,cprj_k,dimlmn,dkinpw,dtset,eig_k,gcg1_k,gs_hamk,&
2064 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,mpw,nband_k,ngfft4,ngfft5,ngfft6,npw_k,&
2065 0 : & nucdip_dirs,occ_k,vectornd_pac)
2066 :
2067 : !Arguments ------------------------------------
2068 : !scalars
2069 : integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpw,nband_k,ngfft4,ngfft5,ngfft6
2070 : integer,intent(in) :: npw_k,nucdip_dirs
2071 : type(dataset_type),intent(in) :: dtset
2072 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
2073 : type(MPI_type), intent(inout) :: mpi_enreg
2074 :
2075 : !arrays
2076 : integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom)
2077 : real(dp),intent(in) :: cg_k(2,mcgk),cg1_k(2,mcgk,3),eig_k(nband_k),dkinpw(mpw,3),occ_k(nband_k)
2078 : real(dp),intent(in) :: vectornd_pac(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,nucdip_dirs)
2079 : real(dp),intent(out) :: gcg1_k(2,mcgk,3)
2080 : type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk)
2081 :
2082 : !Local variables -------------------------
2083 : !scalars
2084 : integer :: adir,berryopt,cplex,iband,ipert,jband,ndat,npwsp
2085 : integer :: optlocal,optnl,opt_gvnlx1,sij_opt,tim_getgh1c,usevnl
2086 : real(dp) :: corrfac,deltae,pertr,perti,pertsize
2087 0 : type(rf_hamiltonian_type) :: rf_hamk
2088 : !arrays
2089 : real(dp) :: hij(2),lambda(1),sij(2)
2090 0 : real(dp),allocatable :: cwavef(:,:),dcg1(:,:),gh1c(:,:)
2091 0 : real(dp),allocatable :: grad_berry(:,:),gs1c(:,:),gvnlx1(:,:)
2092 0 : real(dp),allocatable :: vectornd_pac_idir(:,:,:,:)
2093 0 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
2094 :
2095 : !--------------------------------------------------------------------
2096 :
2097 0 : berryopt=0
2098 0 : ndat=1
2099 0 : optlocal=0
2100 0 : optnl=2
2101 0 : opt_gvnlx1=0
2102 0 : sij_opt=1
2103 0 : tim_getgh1c=0
2104 0 : usevnl = 0
2105 0 : lambda(1) = zero
2106 :
2107 0 : ipert=dtset%natom+1 ! DDK
2108 0 : cplex=1 ! real space 1-order functions on FFT grid are REAL
2109 0 : call rf_hamk%init(cplex,gs_hamk,ipert)
2110 :
2111 0 : npwsp = npw_k*dtset%nspinor
2112 :
2113 0 : ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor))
2114 0 : call pawcprj_alloc(cwaveprj,3,dimlmn)
2115 0 : ABI_MALLOC(cwavef,(2,npwsp))
2116 0 : ABI_MALLOC(gh1c,(2,gs_hamk%npw_kp*gs_hamk%nspinor*ndat))
2117 0 : ABI_MALLOC(gs1c,(2,gs_hamk%npw_kp*gs_hamk%nspinor*ndat))
2118 0 : ABI_MALLOC(gvnlx1,(2,gs_hamk%npw_kp*gs_hamk%nspinor*ndat))
2119 0 : ABI_MALLOC(dcg1,(2,npwsp))
2120 :
2121 0 : if (nucdip_dirs .EQ. 3) then
2122 0 : ABI_MALLOC(vectornd_pac_idir,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc))
2123 : end if
2124 :
2125 0 : gcg1_k = zero
2126 :
2127 0 : do adir = 1, 3
2128 :
2129 0 : call rf_hamk%load_k(dkinpw_k=dkinpw(:,adir))
2130 :
2131 0 : if (nucdip_dirs .EQ. 3) then
2132 0 : vectornd_pac_idir(:,:,:,:)=vectornd_pac(:,:,:,:,adir)
2133 0 : call rf_hamk%load_spin(isppol, vectornd=vectornd_pac_idir)
2134 : end if
2135 :
2136 0 : do iband = 1, nband_k
2137 :
2138 0 : cwavef(1:2,1:npwsp)=cg_k(1:2,(iband-1)*npwsp+1:iband*npwsp)
2139 :
2140 : call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,iband,0,ikpt,0,isppol,dtset%mband,&
2141 0 : & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0)
2142 :
2143 : call getgh1c(berryopt,cwavef,cwaveprj,gh1c,grad_berry,gs1c,gs_hamk,gvnlx1,adir,ipert,&
2144 : & lambda(1),mpi_enreg,ndat,optlocal,optnl,opt_gvnlx1,rf_hamk,sij_opt,&
2145 0 : & tim_getgh1c,usevnl)
2146 :
2147 0 : dcg1=zero
2148 0 : do jband = 1, nband_k
2149 0 : if (jband .EQ. iband) cycle
2150 0 : if(abs(occ_k(jband)).LT.tol8) cycle
2151 0 : deltae = eig_k(jband) - eig_k(iband)
2152 : ! deltae test seems to work best compared to deltapert test
2153 : !if (abs(deltae) .LT. dtset%ggtrcut) cycle
2154 0 : cwavef(1:2,1:npwsp)=cg_k(1:2,(jband-1)*npwsp+1:jband*npwsp)
2155 0 : hij=cg_zdotc(npwsp,cwavef,gh1c)
2156 0 : sij=cg_zdotc(npwsp,cwavef,gs1c)
2157 0 : select case (dtset%orbmag)
2158 : case ( 3 )
2159 0 : lambda(1) = half*(eig_k(jband)+eig_k(iband))
2160 0 : corrfac=-one
2161 : case ( -3 )
2162 0 : lambda(1) = eig_k(iband)
2163 0 : corrfac=-one
2164 : end select
2165 0 : pertr = (hij(1)-lambda(1)*sij(1))/deltae
2166 0 : perti = (hij(2)-lambda(1)*sij(2))/deltae
2167 0 : pertsize=sqrt(pertr*pertr+perti*perti)
2168 0 : if (pertsize .GT. dtset%ggtrcut) cycle
2169 0 : dcg1(1,:) = dcg1(1,:) + pertr*cwavef(1,:) - perti*cwavef(2,:)
2170 0 : dcg1(2,:) = dcg1(2,:) + pertr*cwavef(2,:) + perti*cwavef(1,:)
2171 : end do
2172 : gcg1_k(1:2,(iband-1)*npwsp+1:iband*npwsp,adir) =&
2173 0 : &cg1_k(1:2,(iband-1)*npwsp+1:iband*npwsp,adir)+corrfac*dcg1(1:2,1:npwsp)
2174 : end do
2175 : end do
2176 :
2177 0 : call rf_hamk%free()
2178 0 : if(allocated(vectornd_pac_idir)) then
2179 0 : ABI_FREE(vectornd_pac_idir)
2180 : end if
2181 0 : ABI_FREE(cwavef)
2182 0 : call pawcprj_free(cwaveprj)
2183 0 : ABI_FREE(cwaveprj)
2184 0 : ABI_FREE(gh1c)
2185 0 : ABI_FREE(gs1c)
2186 0 : ABI_FREE(gvnlx1)
2187 0 : ABI_FREE(dcg1)
2188 :
2189 0 : end subroutine para_to_diag
2190 : !!***
2191 :
2192 : !!****f* ABINIT/gauge_treatment
2193 : !! NAME
2194 : !! gauge_treatment
2195 : !!
2196 : !! FUNCTION
2197 : !! convert cg1_k wavefunction to requested gauge and/or projection
2198 : !!
2199 : !! INPUTS
2200 : !! atindx(natom)=index table for atoms (see gstate.f)
2201 : !! cg_k(2,mcgk)=ground state wavefunctions at this k point
2202 : !! cg1_k(2,mcgk,3)=DDK wavefunctions at this k point, all 3 directions
2203 : !! cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k
2204 : !! dimlmn(dtset%natom)=cprj lmn dimensions
2205 : !! dtset <type(dataset_type)>=all input variables for this dataset
2206 : !! gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k
2207 : !! ikpt=current k pt
2208 : !! isppol=current spin polarization
2209 : !! mcgk=dimension of cg_k
2210 : !! mcprjk=dimension of cprj_k
2211 : !! mkmem_rbz=kpts in memory
2212 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
2213 : !! nband_k=bands at this kpt
2214 : !! npw_k=number of planewaves at this kpt
2215 : !! occ_k=band occupations at this kpt
2216 : !!
2217 : !! OUTPUT
2218 : !! gcg1_k(2,mcgk,3)=cg1_k converted to requested gauge and/or projection
2219 : !!
2220 : !! NOTES
2221 : !!
2222 : !! SOURCE
2223 :
2224 675 : subroutine gauge_treatment(atindx,cg_k,cg1_k,cprj_k,dimlmn,dkinpw,dtset,eig_k,gcg1_k,gs_hamk,&
2225 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,mpw,nband_k,ngfft4,ngfft5,ngfft6,npw_k,&
2226 675 : & nucdip_dirs,occ_k,vectornd_pac)
2227 :
2228 : !Arguments ------------------------------------
2229 : !scalars
2230 : integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpw,nband_k,ngfft4,ngfft5,ngfft6
2231 : integer,intent(in) :: npw_k,nucdip_dirs
2232 : type(dataset_type),intent(in) :: dtset
2233 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
2234 : type(MPI_type), intent(inout) :: mpi_enreg
2235 :
2236 : !arrays
2237 : integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom)
2238 : real(dp),intent(in) :: cg_k(2,mcgk),cg1_k(2,mcgk,3),eig_k(nband_k),dkinpw(mpw,3),occ_k(nband_k)
2239 : real(dp),intent(in) :: vectornd_pac(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,nucdip_dirs)
2240 : real(dp),intent(out) :: gcg1_k(2,mcgk,3)
2241 : type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk)
2242 :
2243 : !Local variables -------------------------
2244 : !scalars
2245 :
2246 : !--------------------------------------------------------------------
2247 :
2248 : ! orbmag > 0: cg1_k contains PAW DDK in parallel gauge, which has a ground state part
2249 : ! orbmag < 0: cg1_k contains Berry phase DDK, which is projected onto the
2250 : ! conduction space by construction
2251 :
2252 675 : select case (dtset%orbmag)
2253 : case ( -3 )
2254 : ! Convert Berry DDK to diagonal gauge
2255 : call para_to_diag(atindx,cg_k,cg1_k,cprj_k,dimlmn,dkinpw,dtset,eig_k,gcg1_k,gs_hamk,&
2256 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,mpw,nband_k,ngfft4,ngfft5,ngfft6,npw_k,&
2257 0 : & nucdip_dirs,occ_k,vectornd_pac)
2258 : case ( -2:-1 )
2259 : ! Berry DDK already projected onto conduction space, leave in berry gauge
2260 168608 : gcg1_k(1:2,1:mcgk,1:3) = cg1_k(1:2,1:mcgk,1:3)
2261 : case ( 3 )
2262 : ! Convert PAW DDK to diagonal gauge
2263 : call para_to_diag(atindx,cg_k,cg1_k,cprj_k,dimlmn,dkinpw,dtset,eig_k,gcg1_k,gs_hamk,&
2264 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,mpw,nband_k,ngfft4,ngfft5,ngfft6,npw_k,&
2265 0 : & nucdip_dirs,occ_k,vectornd_pac)
2266 : case default
2267 : ! project cg1_k onto conduction space by removing ground PAW part;
2268 : ! stay in parallel transport gauge
2269 : call make_pcg1(atindx,cg_k,cg1_k,cprj_k,dimlmn,dtset,gs_hamk,&
2270 675 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,occ_k,gcg1_k)
2271 : end select
2272 :
2273 675 : end subroutine gauge_treatment
2274 : !!***
2275 :
2276 : !!****f* ABINIT/make_pcg1
2277 : !! NAME
2278 : !! make_pcg1
2279 : !!
2280 : !! FUNCTION
2281 : !! compute Pc|cg1> from |cg1> and |cg>
2282 : !!
2283 : !! INPUTS
2284 : !! atindx(natom)=index table for atoms (see gstate.f)
2285 : !! cg_k(2,mcgk)=ground state wavefunctions at this k point
2286 : !! cg1_k(2,mcgk,3)=DDK wavefunctions at this k point, all 3 directions
2287 : !! cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k
2288 : !! dimlmn(dtset%natom)=cprj lmn dimensions
2289 : !! dtset <type(dataset_type)>=all input variables for this dataset
2290 : !! gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k
2291 : !! ikpt=current k pt
2292 : !! isppol=current spin polarization
2293 : !! mcgk=dimension of cg_k
2294 : !! mcprjk=dimension of cprj_k
2295 : !! mkmem_rbz=kpts in memory
2296 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
2297 : !! nband_k=bands at this kpt
2298 : !! npw_k=number of planewaves at this kpt
2299 : !! occ_k=band occupations at this kpt
2300 : !!
2301 : !! OUTPUT
2302 : !! pcg1_k(2,mcgk,3)=cg1_k projected on conduction space
2303 : !!
2304 : !! NOTES
2305 : !! see Audouze et al PRB 78, 035105 (2008) Eq. 40
2306 : !!
2307 : !! SOURCE
2308 :
2309 643 : subroutine make_pcg1(atindx,cg_k,cg1_k,cprj_k,dimlmn,dtset,gs_hamk,&
2310 643 : & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,occ_k,pcg1_k)
2311 :
2312 : !Arguments ------------------------------------
2313 : !scalars
2314 : integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k
2315 : type(dataset_type),intent(in) :: dtset
2316 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
2317 : type(MPI_type), intent(inout) :: mpi_enreg
2318 :
2319 : !arrays
2320 : integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom)
2321 : real(dp),intent(in) :: cg_k(2,mcgk),cg1_k(2,mcgk,3),occ_k(nband_k)
2322 : real(dp),intent(out) :: pcg1_k(2,mcgk,3)
2323 : type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk)
2324 :
2325 : !Local variables -------------------------
2326 : !scalars
2327 : integer :: adir,choice,cpopt,iband,jband
2328 : integer :: ndat,nnlout,npwsp,paw_opt,signs,tim_nonlop
2329 : !arrays
2330 : real(dp) :: dotp(2),lambda(1)
2331 643 : real(dp),allocatable :: cwavef(:,:),enlout(:),svectout(:,:)
2332 643 : real(dp),allocatable :: vcg1(:,:),vectout(:,:)
2333 643 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
2334 :
2335 : !--------------------------------------------------------------------
2336 :
2337 643 : choice = 5 ! dS/dk in nonlop
2338 643 : cpopt = 4 ! cprj and derivatives already in memory
2339 643 : paw_opt = 3
2340 643 : signs = 2
2341 643 : tim_nonlop = 0
2342 643 : lambda = zero
2343 643 : nnlout = 0
2344 643 : ndat = 1
2345 :
2346 643 : npwsp = npw_k*dtset%nspinor
2347 :
2348 4626 : ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor))
2349 643 : call pawcprj_alloc(cwaveprj,3,dimlmn)
2350 1929 : ABI_MALLOC(cwavef,(2,npwsp))
2351 1286 : ABI_MALLOC(vectout,(2,npwsp))
2352 1286 : ABI_MALLOC(svectout,(2,npwsp))
2353 1286 : ABI_MALLOC(vcg1,(2,npwsp))
2354 :
2355 7824976 : pcg1_k = zero
2356 :
2357 2572 : do adir = 1, 3
2358 :
2359 14896 : do iband = 1, nband_k
2360 :
2361 7834728 : cwavef(1:2,1:npwsp)=cg_k(1:2,(iband-1)*npwsp+1:iband*npwsp)
2362 :
2363 : call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,iband,0,ikpt,0,isppol,dtset%mband,&
2364 12324 : & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0)
2365 :
2366 : ! compute S^1|u_i^0> where S^1 = \partial S/\partial k_adir, the k derivative of S in
2367 : ! direction adir
2368 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,adir,lambda,mpi_enreg,ndat,&
2369 12324 : & nnlout,paw_opt,signs,svectout,tim_nonlop,cwavef,vectout)
2370 :
2371 : !! form vcg1 = -1/2 \sum |u_j^0><u_j^0|S^1|u_i^0>, the valence band part of cg1
2372 7834728 : vcg1 = zero
2373 98484 : do jband = 1, nband_k
2374 86160 : if(abs(occ_k(jband)).LT.tol8) cycle
2375 44085888 : cwavef(1:2,1:npwsp)=cg_k(1:2,(jband-1)*npwsp+1:jband*npwsp)
2376 66096 : dotp=cg_zdotc(npwsp,cwavef,svectout)
2377 14739360 : vcg1(1,:) = vcg1(1,:) - half*( dotp(1)*cwavef(1,:) - dotp(2)*cwavef(2,:))
2378 14771748 : vcg1(2,:) = vcg1(2,:) - half*( dotp(1)*cwavef(2,:) + dotp(2)*cwavef(1,:))
2379 : end do
2380 :
2381 : ! subtract vcg1 from cg1_k to obtain pcg1, the conduction band part of cg1
2382 : pcg1_k(1:2,(iband-1)*npwsp+1:iband*npwsp,adir) =cg1_k(1:2,(iband-1)*npwsp+1:iband*npwsp,adir)-&
2383 7836657 : & vcg1(1:2,1:npwsp)
2384 :
2385 : end do
2386 : end do
2387 :
2388 643 : ABI_FREE(cwavef)
2389 643 : ABI_FREE(vectout)
2390 643 : ABI_FREE(vcg1)
2391 643 : ABI_FREE(svectout)
2392 643 : call pawcprj_free(cwaveprj)
2393 1926 : ABI_FREE(cwaveprj)
2394 :
2395 643 : end subroutine make_pcg1
2396 : !!***
2397 :
2398 : !!****f* ABINIT/lamb_core
2399 : !! NAME
2400 : !! lamb_core
2401 : !!
2402 : !! FUNCTION
2403 : !! add core electron contribution to the orbital magnetic moment
2404 : !!
2405 : !! INPUTS
2406 : !! atindx(dtset%natom)=index table for atoms (see gstate.f)
2407 : !! dtset <type(dataset_type)>=all input variables for this dataset
2408 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
2409 : !!
2410 : !! OUTPUT
2411 : !! omlamb(2,3)=contribution of Lamb shielding to magnetic moment
2412 : !!
2413 : !! NOTES
2414 : !! lamb shielding of core electrons contributes -m.lambsig to orbital magnetic
2415 : !! moment
2416 : !!
2417 : !! SOURCE
2418 :
2419 17 : subroutine lamb_core(atindx,dtset,omlamb,pawtab,usepaw)
2420 :
2421 : !Arguments ------------------------------------
2422 : !scalars
2423 : integer,intent(in) :: usepaw
2424 : type(dataset_type),intent(in) :: dtset
2425 :
2426 : !arrays
2427 : integer,intent(in) :: atindx(dtset%natom)
2428 : real(dp),intent(out) :: omlamb(3)
2429 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
2430 :
2431 : !Local variables -------------------------
2432 : !scalars
2433 : integer :: adir,iat,iatom,itypat
2434 : real(dp) :: lambsig
2435 :
2436 : !--------------------------------------------------------------------
2437 :
2438 17 : omlamb = zero
2439 47 : do iat=1,dtset%natom
2440 30 : iatom = atindx(iat)
2441 30 : itypat = dtset%typat(iat)
2442 : ! if user input lambsig specifically in the input file, use it
2443 30 : if (abs(dtset%lambsig(itypat)).GT.tol8) then
2444 : lambsig=dtset%lambsig(itypat)
2445 : ! else if paw use the value read in to pawtab structure (which might well be zero)
2446 18 : else if (usepaw .EQ. 1) then
2447 16 : lambsig=pawtab(itypat)%lamb_shielding
2448 : else
2449 : lambsig = zero
2450 : end if
2451 137 : do adir = 1, 3
2452 120 : omlamb(adir) = omlamb(adir) - lambsig*dtset%nucdipmom(adir,iat)
2453 : end do ! end loop over adir
2454 : end do ! end loop over iat
2455 :
2456 17 : end subroutine lamb_core
2457 : !!***
2458 :
2459 : !!****f* ABINIT/nonlop_orbmag_nl
2460 : !! NAME
2461 : !! nonlop_orbmag_nl
2462 : !!
2463 : !! FUNCTION
2464 : !! nonlop-like routine for applying the very specific H1-E_nk*S1 operator
2465 : !! of orbital magnetism
2466 : !!
2467 : !! INPUTS
2468 : !!
2469 : !! OUTPUT
2470 : !! |vectout> = (H1-Enk*S1)|vectin>
2471 : !!
2472 : !! NOTES
2473 : !! See ZTG Eq. 42
2474 : !! computes on-site prefac*\sum_{Rij}<bra|d_bra_dir p_i>aij<d_ket_dir p_j|ket>
2475 : !! dnlbra = 0 if no derivative, dnlbra = adir,bdir,gdir for derivative in *dir direction
2476 : !! dnlket = 0 if no derivative, dnlket = adir,bdir,gdir for derivative in *dir direction
2477 : !!
2478 : !! SOURCE
2479 :
2480 12708 : subroutine nonlop_orbmag_nl(atindx,cwaveprj,dnlbra,dnlket,dterm,dtset,eignk,gs_hamk,npw_k,&
2481 12708 : & pawtab,vectout)
2482 : !Arguments ------------------------------------
2483 : !scalars
2484 : integer,intent(in) :: dnlbra,dnlket,npw_k
2485 : real(dp),intent(in) :: eignk
2486 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
2487 : type(dataset_type),intent(in) :: dtset
2488 : type(dterm_type),intent(in) :: dterm
2489 : !arrays
2490 : integer,intent(in) :: atindx(dtset%natom)
2491 : real(dp),intent(out) :: vectout(2,npw_k*dtset%nspinor)
2492 : type(pawcprj_type),intent(in) :: cwaveprj(dtset%natom,dtset%nspinor)
2493 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
2494 :
2495 : !Local variables -------------------------
2496 : !scalars
2497 : integer :: iat,iatom,il,ilmn,ipw,isp,itypat,jlmn,klmn,npwsp
2498 : real(dp) :: wt
2499 : complex(dp) :: cpj,dij,dij_cpj,proj_i
2500 : ! arrays
2501 : complex(dp),dimension(0:3) :: iexpl=(/cone,j_dpc,-cone,-j_dpc/)
2502 12708 : complex(dp),allocatable :: dij_data(:,:,:)
2503 : !--------------------------------------------------------------------
2504 :
2505 12708 : npwsp = npw_k*dtset%nspinor
2506 12708 : wt = four_pi/SQRT(gs_hamk%ucvol)
2507 :
2508 63540 : ABI_MALLOC(dij_data,(dtset%natom,dterm%lmn2max,dterm%ndij))
2509 2276364 : dij_data = dterm%aij - eignk*dterm%qij
2510 :
2511 8003592 : vectout = zero
2512 31944 : do iat = 1, dtset%natom
2513 19236 : iatom = atindx(iat)
2514 19236 : itypat=dtset%typat(iat)
2515 57324 : do isp = 1, dtset%nspinor
2516 247656 : do ilmn = 1, pawtab(itypat)%lmn_size
2517 :
2518 : dij_cpj = czero
2519 1827360 : do jlmn = 1, pawtab(itypat)%lmn_size
2520 :
2521 1624320 : if (dnlket .NE. 0) then
2522 1624320 : cpj=CMPLX(cwaveprj(iatom,isp)%dcp(1,dnlket,jlmn),cwaveprj(iatom,isp)%dcp(2,dnlket,jlmn))
2523 : else
2524 0 : cpj=CMPLX(cwaveprj(iatom,isp)%cp(1,jlmn),cwaveprj(iatom,isp)%cp(2,jlmn))
2525 : end if
2526 :
2527 1624320 : klmn=MATPACK(ilmn,jlmn)
2528 1624320 : dij = dij_data(iatom,klmn,isp)
2529 : ! see note at top of file near definition of MATPACK macro
2530 1624320 : if (ilmn .GT. jlmn) dij = CONJG(dij)
2531 :
2532 1624320 : dij_cpj = dij_cpj + dij*cpj
2533 :
2534 : ! in ndij = 4 case, isp 1 delivers up-up, isp 2 delivers down-down
2535 1827360 : if (dterm%ndij == 4) then
2536 786432 : if (isp == 1) then
2537 393216 : dij = dij_data(iatom,klmn,3) ! up-down
2538 : ! D^ss'_ij=D^s's_ji^*
2539 393216 : if (ilmn .GT. jlmn) dij = CONJG(dij_data(iatom,klmn,4))
2540 393216 : if (dnlket .NE. 0) then
2541 393216 : cpj=CMPLX(cwaveprj(iatom,2)%dcp(1,dnlket,jlmn),cwaveprj(iatom,2)%dcp(2,dnlket,jlmn))
2542 : else
2543 0 : cpj=CMPLX(cwaveprj(iatom,2)%cp(1,jlmn),cwaveprj(iatom,2)%cp(2,jlmn))
2544 : end if
2545 : else
2546 393216 : dij = dij_data(iatom,klmn,4) ! down-up
2547 : ! D^ss'_ij=D^s's_ji^*
2548 393216 : if (ilmn .GT. jlmn) dij = CONJG(dij_data(iatom,klmn,3))
2549 393216 : if (dnlket .NE. 0) then
2550 393216 : cpj=CMPLX(cwaveprj(iatom,1)%dcp(1,dnlket,jlmn),cwaveprj(iatom,1)%dcp(2,dnlket,jlmn))
2551 : else
2552 0 : cpj=CMPLX(cwaveprj(iatom,1)%cp(1,jlmn),cwaveprj(iatom,1)%cp(2,jlmn))
2553 : end if
2554 : end if
2555 786432 : dij_cpj = dij_cpj + dij*cpj
2556 : end if
2557 :
2558 : end do !jlmn
2559 :
2560 203040 : il = MOD(pawtab(itypat)%indlmn(1,ilmn),4)
2561 32769444 : do ipw = 1, npw_k
2562 : proj_i = CMPLX(gs_hamk%ffnl_k(ipw,1+dnlbra,ilmn,itypat),0.0)*wt*CONJG(iexpl(il))*dij_cpj*&
2563 32541024 : & CMPLX(gs_hamk%ph3d_k(1,ipw,iatom),-gs_hamk%ph3d_k(2,ipw,iatom))
2564 :
2565 32541024 : vectout(1,npw_k*(isp-1)+ipw) = vectout(1,npw_k*(isp-1)+ipw) + REAL(proj_i)
2566 32744064 : vectout(2,npw_k*(isp-1)+ipw) = vectout(2,npw_k*(isp-1)+ipw) + AIMAG(proj_i)
2567 :
2568 : end do
2569 :
2570 : end do !ilmn
2571 : end do ! isp
2572 : end do !iat
2573 :
2574 12708 : ABI_SFREE(dij_data)
2575 :
2576 12708 : end subroutine nonlop_orbmag_nl
2577 : !!***
2578 :
2579 : !!!****f* ABINIT/cprj_test
2580 : !!! NAME
2581 : !!! cprj_test
2582 : !!!
2583 : !!! FUNCTION
2584 : !!! General application of <bra|p>dij<p|ket>
2585 : !!!
2586 : !!! INPUTS
2587 : !!!
2588 : !!! OUTPUT
2589 : !!! complex(dp) nlme
2590 : !!!
2591 : !!! NOTES
2592 : !!! computes on-site prefac*\sum_{Rij}<bra|d_bra_dir p_i>aij<d_ket_dir p_j|ket>
2593 : !!! dnlbra = 0 if no derivative, dnlbra = adir,bdir,gdir for derivative in *dir direction
2594 : !!! dnlket = 0 if no derivative, dnlket = adir,bdir,gdir for derivative in *dir direction
2595 : !!!
2596 : !!! SOURCE
2597 : !
2598 : !subroutine cprj_test(adir,atindx,cwaveprj,dnlbra,dnlket,dterm,dtset,&
2599 : ! & eignk,fermie,gs_hamk,npw_k,oterm,prefac,pawtab,vectin,vectout)
2600 : ! !Arguments ------------------------------------
2601 : ! !scalars
2602 : ! integer,intent(in) :: adir,dnlbra,dnlket,npw_k,oterm
2603 : ! real(dp),intent(in) :: eignk,fermie
2604 : ! complex(dp),intent(in) :: prefac
2605 : ! type(gs_hamiltonian_type),intent(inout) :: gs_hamk
2606 : ! type(dataset_type),intent(in) :: dtset
2607 : ! type(dterm_type),intent(in) :: dterm
2608 : ! !arrays
2609 : ! integer,intent(in) :: atindx(dtset%natom)
2610 : ! real(dp),intent(in),pointer :: vectin(:,:)
2611 : ! real(dp),intent(out) :: vectout(2,npw_k)
2612 : ! type(pawcprj_type),intent(in) :: cwaveprj(dtset%natom,dtset%nspinor)
2613 : ! type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
2614 : !
2615 : ! !Local variables -------------------------
2616 : ! !scalars
2617 : ! integer :: iat,iatom,il,ilmn,ipw,isp,itypat,jlmn,klmn,npwsp
2618 : ! real(dp) :: wt
2619 : ! logical :: il_parity
2620 : ! complex(dp) :: cgfac,my_cpi,cpi,cpj,dij,dij_cpj,phfac,proj_i,scale_fac
2621 : ! ! arrays
2622 : ! complex(dp),dimension(0:3) :: iexpl=(/cone,j_dpc,-cone,-j_dpc/)
2623 : ! complex(dp),allocatable :: dij_data(:,:,:)
2624 : !!--------------------------------------------------------------------
2625 : !
2626 : ! npwsp = npw_k*dtset%nspinor
2627 : ! wt = four_pi/SQRT(gs_hamk%ucvol)
2628 : !
2629 : ! ABI_MALLOC(dij_data,(dtset%natom,dterm%lmn2max,dterm%ndij))
2630 : ! select case (oterm)
2631 : ! !case (inlr)
2632 : ! ! dij_data = dterm%LR(:,:,:,adir)
2633 : ! !case (inbm)
2634 : ! ! dij_data = dterm%BM(:,:,:,adir)
2635 : ! case (innl)
2636 : ! dij_data = dterm%aij - eignk*dterm%qij
2637 : ! case (incc)
2638 : ! dij_data = dterm%aij + (eignk-two*fermie)*dterm%qij
2639 : ! case (invv1)
2640 : ! dij_data = (eignk-fermie)*dterm%qij
2641 : ! case (invv2)
2642 : ! dij_data = (eignk-fermie)*dterm%qij
2643 : ! case DEFAULT
2644 : ! dij_data = czero
2645 : ! end select
2646 : !
2647 : ! !nlme = czero
2648 : ! vectout = zero
2649 : ! do iat = 1, dtset%natom
2650 : ! iatom = atindx(iat)
2651 : ! itypat=dtset%typat(iat)
2652 : ! do isp = 1, dtset%nspinor
2653 : ! do ilmn = 1, pawtab(itypat)%lmn_size
2654 : !
2655 : ! if (dnlbra .NE. 0) then
2656 : ! cpi=CMPLX(cwaveprj(iatom,isp)%dcp(1,dnlbra,ilmn),cwaveprj(iatom,isp)%dcp(2,dnlbra,ilmn))
2657 : ! else
2658 : ! cpi=CMPLX(cwaveprj(iatom,isp)%cp(1,ilmn),cwaveprj(iatom,isp)%cp(2,ilmn))
2659 : ! end if
2660 : !
2661 : ! ! this loop reconstructs cpi and dcpi correctly
2662 : ! !il = MOD(pawtab(itypat)%indlmn(1,ilmn),4)
2663 : ! !my_cpi = czero
2664 : ! !do ipw = 1,npw_k
2665 : ! ! cgfac=CMPLX(vectin(1,ipw),vectin(2,ipw))
2666 : ! ! phfac=CMPLX(gs_hamk%ph3d_k(1,ipw,iatom),gs_hamk%ph3d_k(2,ipw,iatom))
2667 : ! ! my_cpi = my_cpi + cgfac*phfac*gs_hamk%ffnl_k(ipw,1+dnlbra,ilmn,itypat)
2668 : ! !end do
2669 : ! !my_cpi = my_cpi*four_pi*iexpl(il)/SQRT(gs_hamk%ucvol)
2670 : ! !write(std_out,'(a,4es16.8)')'JWZ debug cpi, mycpi : ',&
2671 : ! ! & REAL(cpi),AIMAG(cpi),REAL(my_cpi),AIMAG(my_cpi)
2672 : !
2673 : ! ! |p_i> in recip space
2674 : ! !proj_i = czero
2675 : ! !il = pawtab(itypat)%indlmn(1,ilmn)
2676 : ! !do ipw = 1, npw_k
2677 : ! ! phfac=CMPLX(gs_hamk%ph3d_k(1,ipw,iatom),gs_hamk%ph3d_k(2,ipw,iatom))
2678 : ! ! proj_i(ipw) = CMPLX(gs_hamk%ffnl_k(ipw,1+dnlbra,ilmn,itypat),0.0)*&
2679 : ! ! &(j_dpc**il)*four_pi*CONJG(phfac)/SQRT(gs_hamk%ucvol)
2680 : ! !end do
2681 : !
2682 : ! dij_cpj = czero
2683 : ! do jlmn = 1, pawtab(itypat)%lmn_size
2684 : !
2685 : ! if (dnlket .NE. 0) then
2686 : ! cpj=CMPLX(cwaveprj(iatom,isp)%dcp(1,dnlket,jlmn),cwaveprj(iatom,isp)%dcp(2,dnlket,jlmn))
2687 : ! else
2688 : ! cpj=CMPLX(cwaveprj(iatom,isp)%cp(1,jlmn),cwaveprj(iatom,isp)%cp(2,jlmn))
2689 : ! end if
2690 : !
2691 : ! klmn=MATPACK(ilmn,jlmn)
2692 : ! if ( oterm == inlr ) then
2693 : ! dij = CMPLX(dterm%ekb_LR(2*klmn-1,iatom,1,adir),dterm%ekb_LR(2*klmn,iatom,1,adir))
2694 : ! else if ( oterm == inbm ) then
2695 : ! dij = CMPLX(dterm%ekb_BM(2*klmn-1,iatom,1,adir),dterm%ekb_BM(2*klmn,iatom,1,adir))
2696 : ! else
2697 : ! dij = dij_data(iatom,klmn,1)
2698 : ! end if
2699 : ! ! see note at top of file near definition of MATPACK macro
2700 : ! if (ilmn .GT. jlmn) dij = CONJG(dij)
2701 : !
2702 : ! dij_cpj = dij_cpj + dij*cpj
2703 : !
2704 : ! end do !jlmn
2705 : !
2706 : ! il = MOD(pawtab(itypat)%indlmn(1,ilmn),4)
2707 : ! do ipw = 1, npw_k
2708 : ! proj_i = CMPLX(gs_hamk%ffnl_k(ipw,1+dnlbra,ilmn,itypat),0.0)*wt*CONJG(iexpl(il))*dij_cpj*&
2709 : ! & CMPLX(gs_hamk%ph3d_k(1,ipw,iatom),-gs_hamk%ph3d_k(2,ipw,iatom))
2710 : !
2711 : ! vectout(1,ipw) = vectout(1,ipw) + REAL(proj_i)
2712 : ! vectout(2,ipw) = vectout(2,ipw) + AIMAG(proj_i)
2713 : !
2714 : ! end do
2715 : !
2716 : ! end do !ilmn
2717 : ! end do ! isp
2718 : ! end do !iat
2719 : !
2720 : ! ABI_SFREE(dij_data)
2721 : !
2722 : !end subroutine cprj_test
2723 : !!!***
2724 :
2725 : !!****f* ABINIT/dterm_qij
2726 : !! NAME
2727 : !! dterm_qij
2728 : !!
2729 : !! FUNCTION
2730 : !! Transfer pawtab%sij to dterm, as complex, solely for convenience
2731 : !!
2732 : !! INPUTS
2733 : !! atindx(natom)=index table for atoms (see gstate.f)
2734 : !! dtset <type(dataset_type)>=all input variables for this dataset
2735 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
2736 : !!
2737 : !! OUTPUT
2738 : !!
2739 : !! SIDE EFFECTS
2740 : !! dterm <type(dterm_type)> data related to onsite interactions
2741 : !!
2742 : !! NOTES
2743 : !! Transfer pawtab%sij to dterm, as complex, solely for convenience
2744 : !!
2745 : !! SOURCE
2746 :
2747 16 : subroutine dterm_qij(atindx,dterm,dtset,pawtab)
2748 :
2749 : !Arguments ------------------------------------
2750 : !scalars
2751 : type(dterm_type),intent(inout) :: dterm
2752 : type(dataset_type),intent(in) :: dtset
2753 :
2754 : !arrays
2755 : integer,intent(in) :: atindx(dtset%natom)
2756 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
2757 :
2758 : !Local variables -------------------------
2759 : !scalars
2760 : integer :: iat,iatom,itypat,lmn2_size
2761 :
2762 : !arrays
2763 :
2764 : !--------------------------------------------------------------------
2765 :
2766 3033 : dterm%qij = czero
2767 44 : do iat = 1, dtset%natom
2768 28 : iatom = atindx(iat)
2769 28 : itypat = dtset%typat(iat)
2770 28 : lmn2_size = pawtab(itypat)%lmn2_size
2771 : dterm%qij(iatom,1:lmn2_size,1) = &
2772 1036 : & CMPLX(pawtab(itypat)%sij(1:lmn2_size),zero)
2773 44 : if (dterm%ndij > 1) then
2774 : dterm%qij(iatom,1:lmn2_size,2) = &
2775 370 : & CMPLX(pawtab(itypat)%sij(1:lmn2_size),zero)
2776 : end if
2777 : end do ! iat
2778 :
2779 16 : dterm%has_qij=2
2780 :
2781 16 : end subroutine dterm_qij
2782 : !!***
2783 :
2784 : !!****f* ABINIT/dterm_BM
2785 : !! NAME
2786 : !! dterm_BM
2787 : !!
2788 : !! FUNCTION
2789 : !! Compute onsite <A0.AN>
2790 : !!
2791 : !! INPUTS
2792 : !! atindx(natom)=index table for atoms (see gstate.f)
2793 : !! dtset <type(dataset_type)>=all input variables for this dataset
2794 : !! gntselect((2*my_lmax-1)**2,my_lmax**2*(my_lmax**2+1)/2)=nonzero gaunt integral indices
2795 : !! gprimd(3,3)=reciprocal space lattice vectors
2796 : !! my_lmax=augmented l_max over all psp
2797 : !! pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data
2798 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
2799 : !! realgnt((2*my_lmax-1)**2*(my_lmax)**4)=nonzero gaunt integral values
2800 : !!
2801 : !! OUTPUT
2802 : !!
2803 : !! SIDE EFFECTS
2804 : !! dterm <type(dterm_type)> data related to onsite interactions
2805 : !!
2806 : !! NOTES
2807 : !! ZTG23 Eq. 43
2808 : !! this term is A0.An = \frac{1}{2}(B x r).\alpha^2(m x r) which can be rewritten
2809 : !! as \frac{\alpha^2}{2} [B.(1-\hat{r}\hat{r}).m]/r .
2810 : !!
2811 : !! SOURCE
2812 :
2813 16 : subroutine dterm_BM(atindx,dterm,dtset,gntselect,gprimd,my_lmax,pawrad,pawtab,realgnt)
2814 :
2815 : !Arguments ------------------------------------
2816 : !scalars
2817 : integer,intent(in) :: my_lmax
2818 : type(dterm_type),intent(inout) :: dterm
2819 : type(dataset_type),intent(in) :: dtset
2820 :
2821 : !arrays
2822 : integer,intent(in) :: atindx(dtset%natom)
2823 : integer,intent(in) :: gntselect((2*my_lmax-1)**2,my_lmax**2*(my_lmax**2+1)/2)
2824 : real(dp),intent(in) :: gprimd(3,3),realgnt((2*my_lmax-1)**2*(my_lmax)**4)
2825 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
2826 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
2827 :
2828 : !Local variables -------------------------
2829 : !scalars
2830 : integer :: adir,iat,iatom,itypat,gs1,gs2
2831 : integer :: klmn,klm,kln,mdir,mesh_size,ngnt,pwave_size
2832 : real(dp) :: a2,afact,intg
2833 :
2834 : !arrays
2835 : complex(dp) :: dij_cart(3),dij_red(3)
2836 16 : real(dp),allocatable :: dyadic(:,:,:),ff(:),radint(:)
2837 : !--------------------------------------------------------------------
2838 :
2839 11977 : dterm%ekb_BM = zero
2840 : a2 = FineStructureConstant2
2841 :
2842 44 : do iat = 1, dtset%natom
2843 64 : if(.NOT. ANY(ABS(dtset%nucdipmom(:,iat))>tol8) ) cycle
2844 :
2845 16 : iatom = atindx(iat)
2846 16 : itypat = dtset%typat(iat)
2847 :
2848 16 : mesh_size=pawtab(itypat)%mesh_size
2849 16 : pwave_size=size(pawtab(itypat)%phiphj(:,1))
2850 :
2851 : ! compute angular integrals of S_i (1-rr) S_j
2852 16 : gs1=size(gntselect,1)
2853 16 : gs2=size(gntselect,2)
2854 16 : ngnt=size(realgnt)
2855 48 : ABI_MALLOC(dyadic,(3,3,gs2))
2856 16 : call make_dyadic(one,one,dyadic,gntselect,gs1,gs2,gs2,ngnt,realgnt)
2857 :
2858 : ! compute radial integrals of (ui*uj - tilde{ui}tilde{uj})/r
2859 48 : ABI_MALLOC(radint,(pawtab(itypat)%ij_size))
2860 48 : ABI_MALLOC(ff,(mesh_size))
2861 176 : do kln=1,pawtab(itypat)%ij_size
2862 : ff(2:pwave_size) = &
2863 : & (pawtab(itypat)%phiphj(2:pwave_size,kln)-pawtab(itypat)%tphitphj(2:pwave_size,kln))/&
2864 220730 : & (pawrad(itypat)%rad(2:pwave_size))
2865 160 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
2866 160 : call simp_gen(intg,ff,pawrad(itypat))
2867 176 : radint(kln)=intg
2868 : end do
2869 16 : ABI_FREE(ff)
2870 :
2871 592 : do klmn=1, pawtab(itypat)%lmn2_size
2872 576 : klm = pawtab(itypat)%indklmn(1,klmn)
2873 576 : kln = pawtab(itypat)%indklmn(2,klmn)
2874 :
2875 576 : dij_cart=zero
2876 2304 : do adir = 1, 3 ! B field direction
2877 7488 : do mdir = 1, 3 ! mag dipole direction
2878 5184 : afact=half*a2*radint(kln)*dyadic(adir,mdir,klm)*dtset%nucdipmom(mdir,iat)
2879 6912 : dij_cart(adir)=dij_cart(adir)-CMPLX(afact,zero)
2880 : end do
2881 : end do
2882 :
2883 14400 : dij_red = MATMUL(TRANSPOSE(gprimd),dij_cart)
2884 :
2885 2320 : do adir=1, 3
2886 1728 : dterm%ekb_BM(2*klmn-1,iat,1,adir) = REAL(dij_red(adir))
2887 1728 : dterm%ekb_BM(2*klmn,iat,1,adir) = AIMAG(dij_red(adir))
2888 2304 : if (dterm%ndij > 1) then
2889 540 : dterm%ekb_BM(2*klmn-1,iat,2,adir) = dterm%ekb_BM(2*klmn-1,iat,1,adir)
2890 540 : dterm%ekb_BM(2*klmn,iat,2,adir) = dterm%ekb_BM(2*klmn,iat,1,adir)
2891 : end if
2892 : end do
2893 :
2894 : end do ! end loop over klmn
2895 :
2896 16 : ABI_FREE(dyadic)
2897 44 : ABI_FREE(radint)
2898 : end do ! end loop over iatom
2899 :
2900 16 : dterm%has_BM = 2
2901 :
2902 16 : end subroutine dterm_BM
2903 : !!***
2904 :
2905 : !!****f* ABINIT/dterm_LR
2906 : !! NAME
2907 : !! dterm_LR
2908 : !!
2909 : !! FUNCTION
2910 : !! Compute onsite <L_R/2>
2911 : !!
2912 : !! INPUTS
2913 : !! atindx(natom)=index table for atoms (see gstate.f)
2914 : !! dtset <type(dataset_type)>=all input variables for this dataset
2915 : !! gprimd(3,3)=reciprocal space lattice vectors
2916 : !! pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data
2917 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
2918 : !!
2919 : !! OUTPUT
2920 : !!
2921 : !! SIDE EFFECTS
2922 : !! dterm <type(dterm_type)> data related to onsite interactions
2923 : !!
2924 : !! NOTES
2925 : !! ZTG23 text after Eq 42, the on-site angular momentum
2926 : !!
2927 : !! SOURCE
2928 :
2929 16 : subroutine dterm_LR(atindx,dterm,dtset,gprimd,pawrad,pawtab)
2930 :
2931 : !Arguments ------------------------------------
2932 : !scalars
2933 : type(dterm_type),intent(inout) :: dterm
2934 : type(dataset_type),intent(in) :: dtset
2935 :
2936 : !arrays
2937 : integer,intent(in) :: atindx(dtset%natom)
2938 : real(dp),intent(in) :: gprimd(3,3)
2939 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
2940 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
2941 :
2942 : !Local variables -------------------------
2943 : !scalars
2944 : integer :: adir,iat,iatom,ilmn,il,im,itypat,jlmn,jl,jm
2945 : integer :: klmn,kln,mesh_size,pwave_size
2946 : real(dp) :: intg
2947 : complex(dp) :: orbl_me
2948 : !arrays
2949 : complex(dp) :: dij_cart(3),dij_red(3)
2950 16 : real(dp),allocatable :: ff(:)
2951 : !--------------------------------------------------------------------
2952 :
2953 11977 : dterm%ekb_LR = zero
2954 :
2955 44 : do iat = 1, dtset%natom
2956 28 : iatom = atindx(iat)
2957 28 : itypat = dtset%typat(iat)
2958 28 : mesh_size=pawtab(itypat)%mesh_size
2959 28 : pwave_size=size(pawtab(itypat)%phiphj(:,1))
2960 84 : ABI_MALLOC(ff,(mesh_size))
2961 1036 : do klmn=1, pawtab(itypat)%lmn2_size
2962 :
2963 1008 : ilmn = pawtab(itypat)%indklmn(7,klmn)
2964 1008 : il=pawtab(itypat)%indlmn(1,ilmn)
2965 1008 : im=pawtab(itypat)%indlmn(2,ilmn)
2966 :
2967 1008 : jlmn = pawtab(itypat)%indklmn(8,klmn)
2968 1008 : jl=pawtab(itypat)%indlmn(1,jlmn)
2969 1008 : jm=pawtab(itypat)%indlmn(2,jlmn)
2970 :
2971 1008 : if ( il /= jl ) cycle ! <l'm'|L|lm> = 0 if l' /= l
2972 672 : if ( il == 0 ) cycle ! <00|L|00> = 0
2973 :
2974 588 : kln = pawtab(itypat)%indklmn(2,klmn)
2975 815157 : ff=0
2976 814569 : ff(2:pwave_size) = pawtab(itypat)%phiphj(2:pwave_size,kln)-pawtab(itypat)%tphitphj(2:pwave_size,kln)
2977 588 : call pawrad_deducer0(ff,mesh_size,pawrad(itypat))
2978 588 : call simp_gen(intg,ff,pawrad(itypat))
2979 :
2980 2352 : do adir = 1, 3
2981 : ! compute <L_dir>/2
2982 1764 : call slxyzs(il,im,adir,jl,jm,orbl_me)
2983 2352 : dij_cart(adir) = -half*orbl_me*intg
2984 : end do ! end loop over adir
2985 :
2986 : ! convert to crystal frame
2987 14700 : dij_red = MATMUL(TRANSPOSE(gprimd),dij_cart)
2988 :
2989 2968 : do adir = 1, 3
2990 1764 : dterm%ekb_LR(2*klmn-1,iat,1,adir) = REAL(dij_red(adir))
2991 1764 : dterm%ekb_LR(2*klmn,iat,1,adir) = AIMAG(dij_red(adir))
2992 2352 : if (dterm%ndij > 1) then
2993 630 : dterm%ekb_LR(2*klmn-1,iat,2,adir) = dterm%ekb_LR(2*klmn-1,iat,1,adir)
2994 630 : dterm%ekb_LR(2*klmn,iat,2,adir) = dterm%ekb_LR(2*klmn,iat,1,adir)
2995 : end if
2996 : end do
2997 : end do ! end loop over klmn
2998 44 : ABI_FREE(ff)
2999 : end do ! end loop over itypat
3000 :
3001 16 : dterm%has_LR = 2
3002 :
3003 16 : end subroutine dterm_LR
3004 : !!***
3005 :
3006 : !!****f* ABINIT/orbmag_output
3007 : !! NAME
3008 : !! orbmag_output
3009 : !!
3010 : !! FUNCTION
3011 : !! Only printing. This routine outputs orbmag terms to the normal abinit output file
3012 : !!
3013 : !! INPUTS
3014 : !! dtset <type(dataset_type)>=all input variables for this dataset
3015 : !! omlamb(3)=Lamb shielding
3016 : !!
3017 : !! OUTPUT
3018 : !!
3019 : !! SIDE EFFECTS
3020 : !!
3021 : !! SOURCE
3022 :
3023 17 : subroutine orbmag_output(self,dtset,omlamb)
3024 :
3025 : !Arguments ------------------------------------
3026 : !scalars
3027 : class(orbmag_mesh_type),intent(inout),target :: self
3028 : type(dataset_type),intent(in) :: dtset
3029 :
3030 : !arrays
3031 : real(dp),intent(in) :: omlamb(3)
3032 :
3033 : !Local variables -------------------------
3034 : !scalars
3035 : integer :: adir,iband,isppol,iterms
3036 : character(len=500) :: message
3037 :
3038 : !arrays
3039 17 : real(dp) :: berry_bb(dtset%mband,3),berry_total(3),orbmag_bb(dtset%mband,3),orbmag_total(3)
3040 :
3041 : ! ***********************************************************************
3042 :
3043 332 : orbmag_bb=zero;orbmag_total=zero
3044 119 : do iterms = 1,orbmag_nterms
3045 408 : orbmag_total(1:3)=orbmag_total(1:3) + self%orbmag_trace(1:3,iterms)
3046 227 : do isppol = 1, dtset%nsppol
3047 762 : do iband=1, dtset%mband
3048 2316 : orbmag_bb(iband,1:3) = orbmag_bb(iband,1:3) + self%orbmag_terms(iband,isppol,1:3,iterms)
3049 : end do ! iband
3050 : end do ! isppol
3051 : end do
3052 :
3053 68 : orbmag_total=orbmag_total+omlamb
3054 :
3055 332 : berry_bb=zero;berry_total=zero
3056 51 : do iterms = 1,chern_nterms
3057 136 : berry_total(1:3)=berry_total(1:3) + self%chern_trace(1:3,iterms)
3058 87 : do isppol = 1, dtset%nsppol
3059 254 : do iband=1, dtset%mband
3060 772 : berry_bb(iband,1:3) = berry_bb(iband,1:3) + self%chern_terms(iband,isppol,1:3,iterms)
3061 : end do ! iband
3062 : end do ! isppol
3063 : end do
3064 :
3065 17 : write(message,'(a,a,a)')ch10,'====================================================',ch10
3066 17 : call wrtout(ab_out,message,'COLL')
3067 :
3068 17 : write(message,'(a,a)')' Orbital magnetic moment computed with DFPT derivative wavefunctions ',ch10
3069 17 : call wrtout(ab_out,message,'COLL')
3070 :
3071 17 : write(message,'(a)')' Orbital magnetic moment, Cartesian directions : '
3072 17 : call wrtout(ab_out,message,'COLL')
3073 17 : write(message,'(3es16.8)') (orbmag_total(adir),adir=1,3)
3074 17 : call wrtout(ab_out,message,'COLL')
3075 17 : write(message,'(a)')ch10
3076 17 : call wrtout(ab_out,message,'COLL')
3077 17 : write(message,'(a)')' Chern vector, Cartesian directions : '
3078 17 : call wrtout(ab_out,message,'COLL')
3079 17 : write(message,'(3es16.8)') (berry_total(adir),adir=1,3)
3080 17 : call wrtout(ab_out,message,'COLL')
3081 :
3082 17 : if(abs(dtset%orbmag) .GE. 2) then
3083 17 : write(message,'(a)')ch10
3084 17 : call wrtout(ab_out,message,'COLL')
3085 17 : write(message,'(a)')' Orbital magnetic moment, term-by-term breakdown : '
3086 17 : call wrtout(ab_out,message,'COLL')
3087 68 : write(message,'(a,3es16.8)') ' rho(1) CC : ',(self%orbmag_trace(adir,incc),adir=1,3)
3088 17 : call wrtout(ab_out,message,'COLL')
3089 68 : write(message,'(a,3es16.8)') ' rho(1) VV1 : ',(self%orbmag_trace(adir,invv1),adir=1,3)
3090 17 : call wrtout(ab_out,message,'COLL')
3091 68 : write(message,'(a,3es16.8)') ' rho(1) VV2 : ',(self%orbmag_trace(adir,invv2),adir=1,3)
3092 17 : call wrtout(ab_out,message,'COLL')
3093 68 : write(message,'(a,3es16.8)') ' rho(0) NL : ',(self%orbmag_trace(adir,innl),adir=1,3)
3094 17 : call wrtout(ab_out,message,'COLL')
3095 68 : write(message,'(a,3es16.8)') ' <L_R> : ',(self%orbmag_trace(adir,inlr),adir=1,3)
3096 17 : call wrtout(ab_out,message,'COLL')
3097 68 : write(message,'(a,3es16.8)') ' <A0.An> : ',(self%orbmag_trace(adir,inbm),adir=1,3)
3098 17 : call wrtout(ab_out,message,'COLL')
3099 17 : write(message,'(a,3es16.8)') ' Lamb (core) : ',(omlamb(adir),adir=1,3)
3100 17 : call wrtout(ab_out,message,'COLL')
3101 17 : write(message,'(a)')ch10
3102 17 : call wrtout(ab_out,message,'COLL')
3103 17 : write(message,'(a)')' Chern vector, term-by-term breakdown : '
3104 17 : call wrtout(ab_out,message,'COLL')
3105 68 : write(message,'(a,3es16.8)') ' Chern CC : ',(self%chern_trace(adir,ibcc),adir=1,3)
3106 17 : call wrtout(ab_out,message,'COLL')
3107 68 : write(message,'(a,3es16.8)') ' Chern VV1 : ',(self%chern_trace(adir,ibvv1),adir=1,3)
3108 17 : call wrtout(ab_out,message,'COLL')
3109 : end if
3110 :
3111 17 : if(abs(dtset%orbmag) .EQ. 4) then
3112 1 : write(message,'(a)')ch10
3113 1 : call wrtout(ab_out,message,'COLL')
3114 1 : write(message,'(a)')' Integrals of densities on real space mesh, should be close to values above : '
3115 1 : call wrtout(ab_out,message,'COLL')
3116 4 : write(message,'(a,3es16.8)') ' rho(1) CC : ',(self%rmesh_int(adir,incc),adir=1,3)
3117 1 : call wrtout(ab_out,message,'COLL')
3118 4 : write(message,'(a,3es16.8)') ' rho(1) VV1 : ',(self%rmesh_int(adir,invv1),adir=1,3)
3119 1 : call wrtout(ab_out,message,'COLL')
3120 4 : write(message,'(a,3es16.8)') ' rho(1) VV2 : ',(self%rmesh_int(adir,invv2),adir=1,3)
3121 1 : call wrtout(ab_out,message,'COLL')
3122 4 : write(message,'(a,3es16.8)') ' rho(0) NL : ',(self%rmesh_int(adir,innl),adir=1,3)
3123 1 : call wrtout(ab_out,message,'COLL')
3124 4 : write(message,'(a,3es16.8)') ' <L_R> : ',(self%rmesh_int(adir,inlr),adir=1,3)
3125 1 : call wrtout(ab_out,message,'COLL')
3126 4 : write(message,'(a,3es16.8)') ' <A0.An> : ',(self%rmesh_int(adir,inbm),adir=1,3)
3127 1 : call wrtout(ab_out,message,'COLL')
3128 : end if
3129 :
3130 :
3131 17 : write(message,'(a,a,a)')ch10,'====================================================',ch10
3132 17 : call wrtout(ab_out,message,'COLL')
3133 :
3134 17 : end subroutine orbmag_output
3135 : !!***
3136 :
3137 : !!****f* ABINIT/dterm_free
3138 : !! NAME
3139 : !! dterm_free
3140 : !!
3141 : !! FUNCTION
3142 : !! free space in dterm_type
3143 : !!
3144 : !! SOURCE
3145 :
3146 16 : subroutine dterm_free(self)
3147 :
3148 : !Arguments ------------------------------------
3149 : !scalars
3150 : class(dterm_type),intent(inout),target :: self
3151 : !--------------------------------------------------------------------
3152 :
3153 16 : ABI_SFREE(self%aij)
3154 16 : self%has_aij=0
3155 :
3156 16 : ABI_SFREE(self%qij)
3157 16 : self%has_qij=0
3158 :
3159 16 : ABI_SFREE(self%ekb_LR)
3160 16 : self%has_LR=0
3161 :
3162 16 : ABI_SFREE(self%ekb_BM)
3163 16 : self%has_BM=0
3164 :
3165 16 : end subroutine dterm_free
3166 : !!***
3167 :
3168 : !!****f* ABINIT/dterm_init
3169 : !! NAME
3170 : !! dterm_init
3171 : !!
3172 : !! FUNCTION
3173 : !! allocate space in dterm_type
3174 : !!
3175 : !! INPUTS
3176 : !! lmnmax=max value of lmn over all psps
3177 : !! lmn2max=max value of lmn2 over all psps
3178 : !! natom=number of atoms in cell
3179 : !! ndij=spin channels in dij
3180 : !!
3181 : !! OUTPUT
3182 : !!
3183 : !! SIDE EFFECTS
3184 : !! dterm <type(dterm_type)> data related to onsite interactions
3185 : !!
3186 : !! SOURCE
3187 :
3188 16 : subroutine dterm_init(self,lmnmax,lmn2max,natom,ndij)
3189 :
3190 : !Arguments ------------------------------------
3191 : !scalars
3192 : class(dterm_type),intent(inout),target :: self
3193 : integer,intent(in) :: lmnmax,lmn2max,natom,ndij
3194 : !--------------------------------------------------------------------
3195 :
3196 16 : self%lmnmax = lmnmax
3197 16 : self%lmn2max = lmn2max
3198 16 : self%natom = natom
3199 16 : self%ndij = ndij
3200 :
3201 80 : ABI_REMALLOC(self%aij,(natom,lmn2max,ndij))
3202 16 : self%has_aij=1
3203 :
3204 64 : ABI_REMALLOC(self%qij,(natom,lmn2max,ndij))
3205 16 : self%has_qij=1
3206 :
3207 96 : ABI_REMALLOC(self%ekb_LR,(2*lmn2max,natom,ndij,3))
3208 16 : self%has_LR=1
3209 :
3210 64 : ABI_REMALLOC(self%ekb_BM,(2*lmn2max,natom,ndij,3))
3211 16 : self%has_BM=1
3212 :
3213 16 : end subroutine dterm_init
3214 : !!***
3215 :
3216 : !!****f* ABINIT/orbmag_init
3217 : !! NAME
3218 : !! orbmag_init
3219 : !!
3220 : !! FUNCTION
3221 : !! allocate space in orbmag_mesh_type
3222 : !!
3223 : !! INPUTS
3224 : !! dtset <type(dataset_type)>=all input variables for this dataset
3225 : !!
3226 : !! OUTPUT
3227 : !!
3228 : !! SIDE EFFECTS
3229 : !! orbmag_mesh <type(orbmag_mesh_type)> data related to orbmag terms on kpt mesh
3230 : !!
3231 : !! SOURCE
3232 :
3233 17 : subroutine orbmag_init(self,dtset)
3234 :
3235 : !Arguments ------------------------------------
3236 : !scalars
3237 : class(orbmag_mesh_type),intent(inout),target :: self
3238 : type(dataset_type),intent(in) :: dtset
3239 : !--------------------------------------------------------------------
3240 :
3241 17 : self%mband = dtset%mband
3242 17 : self%nkpt = dtset%nkpt
3243 17 : self%nsppol = dtset%nsppol
3244 17 : self%natom = dtset%natom
3245 17 : self%ntypat = dtset%ntypat
3246 17 : self%n4=dtset%ngfft(4)
3247 17 : self%n5=dtset%ngfft(5)
3248 17 : self%n6=dtset%ngfft(6)
3249 17 : self%chern_nterms = chern_nterms
3250 17 : self%orbmag_nterms = orbmag_nterms
3251 :
3252 51 : ABI_REMALLOC(self%lambsig,(self%ntypat))
3253 51 : ABI_REMALLOC(self%nucdipmom,(3,self%natom))
3254 119 : ABI_REMALLOC(self%cmesh,(self%mband,self%nkpt,self%nsppol,3,chern_nterms))
3255 30687 : self%cmesh=zero
3256 102 : ABI_REMALLOC(self%chern_terms,(dtset%mband,dtset%nsppol,3,chern_nterms))
3257 813 : self%chern_terms=zero
3258 17 : ABI_REMALLOC(self%chern_trace,(3,chern_nterms))
3259 153 : self%chern_trace=zero
3260 119 : ABI_REMALLOC(self%omesh,(self%mband,self%nkpt,self%nsppol,3,orbmag_nterms))
3261 92027 : self%omesh=zero
3262 85 : ABI_REMALLOC(self%orbmag_terms,(dtset%mband,dtset%nsppol,3,orbmag_nterms))
3263 2405 : self%orbmag_terms=zero
3264 17 : ABI_REMALLOC(self%orbmag_trace,(3,orbmag_nterms))
3265 425 : self%orbmag_trace=zero
3266 17 : if (dtset%orbmag .EQ. 4) then
3267 7 : ABI_REMALLOC(self%rmesh,(self%n4,self%n5,self%n6,3,orbmag_nterms))
3268 1677535 : self%rmesh=zero
3269 1 : ABI_REMALLOC(self%rmesh_int,(3,orbmag_nterms))
3270 25 : self%rmesh_int=zero
3271 : end if
3272 :
3273 17 : end subroutine orbmag_init
3274 : !!***
3275 :
3276 : !!****f* ABINIT/orbmag_free
3277 : !! NAME
3278 : !! orbmag_free
3279 : !!
3280 : !! FUNCTION
3281 : !! free space in orbmag_mesh_type
3282 : !!
3283 : !! SIDE EFFECTS
3284 : !! orbmag_mesh <type(orbmag_mesh_type)> data related to orbmag terms on kpt mesh
3285 : !!
3286 : !! SOURCE
3287 :
3288 17 : subroutine orbmag_free(self)
3289 :
3290 : !Arguments ------------------------------------
3291 : !scalars
3292 : class(orbmag_mesh_type),intent(inout),target :: self
3293 : !--------------------------------------------------------------------
3294 :
3295 17 : ABI_SFREE(self%lambsig)
3296 17 : ABI_SFREE(self%nucdipmom)
3297 17 : ABI_SFREE(self%cmesh)
3298 17 : ABI_SFREE(self%chern_terms)
3299 17 : ABI_SFREE(self%chern_trace)
3300 17 : ABI_SFREE(self%omesh)
3301 17 : ABI_SFREE(self%orbmag_terms)
3302 17 : ABI_SFREE(self%orbmag_trace)
3303 17 : ABI_SFREE(self%rmesh)
3304 17 : ABI_SFREE(self%rmesh_int)
3305 :
3306 17 : end subroutine orbmag_free
3307 : !!***
3308 :
3309 : !!****f* ABINIT/dterm_aij
3310 : !! NAME
3311 : !! dterm_aij
3312 : !!
3313 : !! FUNCTION
3314 : !! transfer paw_ij to dterm%aij in more convenient format
3315 : !!
3316 : !! INPUTS
3317 : !! atindx(natom)=index table for atoms (see gstate.f)
3318 : !! dtset <type(dataset_type)>=all input variables for this dataset
3319 : !! paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
3320 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
3321 : !!
3322 : !! OUTPUT
3323 : !!
3324 : !! SIDE EFFECTS
3325 : !! dterm <type(dterm_type)> data related to onsite interactions
3326 : !!
3327 : !! SOURCE
3328 :
3329 16 : subroutine dterm_aij(atindx,dterm,dtset,paw_ij,pawtab)
3330 :
3331 : !Arguments ------------------------------------
3332 : !scalars
3333 : type(dterm_type),intent(inout) :: dterm
3334 : type(dataset_type),intent(in) :: dtset
3335 :
3336 : !arrays
3337 : integer,intent(in) :: atindx(dtset%natom)
3338 : type(paw_ij_type),intent(inout) :: paw_ij(dtset%natom)
3339 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
3340 :
3341 : !Local variables -------------------------
3342 : !scalars
3343 : integer :: iat,iatom,idij,itypat,klmn
3344 : !arrays
3345 : !--------------------------------------------------------------------
3346 :
3347 3033 : dterm%aij = czero
3348 :
3349 : ! note that paw_ij has atoms ordered by input, while
3350 : ! we would like to order by groups of atoms with atindx
3351 44 : do iat=1,dtset%natom
3352 28 : iatom=atindx(iat)
3353 28 : itypat=dtset%typat(iat)
3354 1052 : do klmn=1,pawtab(itypat)%lmn2_size
3355 2980 : do idij = 1, dterm%ndij
3356 2952 : if (paw_ij(iat)%cplex_dij .EQ. 2) then
3357 : dterm%aij(iatom,klmn,idij) = &
3358 1944 : & CMPLX(paw_ij(iat)%dij(2*klmn-1,idij),paw_ij(iat)%dij(2*klmn,idij))
3359 : else
3360 0 : dterm%aij(iatom,klmn,idij) = CMPLX(paw_ij(iat)%dij(klmn,idij),zero)
3361 : end if
3362 : end do ! idij
3363 : end do ! klmn
3364 : end do ! iat
3365 :
3366 16 : dterm%has_aij = 2
3367 :
3368 16 : end subroutine dterm_aij
3369 :
3370 : !!***
3371 :
3372 : !!****f* ABINIT/make_d
3373 : !! NAME
3374 : !! make_d
3375 : !!
3376 : !! FUNCTION
3377 : !! this is a driver to compute different onsite terms, in convenient (complex) format
3378 : !!
3379 : !! INPUTS
3380 : !! atindx(natom)=index table for atoms (see gstate.f)
3381 : !! dtset <type(dataset_type)>=all input variables for this dataset
3382 : !! gprimd(3,3)=reciprocal space lattice vectors
3383 : !! paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
3384 : !! pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data
3385 : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
3386 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
3387 : !!
3388 : !! OUTPUT
3389 : !!
3390 : !! SIDE EFFECTS
3391 : !! dterm <type(dterm_type)> data related to onsite interactions
3392 : !!
3393 : !! SOURCE
3394 :
3395 16 : subroutine make_d(atindx,dterm,dtset,gprimd,paw_ij,pawrad,pawtab,psps)
3396 :
3397 : !Arguments ------------------------------------
3398 : !scalars
3399 : type(dterm_type),intent(inout) :: dterm
3400 : type(dataset_type),intent(in) :: dtset
3401 : type(pseudopotential_type), intent(in) :: psps
3402 :
3403 : !arrays
3404 : integer,intent(in) :: atindx(dtset%natom)
3405 : real(dp),intent(in) :: gprimd(3,3)
3406 : type(paw_ij_type),intent(inout) :: paw_ij(dtset%natom)
3407 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
3408 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
3409 :
3410 : !Local variables -------------------------
3411 : !scalars
3412 : integer :: my_lmax,ngnt
3413 :
3414 : !arrays
3415 : integer,allocatable :: gntselect(:,:)
3416 : real(dp),allocatable :: realgnt(:)
3417 : !--------------------------------------------------------------------
3418 :
3419 : ! make Gaunt integrals
3420 16 : my_lmax = psps%mpsang + 1
3421 48 : ABI_MALLOC(realgnt,((2*my_lmax-1)**2*(my_lmax)**4))
3422 64 : ABI_MALLOC(gntselect,((2*my_lmax-1)**2,my_lmax**2*(my_lmax**2+1)/2))
3423 16 : call realgaunt(my_lmax,ngnt,gntselect,realgnt)
3424 :
3425 : ! generate d terms
3426 :
3427 : ! normal PAW sij overlap, in complex form because it's convenient
3428 16 : call dterm_qij(atindx,dterm,dtset,pawtab)
3429 :
3430 : ! onsite angular momentum expectation values
3431 16 : call dterm_LR(atindx,dterm,dtset,gprimd,pawrad,pawtab)
3432 :
3433 : ! onsite <A_0.A_n> interaction between magnetic field and nuclear dipole
3434 : !call dterm_BM(atindx,dterm,dtset,gntselect,gprimd,my_lmax,pawrad,pawtab,realgnt)
3435 16 : call dterm_BM(atindx,dterm,dtset,gntselect,gprimd,my_lmax,pawrad,pawtab,realgnt)
3436 :
3437 : ! transfers paw_ij to dterm%aij because it's convenient
3438 16 : call dterm_aij(atindx,dterm,dtset,paw_ij,pawtab)
3439 :
3440 16 : ABI_FREE(realgnt)
3441 16 : ABI_FREE(gntselect)
3442 :
3443 16 : end subroutine make_d
3444 : !!***
3445 :
3446 : !!****f* ABINIT/me_proj_mesh
3447 : !! NAME
3448 : !! me_proj_mesh
3449 : !!
3450 : !! FUNCTION
3451 : !! compute realspace mesh for matrix element
3452 : !!
3453 : !! INPUTS
3454 : !!
3455 : !! OUTPUT
3456 : !!
3457 : !! SIDE EFFECTS
3458 : !!
3459 : !! CHILDREN
3460 : !!
3461 : !! SOURCE
3462 :
3463 84 : subroutine me_proj_mesh(bra,fofr,gs_hamk,ket,mpi_enreg,n4,n5,n6,ndat,npw_k,scalar_factor)
3464 :
3465 : !Arguments ------------------------------------
3466 : !scalars
3467 : integer,intent(in) :: n4,n5,n6,ndat,npw_k
3468 : complex(dp),intent(in) :: scalar_factor
3469 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
3470 : type(MPI_type), intent(inout) :: mpi_enreg
3471 : !arrays
3472 : real(dp),intent(out) :: fofr(2,n4,n5,n6*ndat)
3473 : real(dp),intent(in),pointer :: bra(:,:),ket(:,:)
3474 :
3475 : !Local variables -------------------------
3476 : !scalars
3477 : integer :: fourwf_cplex,fourwf_option,tim_fourwf
3478 : real(dp) :: density_norm,weight_i,weight_r
3479 : !arrays
3480 84 : real(dp),allocatable :: fbra(:,:,:,:),fket(:,:,:,:),denpot(:,:),fofgout(:,:)
3481 :
3482 : !--------------------------------------------------------------------
3483 :
3484 420 : ABI_MALLOC(fket,(2,n4,n5,n6*ndat))
3485 336 : ABI_MALLOC(fbra,(2,n4,n5,n6*ndat))
3486 84 : fourwf_cplex = 1
3487 84 : fourwf_option = 0
3488 84 : tim_fourwf = 1
3489 : call fourwf(fourwf_cplex,denpot,ket,fofgout,fket,gs_hamk%gbound_k,&
3490 : & gs_hamk%gbound_k,gs_hamk%istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,&
3491 : & gs_hamk%mgfft,mpi_enreg,ndat,gs_hamk%ngfft,npw_k,npw_k,&
3492 84 : & n4,n5,n6,fourwf_option,tim_fourwf,weight_r,weight_i)
3493 : call fourwf(fourwf_cplex,denpot,bra,fofgout,fbra,gs_hamk%gbound_k,&
3494 : & gs_hamk%gbound_k,gs_hamk%istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,&
3495 : & gs_hamk%mgfft,mpi_enreg,ndat,gs_hamk%ngfft,npw_k,npw_k,&
3496 84 : & n4,n5,n6,fourwf_option,tim_fourwf,weight_r,weight_i)
3497 :
3498 : ! now form scalar*conj(bra)*ket
3499 : fofr(1,:,:,:) = &
3500 : & REAL(scalar_factor)*( fbra(1,:,:,:)*fket(1,:,:,:)+fbra(2,:,:,:)*fket(2,:,:,:)) - &
3501 7828464 : & AIMAG(scalar_factor)*(-fbra(2,:,:,:)*fket(1,:,:,:)+fbra(1,:,:,:)*fket(2,:,:,:))
3502 : fofr(2,:,:,:) = &
3503 : & AIMAG(scalar_factor)*( fbra(1,:,:,:)*fket(1,:,:,:)+fbra(2,:,:,:)*fket(2,:,:,:)) + &
3504 7828464 : & REAL(scalar_factor)*(-fbra(2,:,:,:)*fket(1,:,:,:)+fbra(1,:,:,:)*fket(2,:,:,:))
3505 :
3506 84 : density_norm=gs_hamk%ucvol/gs_hamk%nfft
3507 23137464 : fofr = fofr*density_norm
3508 :
3509 84 : ABI_SFREE(fket)
3510 84 : ABI_SFREE(fbra)
3511 :
3512 84 : end subroutine me_proj_mesh
3513 : !!***
3514 :
3515 : !!****f* ABINIT/local_fermie
3516 : !! NAME
3517 : !! local_fermie
3518 : !!
3519 : !! FUNCTION
3520 : !! estimate Fermi energy as max value of all occupied bands/kpts
3521 : !!
3522 : !! INPUTS
3523 : !! dtset <type(dataset_type)>=all input variables for this dataset
3524 : !! eigen0(dtset%mband*dtset%nkpt*dtset%nsppol)=ground state eigenvalues at each band and kpt
3525 : !! mpi_enreg<type(MPI_type)>=information about MPI parallelization
3526 : !! occ(dtset%mband*dtset%nkpt*dtset%nsppol)=occup number for each band (often 2) at each k point
3527 : !!
3528 : !! OUTPUT
3529 : !! fermie=maximum energy (real(dp)) found over all occupied input bands
3530 : !!
3531 : !! CHILDREN
3532 : !!
3533 : !! SOURCE
3534 :
3535 17 : subroutine local_fermie(dtset,ebands_k,fermie,mpi_enreg)
3536 :
3537 : !Arguments ------------------------------------
3538 : !scalars
3539 : real(dp),intent(out) :: fermie
3540 : type(dataset_type),intent(in) :: dtset
3541 : type(ebands_t) :: ebands_k
3542 : type(MPI_type), intent(inout) :: mpi_enreg
3543 :
3544 : !arrays
3545 :
3546 : !Local variables -------------------------
3547 : !scalars
3548 : integer :: ierr,ikpt,isppol,me
3549 : integer :: nband_k,nn,nproc,spaceComm
3550 : real(dp) :: fermie_proc
3551 :
3552 : !arrays
3553 17 : real(dp),allocatable :: eig_k(:),occ_k(:)
3554 :
3555 : !--------------------------------------------------------------------
3556 :
3557 17 : spaceComm=mpi_enreg%comm_cell
3558 17 : nproc=xmpi_comm_size(spaceComm)
3559 17 : me = mpi_enreg%me_kpt
3560 :
3561 17 : fermie_proc = -1.0D99
3562 35 : do isppol = 1, dtset%nsppol
3563 742 : do ikpt = 1, dtset%nkpt
3564 :
3565 707 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
3566 :
3567 : ! if the current kpt is not on the current processor, cycle
3568 707 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
3569 :
3570 2121 : ABI_MALLOC(occ_k,(nband_k))
3571 5071 : occ_k(:)=ebands_k%occ(1:nband_k,ikpt,isppol)
3572 :
3573 1414 : ABI_MALLOC(eig_k,(nband_k))
3574 5071 : eig_k(:)=ebands_k%eig(1:nband_k,ikpt,isppol)
3575 :
3576 5071 : do nn = 1, nband_k
3577 5071 : if ( (abs(occ_k(nn)).GT.tol8) .AND. (eig_k(nn).GT.fermie_proc) ) then
3578 182 : fermie_proc = eig_k(nn)
3579 : end if
3580 : end do ! nn
3581 :
3582 707 : ABI_FREE(occ_k)
3583 725 : ABI_FREE(eig_k)
3584 :
3585 : end do ! end loop over kpts
3586 : end do ! end loop over isppol
3587 :
3588 17 : call xmpi_max(fermie_proc,fermie,spaceComm,ierr)
3589 :
3590 17 : end subroutine local_fermie
3591 : !!***
3592 :
3593 : !!****f* m_orbmag/orbmag_ncwrite
3594 : !! NAME
3595 : !! orbmag_ncwrite
3596 : !!
3597 : !! FUNCTION
3598 : !! Write orbmag_mesh contributions to netcdf file.
3599 : !!
3600 : !! INPUTS
3601 : !! crystal<crystal_t>=Object defining the unit cell and its symmetries.
3602 : !! dtset<dtset_type>=Dataset type
3603 : !! ebands<ebands_t>=Band structure data.
3604 : !! hdr<hdr_t>=Abinit header
3605 : !! orbmag_mesh<orbmag_mesh_type>=orbmag_mesh data
3606 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
3607 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
3608 : !! ncid=NC file handle.
3609 : !!
3610 : !! OUTPUT
3611 : !! Only writing
3612 : !!
3613 : !! SOURCE
3614 :
3615 34 : subroutine orbmag_ncwrite(crystal,dtset,ebands,hdr,ncid,orbmag_mesh)
3616 :
3617 : !Arguments ------------------------------------
3618 : !scalars
3619 : integer,intent(in) :: ncid
3620 : type(crystal_t),intent(in) :: crystal
3621 : type(dataset_type),intent(in) :: dtset
3622 : type(ebands_t),intent(in) :: ebands
3623 : type(hdr_type),intent(in) :: hdr
3624 : type(orbmag_mesh_type),intent(in) :: orbmag_mesh
3625 : !arrays
3626 :
3627 : !Local variables-------------------------------
3628 : !scalars
3629 : integer :: ncerr,fform
3630 : real(dp) :: cpu,wall,gflops
3631 : logical :: has_ormesh
3632 : character(len=500) :: msg
3633 : !arrays
3634 : !*************************************************************************
3635 :
3636 17 : call cwtime(cpu, wall, gflops, "start")
3637 :
3638 17 : has_ormesh = (dtset%orbmag .EQ. 4)
3639 :
3640 17 : fform = fform_from_ext("ORBMAG.nc")
3641 17 : ABI_CHECK(fform /= 0, "Cannot find fform associated to ORBMAG.nc")
3642 :
3643 : ! Write header, crystal structure and band energies.
3644 17 : NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.))
3645 17 : NCF_CHECK(crystal%ncwrite(ncid))
3646 17 : NCF_CHECK(ebands%ncwrite(ncid))
3647 :
3648 : !! Add orbmag-mesh-specific quantities
3649 : ncerr = nctk_def_dims(ncid, [ &
3650 : nctkdim_t("ntypat", dtset%ntypat),&
3651 : nctkdim_t("mband", orbmag_mesh%mband),&
3652 : nctkdim_t("nkpt", orbmag_mesh%nkpt),&
3653 : nctkdim_t("nsppol", orbmag_mesh%nsppol),&
3654 : nctkdim_t("chern_nterms", orbmag_mesh%chern_nterms),&
3655 : nctkdim_t("orbmag_nterms", orbmag_mesh%orbmag_nterms),&
3656 : nctkdim_t("ndir",3),&
3657 153 : nctkdim_t("natom",dtset%natom)],defmode=.True.)
3658 17 : NCF_CHECK(ncerr)
3659 :
3660 : !! add orbmag_rmesh_cplex,n4,n5,n6 only if orbmag_rmesh will be output
3661 17 : if (has_ormesh) then
3662 : ncerr = nctk_def_dims(ncid, [ &
3663 : nctkdim_t("n4", orbmag_mesh%n4),&
3664 : nctkdim_t("n5", orbmag_mesh%n5),&
3665 4 : nctkdim_t("n6", orbmag_mesh%n6)],defmode=.True.)
3666 1 : NCF_CHECK(ncerr)
3667 : endif
3668 :
3669 : ncerr = nctk_def_arrays(ncid, [&
3670 : nctkarr_t("chern_mesh", "dp", "mband, nkpt, nsppol, ndir, chern_nterms"),&
3671 : nctkarr_t("orbmag_mesh", "dp", "mband, nkpt, nsppol, ndir, orbmag_nterms"),&
3672 : nctkarr_t("lambsig", "dp", "ntypat"),&
3673 85 : nctkarr_t("nucdipmom", "dp", "ndir, natom")])
3674 17 : NCF_CHECK(ncerr)
3675 :
3676 : !! orbmag_rmesh dimensions, only if output
3677 17 : if (has_ormesh) then
3678 : ncerr = nctk_def_arrays(ncid, [&
3679 : nctkarr_t("orbmag_rmesh", "dp", "n4,n5,n6,ndir,orbmag_nterms"),&
3680 3 : nctkarr_t("orbmag_rmesh_int","dp","ndir,orbmag_nterms")])
3681 1 : NCF_CHECK(ncerr)
3682 : endif
3683 :
3684 :
3685 17 : NCF_CHECK(nctk_set_datamode(ncid))
3686 :
3687 17 : NCF_CHECK(nf90_put_var(ncid, vid("chern_mesh"), orbmag_mesh%cmesh))
3688 17 : NCF_CHECK(nf90_put_var(ncid, vid("orbmag_mesh"), orbmag_mesh%omesh))
3689 17 : NCF_CHECK(nf90_put_var(ncid, vid("nucdipmom"), orbmag_mesh%nucdipmom))
3690 17 : NCF_CHECK(nf90_put_var(ncid, vid("lambsig"), orbmag_mesh%lambsig))
3691 17 : if ( has_ormesh ) then
3692 1 : NCF_CHECK(nf90_put_var(ncid, vid("orbmag_rmesh"), orbmag_mesh%rmesh))
3693 1 : NCF_CHECK(nf90_put_var(ncid, vid("orbmag_rmesh_int"), orbmag_mesh%rmesh_int))
3694 : end if
3695 :
3696 17 : call cwtime(cpu,wall,gflops,"stop")
3697 17 : write(msg,'(2(a,f8.2),a)')" orbmag_ncwrite: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
3698 17 : call wrtout(std_out,msg,"PERS")
3699 :
3700 : contains
3701 70 : integer function vid(vname)
3702 : character(len=*),intent(in) :: vname
3703 70 : vid = nctk_idname(ncid, vname)
3704 : end function vid
3705 :
3706 : end subroutine orbmag_ncwrite
3707 : !!***
3708 :
3709 0 : end module m_orbmag
|
