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