Branch data 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 : : !!***
|
