Line data Source code
1 : !{\src2tex{textfont=tt}}
2 : !!****m* ABINIT/m_ddb_magpen
3 : !! NAME
4 : !! m_ddb_magpen
5 : !!
6 : !! FUNCTION
7 : !! Convert second- and -third (Berry curvatures) order total energy derivatives
8 : !! calculated with the magnetic penalty (constrained-B functional) into the
9 : !! corresponding quantities of different magnetic functionals:
10 : !! --constrained-M (fixed-spin)
11 : !! --constrained-H (relaxed spin)
12 : !! Calculate and write the ensuing clamped-ion susceptibilities. Lattice-mediated
13 : !! contributions are incorporated in m_ddb_omega_interpol.
14 : !!
15 : !! COPYRIGHT
16 : !! Copyright (C) 2023 ABINIT group (MR and MS)
17 : !! This file is distributed under the terms of the
18 : !! GNU General Public License, see ~abinit/COPYING
19 : !! or http://www.gnu.org/copyleft/gpl.txt .
20 : !!
21 : !! NOTES
22 : !!
23 : !! PARENTS
24 : !!
25 : !! CHILDREN
26 : !!
27 : !! SOURCE
28 :
29 : #if defined HAVE_CONFIG_H
30 : #include "config.h"
31 : #endif
32 :
33 : #include "abi_common.h"
34 :
35 : module m_ddb_magpen
36 :
37 : use defs_basis
38 : use m_abicore
39 : use m_profiling_abi
40 : use m_errors
41 : use m_ddb
42 : use m_fstrings, only : itoa, sjoin
43 : use m_macroave, only : POLINT
44 : use m_io_tools, only : open_file
45 : use m_cgtools, only : fxphas_seq
46 : use m_dynmat, only : pheigvec_normalize
47 : use m_numeric_tools, only : polcoe
48 :
49 : implicit none
50 :
51 : public :: ddb_magpen ! Convert the derivatives calculated with the magnetic penalty into the physically relevant ones.
52 : public :: local_spinsus ! Treat local spin susceptibility (2nd-order magnetic derivatives)
53 : public :: magmom ! Treat first-order magnetic moments (2nd-order mixed derivatives)
54 : public :: mp_d2etot ! Treat 2nd-order nonmagnetic derivatives
55 : public :: asrw0 ! Apply the ASR correction calculated at w=0 at any value of w
56 :
57 : private
58 :
59 : ! *************************************************************************
60 :
61 : contains
62 : !!***
63 :
64 : !!****f* m_ddb_magpen/ddb_magpen
65 : !! NAME
66 : !! ddb_magpen
67 : !!
68 : !! FUNCTION
69 : !! Convert the second- and possibly third-order derivatives calculated with the magnetic penalty
70 : !! into physically relevant quantities.
71 : !!
72 : !! INPUTS
73 : !! ddb (INOUT) = ddb block datastructure
74 : !! ddb_lw (INOUT) = ddb_lw block datastructure
75 : !! magpen = amplitude (in Ha) of the applied magnetic penalty
76 : !! mpatpol(2) = Atoms on which the magnetic penalty has been applied
77 : !! mpdir(3) = Directions along which the spin-degrees of freedom have been stiffened
78 : !! mpert = maximum number of ipert
79 : !! mpopt = 1 calculate the frozen-magnetic second-order quantities
80 : !! 2 calculate the spin-relaxed second-order quantities
81 : !! natom= number of atoms in unit cell
82 : !! ntypat= number of atom types
83 : !! timdisp= 1 calculate the third-order Berry curvatures
84 : !! ucvol= unit cell volume
85 : !!
86 : !! OUTPUT
87 : !! ddb= ddb%val updated with the corrected second-order derivatives
88 : !!
89 : !! SOURCE
90 :
91 3 : subroutine ddb_magpen(ddb,ddb_lw,&
92 : & magpen,mpatpol,mpdir,mpert,mpopt,natom, &
93 3 : & prtvol,rftyp,ucvol,timdisp,xred)
94 :
95 : !Arguments -------------------------------
96 : !scalars
97 : integer,intent(in) :: mpert,mpopt,natom,prtvol,rftyp,timdisp
98 : real(dp),intent(in) :: magpen,ucvol
99 : !arrays
100 : type(ddb_type),intent(inout) :: ddb,ddb_lw
101 : integer,intent(in) :: mpatpol(2),mpdir(3)
102 : real(dp),intent(in) :: xred(3,natom)
103 :
104 : !Local variables -------------------------
105 : !scalars
106 : integer :: iblok,ii,jblok,kblok,lblok,nblok,ndim
107 : integer :: nmat,nmdir,optgb,prtopt
108 : character(len=500) :: msg
109 : logical :: qeq0
110 : !arrays
111 : integer :: rfelfd(4),rfphon(4),rfstrs(4),rfmagn(4),rffreq(4)
112 : real(dp) :: omega(3),qphnrm(3),qphon(3,3),qpt(3)
113 3 : complex(dpc), allocatable :: barmagsus(:,:),invbarmagsus(:,:)
114 3 : complex(dpc), allocatable :: invmagsus(:,:), magsus(:,:), invhmat(:,:)
115 3 : complex(dpc), allocatable :: barmmom(:,:),barmmom_tr(:,:),mmom(:,:),mmom_tr(:,:)
116 3 : complex(dpc), allocatable :: zfield(:,:), zfield_tr(:,:)
117 3 : complex(dpc), allocatable :: bc_barmagsus(:,:),bc_ss(:,:),bc_sp(:,:)
118 3 : complex(dpc), allocatable :: ifcmat(:,:),ifcmat_fm(:,:),zeff(:,:),zeff_tr(:,:)
119 3 : complex(dpc), allocatable :: fmzeff(:,:),fmzeff_tr(:,:)
120 3 : complex(dpc), allocatable :: lm_epsilon(:,:),dum_phongreen(:,:)
121 :
122 : ! *********************************************************************
123 246 : write(msg, '(2a,(80a),6a)' ) ch10,('=',ii=1,80),ch10,ch10,&
124 3 : ' Constrained DFPT section ',ch10,ch10,&
125 249 : ' (Quantities are in atomic units and along Cartesian directions)'
126 9 : call wrtout([std_out, ab_out], msg)
127 :
128 3 : prtopt=1
129 3 : if (magpen<zero) then
130 : nmat= 1
131 3 : else if (magpen>zero) then
132 3 : nmat= mpatpol(2) - mpatpol(1) + 1
133 : end if
134 12 : nmdir=sum(mpdir(:))
135 3 : ndim=nmat*nmdir
136 3 : optgb=1
137 12 : ABI_MALLOC(barmagsus,(ndim,ndim))
138 9 : ABI_MALLOC(magsus,(ndim,ndim))
139 9 : ABI_MALLOC(invbarmagsus,(ndim,ndim))
140 9 : ABI_MALLOC(invmagsus,(ndim,ndim))
141 9 : ABI_MALLOC(invhmat,(ndim,ndim))
142 12 : ABI_MALLOC(barmmom,(ndim,(natom+5)*3))
143 9 : ABI_MALLOC(barmmom_tr,((natom+5)*3,ndim))
144 9 : ABI_MALLOC(mmom,(ndim,(natom+5)*3))
145 9 : ABI_MALLOC(mmom_tr,((natom+5)*3,ndim))
146 9 : ABI_MALLOC(zfield,(ndim,(natom+5)*3))
147 9 : ABI_MALLOC(zfield_tr,((natom+5)*3,ndim))
148 :
149 12 : ABI_MALLOC(ddb%val_fs,(2,ddb%msize,ddb%nblok))
150 : ! if (mpopt==2) ABI_MALLOC(ddb%val_rs,(2,ddb%msize,ddb%nblok))
151 9 : ABI_MALLOC(ddb%val_rs,(2,ddb%msize,ddb%nblok))
152 :
153 15 : ABI_MALLOC(ifcmat,(3*natom,3*natom))
154 9 : ABI_MALLOC(ifcmat_fm,(3*natom,3*natom))
155 9 : ABI_MALLOC(fmzeff,(3,3*natom))
156 9 : ABI_MALLOC(fmzeff_tr,(3*natom,3))
157 6 : ABI_MALLOC(zeff,(3,3*natom))
158 6 : ABI_MALLOC(zeff_tr,(3*natom,3))
159 3 : ABI_MALLOC(lm_epsilon,(3,3))
160 9 : ABI_MALLOC(dum_phongreen,(3*natom,3*natom))
161 :
162 3 : nblok=ddb%nblok
163 7 : do kblok=1,nblok
164 :
165 : ! Look for the local spin-susceptibility block in the DDB
166 4 : omega=zero
167 4 : qphon=zero
168 16 : qphon(:,1)=ddb%qpt(1:3,kblok)
169 16 : qpt=qphon(:,1)
170 16 : qeq0=(sqrt(sum(qphon(:,1)**2))<tol8)
171 16 : qphnrm(:)=ddb%nrm(1,kblok)
172 4 : omega(1)=ddb%omega(1,kblok)
173 12 : rfphon(1:2)=0
174 12 : rfelfd(1:2)=0
175 12 : rfstrs(1:2)=0
176 4 : if (magpen<zero) then
177 0 : rfmagn(1:2)= 1
178 4 : else if (magpen>zero) then
179 12 : rfmagn(1:2)= 2
180 : end if
181 :
182 328 : write(msg, '(1a,(80a),2a,3f16.8,2a,f16.8,a)' ) ch10,('-',ii=1,80),ch10, &
183 4 : ' q point ', qphon(:,1),ch10,&
184 332 : ' frequency', ddb%omega(1,kblok), ch10
185 12 : call wrtout([std_out, ab_out], msg)
186 :
187 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, &
188 4 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn)
189 :
190 : ! Calculate and write the spin-susceptibility matrices
191 4 : if (iblok /= 0) then
192 4 : if (prtvol>1) then
193 4 : if (magpen<zero) then
194 0 : write(msg, '(2a)' ) ' Spin susceptibility (Uniform Zeeman) ',ch10
195 0 : call wrtout([std_out, ab_out], msg)
196 4 : else if (magpen>zero) then
197 4 : write(msg, '(2a)' ) ' Spin susceptibility (Local Zeeman) ',ch10
198 12 : call wrtout([std_out, ab_out], msg)
199 : end if
200 : end if
201 :
202 : call local_spinsus(barmagsus,ddb,iblok,invbarmagsus,invmagsus,invhmat,magpen,magsus,&
203 4 : & mpatpol,mpdir,mpert,natom,ndim,nmdir,prtopt,prtvol,qpt,xred)
204 :
205 : end if
206 :
207 : ! Calculate and write the induced magnetic moments
208 4 : if (prtvol>1) then
209 328 : write(msg, '(2a,(80a),4a)' ) ch10,('-',ii=1,80),ch10,ch10,&
210 332 : ' First-order magnetic moments ',ch10
211 12 : call wrtout([std_out, ab_out], msg)
212 : end if
213 :
214 : ! First atomic-displacement
215 4 : rfphon(2)=1
216 12 : rfelfd(1:2)=0
217 12 : rfstrs(1:2)=0
218 4 : rfmagn(:)=0
219 4 : if (magpen<zero) then
220 0 : rfmagn(1)= 1
221 4 : else if (magpen>zero) then
222 4 : rfmagn(1)= 2
223 : end if
224 :
225 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, &
226 4 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn)
227 :
228 : ! Then electric field
229 : ! Look for the induced magnetic moments block in the DDB
230 4 : jblok=0
231 4 : if (qeq0) then
232 3 : rfphon(2)=0
233 3 : rfelfd(2)=2
234 9 : rfstrs(1:2)=0
235 3 : rfmagn(:)=0
236 3 : if (magpen<zero) then
237 0 : rfmagn(1)= 1
238 3 : else if (magpen>zero) then
239 3 : rfmagn(1)= 2
240 : end if
241 :
242 : call ddb%get_block(jblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, &
243 3 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn)
244 : end if
245 :
246 : ! Then macroscopic Zeeman field
247 : ! Look for the induced magnetic moments block in the DDB
248 : lblok=0
249 4 : rfphon(:)=0
250 4 : rfelfd(:)=0
251 12 : rfstrs(1:2)=0
252 4 : rfmagn(2)=1
253 4 : if (magpen<zero) then
254 0 : rfmagn(1)= 1
255 4 : else if (magpen>zero) then
256 4 : rfmagn(1)= 2
257 : end if
258 :
259 : call ddb%get_block(lblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, &
260 4 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn)
261 :
262 4 : if (iblok /= 0 .or. jblok /=0 .or. lblok/=0) then
263 : call magmom(barmmom,barmmom_tr,ddb,invbarmagsus,invhmat,iblok,jblok,lblok,magpen,magsus,mmom,mmom_tr,&
264 4 : & mpatpol,mpdir,mpert,natom,ndim,nmdir,prtopt,prtvol,qpt,xred,zfield,zfield_tr)
265 : end if
266 :
267 : !Now calculate the non-magnetic second-order quantities
268 328 : write(msg, '(2a,(80a),4a)' ) ch10,('-',ii=1,80),ch10,ch10,&
269 332 : ' Second-order linear-response tensors ',ch10
270 12 : call wrtout([std_out, ab_out], msg)
271 :
272 : !Convert ddb%val to second-order energies
273 4 : call ddb%to_d2etot(ddb%val,kblok,0,qeq0,qphon,qphnrm,ucvol,optgb,omega=omega)
274 :
275 : !Convert second-order derivatives to diferent magnetic boundary conditions
276 : call mp_d2etot(barmagsus,ddb,kblok,magsus,mpert,mpopt,natom, &
277 4 : & ndim,zfield,zfield_tr)
278 :
279 : !Convert second-order energies to the physical quantities of ddb%val
280 4 : call ddb%to_d2etot(ddb%val,kblok,1,qeq0,qphon,qphnrm,ucvol,optgb,omega=omega)
281 4 : call ddb%to_d2etot(ddb%val_fs,kblok,1,qeq0,qphon,qphnrm,ucvol,optgb,omega=omega)
282 4 : if (mpopt==2) call ddb%to_d2etot(ddb%val_rs,kblok,1,qeq0,qphon,qphnrm,ucvol,optgb,omega=omega)
283 :
284 : !Print the physical quantities in the new magnetic boundary conditions
285 : if (prtopt==1) then
286 4 : if (mpopt==1) then
287 0 : call mp_d2etot_print(ddb,ddb%val_fs,kblok,mpert,natom,nblok,1,omega,qeq0,qphnrm,qphon,ddb%msize)
288 4 : else if (mpopt==2) then
289 4 : call mp_d2etot_print(ddb,ddb%val_rs,kblok,mpert,natom,nblok,2,omega,qeq0,qphnrm,qphon,ddb%msize)
290 : end if
291 : end if
292 :
293 4 : rfmagn(:)=0
294 4 : rfelfd(:)=0
295 7 : rfphon(:)=0
296 :
297 : end do
298 :
299 : ! BERRY CURVATURES
300 3 : if (timdisp==1) then
301 :
302 8 : ABI_MALLOC(ddb_lw%val_fs,(2,ddb_lw%msize,ddb_lw%nblok))
303 :
304 6 : ABI_MALLOC(bc_barmagsus,(ndim,ndim))
305 6 : ABI_MALLOC(bc_ss,(ndim,ndim))
306 6 : ABI_MALLOC(bc_sp,(ndim,(natom+5)*3))
307 :
308 164 : write(msg, '(2a,(80a),4a)' ) ch10,('=',ii=1,80),ch10,ch10,&
309 166 : ' Frequency-derivatives (Berry curvatures) constrained DFPT section ',ch10
310 6 : call wrtout([std_out, ab_out], msg)
311 :
312 2 : rffreq(:)=0
313 2 : nblok=ddb_lw%nblok
314 6 : do kblok=1,nblok
315 :
316 4 : if (ddb_lw%typ(kblok)/=33) cycle
317 :
318 : !Berry curvature of the penalized spin-susceptibility
319 2 : qphon=zero
320 8 : qphon(:,1)=ddb_lw%qpt(1:3,kblok)
321 8 : qpt= qphon(:,1)
322 8 : qphnrm(:)=ddb_lw%nrm(1,kblok)
323 8 : omega(:)=ddb_lw%omega(:,kblok)
324 8 : rfphon(1:3)=0
325 8 : rfelfd(1:3)=0
326 8 : rfstrs(1:3)=0
327 2 : rffreq(3)=1
328 2 : iblok=0
329 2 : if (magpen<zero) then
330 0 : rfmagn(1:2)= 1
331 2 : else if (magpen>zero) then
332 6 : rfmagn(1:2)= 2
333 : end if
334 :
335 164 : write(msg, '(a,(80a),a,3(a,3f16.8,a),3(a,f16.8,a))' ) ch10,('-',ii=1,80),ch10, &
336 2 : ' q point 1 ', qphon(:,1),ch10,&
337 2 : ' q point 2 ', qphon(:,2),ch10,&
338 2 : ' q point 3 ', qphon(:,3),ch10,&
339 2 : ' frequency 1', ddb_lw%omega(1,kblok), ch10,&
340 2 : ' frequency 2', ddb_lw%omega(2,kblok), ch10,&
341 166 : ' frequency 3', ddb_lw%omega(3,kblok), ch10
342 6 : call wrtout([std_out, ab_out], msg)
343 :
344 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, 33, &
345 2 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn,rffreq=rffreq)
346 :
347 2 : if (iblok /= 0) then
348 : call berrycurv_ss(bc_barmagsus,bc_ss,ddb_lw,iblok,invbarmagsus,mpatpol,mpdir,mpert,&
349 2 : & natom,ndim,nmdir,prtvol,qpt,xred)
350 : end if
351 :
352 : !Berry curvature of the induced Zeeman fields
353 :
354 : !First atomic-displacement
355 : iblok=0
356 2 : rfphon(2)=1
357 2 : rfmagn(:)=0
358 2 : if (magpen<zero) then
359 0 : rfmagn(1)= 1
360 2 : else if (magpen>zero) then
361 2 : rfmagn(1)= 2
362 : end if
363 :
364 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, 33, &
365 2 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn,rffreq=rffreq)
366 :
367 : ! Then electric field
368 2 : jblok=0
369 2 : if (qeq0) then
370 1 : rfphon(2)=0
371 1 : rfelfd(2)=2
372 :
373 : call ddb_lw%get_block(jblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, 33, &
374 1 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn,rffreq=rffreq)
375 : end if
376 :
377 : ! Then macroscopic Zeeman field
378 2 : lblok=0
379 2 : if (qeq0) then
380 1 : rfelfd(2)=0
381 1 : rfmagn(2)=1
382 :
383 : call ddb_lw%get_block(lblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, 33, &
384 1 : & mpatpol=mpatpol,mpdir=mpdir,omega=omega,rfmagn=rfmagn,rffreq=rffreq)
385 : end if
386 :
387 2 : if (iblok /= 0 .or. jblok /=0 .or. lblok /= 0) then
388 : call berrycurv_sp(barmmom,bc_sp,bc_ss,ddb_lw,iblok,invbarmagsus,jblok,lblok, &
389 2 : & mpatpol,mpdir,mpert,natom,ndim,nmdir,prtvol,qpt,xred)
390 : end if
391 :
392 : !Berry curvature of other second-order quantites
393 : call berrycurv_pp(barmagsus,bc_barmagsus,bc_sp,ddb_lw,kblok, &
394 2 : & mpert,natom,ndim,qeq0,ucvol,zfield)
395 :
396 : !Print them
397 : call mp_d3etot_print(ddb_lw,ddb_lw%val_fs,kblok,mpert,natom,nblok,1,omega,&
398 6 : & prtvol,qeq0,qphnrm,qphon,ucvol,ddb_lw%msize)
399 :
400 : end do
401 2 : ABI_FREE(bc_ss)
402 2 : ABI_FREE(bc_sp)
403 2 : ABI_FREE(bc_barmagsus)
404 : end if
405 :
406 : !Deallocations
407 3 : ABI_FREE(barmagsus)
408 3 : ABI_FREE(barmmom)
409 3 : ABI_FREE(barmmom_tr)
410 3 : ABI_FREE(magsus)
411 3 : ABI_FREE(invbarmagsus)
412 3 : ABI_FREE(invmagsus)
413 3 : ABI_FREE(mmom)
414 3 : ABI_FREE(mmom_tr)
415 3 : ABI_FREE(zfield)
416 3 : ABI_FREE(zfield_tr)
417 3 : ABI_FREE(ifcmat)
418 3 : ABI_FREE(ifcmat_fm)
419 3 : ABI_FREE(zeff)
420 3 : ABI_FREE(zeff_tr)
421 3 : ABI_FREE(fmzeff)
422 3 : ABI_FREE(fmzeff_tr)
423 3 : ABI_FREE(lm_epsilon)
424 3 : ABI_FREE(dum_phongreen)
425 3 : ABI_FREE(invhmat)
426 :
427 3 : end subroutine ddb_magpen
428 : !!***
429 :
430 : !!****f* m_ddb_magpen/local_spinsus
431 : !! NAME
432 : !! local_spinsus
433 : !!
434 : !! FUNCTION
435 : !! Calculate the spin-susceptibility matrix and its inverse
436 : !!
437 : !! INPUTS
438 : !! ddb= Second-order derivative arrais
439 : !! iblok= index of the current block
440 : !! magpen = amplitude (in Ha) of the applied magnetic penalty
441 : !! mpatpol(2) = Atoms on which the magnetic penalty has been applied
442 : !! mpdir(3) = Directions along which the spin-degrees of freedom have been stiffened
443 : !! mpert =maximum number of ipert
444 : !! natom= number of atoms in unit cell
445 : !! nblok= number of blocks in the DDB
446 : !! nmdir= number of directions along which the magnetic penalty was applied
447 : !! ndim= dimension of the square susceptibilities
448 : !! prtvol= control the volume of information written on output
449 : !! fs2rs= (optional) if 1, the routine starts from a precalculated blkval_fs
450 : !! blkval_fs(2,3,mpert,3,mpert)= fixed-spin 2nd-order derivatives
451 : !! xred(3,natom)= reduced atomic coordinates
452 : !!
453 : !! OUTPUT
454 : !! barmagsus(ndim,ndim)= Penalized spin-sussceptibility tensor
455 : !! invbarmagsus(ndim,ndim)= Inverse of the penalized spin-sussceptibility tensor
456 : !! magsus(ndim,ndim)= Spin-sussceptibility tensor
457 : !! invmagsus(ndim,ndim)= Inverse of spin-sussceptibility tensor
458 : !!
459 : !! SOURCE
460 :
461 68 : subroutine local_spinsus(barmagsus,ddb,iblok,invbarmagsus,invmagsus,invhmat,magpen,magsus, &
462 34 : & mpatpol,mpdir,mpert,natom,ndim,nmdir,prtopt,prtvol,qphon,xred, &
463 30 : & fs2rs,blkval_fs) !optional
464 :
465 : !Arguments -------------------------------
466 : !scalars
467 : integer,intent(in) :: iblok,mpert,natom,ndim,nmdir
468 : integer,intent(in) :: prtopt,prtvol
469 : integer,intent(in),optional :: fs2rs
470 : real(dp),intent(in) :: magpen
471 : !arrays
472 : type(ddb_type),intent(inout) :: ddb
473 : integer,intent(in) :: mpatpol(2),mpdir(3)
474 : real(dp),intent(in),optional :: blkval_fs(2,3,mpert,3,mpert,1)
475 : real(dp),intent(in) :: qphon(3),xred(3,natom)
476 : complex(dpc),intent(out) :: barmagsus(ndim,ndim)
477 : complex(dpc),intent(out) :: invbarmagsus(ndim,ndim)
478 : complex(dpc),intent(out) :: invhmat(ndim,ndim)
479 : complex(dpc),intent(out) :: magsus(ndim,ndim)
480 : complex(dpc),intent(out) :: invmagsus(ndim,ndim)
481 :
482 : !Local variables -------------------------
483 : !scalars
484 : integer :: fs2rs_
485 : integer :: iat1,iat2,icol,idir1,idir2,index,info,ipert1,ipert2,irow,lwork
486 : integer :: ipert1_red,ipert2_red,idir1_red,idir2_red
487 : !real(dp) :: fac
488 : character(len=1000) :: msg
489 : !arrays
490 68 : complex(dpc) :: idty(ndim,ndim)
491 34 : integer :: indexat(ndim),indexdir(ndim)
492 34 : integer, allocatable :: ipiv(:)
493 34 : complex(dpc),allocatable :: work(:),work1(:,:),work2(:,:)
494 : character(len=1) :: cart(3)=(/'x','y','z'/)
495 :
496 : ! *********************************************************************
497 :
498 : !If fixed-spin case has been precalculated do less stuff
499 34 : fs2rs_=0; if (present(fs2rs)) fs2rs_=fs2rs
500 30 : if (fs2rs_==1) then
501 30 : if (.not.present(blkval_fs)) then
502 0 : write(msg, '(3a)' )' No fixed-spin array has been passed to local_spinsus', &
503 0 : & ' but fs2rs=1 ',ch10
504 0 : ABI_ERROR(msg)
505 : end if
506 : end if
507 :
508 : !Extract the penalized susceptibility
509 714 : idty=(zero,zero)
510 34 : ipert2_red= 0
511 102 : do iat2= mpatpol(1), mpatpol(2)
512 68 : ipert2= natom + 11 + iat2
513 68 : ipert2_red= ipert2_red + 1
514 68 : idir2_red= 0
515 306 : do idir2= 1, 3
516 204 : if (mpdir(idir2)==0) cycle
517 136 : idir2_red= idir2_red + 1
518 136 : icol=idir2_red+(ipert2_red-1)*nmdir
519 136 : indexat(icol)=iat2
520 136 : indexdir(icol)=idir2
521 136 : idty(icol,icol)=(one,zero)
522 136 : ipert1_red=0
523 476 : do iat1= mpatpol(1), mpatpol(2)
524 272 : ipert1= natom + 11 + iat1
525 272 : ipert1_red= ipert1_red + 1
526 272 : idir1_red= 0
527 1292 : do idir1= 1, 3
528 816 : if (mpdir(idir1)==0) cycle
529 544 : idir1_red=idir1_red+1
530 544 : irow=idir1_red+(ipert1_red-1)*nmdir
531 544 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
532 :
533 : !Apply a phase factor to adopt the same Gonze&Lee convention as for atomic
534 : !displacements at finite q (See M. Stengel PRB 2013).
535 816 : if (fs2rs_==0) then
536 : barmagsus(irow,icol)= &
537 : & cmplx(ddb%val(1,index,iblok),ddb%val(2,index,iblok),kind=dpc) &
538 256 : & * exp(two_pi*(0.d0,1.d0)*dot_product(qphon,xred(:,iat1)-xred(:,iat2)))
539 480 : else if (fs2rs_==1) then
540 : invmagsus(irow,icol)= &
541 : & cmplx(blkval_fs(1,idir1,ipert1,idir2,ipert2,iblok), &
542 480 : & blkval_fs(2,idir1,ipert1,idir2,ipert2,iblok),kind=dpc)
543 : end if
544 :
545 : end do
546 : end do
547 : end do
548 : end do
549 :
550 : !Use magsus to store the intermediate array
551 714 : magsus=idty-magpen*barmagsus
552 :
553 : !Invert the arrays
554 136 : ABI_MALLOC(work1,(ndim,ndim))
555 102 : ABI_MALLOC(work2,(ndim,ndim))
556 34 : if (fs2rs_==0) then
557 88 : work1=barmagsus
558 30 : else if (fs2rs_==1) then
559 660 : work1=invmagsus
560 : end if
561 748 : work2=magsus
562 :
563 102 : ABI_MALLOC(ipiv,(ndim))
564 34 : call zgetrf( ndim, ndim, work1, ndim, ipiv, info )
565 34 : ABI_CHECK(info == 0, sjoin('zgetrf returned:', itoa(info)))
566 :
567 34 : ABI_MALLOC(work,(2))
568 34 : call zgetri( ndim, work1, ndim, ipiv, work, -1, info )
569 34 : ABI_CHECK(info == 0, sjoin('zgetri returned:', itoa(info)))
570 34 : lwork=int(work(1))
571 :
572 102 : ABI_REMALLOC(work,(lwork))
573 34 : call zgetri( ndim, work1, ndim, ipiv, work, lwork, info )
574 34 : ABI_CHECK(info == 0, sjoin('zgetri returned:', itoa(info)))
575 :
576 34 : if (fs2rs_==1) then
577 630 : magsus=work1
578 30 : ABI_SFREE(ipiv)
579 30 : ABI_SFREE(work1)
580 30 : ABI_SFREE(work2)
581 30 : ABI_SFREE(work)
582 30 : return
583 : end if
584 :
585 4 : call zgetrf( ndim, ndim, work2, ndim, ipiv, info )
586 4 : ABI_CHECK(info == 0, sjoin('zgetrf returned:', itoa(info)))
587 :
588 4 : call zgetri( ndim, work2, ndim, ipiv, work, -1, info )
589 4 : ABI_CHECK(info == 0, sjoin('zgetri returned:', itoa(info)))
590 4 : lwork=int(work(1))
591 :
592 12 : ABI_REMALLOC(work,(lwork))
593 4 : call zgetri( ndim, work2, ndim, ipiv, work, lwork, info )
594 4 : ABI_CHECK(info == 0, sjoin('zgetri returned:', itoa(info)))
595 4 : ABI_FREE(work)
596 :
597 : !Write the results in meaningfull arrays
598 84 : invbarmagsus=work1
599 :
600 : !At last, calculate the susceptibility and its inverse
601 84 : invhmat=work2
602 420 : magsus=matmul(work2,barmagsus)
603 84 : invmagsus=invbarmagsus-magpen*idty
604 :
605 4 : ABI_FREE(ipiv)
606 4 : ABI_FREE(work1)
607 4 : ABI_FREE(work2)
608 :
609 4 : if (prtopt==1.and.prtvol>1) then
610 :
611 : !Write results in output
612 4 : if (magpen > zero) then
613 12 : call wrtout([ab_out,std_out], ' Local spin susceptibility ')
614 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
615 20 : do irow=1, ndim
616 84 : do icol=1, ndim
617 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
618 64 : & indexat(irow), cart(indexdir(irow)), indexat(icol), cart(indexdir(icol)), &
619 128 : & real(magsus(irow,icol)), aimag(magsus(irow,icol))
620 208 : call wrtout([ab_out,std_out], msg)
621 : end do
622 : end do
623 12 : call wrtout([ab_out,std_out], ' ')
624 :
625 12 : call wrtout([ab_out,std_out], ' Inverse of local spin susceptibility ')
626 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
627 20 : do irow=1, ndim
628 84 : do icol=1, ndim
629 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
630 64 : & indexat(irow), cart(indexdir(irow)), indexat(icol), cart(indexdir(icol)), &
631 128 : & real(invmagsus(irow,icol)), aimag(invmagsus(irow,icol))
632 208 : call wrtout([ab_out,std_out], msg)
633 : end do
634 : end do
635 12 : call wrtout([ab_out,std_out], ' ')
636 :
637 4 : if (prtvol > 2) then
638 12 : call wrtout([ab_out,std_out], ' Penalized local spin susceptibility ')
639 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
640 20 : do irow=1, ndim
641 84 : do icol=1, ndim
642 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
643 64 : & indexat(irow), cart(indexdir(irow)), indexat(icol), cart(indexdir(icol)), &
644 128 : & real(barmagsus(irow,icol)), aimag(barmagsus(irow,icol))
645 208 : call wrtout([ab_out,std_out], msg)
646 : end do
647 : end do
648 12 : call wrtout([ab_out,std_out], ' ')
649 :
650 12 : call wrtout([ab_out,std_out], ' Inverse of penalized local spin susceptibility ')
651 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
652 20 : do irow=1, ndim
653 84 : do icol=1, ndim
654 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
655 64 : & indexat(irow), cart(indexdir(irow)), indexat(icol), cart(indexdir(icol)), &
656 128 : & real(invbarmagsus(irow,icol)), aimag(invbarmagsus(irow,icol))
657 208 : call wrtout([ab_out,std_out], msg)
658 : end do
659 : end do
660 12 : call wrtout([ab_out,std_out], ' ')
661 :
662 12 : call wrtout([ab_out,std_out], ' Inverse of H matrix (I-alpha barchi)^{-1} ')
663 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
664 20 : do irow=1, ndim
665 84 : do icol=1, ndim
666 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
667 64 : & indexat(irow), cart(indexdir(irow)), indexat(icol), cart(indexdir(icol)), &
668 128 : & real(invhmat(irow,icol)), aimag(invhmat(irow,icol))
669 208 : call wrtout([ab_out,std_out], msg)
670 : end do
671 : end do
672 12 : call wrtout([ab_out,std_out], ' ')
673 : end if
674 :
675 : end if !magpen>zero
676 :
677 : end if !prtopt
678 :
679 : !Store the FS and RS flavors on the DDB array
680 : ipert2_red= 0
681 12 : do iat2= mpatpol(1), mpatpol(2)
682 8 : ipert2= natom + 11 + iat2
683 8 : ipert2_red= ipert2_red + 1
684 8 : idir2_red= 0
685 36 : do idir2= 1, 3
686 24 : if (mpdir(idir2)==0) cycle
687 16 : idir2_red= idir2_red + 1
688 16 : icol=idir2_red+(ipert2_red-1)*nmdir
689 16 : ipert1_red=0
690 56 : do iat1= mpatpol(1), mpatpol(2)
691 32 : ipert1= natom + 11 + iat1
692 32 : ipert1_red= ipert1_red + 1
693 32 : idir1_red= 0
694 152 : do idir1= 1, 3
695 96 : if (mpdir(idir1)==0) cycle
696 64 : idir1_red=idir1_red+1
697 64 : irow=idir1_red+(ipert1_red-1)*nmdir
698 64 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
699 :
700 64 : ddb%val_fs(1,index,iblok)=real(invmagsus(irow,icol))
701 64 : ddb%val_fs(2,index,iblok)=aimag(invmagsus(irow,icol))
702 :
703 64 : ddb%val_rs(1,index,iblok)=real(magsus(irow,icol))
704 128 : ddb%val_rs(2,index,iblok)=aimag(magsus(irow,icol))
705 :
706 : end do
707 : end do
708 : end do
709 : end do
710 :
711 64 : end subroutine local_spinsus
712 : !!***
713 :
714 : !!****f* m_ddb_magpen/magmom
715 : !! NAME
716 : !! magmom
717 : !!
718 : !! FUNCTION
719 : !! Calculate the first-order magnetic moments and the constrained
720 : !! Zeeman fields
721 : !!
722 : !! INPUTS
723 : !! ddb= Second-order derivative arrais
724 : !! invbarmagsus(ndim,ndim)= Inverse of the penalized spin-sussceptibility tensor
725 : !! iblok= index of the atomic displacement block
726 : !! jblok= index of the electric field block
727 : !! magsus(ndim,ndim)= Spin-sussceptibility tensor
728 : !! magpen = amplitude (in Ha) of the applied magnetic penalty
729 : !! mpatpol(2) = Atoms on which the magnetic penalty has been applied
730 : !! mpdir(3) = Directions along which the spin-degrees of freedom have been stiffened
731 : !! mpert =maximum number of ipert
732 : !! natom= number of atoms in unit cell
733 : !! nblok= number of blocks in the DDB
734 : !! nmdir= number of directions along which the magnetic penalty was applied
735 : !! ndim= dimension of the square susceptibilities
736 : !! prtvol= control the volume of information written on output
737 : !! fs2rs= (optional) if 1, the routine starts from a precalculated blkval_fs
738 : !! blkval_fs(2,3,mpert,3,mpert)= fixed-spin 2nd-order derivatives
739 : !!
740 : !! OUTPUT
741 : !! barmmom(ndim,(natom+2)*3)= penalized first order magnetic moments on the atoms and
742 : !! directions of the penalty induced by atomic displacements and/or electric fields.
743 : !! mmom(ndim,(natom+2)*3)= first order magnetic moments on the atoms and
744 : !! directions of the penalty induced by atomic displacements and/or electric fields.
745 : !! zfield(ndim,(natom+2)*3)= Zeeman fields at constrained magnetic moments.
746 : !!
747 : !! SOURCE
748 :
749 68 : subroutine magmom(barmmom,barmmom_tr,ddb,invbarmagsus,invhmat,iblok,jblok,lblok,magpen,magsus,mmom,mmom_tr,&
750 34 : & mpatpol,mpdir,mpert,natom,ndim,nmdir,prtopt,prtvol,qphon,xred,zfield,zfield_tr, &
751 30 : & fs2rs,blkval_fs) !optional
752 :
753 : !Arguments -------------------------------
754 : !scalars
755 : integer,intent(in) :: iblok,jblok,lblok,mpert,natom,ndim,nmdir,prtopt,prtvol
756 : integer,intent(in),optional :: fs2rs
757 : real(dp),intent(in) :: magpen
758 : !arrays
759 : type(ddb_type),intent(inout) :: ddb
760 : real(dp),intent(in) :: qphon(3),xred(3,natom)
761 : integer,intent(in) :: mpatpol(2),mpdir(3)
762 : real(dp),intent(in),optional :: blkval_fs(2,3,mpert,3,mpert,1)
763 : complex(dpc),intent(out) :: barmmom(ndim,(natom+5)*3)
764 : complex(dpc),intent(out) :: barmmom_tr((natom+5)*3,ndim)
765 : complex(dpc),intent(in) :: invbarmagsus(ndim,ndim)
766 : complex(dpc),intent(in) :: invhmat(ndim,ndim)
767 : complex(dpc),intent(in) :: magsus(ndim,ndim)
768 : complex(dpc),intent(out) :: mmom(ndim,(natom+5)*3)
769 : complex(dpc),intent(out) :: mmom_tr((natom+5)*3,ndim)
770 : complex(dpc),intent(out) :: zfield(ndim,(natom+5)*3)
771 : complex(dpc),intent(out) :: zfield_tr((natom+5)*3,ndim)
772 : !Local variables -------------------------
773 : !scalars
774 : integer :: fs2rs_
775 : integer :: iat1,icol,idir1,idir2,index,ipert1,ipert2,irow
776 : integer :: ipert1_red,idir1_red,jndex,zblok
777 : complex(dpc) :: qfac, qfac_tr
778 : character(len=1000) :: msg
779 : !arrays
780 68 : integer :: indexat1(ndim),indexdir1(ndim)
781 68 : integer :: indexat2((natom+5)*3),indexdir2((natom+5)*3)
782 34 : complex(dpc) :: mmom_alt(ndim,(natom+5)*3)
783 : character(len=1) :: cart(3)=(/'x','y','z'/)
784 :
785 : ! *********************************************************************
786 :
787 : !If fixed-spin case has been precalculated do less stuff
788 34 : fs2rs_=0; if (present(fs2rs)) fs2rs_=fs2rs
789 30 : if (fs2rs_==1) then
790 30 : if (.not.present(blkval_fs)) then
791 0 : write(msg, '(3a)' )' No fixed-spin array has been passed to magmom', &
792 0 : & ' but fs2rs=1 ',ch10
793 0 : ABI_ERROR(msg)
794 : end if
795 : end if
796 :
797 : !Extract the penalized moments
798 408 : do ipert2=1,natom+5
799 : !exclude strain perturbation
800 374 : if (ipert2==natom+3.or.ipert2==natom+4) cycle
801 1258 : do idir2=1,3
802 918 : icol=idir2+(ipert2-1)*3
803 918 : indexat2(icol)=ipert2
804 918 : indexdir2(icol)=idir2
805 :
806 918 : ipert1_red= 0
807 3128 : do iat1= mpatpol(1), mpatpol(2)
808 1836 : ipert1= natom + 11 + iat1
809 1836 : ipert1_red= ipert1_red + 1
810 1836 : idir1_red= 0
811 :
812 : !Apply a phase factor to adopt the same Gonze&Lee convention as for atomic
813 : !displacements at finite q (See M. Stengel PRB 2013).
814 7344 : qfac= exp(two_pi*(0.d0,1.d0)* dot_product(qphon,xred(:,iat1)))
815 7344 : qfac_tr= exp(-two_pi*(0.d0,1.d0)* dot_product(qphon,xred(:,iat1)))
816 8262 : do idir1= 1, 3
817 5508 : if (mpdir(idir1)==0) cycle
818 3672 : idir1_red= idir1_red + 1
819 3672 : irow=idir1_red+(ipert1_red-1)*nmdir
820 3672 : indexat1(irow)=iat1
821 3672 : indexdir1(irow)=idir1
822 3672 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
823 3672 : jndex= idir2 + 3*((ipert2-1)+mpert*((idir1-1)+3*(ipert1-1)))
824 :
825 5508 : if (fs2rs_==0) then
826 432 : if (iblok /=0 .and. ipert2 <= natom) then
827 288 : barmmom(irow,icol)= cmplx(ddb%val(1,index,iblok),ddb%val(2,index,iblok),kind=dpc) * qfac
828 288 : barmmom_tr(icol,irow)= cmplx(ddb%val(1,jndex,iblok),ddb%val(2,jndex,iblok),kind=dpc) * qfac_tr
829 144 : else if (jblok /=0 .and. ipert2 == natom+2) then
830 36 : barmmom(irow,icol)= cmplx(ddb%val(1,index,jblok),ddb%val(2,index,jblok),kind=dpc)
831 36 : barmmom_tr(icol,irow)= cmplx(ddb%val(1,jndex,jblok),ddb%val(2,jndex,jblok),kind=dpc)
832 108 : else if (lblok /=0 .and. ipert2 == natom+5) then
833 48 : barmmom(irow,icol)= cmplx(ddb%val(1,index,lblok),ddb%val(2,index,lblok),kind=dpc) * qfac
834 48 : barmmom_tr(icol,irow)= cmplx(ddb%val(1,jndex,lblok),ddb%val(2,jndex,lblok),kind=dpc) * qfac_tr
835 : end if
836 3240 : else if (fs2rs_==1) then
837 : zfield(irow,icol)= &
838 : & cmplx(blkval_fs(1,idir1,ipert1,idir2,ipert2,iblok), &
839 3240 : & blkval_fs(2,idir1,ipert1,idir2,ipert2,iblok),kind=dpc)
840 : zfield_tr(icol,irow)= &
841 : & cmplx(blkval_fs(1,idir2,ipert2,idir1,ipert1,iblok), &
842 3240 : & blkval_fs(2,idir2,ipert2,idir1,ipert1,iblok),kind=dpc)
843 : end if
844 :
845 : end do
846 : end do
847 : end do
848 : end do
849 :
850 : !Compute the Zeeman fields
851 34 : if (fs2rs_==0) then
852 4112 : zfield=-matmul(invbarmagsus,barmmom)
853 3300 : zfield_tr=-matmul(barmmom_tr,invbarmagsus)
854 : end if
855 :
856 : !Compute the moments
857 : if (fs2rs_==0) then
858 4112 : mmom=-matmul(magsus,zfield)
859 3436 : mmom_alt=matmul(invhmat,barmmom)
860 2740 : mmom_tr=matmul(barmmom_tr,invhmat)
861 30 : else if (fs2rs_==1) then
862 30840 : mmom=-matmul(magsus,zfield)
863 24750 : mmom_tr=-matmul(zfield_tr,magsus)
864 : return
865 : end if
866 :
867 4 : if (prtopt==1.and.prtvol>1) then
868 :
869 : !Write the results
870 4 : if (magpen > zero) then
871 :
872 : !Atomic displacements
873 4 : if (iblok /= 0) then
874 12 : call wrtout([ab_out,std_out], ' Local Zeeman fields induced by atomic displacements (at constrained magnetic moments)')
875 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
876 20 : do irow=1, ndim
877 308 : do icol=1, natom*3
878 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
879 288 : & indexat1(irow), cart(indexdir1(irow)), indexat2(icol), cart(indexdir2(icol)), &
880 576 : & real(zfield(irow,icol)), aimag(zfield(irow,icol))
881 880 : call wrtout([ab_out,std_out], msg)
882 : end do
883 : end do
884 12 : call wrtout([ab_out,std_out], ' ')
885 :
886 12 : call wrtout([ab_out,std_out], ' Local magnetic moments induced by atomic displacements (from induced Zeeman fields)')
887 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
888 20 : do irow=1, ndim
889 308 : do icol=1, natom*3
890 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
891 288 : & indexat1(irow), cart(indexdir1(irow)), indexat2(icol), cart(indexdir2(icol)), &
892 576 : & real(mmom(irow,icol)), aimag(mmom(irow,icol))
893 880 : call wrtout([ab_out,std_out], msg)
894 : end do
895 : end do
896 12 : call wrtout([ab_out,std_out], ' ')
897 :
898 12 : call wrtout([ab_out,std_out], ' Local magnetic moments induced by atomic displacements (from induced penalized moments)')
899 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
900 20 : do irow=1, ndim
901 308 : do icol=1, natom*3
902 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
903 288 : & indexat1(irow), cart(indexdir1(irow)), indexat2(icol), cart(indexdir2(icol)), &
904 576 : & real(mmom_alt(irow,icol)), aimag(mmom_alt(irow,icol))
905 880 : call wrtout([ab_out,std_out], msg)
906 : end do
907 : end do
908 12 : call wrtout([ab_out,std_out], ' ')
909 : end if
910 :
911 : !Electric field
912 4 : if (jblok /= 0) then
913 9 : call wrtout([ab_out,std_out], ' Local Zeeman fields induced by electric field (at constrained magnetic moments)')
914 9 : call wrtout([ab_out,std_out], ' atom1 dir E-dir Real Imag')
915 15 : do irow=1, ndim
916 51 : do icol=(natom+1)*3+1, (natom+2)*3
917 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
918 36 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
919 72 : & real(zfield(irow,icol)), aimag(zfield(irow,icol))
920 120 : call wrtout([ab_out,std_out], msg)
921 : end do
922 : end do
923 9 : call wrtout([ab_out,std_out], ' ')
924 9 : call wrtout([ab_out,std_out], ' Local magnetic moments induced by electric field (from induced Zeeman fields)')
925 9 : call wrtout([ab_out,std_out], ' atom1 dir E-dir Real Imag')
926 15 : do irow=1, ndim
927 51 : do icol=(natom+1)*3+1, (natom+2)*3
928 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
929 36 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
930 72 : & real(mmom(irow,icol)), aimag(mmom(irow,icol))
931 120 : call wrtout([ab_out,std_out], msg)
932 : end do
933 : end do
934 9 : call wrtout([ab_out,std_out], ' ')
935 9 : call wrtout([ab_out,std_out], ' Local magnetic moments induced by electric field (from induced penalized moments)')
936 9 : call wrtout([ab_out,std_out], ' atom1 dir E-dir Real Imag')
937 15 : do irow=1, ndim
938 51 : do icol=(natom+1)*3+1, (natom+2)*3
939 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
940 36 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
941 72 : & real(mmom_alt(irow,icol)), aimag(mmom_alt(irow,icol))
942 120 : call wrtout([ab_out,std_out], msg)
943 : end do
944 : end do
945 9 : call wrtout([ab_out,std_out], ' ')
946 : end if
947 :
948 : !Macroscopic Zeeman
949 4 : if (lblok /= 0) then
950 12 : call wrtout([ab_out,std_out], ' Local Zeeman fields induced by macroscopic Zeeman field (at constrained magnetic moments)')
951 12 : call wrtout([ab_out,std_out], ' atom1 dir B-dir Real Imag')
952 20 : do irow=1, ndim
953 68 : do icol=(natom+4)*3+1, (natom+5)*3
954 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
955 48 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
956 96 : & real(zfield(irow,icol)), aimag(zfield(irow,icol))
957 160 : call wrtout([ab_out,std_out], msg)
958 : end do
959 : end do
960 12 : call wrtout([ab_out,std_out], ' ')
961 12 : call wrtout([ab_out,std_out], ' Local magnetic moments induced by macroscopic Zeeman field (from induced Zeeman fields)')
962 12 : call wrtout([ab_out,std_out], ' atom1 dir B-dir Real Imag')
963 20 : do irow=1, ndim
964 68 : do icol=(natom+4)*3+1, (natom+5)*3
965 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
966 48 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
967 96 : & real(mmom(irow,icol)), aimag(mmom(irow,icol))
968 160 : call wrtout([ab_out,std_out], msg)
969 : end do
970 : end do
971 12 : call wrtout([ab_out,std_out], ' ')
972 12 : call wrtout([ab_out,std_out], ' Local magnetic moments induced by macroscopic Zeeman field (from induced penalized moments)')
973 12 : call wrtout([ab_out,std_out], ' atom1 dir B-dir Real Imag')
974 20 : do irow=1, ndim
975 68 : do icol=(natom+4)*3+1, (natom+5)*3
976 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
977 48 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
978 96 : & real(mmom_alt(irow,icol)), aimag(mmom_alt(irow,icol))
979 160 : call wrtout([ab_out,std_out], msg)
980 : end do
981 : end do
982 12 : call wrtout([ab_out,std_out], ' ')
983 : end if
984 :
985 4 : if (prtvol > 2) then
986 4 : if (iblok /= 0) then
987 12 : call wrtout([ab_out,std_out], ' Penalized local magnetic moments induced by atomic displacements ')
988 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
989 20 : do irow=1, ndim
990 308 : do icol=1, natom*3
991 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
992 288 : & indexat1(irow), cart(indexdir1(irow)), indexat2(icol), cart(indexdir2(icol)), &
993 576 : & real(barmmom(irow,icol)), aimag(barmmom(irow,icol))
994 880 : call wrtout([ab_out,std_out], msg)
995 : end do
996 : end do
997 12 : call wrtout([ab_out,std_out], ' ')
998 : end if
999 4 : if (jblok /= 0) then
1000 9 : call wrtout([ab_out,std_out], ' Penalized local magnetic moments induced by electric field ')
1001 9 : call wrtout([ab_out,std_out], ' atom1 dir E-dir Real Imag')
1002 15 : do irow=1, ndim
1003 51 : do icol=(natom+1)*3+1, (natom+2)*3
1004 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
1005 36 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
1006 72 : & real(barmmom(irow,icol)), aimag(barmmom(irow,icol))
1007 120 : call wrtout([ab_out,std_out], msg)
1008 : end do
1009 : end do
1010 9 : call wrtout([ab_out,std_out], ' ')
1011 : end if
1012 4 : if (lblok /= 0) then
1013 12 : call wrtout([ab_out,std_out], ' Penalized local magnetic moments induced by macroscopic Zeeman field ')
1014 12 : call wrtout([ab_out,std_out], ' atom1 dir B-dir Real Imag')
1015 20 : do irow=1, ndim
1016 68 : do icol=(natom+4)*3+1, (natom+5)*3
1017 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
1018 48 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
1019 96 : & real(barmmom(irow,icol)), aimag(barmmom(irow,icol))
1020 160 : call wrtout([ab_out,std_out], msg)
1021 : end do
1022 : end do
1023 12 : call wrtout([ab_out,std_out], ' ')
1024 : end if
1025 : end if
1026 : end if
1027 : end if
1028 :
1029 : !Store the FS and RS flavors on the DDB array
1030 48 : do ipert2=1,natom+5
1031 180 : do idir2=1,3
1032 132 : icol=idir2+(ipert2-1)*3
1033 132 : indexat2(icol)=ipert2
1034 132 : indexdir2(icol)=idir2
1035 :
1036 132 : ipert1_red= 0
1037 440 : do iat1= mpatpol(1), mpatpol(2)
1038 264 : ipert1= natom + 11 + iat1
1039 264 : ipert1_red= ipert1_red + 1
1040 264 : idir1_red= 0
1041 1188 : do idir1= 1, 3
1042 792 : if (mpdir(idir1)==0) cycle
1043 528 : idir1_red= idir1_red + 1
1044 528 : irow=idir1_red+(ipert1_red-1)*nmdir
1045 528 : indexat1(irow)=iat1
1046 528 : indexdir1(irow)=idir1
1047 528 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1048 528 : jndex= idir2 + 3*((ipert2-1)+mpert*((idir1-1)+3*(ipert1-1)))
1049 :
1050 528 : zblok= 0
1051 528 : if (iblok /=0 .and. ipert2 <= natom) zblok= iblok
1052 528 : if (jblok /=0 .and. ipert2 == natom+2) zblok= jblok
1053 528 : if (lblok /=0 .and. ipert2 == natom+5) zblok= lblok
1054 792 : if (zblok /= 0) then
1055 372 : ddb%val_fs(1,index,zblok)=real(zfield(irow,icol))
1056 372 : ddb%val_fs(2,index,zblok)=aimag(zfield(irow,icol))
1057 372 : ddb%val_fs(1,jndex,zblok)=real(zfield_tr(icol,irow))
1058 372 : ddb%val_fs(2,jndex,zblok)=aimag(zfield_tr(icol,irow))
1059 :
1060 372 : ddb%val_rs(1,index,zblok)=real(mmom(irow,icol))
1061 372 : ddb%val_rs(2,index,zblok)=aimag(mmom(irow,icol))
1062 372 : ddb%val_rs(1,jndex,zblok)=real(mmom_tr(icol,irow))
1063 372 : ddb%val_rs(2,jndex,zblok)=aimag(mmom_tr(icol,irow))
1064 : end if
1065 :
1066 : end do
1067 : end do
1068 : end do
1069 : end do
1070 :
1071 34 : end subroutine magmom
1072 : !!***
1073 :
1074 : !!****f* m_ddb_magpen/mp_d2etot
1075 : !! NAME
1076 : !! mp_d2etot
1077 : !!
1078 : !! FUNCTION
1079 : !! Calculate the different magnetic flavors (see mpopt below) of the second-
1080 : !! order derivatives of total energy
1081 : !!
1082 : !! INPUTS
1083 : !! barmagsus(ndim,ndim)= Penalized spin-sussceptibility tensor (\bar{\chi})
1084 : !! (equal to barmom^{\dagger} in the nondissipative regime)
1085 : !! ddb= the ddb object
1086 : !! iblok= index of the IFCs block
1087 : !! invhmat(ndim,ndim)= (I-\alpha\bar{\chi})^-1 matrix
1088 : !! magsus(ndim,ndim)= Spin-sussceptibility tensor
1089 : !! magpen= magnetic penalty amplitude
1090 : !! mpert =maximum number of ipert
1091 : !! mpopt = 1 calculate the frozen-spin second-order quantities
1092 : !! 2 calculate the relaxed-spin second-order quantities
1093 : !! natom= number of atoms in unit cell
1094 : !! nblok= number of blocks in the DDB
1095 : !! ndim= number of local magnetic degres of freedom
1096 : !! qphon= momentum wave-vector
1097 : !! xred(3,natom)= reduced atomic coordinates
1098 : !! zfield(ndim,(natom+2)*3)= First-order induced Zeeman fields
1099 : !! zfield_tr(natom+2)*3,ndim)= Linear-responses to external Zeeman fields
1100 : !! (equal to zfield^{\dagger} in the nondissipative regime)
1101 : !! fs2rs= (optional) if 1, the routine starts from a precalculated blkval_fs
1102 : !! blkval_fs(2,3,mpert,3,mpert,1)= fixed-spin 2nd-order derivatives
1103 : !!
1104 : !! OUTPUT
1105 : !! ddb%val_fs(2,msize,nblok)= second-order derivatives at fixed spin.
1106 : !! ddb%val_rs(2,msize,nblok)= second-order derivatives at relaxed spin.
1107 : !! blkval_rs(2,3,mpert,3,mpert,1)= (optional) relaxed-spin 2nd-order derivatives
1108 : !!
1109 : !! SOURCE
1110 :
1111 68 : subroutine mp_d2etot(barmagsus,ddb,&
1112 34 : & iblok,magsus,mpert,mpopt,&
1113 34 : & natom,ndim,zfield,zfield_tr,&
1114 60 : & fs2rs,blkval_fs,blkval_rs) !optional
1115 :
1116 : !Arguments -------------------------------
1117 : !scalars
1118 : integer,intent(in) :: iblok,mpert,mpopt,natom,ndim
1119 : integer,intent(in),optional :: fs2rs
1120 : ! real(dp),intent(in) :: magpen
1121 : !arrays
1122 : type(ddb_type),intent(inout) :: ddb
1123 : real(dp),intent(in),optional :: blkval_fs(2,3,mpert,3,mpert,1)
1124 : real(dp),intent(out),optional :: blkval_rs(2,3,mpert,3,mpert,1)
1125 : complex(dpc),intent(in) :: barmagsus(ndim,ndim)
1126 : ! complex(dpc),intent(in) :: invhmat(ndim,ndim)
1127 : complex(dpc),intent(in) :: magsus(ndim,ndim)
1128 : complex(dpc),intent(in) :: zfield(ndim,(natom+5)*3)
1129 : complex(dpc),intent(in) :: zfield_tr((natom+5)*3,ndim)
1130 : !Local variables -------------------------
1131 : !scalars
1132 : integer :: fs2rs_
1133 : integer :: idir1,idir2,ipert1,ipert2,index,irow,icol
1134 : complex(dpc) :: val_ps,val_fs, val_rs
1135 : character(len=1000) :: msg
1136 : !arrays
1137 :
1138 : ! *********************************************************************
1139 :
1140 : !If fixed-spin case has been precalculated do less stuff
1141 34 : fs2rs_=0; if (present(fs2rs)) fs2rs_=fs2rs
1142 30 : if (fs2rs_==1) then
1143 30 : if (.not.present(blkval_fs)) then
1144 0 : write(msg, '(3a)' )' No fixed-spin array has been passed to mp_d2etot', &
1145 0 : & ' but fs2rs=1 ',ch10
1146 0 : ABI_ERROR(msg)
1147 : end if
1148 30 : if (.not.present(blkval_rs)) then
1149 0 : write(msg, '(3a)' )' No relaxed-spin array has been passed to mp_d2etot', &
1150 0 : & ' but fs2rs=1 ',ch10
1151 0 : ABI_ERROR(msg)
1152 : end if
1153 : end if
1154 :
1155 : !Extract the penalized/constrained quantities
1156 408 : do ipert2= 1, natom+5
1157 1530 : do idir2= 1, 3
1158 1122 : icol= (ipert2-1)*3 + idir2
1159 13838 : do ipert1= 1, natom+5
1160 50490 : do idir1= 1, 3
1161 37026 : irow= (ipert1-1)*3 + idir1
1162 37026 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1163 :
1164 37026 : if (fs2rs_==0) then
1165 : !Extract the penalized second-order derivatives
1166 4356 : val_ps= cmplx(ddb%val(1,index,iblok),ddb%val(2,index,iblok),kind=dpc)
1167 : !Calculate the fixed-spin flavor
1168 : val_fs= val_ps + &
1169 139392 : & sum( zfield_tr(irow,:) * matmul( barmagsus,zfield(:,icol) ) )
1170 4356 : ddb%val_fs(1,index,iblok)= real(val_fs)
1171 4356 : ddb%val_fs(2,index,iblok)= aimag(val_fs)
1172 32670 : else if (fs2rs_==1) then
1173 : val_fs= &
1174 : & cmplx(blkval_fs(1,idir1,ipert1,idir2,ipert2,iblok), &
1175 32670 : & blkval_fs(2,idir1,ipert1,idir2,ipert2,iblok),kind=dpc)
1176 : end if
1177 :
1178 49368 : if (mpopt==2) then
1179 : !Calculate the relaxed-spin flavor
1180 : val_rs= val_fs - &
1181 1184832 : & sum( zfield_tr(irow,:) * matmul( magsus,zfield(:,icol) ) )
1182 37026 : if (fs2rs_==0) then
1183 4356 : ddb%val_rs(1,index,iblok)= real(val_rs)
1184 4356 : ddb%val_rs(2,index,iblok)= aimag(val_rs)
1185 32670 : else if (fs2rs_==1) then
1186 32670 : blkval_rs(1,idir1,ipert1,idir2,ipert2,iblok)=real(val_rs)
1187 32670 : blkval_rs(2,idir1,ipert1,idir2,ipert2,iblok)=aimag(val_rs)
1188 : end if
1189 : end if
1190 :
1191 : end do
1192 : end do
1193 : end do
1194 : end do
1195 :
1196 64 : end subroutine mp_d2etot
1197 : !!***
1198 :
1199 : !!****f* m_ddb_magpen/mp_d2etot_print
1200 : !! NAME
1201 : !! mp_d2etot_print
1202 : !!
1203 : !! FUNCTION
1204 : !! Write on output file the fixed- and relaxed-spin susceptibilities
1205 : !!
1206 : !! INPUTS
1207 : !! blkval= 2nd-order susceptibilities matrix
1208 : !! kblok= index of the current block
1209 : !! opt= 1 write the frozen-spin second-order quantities
1210 : !! 2 write the relaxed-spin second-order quantities
1211 : !! omega= frequency of the perturbation
1212 : !! qphon= momentum wave-vector
1213 : !! ucvol= unit-cell volume
1214 : !!
1215 : !! OUTPUT
1216 : !!
1217 : !! SOURCE
1218 :
1219 16 : subroutine mp_d2etot_print(ddb,blkval,kblok,mpert,natom,nblok,opt,omega,qeq0,qphnrm,qphon,msize)
1220 :
1221 : !Arguments -------------------------------
1222 : !scalars
1223 : class(ddb_type),intent(in) :: ddb
1224 : integer,intent(in) :: kblok,mpert,natom,nblok,opt,msize
1225 : logical,intent(in) :: qeq0
1226 : !arrays
1227 : real(dp),intent(in) :: omega(3)
1228 : real(dp),intent(in) :: blkval(2,msize,nblok)
1229 : real(dp),intent(inout) :: qphnrm(3),qphon(3,3)
1230 :
1231 : !Local variables -------------------------
1232 : !scalars
1233 : integer :: iblok,idir1,idir2,ipert1,ipert2,irow,icol
1234 : integer :: rftyp,index
1235 : character(len=1000) :: msg
1236 : !arrays
1237 : integer :: rfelfd(4),rfmagn(4),rfphon(4),rfstrs(4)
1238 : real(dp) :: val(2)
1239 : character(len=1) :: cart(3)=(/'x','y','z'/)
1240 :
1241 : ! *********************************************************************
1242 :
1243 4 : rfelfd(:)=0
1244 4 : rfphon(:)=0
1245 4 : rfstrs(:)=0
1246 4 : rfmagn(:)=0
1247 4 : rftyp = 1
1248 :
1249 : !IFCs
1250 12 : rfphon(1:2)=1
1251 4 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega)
1252 4 : if (iblok/=0.and.iblok==kblok) then
1253 4 : if (opt==1) then
1254 0 : call wrtout([ab_out,std_out], ' Frozen-spin interatomic force constants')
1255 4 : else if (opt==2) then
1256 12 : call wrtout([ab_out,std_out], ' Relaxed-spin interatomic force constants')
1257 : end if
1258 12 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
1259 28 : do ipert1= 1, natom
1260 100 : do idir1= 1, 3
1261 72 : irow=( ipert1-1)*3 + idir1
1262 504 : do ipert2= 1, natom
1263 1800 : do idir2= 1, 3
1264 1296 : icol=( ipert2-1)*3 + idir2
1265 1296 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1266 3888 : val(:)=blkval(:,index,kblok)
1267 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)') &
1268 1296 : & ipert1, cart(idir1), ipert2, cart(idir2), val(1), val(2)
1269 4320 : call wrtout([ab_out,std_out], msg)
1270 : end do
1271 : end do
1272 240 : call wrtout([ab_out,std_out], ' ')
1273 : end do
1274 : end do
1275 : end if
1276 :
1277 4 : if (qeq0) then
1278 :
1279 : !Born charges
1280 9 : rfphon(1:2)=1
1281 9 : rfelfd(1:2)=2
1282 3 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega)
1283 3 : if (iblok/=0.and.iblok==kblok) then
1284 3 : if (opt==1) then
1285 0 : call wrtout([ab_out,std_out], ' Frozen-spin Born effective charges')
1286 3 : else if (opt==2) then
1287 9 : call wrtout([ab_out,std_out], ' Relaxed-spin Born effective charges')
1288 : end if
1289 9 : call wrtout([ab_out,std_out], ' E-dir atom dir Real Imag')
1290 3 : ipert1= natom + 2
1291 12 : do idir1= 1, 3
1292 63 : do ipert2= 1, natom
1293 225 : do idir2= 1, 3
1294 162 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1295 486 : val(:)=blkval(:,index,kblok)
1296 : write(msg,'(3x,a2,7x,i3,4x,a2,2x,2es18.9)') &
1297 162 : & cart(idir1), ipert2, cart(idir2), val(1), val(2)
1298 540 : call wrtout([ab_out,std_out], msg)
1299 : end do
1300 : end do
1301 30 : call wrtout([ab_out,std_out], ' ')
1302 : end do
1303 : end if
1304 :
1305 : !Dielectric tensor
1306 : iblok=0
1307 3 : rfphon(:)=0
1308 9 : rfelfd(1:2)=2
1309 3 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega)
1310 3 : if (iblok/=0.and.iblok==kblok) then
1311 3 : if (opt==1) then
1312 0 : call wrtout([ab_out,std_out], ' Frozen-spin clamped-ion dielectric tensor')
1313 3 : else if (opt==2) then
1314 9 : call wrtout([ab_out,std_out], ' Relaxed-spin clamped-ion dielectric tensor')
1315 : end if
1316 9 : call wrtout([ab_out,std_out], ' dir dir Real Imag')
1317 3 : ipert1= ddb%natom + 2
1318 3 : ipert2= ddb%natom + 2
1319 12 : do idir2= 1, 3
1320 36 : do idir1= 1, 3
1321 27 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1322 81 : val(:)=blkval(:,index,kblok)
1323 27 : write(msg,'(2x,a2,3x,a2,2x,2es18.9)' ) cart(idir1), cart(idir2), &
1324 54 : & val(1), val(2)
1325 90 : call wrtout([ab_out,std_out], msg)
1326 : end do
1327 30 : call wrtout([ab_out,std_out], ' ')
1328 : end do
1329 : end if
1330 :
1331 : !Magnetoelectric susceptibility
1332 : iblok=0
1333 3 : rfphon(:)=0
1334 3 : rfelfd(1)=0
1335 3 : rfelfd(2)=2
1336 3 : rfmagn(1)=1
1337 : rfmagn(2)=0
1338 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
1339 3 : & rfmagn=rfmagn)
1340 3 : if (iblok/=0.and.iblok==kblok) then
1341 3 : if (opt==1) then
1342 0 : call wrtout([ab_out,std_out], ' Frozen-spin clamped-ion magnetoelectric susceptibility')
1343 3 : else if (opt==2) then
1344 9 : call wrtout([ab_out,std_out], ' Relaxed-spin clamped-ion magnetoelectric susceptibility')
1345 : end if
1346 9 : call wrtout([ab_out,std_out], ' M-dir E-dir Real Imag')
1347 3 : ipert1= ddb%natom + 5
1348 3 : ipert2= ddb%natom + 2
1349 12 : do idir2= 1, 3
1350 36 : do idir1= 1, 3
1351 27 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1352 81 : val(:)=blkval(:,index,kblok)
1353 27 : write(msg,'(2x,a2,3x,a2,2x,2es18.9)' ) cart(idir1), cart(idir2), &
1354 54 : & val(1), val(2)
1355 90 : call wrtout([ab_out,std_out], msg)
1356 : end do
1357 30 : call wrtout([ab_out,std_out], ' ')
1358 : end do
1359 9 : call wrtout([ab_out,std_out], ' P-dir H-dir Real Imag')
1360 3 : ipert1= ddb%natom + 2
1361 3 : ipert2= ddb%natom + 5
1362 12 : do idir2= 1, 3
1363 36 : do idir1= 1, 3
1364 27 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1365 81 : val(:)=blkval(:,index,kblok)
1366 27 : write(msg,'(2x,a2,3x,a2,2x,2es18.9)' ) cart(idir1), cart(idir2), &
1367 54 : & val(1), val(2)
1368 90 : call wrtout([ab_out,std_out], msg)
1369 : end do
1370 30 : call wrtout([ab_out,std_out], ' ')
1371 : end do
1372 : end if
1373 :
1374 : end if
1375 :
1376 : !Magnetic susceptibility
1377 : iblok=0
1378 4 : rfphon(:)=0
1379 4 : rfelfd(:)=0
1380 4 : rfmagn(1)=1
1381 4 : rfmagn(2)=1
1382 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
1383 4 : & rfmagn=rfmagn)
1384 4 : if (iblok/=0.and.iblok==kblok) then
1385 4 : if (opt==1) then
1386 0 : call wrtout([ab_out,std_out], ' Frozen-spin clamped-ion magnetic susceptibility')
1387 4 : else if (opt==2) then
1388 12 : call wrtout([ab_out,std_out], ' Relaxed-spin clamped-ion magnetic susceptibility')
1389 : end if
1390 12 : call wrtout([ab_out,std_out], ' dir dir Real Imag')
1391 4 : ipert1= ddb%natom + 5
1392 4 : ipert2= ddb%natom + 5
1393 16 : do idir2= 1, 3
1394 48 : do idir1= 1, 3
1395 36 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1396 108 : val(:)=blkval(:,index,kblok)
1397 36 : write(msg,'(2x,a2,3x,a2,2x,2es18.9)' ) cart(idir1), cart(idir2), &
1398 72 : & val(1), val(2)
1399 120 : call wrtout([ab_out,std_out], msg)
1400 : end do
1401 40 : call wrtout([ab_out,std_out], ' ')
1402 : end do
1403 : end if
1404 :
1405 : !Magnetic Born effective charges
1406 : iblok=0
1407 4 : rfphon(1)=1
1408 4 : rfelfd(:)=0
1409 4 : rfmagn(1)=0
1410 : rfmagn(2)=1
1411 : call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
1412 4 : & rfmagn=rfmagn)
1413 4 : if (iblok/=0.and.iblok==kblok) then
1414 4 : if (opt==1) then
1415 0 : call wrtout([ab_out,std_out], ' Frozen-spin magnetic Born effective charges')
1416 4 : else if (opt==2) then
1417 12 : call wrtout([ab_out,std_out], ' Relaxed-spin magnetic Born effective charges')
1418 : end if
1419 12 : call wrtout([ab_out,std_out], ' B-dir atom dir Real Imag')
1420 :
1421 4 : ipert1= ddb%natom + 5
1422 16 : do idir1= 1, 3
1423 84 : do ipert2= 1, natom
1424 300 : do idir2= 1, 3
1425 216 : index= idir1 + 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1)))
1426 648 : val(:)=blkval(:,index,kblok)
1427 : write(msg,'(3x,a2,7x,i3,4x,a2,2x,2es18.9)') &
1428 216 : & cart(idir1), ipert2, cart(idir2), val(1), val(2)
1429 720 : call wrtout([ab_out,std_out], msg)
1430 : end do
1431 : end do
1432 40 : call wrtout([ab_out,std_out], ' ')
1433 : end do
1434 : end if
1435 :
1436 4 : end subroutine mp_d2etot_print
1437 : !!***
1438 :
1439 : !!****f* m_ddb_magpen/asrw0
1440 : !! NAME
1441 : !! asrw0
1442 : !!
1443 : !! FUNCTION
1444 : !! Impose the Acoustic Sum Rule from the w=0 IFCs
1445 : !!
1446 : !! INPUTS
1447 : !! natom= namber of atoms
1448 : !! option= if 0, this is the w=0 case, calculate delta_asrw0
1449 : ! if 1, use the previously calculated delta_asrw0
1450 : !!
1451 : !! OUTPUT
1452 : !! ifcmat(3*natom,3*natom)= IFC matrix after ASR has been applied.
1453 : !! delta_asrw0(3*natom,3)= Amount to remove in order to enforce ASR.
1454 : !!
1455 : !! SOURCE
1456 :
1457 0 : subroutine asrw0(delta_asrw0,ifcmat,natom,option)
1458 :
1459 : !Arguments -------------------------------
1460 : !scalars
1461 : integer,intent(in) :: natom,option
1462 : !arrays
1463 : real(dp),intent(inout) :: delta_asrw0(3*natom,3)
1464 : complex(dpc),intent(inout) :: ifcmat(3*natom,3*natom)
1465 : !Local variables -------------------------
1466 : !scalars
1467 : integer :: icol,idir1,idir2,ipert1,ipert2,irow
1468 : !arrays
1469 :
1470 : ! *********************************************************************
1471 :
1472 : !Calculate the ASR correction
1473 0 : if (option == 0) then
1474 0 : delta_asrw0= zero
1475 0 : do idir1= 1, 3
1476 0 : do ipert1= 1, natom
1477 0 : irow= (ipert1-1)*3 + idir1
1478 0 : do idir2= 1, 3
1479 0 : do ipert2= 1, natom
1480 0 : icol= (ipert2-1)*3 + idir2
1481 : delta_asrw0(irow,idir2)=delta_asrw0(irow,idir2) + &
1482 0 : & real(ifcmat(irow,icol))
1483 : end do
1484 : end do
1485 : end do
1486 : end do
1487 : end if
1488 :
1489 : !Apply the ASR
1490 0 : do idir1= 1, 3
1491 0 : do ipert1= 1, natom
1492 0 : irow= (ipert1-1)*3 + idir1
1493 0 : do idir2= 1, 3
1494 0 : icol= (ipert1-1)*3 + idir2
1495 0 : ifcmat(irow,icol)= ifcmat(irow,icol) - delta_asrw0(irow,idir2)
1496 : end do
1497 : end do
1498 : end do
1499 :
1500 0 : end subroutine asrw0
1501 : !!***
1502 :
1503 : !!****f* m_ddb_magpen/berrycurv_ss
1504 : !! NAME
1505 : !! berrycurv_ss
1506 : !!
1507 : !! FUNCTION
1508 : !! Calculate the Berry curvature of the inverse magnetic susceptibility
1509 : !! (this is equivalent to the G^(ss) matrix of S.Ren et al.)
1510 : !!
1511 : !! INPUTS
1512 : !! blkval(2,3*mpert*3*mpert*3*mpert,nblok)= Third-order derivative matrices
1513 : !! In our case, the nblok is restricted to iblok
1514 : !! iblok= index of the current block
1515 : !! invbarmagsus(ndim,ndim)= Inverse of the penalized spin-sussceptibility tensor
1516 : !! mpatpol(2) = Atoms on which the magnetic penalty has been applied
1517 : !! mpdir(3) = Directions along which the spin-degrees of freedom have been stiffened
1518 : !! mpert =maximum number of ipert
1519 : !! natom= number of atoms in unit cell
1520 : !! nblok= number of blocks in the DDB
1521 : !! nmdir= number of directions along which the magnetic penalty was applied
1522 : !! ndim= dimension of the square susceptibilities
1523 : !! prtvol= control the volume of information written on output
1524 : !!
1525 : !! OUTPUT
1526 : !! bc_ss(ndim,ndim)= Berry-curvature of the inverse local-spin susceptibility
1527 : !!
1528 : !! SOURCE
1529 :
1530 2 : subroutine berrycurv_ss(bc_barmagsus,bc_ss,ddb_lw,iblok,invbarmagsus,&
1531 2 : & mpatpol,mpdir,mpert,natom,ndim,nmdir,prtvol,qphon,xred)
1532 :
1533 : !Arguments -------------------------------
1534 : !scalars
1535 : integer,intent(in) :: iblok,mpert,natom,ndim,nmdir,prtvol
1536 : !arrays
1537 : type(ddb_type),intent(inout) :: ddb_lw
1538 : integer,intent(in) :: mpatpol(2),mpdir(3)
1539 : real(dp),intent(in) :: qphon(3),xred(3,natom)
1540 : complex(dpc),intent(in) :: invbarmagsus(ndim,ndim)
1541 : complex(dpc),intent(out) :: bc_ss(ndim,ndim)
1542 : complex(dpc),intent(out) :: bc_barmagsus(ndim,ndim)
1543 : !Local variables -------------------------
1544 : !scalars
1545 : integer :: iat1,iat2,icol,idir1,idir2,idir3,index
1546 : integer :: ipert1,ipert2,ipert3,irow
1547 : integer :: ipert1_red,ipert2_red,idir1_red,idir2_red
1548 : complex(dpc), parameter :: ione=(0.d0,1.d0)
1549 : character(len=1000) :: msg
1550 : !arrays
1551 4 : complex(dpc) :: idty(ndim,ndim)
1552 2 : integer :: indexat(ndim),indexdir(ndim)
1553 : character(len=1) :: cart(3)=(/'x','y','z'/)
1554 :
1555 : ! *********************************************************************
1556 :
1557 : !Extract the Berry-curvature of the penalized susceptibility
1558 42 : idty=(zero,zero)
1559 2 : ipert3= natom + 9
1560 2 : idir3= 1
1561 2 : ipert2_red= 0
1562 6 : do iat2= mpatpol(1), mpatpol(2)
1563 4 : ipert2= natom + 11 + iat2
1564 4 : ipert2_red= ipert2_red + 1
1565 4 : idir2_red= 0
1566 18 : do idir2= 1, 3
1567 12 : if (mpdir(idir2)==0) cycle
1568 8 : idir2_red= idir2_red + 1
1569 8 : icol=idir2_red+(ipert2_red-1)*nmdir
1570 8 : indexat(icol)=iat2
1571 8 : indexdir(icol)=idir2
1572 8 : idty(icol,icol)=(one,zero)
1573 8 : ipert1_red=0
1574 28 : do iat1= mpatpol(1), mpatpol(2)
1575 16 : ipert1= natom + 11 + iat1
1576 16 : ipert1_red= ipert1_red + 1
1577 16 : idir1_red= 0
1578 76 : do idir1= 1, 3
1579 48 : if (mpdir(idir1)==0) cycle
1580 32 : idir1_red=idir1_red+1
1581 32 : irow=idir1_red+(ipert1_red-1)*nmdir
1582 : index = idir1 + &
1583 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
1584 32 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
1585 :
1586 : !Apply a phase factor to adopt the same Gonze&Lee convention as for atomic
1587 : !displacements at finite q (See M. Stengel PRB 2013).
1588 : bc_barmagsus(irow,icol)= -one* & !To convert from d3etot to local susc.
1589 : & cmplx(ddb_lw%val(1,index,iblok),ddb_lw%val(2,index,iblok),kind=dpc) &
1590 160 : & * exp(two_pi*(0.d0,1.d0)*dot_product(qphon,xred(:,iat1)-xred(:,iat2)))
1591 :
1592 : end do
1593 : end do
1594 : end do
1595 : end do
1596 :
1597 : !Calculate the Berry-curvature of the inverse magnetic susceptibility
1598 600 : bc_ss=-matmul(invbarmagsus,matmul(bc_barmagsus,invbarmagsus))
1599 :
1600 2 : if (prtvol>1) then
1601 6 : call wrtout([ab_out,std_out], ' Fixed-spin Berry curvature of the inverse spin susceptibility ')
1602 6 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
1603 10 : do irow=1, ndim
1604 42 : do icol=1, ndim
1605 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
1606 32 : & indexat(irow), cart(indexdir(irow)), indexat(icol), cart(indexdir(icol)), &
1607 64 : & real(bc_ss(irow,icol)), aimag(bc_ss(irow,icol))
1608 104 : call wrtout([ab_out,std_out], msg)
1609 : end do
1610 : end do
1611 6 : call wrtout([ab_out,std_out], ' ')
1612 : end if
1613 :
1614 : !Store the FM flavor in the DDB object
1615 : ipert2_red= 0
1616 6 : do iat2= mpatpol(1), mpatpol(2)
1617 4 : ipert2= natom + 11 + iat2
1618 4 : ipert2_red= ipert2_red + 1
1619 4 : idir2_red= 0
1620 18 : do idir2= 1, 3
1621 12 : if (mpdir(idir2)==0) cycle
1622 8 : idir2_red= idir2_red + 1
1623 8 : icol=idir2_red+(ipert2_red-1)*nmdir
1624 8 : ipert1_red=0
1625 28 : do iat1= mpatpol(1), mpatpol(2)
1626 16 : ipert1= natom + 11 + iat1
1627 16 : ipert1_red= ipert1_red + 1
1628 16 : idir1_red= 0
1629 76 : do idir1= 1, 3
1630 48 : if (mpdir(idir1)==0) cycle
1631 32 : idir1_red=idir1_red+1
1632 32 : irow=idir1_red+(ipert1_red-1)*nmdir
1633 : index = idir1 + &
1634 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
1635 32 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
1636 :
1637 32 : ddb_lw%val_fs(1,index,iblok)= real(bc_ss(irow,icol))
1638 64 : ddb_lw%val_fs(2,index,iblok)= aimag(bc_ss(irow,icol))
1639 : end do
1640 : end do
1641 : end do
1642 : end do
1643 :
1644 2 : end subroutine berrycurv_ss
1645 : !!***
1646 :
1647 : !!****f* m_ddb_magpen/berrycurv_sp
1648 : !! NAME
1649 : !! berrycurv_sp
1650 : !!
1651 : !! FUNCTION
1652 : !! Calculate the Berry curvature of the spin-phonon Hessian
1653 : !! (equivalent to the Berry curvature of the induced Zeeman fields at constrained
1654 : !! magnetic moments)
1655 : !!
1656 : !! INPUTS
1657 : !! ddb_lw= Third-order derivative ddb
1658 : !! In our case, the nblok is restricted to iblok
1659 : !! iblok= index of the current block
1660 : !! invbarmagsus(ndim,ndim)= Inverse of the penalized spin-sussceptibility tensor
1661 : !! mpatpol(2) = Atoms on which the magnetic penalty has been applied
1662 : !! mpdir(3) = Directions along which the spin-degrees of freedom have been stiffened
1663 : !! mpert =maximum number of ipert
1664 : !! natom= number of atoms in unit cell
1665 : !! nblok= number of blocks in the DDB
1666 : !! nmdir= number of directions along which the magnetic penalty was applied
1667 : !! ndim= dimension of the square susceptibilities
1668 : !! prtvol= control the volume of information written on output
1669 : !!
1670 : !! OUTPUT
1671 : !!
1672 : !! SOURCE
1673 :
1674 2 : subroutine berrycurv_sp(barmmom,bc_sp,bc_ss,ddb_lw,iblok,invbarmagsus,&
1675 2 : & jblok,lblok,mpatpol,mpdir,mpert,natom,ndim,nmdir,prtvol,qphon,xred)
1676 :
1677 : !Arguments -------------------------------
1678 : !scalars
1679 : integer,intent(in) :: iblok,jblok,lblok,mpert,natom,ndim,nmdir,prtvol
1680 : !arrays
1681 : type(ddb_type),intent(inout) :: ddb_lw
1682 : integer,intent(in) :: mpatpol(2),mpdir(3)
1683 : real(dp),intent(in) :: qphon(3),xred(3,natom)
1684 : complex(dpc),intent(in) :: barmmom(ndim,(natom+5)*3)
1685 : complex(dpc),intent(in) :: bc_ss(ndim,ndim)
1686 : complex(dpc),intent(in) :: invbarmagsus(ndim,ndim)
1687 : complex(dpc),intent(out) :: bc_sp(ndim,(natom+5)*3)
1688 : !Local variables -------------------------
1689 : !scalars
1690 : integer :: iat1,icol,idir1,idir2,idir3,index,ipert1,ipert2,ipert3,irow
1691 : integer :: ipert1_red,idir1_red,jndex
1692 : complex(dpc), parameter :: ione=(0.d0,1.d0)
1693 : complex(dpc) :: qfac
1694 : character(len=1000) :: msg
1695 : !arrays
1696 4 : integer :: indexat1(ndim),indexdir1(ndim)
1697 4 : integer :: indexat2((natom+5)*3),indexdir2((natom+5)*3)
1698 2 : complex(dpc) :: bc_barsp(ndim,(natom+5)*3)
1699 : ! complex(dpc) :: bc_ps((natom+2)*3,ndim)
1700 : character(len=1) :: cart(3)=(/'x','y','z'/)
1701 :
1702 : ! *********************************************************************
1703 :
1704 : !Extract the berry curvature of the penalized moments
1705 332 : bc_barsp=(zero,zero)
1706 : ipert3= natom + 9
1707 : idir3= 1
1708 24 : do ipert2=1,natom+5
1709 90 : do idir2=1,3
1710 66 : icol=idir2+(ipert2-1)*3
1711 66 : indexat2(icol)=ipert2
1712 66 : indexdir2(icol)=idir2
1713 :
1714 66 : ipert1_red= 0
1715 220 : do iat1= mpatpol(1), mpatpol(2)
1716 132 : ipert1= natom + 11 + iat1
1717 132 : ipert1_red= ipert1_red + 1
1718 132 : idir1_red= 0
1719 :
1720 : !Apply a phase factor to adopt the same Gonze&Lee convention as for atomic
1721 : !displacements at finite q (See M. Stengel PRB 2013).
1722 528 : qfac= exp(two_pi*(0.d0,1.d0)* dot_product(qphon,xred(:,iat1)))
1723 594 : do idir1= 1, 3
1724 396 : if (mpdir(idir1)==0) cycle
1725 264 : idir1_red= idir1_red + 1
1726 264 : irow=idir1_red+(ipert1_red-1)*nmdir
1727 264 : indexat1(irow)=iat1
1728 264 : indexdir1(irow)=idir1
1729 : index = idir1 + &
1730 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
1731 264 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
1732 :
1733 396 : if (iblok /=0 .and. ipert2 <= natom) then
1734 : bc_barsp(irow,icol)= -one* & !To go from d3etot to induced local field
1735 144 : & cmplx(ddb_lw%val(1,index,iblok),ddb_lw%val(2,index,iblok),kind=dpc) * qfac
1736 120 : else if (jblok /=0 .and. ipert2 == natom+2) then
1737 : bc_barsp(irow,icol)= -one* &
1738 12 : & cmplx(ddb_lw%val(1,index,jblok),ddb_lw%val(2,index,jblok),kind=dpc)
1739 108 : else if (lblok /=0 .and. ipert2 == natom+5) then
1740 : bc_barsp(irow,icol)= -one* &
1741 12 : & cmplx(ddb_lw%val(1,index,lblok),ddb_lw%val(2,index,lblok),kind=dpc) * qfac
1742 : end if
1743 :
1744 : end do
1745 : end do
1746 : end do
1747 : end do
1748 :
1749 : !Calculate the Berry curvature of the induced Zeeman fields
1750 3774 : bc_sp= -matmul(bc_ss,barmmom) - matmul(invbarmagsus,bc_barsp)
1751 :
1752 2 : if (prtvol > 1) then
1753 2 : if (iblok /= 0) then
1754 6 : call wrtout([ab_out,std_out], ' Fixed-spin Berry curvature of the Zeeman fields induced by atomic displacements')
1755 6 : call wrtout([ab_out,std_out], ' atom1 dir E-dir Real Imag')
1756 10 : do irow=1, ndim
1757 154 : do icol=1, natom*3
1758 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)' ) &
1759 144 : & indexat1(irow), cart(indexdir1(irow)), indexat2(icol), cart(indexdir2(icol)), &
1760 288 : & real(bc_sp(irow,icol)), aimag(bc_sp(irow,icol))
1761 440 : call wrtout([ab_out,std_out], msg)
1762 : end do
1763 : end do
1764 6 : call wrtout([ab_out,std_out], ' ')
1765 : end if
1766 2 : if (jblok /= 0) then
1767 3 : call wrtout([ab_out,std_out], ' Fixed-spin Berry curvature of the Zeeman fields induced by electric field')
1768 3 : call wrtout([ab_out,std_out], ' atom1 dir E-dir Real Imag')
1769 5 : do irow=1, ndim
1770 17 : do icol=(natom+1)*3+1, (natom+2)*3
1771 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
1772 12 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
1773 24 : & real(bc_sp(irow,icol)), aimag(bc_sp(irow,icol))
1774 40 : call wrtout([ab_out,std_out], msg)
1775 : end do
1776 : end do
1777 3 : call wrtout([ab_out,std_out], ' ')
1778 : end if
1779 2 : if (lblok /= 0) then
1780 3 : call wrtout([ab_out,std_out], ' Fixed-spin Berry curvature of the Zeeman fields induced by macroscopic Zeeman field')
1781 3 : call wrtout([ab_out,std_out], ' atom1 dir B-dir Real Imag')
1782 5 : do irow=1, ndim
1783 17 : do icol=(natom+4)*3+1, (natom+5)*3
1784 : write(msg,'(i4,4x,a2,4x,a2,6x,2es18.9)' ) &
1785 12 : & indexat1(irow), cart(indexdir1(irow)), cart(indexdir2(icol)), &
1786 24 : & real(bc_sp(irow,icol)), aimag(bc_sp(irow,icol))
1787 40 : call wrtout([ab_out,std_out], msg)
1788 : end do
1789 : end do
1790 3 : call wrtout([ab_out,std_out], ' ')
1791 : end if
1792 : end if
1793 :
1794 : !Store the FM flavor in the DDB object
1795 24 : do ipert2=1,natom+5
1796 90 : do idir2=1,3
1797 66 : icol=idir2+(ipert2-1)*3
1798 66 : indexat2(icol)=ipert2
1799 66 : indexdir2(icol)=idir2
1800 :
1801 66 : ipert1_red= 0
1802 220 : do iat1= mpatpol(1), mpatpol(2)
1803 132 : ipert1= natom + 11 + iat1
1804 132 : ipert1_red= ipert1_red + 1
1805 132 : idir1_red= 0
1806 594 : do idir1= 1, 3
1807 396 : if (mpdir(idir1)==0) cycle
1808 264 : idir1_red= idir1_red + 1
1809 264 : irow=idir1_red+(ipert1_red-1)*nmdir
1810 : index = idir1 + &
1811 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
1812 264 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
1813 : jndex = idir2 + &
1814 : & 3*((ipert2 - 1) + mpert*((idir1 - 1) + &
1815 264 : & 3*((ipert1 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
1816 :
1817 396 : if (iblok /=0 .and. ipert2 <= natom) then
1818 144 : ddb_lw%val_fs(1,index,iblok)=real(bc_sp(irow,icol))
1819 144 : ddb_lw%val_fs(2,index,iblok)=aimag(bc_sp(irow,icol))
1820 144 : ddb_lw%val_fs(1,jndex,iblok)=real(bc_sp(irow,icol))
1821 144 : ddb_lw%val_fs(2,jndex,iblok)=-aimag(bc_sp(irow,icol))
1822 120 : else if (jblok /=0 .and. ipert2 == natom+2) then
1823 12 : ddb_lw%val_fs(1,index,jblok)=real(bc_sp(irow,icol))
1824 12 : ddb_lw%val_fs(2,index,jblok)=aimag(bc_sp(irow,icol))
1825 12 : ddb_lw%val_fs(1,jndex,jblok)=real(bc_sp(irow,icol))
1826 12 : ddb_lw%val_fs(2,jndex,jblok)=-aimag(bc_sp(irow,icol))
1827 108 : else if (lblok /=0 .and. ipert2 == natom+5) then
1828 12 : ddb_lw%val_fs(1,index,lblok)=real(bc_sp(irow,icol))
1829 12 : ddb_lw%val_fs(2,index,lblok)=aimag(bc_sp(irow,icol))
1830 12 : ddb_lw%val_fs(1,jndex,lblok)=real(bc_sp(irow,icol))
1831 12 : ddb_lw%val_fs(2,jndex,lblok)=-aimag(bc_sp(irow,icol))
1832 : end if
1833 :
1834 : end do
1835 : end do
1836 : end do
1837 : end do
1838 :
1839 2 : end subroutine berrycurv_sp
1840 : !!***
1841 :
1842 : !!****f* m_ddb_magpen/berrycurv_pp
1843 : !! NAME
1844 : !! berrycurv_pp
1845 : !!
1846 : !! FUNCTION
1847 : !! Calculate the Berry curvature of the phonon-phonon Hessian
1848 : !! (at constrained magnetic moments)
1849 : !!
1850 : !! INPUTS
1851 : !! blkval(2,3*mpert*3*mpert*3*mpert,nblok)= Third-order derivative matrices
1852 : !! In our case, the nblok is restricted to iblok
1853 : !! kblok= index of the current block
1854 : !! invbarmagsus(ndim,ndim)= Inverse of the penalized spin-sussceptibility tensor
1855 : !! mpatpol(2) = Atoms on which the magnetic penalty has been applied
1856 : !! mpdir(3) = Directions along which the spin-degrees of freedom have been stiffened
1857 : !! mpert =maximum number of ipert
1858 : !! natom= number of atoms in unit cell
1859 : !! nblok= number of blocks in the DDB
1860 : !! nmdir= number of directions along which the magnetic penalty was applied
1861 : !! ndim= dimension of the square susceptibilities
1862 : !! prtvol= control the volume of information written on output
1863 : !!
1864 : !! OUTPUT
1865 : !!
1866 : !! SOURCE
1867 :
1868 2 : subroutine berrycurv_pp(barmagsus,bc_barmagsus,bc_sp,ddb_lw,kblok,&
1869 2 : & mpert,natom,ndim,qeq0,ucvol,zfield)
1870 :
1871 : !Arguments -------------------------------
1872 : !scalars
1873 : integer,intent(in) :: kblok,mpert,natom,ndim
1874 : real(dp),intent(in) :: ucvol
1875 : logical,intent(in) :: qeq0
1876 : !arrays
1877 : type(ddb_type),intent(inout) :: ddb_lw
1878 : complex(dpc),intent(in) :: barmagsus(ndim,ndim)
1879 : complex(dpc),intent(in) :: bc_barmagsus(ndim,ndim)
1880 : complex(dpc),intent(in) :: bc_sp(ndim,(natom+5)*3)
1881 : complex(dpc),intent(in) :: zfield(ndim,(natom+5)*3)
1882 : !Local variables -------------------------
1883 : !scalars
1884 : integer :: icol,idir1,idir2,idir3,index,ipert1,ipert2,ipert3,irow
1885 : complex(dpc), parameter :: ione=(0.d0,1.d0)
1886 : complex(dpc) :: cval
1887 : !arrays
1888 4 : complex(dpc) :: bc_barpp((natom+5)*3,(natom+5)*3), bc_pp((natom+5)*3,(natom+5)*3)
1889 4 : complex(dpc) :: term((natom+5)*3,(natom+5)*3,3)
1890 :
1891 : ! *********************************************************************
1892 :
1893 : !Extract the frequency derivative of the penalized 2nd order nonmagnetic quantities
1894 2 : ipert3= natom + 9
1895 2 : idir3= 1
1896 24 : do ipert2= 1, natom+5
1897 90 : do idir2= 1, 3
1898 66 : icol=( ipert2-1)*3 + idir2
1899 814 : do ipert1= 1, natom+5
1900 2970 : do idir1= 1, 3
1901 2178 : irow=( ipert1-1)*3 + idir1
1902 : index = idir1 + &
1903 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
1904 2178 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
1905 :
1906 2904 : bc_barpp(irow,icol)= cmplx(ddb_lw%val(1,index,kblok),ddb_lw%val(2,index,kblok),kind=dpc)
1907 :
1908 : end do
1909 : end do
1910 : end do
1911 : end do
1912 :
1913 : !Calculate the different terms entering the Berry curvature
1914 23646 : term(:,:,1)=matmul(transpose(conjg(bc_sp)),matmul(barmagsus,zfield))
1915 23644 : term(:,:,2)=matmul(transpose(conjg(zfield)),matmul(bc_barmagsus,zfield))
1916 23644 : term(:,:,3)=matmul(transpose(conjg(zfield)),matmul(barmagsus,bc_sp))
1917 :
1918 2246 : bc_pp(:,:)= bc_barpp(:,:) + term(:,:,1) + term(:,:,2) + term(:,:,3)
1919 :
1920 : !Apply factors to convert derivatives of energy into susceptibilities
1921 : !Born charges
1922 2 : if (qeq0) then
1923 7 : ipert1= natom + 2
1924 7 : do ipert2= 1, natom
1925 25 : do idir2= 1, 3
1926 18 : icol=(ipert2-1)*3 + idir2
1927 78 : do idir1= 1, 3
1928 54 : irow=(ipert1-1)*3 + idir1
1929 54 : cval=bc_pp(irow,icol)
1930 54 : bc_pp(irow,icol)=-one*cval
1931 54 : cval=bc_pp(icol,irow)
1932 72 : bc_pp(icol,irow)=-one*cval
1933 : end do
1934 : end do
1935 : end do
1936 : end if
1937 :
1938 : !MR: Not sure why this is different from the electrical charges, | --------------------------------------------------------------------------------
1939 : !but is needed to reproduce freqflag=3 numbers
1940 : !Magnetic charges induced by atomic displacement
1941 : ! ipert1= natom + 5
1942 : ! do ipert2= 1, natom
1943 : ! do idir2= 1, 3
1944 : ! icol=( ipert2-1)*3 + idir2
1945 : ! do idir1= 1, 3
1946 : ! irow=( ipert1-1)*3 + idir1
1947 : ! cval=bc_pp(irow,icol)
1948 : ! bc_pp(irow,icol)=-cval
1949 : ! cval=bc_pp(icol,irow)
1950 : ! bc_pp(icol,irow)=-cval
1951 : ! end do
1952 : ! end do
1953 : ! end do
1954 :
1955 : !Dielectric tensor
1956 2 : if (qeq0) then
1957 4 : ipert1= natom + 2
1958 4 : ipert2= natom + 2
1959 4 : do idir2= 1, 3
1960 3 : icol=( ipert2-1)*3 + idir2
1961 13 : do idir1= 1, 3
1962 9 : irow=( ipert1-1)*3 + idir1
1963 9 : cval=bc_pp(irow,icol)
1964 12 : bc_pp(irow,icol)=-four_pi/ucvol*cval
1965 : end do
1966 : end do
1967 : end if
1968 :
1969 : !Magnetoelectric susceptibility
1970 2 : if (qeq0) then
1971 4 : ipert1= natom + 5
1972 4 : ipert2= natom + 2
1973 4 : do idir2= 1, 3
1974 3 : icol=( ipert2-1)*3 + idir2
1975 13 : do idir1= 1, 3
1976 9 : irow=( ipert1-1)*3 + idir1
1977 9 : cval=bc_pp(irow,icol)
1978 9 : bc_pp(irow,icol)=-cval/ucvol
1979 9 : cval=bc_pp(icol,irow)
1980 12 : bc_pp(icol,irow)=-cval/ucvol
1981 : end do
1982 : end do
1983 : end if
1984 :
1985 : !Magnetic susceptibility
1986 8 : ipert1= natom + 5
1987 8 : ipert2= natom + 5
1988 8 : do idir2= 1, 3
1989 6 : icol=( ipert2-1)*3 + idir2
1990 26 : do idir1= 1, 3
1991 18 : irow=( ipert1-1)*3 + idir1
1992 18 : cval=bc_pp(irow,icol)
1993 24 : bc_pp(irow,icol)=-cval/ucvol
1994 : end do
1995 : end do
1996 :
1997 : !Store the FM flavor in the DDB object
1998 24 : do ipert2= 1, natom+5
1999 90 : do idir2= 1, 3
2000 66 : icol=( ipert2-1)*3 + idir2
2001 814 : do ipert1= 1, natom+5
2002 2970 : do idir1= 1, 3
2003 2178 : irow=( ipert1-1)*3 + idir1
2004 : index = idir1 + &
2005 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
2006 2178 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
2007 2178 : ddb_lw%val_fs(1,index,kblok)=real(bc_pp(irow,icol))
2008 2904 : ddb_lw%val_fs(2,index,kblok)=aimag(bc_pp(irow,icol))
2009 : end do
2010 : end do
2011 : end do
2012 : end do
2013 :
2014 2 : end subroutine berrycurv_pp
2015 : !!***
2016 :
2017 : !!****f* m_ddb_magpen/mp_d3etot_print
2018 : !! NAME
2019 : !! mp_d3etot_print
2020 : !!
2021 : !! FUNCTION
2022 : !! Write on output file the fixed- and relaxed-spin susceptibilities
2023 : !! Only macroscopic quantities are printed if prtvol=1
2024 : !!
2025 : !! INPUTS
2026 : !! blkval= 2nd-order susceptibilities matrix
2027 : !! kblok= index of the current block
2028 : !! opt= 1 write the frozen-spin second-order quantities
2029 : !! 2 write the relaxed-spin second-order quantities
2030 : !! omega= frequency of the perturbation
2031 : !! qphon= momentum wave-vector
2032 : !! ucvol= unit-cell volume
2033 : !!
2034 : !! OUTPUT
2035 : !!
2036 : !! SOURCE
2037 :
2038 6 : subroutine mp_d3etot_print(ddb_lw,blkval,kblok,mpert,natom,nblok,opt,omega,prtvol,qeq0,qphnrm,qphon,ucvol,msize)
2039 :
2040 : !Arguments -------------------------------
2041 : !scalars
2042 : class(ddb_type),intent(in) :: ddb_lw
2043 : integer,intent(in) :: kblok,mpert,natom,nblok,opt,prtvol,msize
2044 : logical,intent(in) :: qeq0
2045 : real(dp),intent(in) :: ucvol
2046 : !arrays
2047 : real(dp),intent(in) :: omega(3)
2048 : real(dp),intent(in) :: blkval(2,msize,nblok)
2049 : real(dp),intent(inout) :: qphnrm(3),qphon(3,3)
2050 :
2051 : !Local variables -------------------------
2052 : !scalars
2053 : integer :: iblok,idir1,idir2,idir3,ipert1,ipert2,ipert3,irow,icol
2054 : integer :: rftyp,index,tmp
2055 : character(len=1000) :: msg
2056 : !arrays
2057 : integer :: rfelfd(4),rfmagn(4),rfphon(4),rfstrs(4),rffreq(4)
2058 : real(dp) :: val(2)
2059 : character(len=1) :: cart(3)=(/'x','y','z'/)
2060 :
2061 : ! *********************************************************************
2062 :
2063 2 : rfelfd(:)=0
2064 2 : rfphon(:)=0
2065 2 : rfstrs(:)=0
2066 2 : rfmagn(:)=0
2067 2 : rffreq(:)=0
2068 2 : rffreq(3)=1
2069 2 : rftyp = 33
2070 :
2071 2 : ipert3= natom + 9
2072 2 : idir3= 1
2073 2 : tmp=opt !
2074 : !IFCs
2075 2 : if (prtvol>1) then
2076 6 : rfphon(1:2)=1
2077 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
2078 2 : & rffreq=rffreq)
2079 2 : if (iblok/=0.and.iblok==kblok) then
2080 : ! if (opt==1) then
2081 6 : call wrtout([ab_out,std_out], ' Frozen-spin Berry curvature of interatomic force constants')
2082 : ! else if (opt==2) then
2083 : ! call wrtout([ab_out,std_out], ' Relaxed-spin Berry curvature of interatomic force constants')
2084 : ! end if
2085 6 : call wrtout([ab_out,std_out], ' atom1 dir atom2 dir Real Imag')
2086 14 : do ipert1= 1, natom
2087 50 : do idir1= 1, 3
2088 36 : irow=( ipert1-1)*3 + idir1
2089 252 : do ipert2= 1, natom
2090 900 : do idir2= 1, 3
2091 648 : icol=( ipert2-1)*3 + idir2
2092 : index = idir1 + &
2093 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
2094 648 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
2095 1944 : val(:)=blkval(:,index,kblok)
2096 : write(msg,'(2(i4,4x,a2,2x),2x,2es18.9)') &
2097 648 : & ipert1, cart(idir1), ipert2, cart(idir2), val(1), val(2)
2098 2160 : call wrtout([ab_out,std_out], msg)
2099 : end do
2100 : end do
2101 120 : call wrtout([ab_out,std_out], ' ')
2102 : end do
2103 : end do
2104 : end if
2105 : end if
2106 :
2107 2 : if (qeq0) then
2108 :
2109 : !Born charges
2110 3 : rfphon(1:2)=1
2111 3 : rfelfd(1:2)=2
2112 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
2113 1 : & rffreq=rffreq)
2114 1 : if (iblok/=0.and.iblok==kblok) then
2115 : ! if (opt==1) then
2116 3 : call wrtout([ab_out,std_out], ' Frozen-spin Berry curvature of Born effective charges')
2117 : ! else if (opt==2) then
2118 : ! call wrtout([ab_out,std_out], ' Relaxed-spin Berry curvature of Born effective charges')
2119 : ! end if
2120 3 : call wrtout([ab_out,std_out], ' E-dir atom dir Real Imag')
2121 1 : ipert1= natom + 2
2122 4 : do idir1= 1, 3
2123 21 : do ipert2= 1, natom
2124 75 : do idir2= 1, 3
2125 : index = idir1 + &
2126 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
2127 54 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
2128 162 : val(:)=blkval(:,index,kblok)
2129 : write(msg,'(3x,a2,7x,i3,4x,a2,2x,2es18.9)') &
2130 54 : & cart(idir1), ipert2, cart(idir2), val(1), val(2)
2131 180 : call wrtout([ab_out,std_out], msg)
2132 : end do
2133 : end do
2134 10 : call wrtout([ab_out,std_out], ' ')
2135 : end do
2136 : end if
2137 :
2138 : !Dielectric tensor
2139 : iblok=0
2140 1 : rfphon(:)=0
2141 3 : rfelfd(1:2)=2
2142 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
2143 1 : & rffreq=rffreq)
2144 1 : if (iblok/=0.and.iblok==kblok) then
2145 : ! if (opt==1) then
2146 3 : call wrtout([ab_out,std_out], ' Frozen-spin Berry curvature of clamped-ion dielectric tensor')
2147 : ! else if (opt==2) then
2148 : ! call wrtout([ab_out,std_out], ' Relaxed-spin Berry curvature of clamped-ion dielectric tensor')
2149 : ! end if
2150 3 : call wrtout([ab_out,std_out], ' dir dir Real Imag')
2151 1 : ipert1= ddb_lw%natom + 2
2152 1 : ipert2= ddb_lw%natom + 2
2153 4 : do idir2= 1, 3
2154 12 : do idir1= 1, 3
2155 : index = idir1 + &
2156 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
2157 9 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
2158 27 : val(:)=blkval(:,index,kblok)
2159 9 : write(msg,'(2x,a2,3x,a2,2x,2es18.9)' ) cart(idir1), cart(idir2), &
2160 18 : & val(1), val(2)
2161 30 : call wrtout([ab_out,std_out], msg)
2162 : end do
2163 10 : call wrtout([ab_out,std_out], ' ')
2164 : end do
2165 : end if
2166 :
2167 : !Magnetoelectric susceptibility
2168 : iblok=0
2169 1 : rfphon(:)=0
2170 1 : rfelfd(1)=0
2171 1 : rfelfd(2)=2
2172 1 : rfmagn(1)=1
2173 : rfmagn(2)=0
2174 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
2175 1 : & rfmagn=rfmagn,rffreq=rffreq)
2176 1 : if (iblok/=0.and.iblok==kblok) then
2177 : ! if (opt==1) then
2178 3 : call wrtout([ab_out,std_out], ' Frozen-spin Berry curvature of clamped-ion magnetoelectric susceptibility')
2179 : ! else if (opt==2) then
2180 : ! call wrtout([ab_out,std_out], ' Relaxed-spin Berry curvature of clamped-ion magnetoelectric susceptibility')
2181 : ! end if
2182 3 : call wrtout([ab_out,std_out], ' M-dir E-dir Real Imag')
2183 1 : ipert1= ddb_lw%natom + 5
2184 1 : ipert2= ddb_lw%natom + 2
2185 4 : do idir2= 1, 3
2186 12 : do idir1= 1, 3
2187 : index = idir1 + &
2188 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
2189 9 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
2190 27 : val(:)=blkval(:,index,kblok)/ucvol
2191 9 : write(msg,'(2x,a2,3x,a2,2x,2es18.9)' ) cart(idir1), cart(idir2), &
2192 18 : & val(1), val(2)
2193 30 : call wrtout([ab_out,std_out], msg)
2194 : end do
2195 10 : call wrtout([ab_out,std_out], ' ')
2196 : end do
2197 : end if
2198 :
2199 : end if
2200 :
2201 : !Magnetic susceptibility
2202 : iblok=0
2203 2 : rfphon(:)=0
2204 2 : rfelfd(:)=0
2205 2 : rfmagn(1)=1
2206 2 : rfmagn(2)=1
2207 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
2208 2 : & rfmagn=rfmagn,rffreq=rffreq)
2209 2 : if (iblok/=0.and.iblok==kblok) then
2210 : ! if (opt==1) then
2211 6 : call wrtout([ab_out,std_out], ' Frozen-spin Berry curvature of clamped-ion magnetic susceptibility')
2212 : ! else if (opt==2) then
2213 : ! call wrtout([ab_out,std_out], ' Relaxed-spin Berry curvature of clamped-ion magnetic susceptibility')
2214 : ! end if
2215 6 : call wrtout([ab_out,std_out], ' dir dir Real Imag')
2216 2 : ipert1= ddb_lw%natom + 5
2217 2 : ipert2= ddb_lw%natom + 5
2218 8 : do idir2= 1, 3
2219 24 : do idir1= 1, 3
2220 : index = idir1 + &
2221 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
2222 18 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
2223 54 : val(:)=blkval(:,index,kblok)/ucvol
2224 18 : write(msg,'(2x,a2,3x,a2,2x,2es18.9)' ) cart(idir1), cart(idir2), &
2225 36 : & val(1), val(2)
2226 60 : call wrtout([ab_out,std_out], msg)
2227 : end do
2228 20 : call wrtout([ab_out,std_out], ' ')
2229 : end do
2230 : end if
2231 :
2232 : !Magnetic Born effective charges
2233 : iblok=0
2234 2 : rfphon(1)=1
2235 2 : rfelfd(:)=0
2236 2 : rfmagn(1)=0
2237 : rfmagn(2)=1
2238 : call ddb_lw%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, omega=omega, &
2239 2 : & rfmagn=rfmagn,rffreq=rffreq)
2240 2 : if (iblok/=0.and.iblok==kblok) then
2241 : ! if (opt==1) then
2242 6 : call wrtout([ab_out,std_out], ' Frozen-spin Berry curvature of magnetic Born effective charges')
2243 : ! else if (opt==2) then
2244 : ! call wrtout([ab_out,std_out], ' Relaxed-spin Berry curvature of magnetic Born effective charges')
2245 : ! end if
2246 6 : call wrtout([ab_out,std_out], ' atom dir B-dir Real Imag')
2247 :
2248 2 : ipert2= ddb_lw%natom + 5
2249 14 : do ipert1= 1, natom
2250 48 : do idir1= 1, 3
2251 156 : do idir2= 1, 3
2252 : index = idir1 + &
2253 : & 3*((ipert1 - 1) + mpert*((idir2 - 1) + &
2254 108 : & 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
2255 324 : val(:)=blkval(:,index,kblok)
2256 : write(msg,'(i3,4x,a2,7x,a2,2x,2es18.9)') &
2257 108 : & ipert1, cart(idir1), cart(idir2), val(1), val(2)
2258 360 : call wrtout([ab_out,std_out], msg)
2259 : end do
2260 : end do
2261 38 : call wrtout([ab_out,std_out], ' ')
2262 : end do
2263 : end if
2264 :
2265 2 : end subroutine mp_d3etot_print
2266 : !!***
2267 :
2268 41464 : end module m_ddb_magpen
2269 : !!***
|
