Branch data Line data Source code
1 : : !{\src2tex{textfont=tt}}
2 : : !!****m* ABINIT/m_ddb_omega_interpol
3 : : !! NAME
4 : : !! m_ddb_omega_interpol
5 : : !!
6 : : !! FUNCTION
7 : : !! Interpolate the nonadiabatic second-order susceptibilities
8 : : !! onto a fine frequency grid and incorporate the lattice-mediated contributions.
9 : : !!
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2024 ABINIT group (MR and MS)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! NOTES
17 : : !!
18 : : !! PARENTS
19 : : !!
20 : : !! CHILDREN
21 : : !!
22 : : !! SOURCE
23 : :
24 : : #if defined HAVE_CONFIG_H
25 : : #include "config.h"
26 : : #endif
27 : :
28 : : #include "abi_common.h"
29 : :
30 : : module m_ddb_omega_interpol
31 : :
32 : : use defs_basis
33 : : use m_abicore
34 : : use m_profiling_abi
35 : : use m_errors
36 : : use m_ddb
37 : : use m_ddb_magpen, only : local_spinsus, magmom, mp_d2etot, asrw0
38 : : use m_fstrings, only : itoa, sjoin
39 : : use m_macroave, only : POLINT
40 : : use m_io_tools, only : open_file
41 : : use m_cgtools, only : fxphas_seq
42 : : use m_dynmat, only : pheigvec_normalize,phdispl_from_eigvec
43 : : use m_numeric_tools, only : polcoe
44 : : use m_ddb_diel, only : alignph
45 : :
46 : : implicit none
47 : :
48 : : public :: ddb_omega_interpol ! Perform an interpolation of the derivatives calculated at constrained magnetic moments
49 : : ! and later on convert them into the physically relevant (spin-relaxed) ones at each value of interpolated omega.
50 : :
51 : : private
52 : :
53 : : ! *************************************************************************
54 : :
55 : : contains
56 : : !!***
57 : :
58 : : !!****f* m_ddb_omega_interpol/ddb_omega_interpol
59 : : !! NAME
60 : : !! ddb_omega_interpol
61 : : !!
62 : : !! FUNCTION
63 : : !! Interpolate over frequency the secon-order derivatives calculated at
64 : : !! constrained magnetic moments, latter on convert them into susceptibilities
65 : : !! calculated at fixed/relaxed ions/spins.
66 : : !!
67 : : !! INPUTS
68 : : !! ddb (INOUT) = ddb block datastructure
69 : : !! magpen = amplitude (in Ha) of the applied magnetic penalty
70 : : !! mpatpol(2) = Atoms on which the magnetic penalty has been applied
71 : : !! mpdir(3) = Directions along which the spin-degrees of freedom have been stiffened
72 : : !! mpert = maximum number of ipert
73 : : !! mpopt = 1 calculate the frozen-magnetic second-order quantities
74 : : !! 2 calculate the spin-relaxed second-order quantities
75 : : !! natom= number of atoms in unit cell
76 : : !! ntypat= number of atom types
77 : : !! ucvol= unit cell volume
78 : : !!
79 : : !! OUTPUT
80 : : !! ddb= ddb%val updated with the corrected second-order derivatives
81 : : !!
82 : : !! SOURCE
83 : :
84 : 3 : subroutine ddb_omega_interpol(amu,ddb,ddb_lw, &
85 : : & eta,outfilename_radix,magpen,mpatpol,mpdir,mpert,mpopt,natom, &
86 : 3 : & nomega,ntypat,omegaflag,omegamax,omegamin,prtvol,typat,ucvol,xred)
87 : :
88 : : !Arguments -------------------------------
89 : : !scalars
90 : : integer,intent(in) :: mpert,mpopt,natom,nomega,ntypat,omegaflag,prtvol
91 : : real(dp),intent(in) :: eta,magpen,omegamax,omegamin,ucvol
92 : : character(len=*),intent(in) :: outfilename_radix
93 : : !arrays
94 : : type(ddb_type),intent(inout) :: ddb,ddb_lw
95 : : integer,intent(in) :: mpatpol(2),mpdir(3),typat(natom)
96 : : real(dp),intent(in) :: amu(ntypat)
97 : : real(dp),intent(in) :: xred(3,natom)
98 : :
99 : : !Local variables -------------------------
100 : : !scalars
101 : : integer :: alpha_unit,alpha_hc_unit,diel_unit,fs2rs,i,iblok,ifound
102 : : integer :: ii,iishift,imode,iw,j,jw,ipert3,idir3
103 : : integer :: mmag_unit,mmom_unit,nblok,ndim
104 : : integer :: nmat,nmdir,nwcalc,optgb,phon_unit,prtopt
105 : : integer :: locmagsus_unit,zeff_unit,zfield_unit
106 : : real(dp) :: mcell,omegastp,convfac
107 : : character(len=5000) :: msg,pfmt
108 : : character(len=fnlen) :: alpha_filename,alpha_hc_filename
109 : : character(len=fnlen) :: diel_filename,mmag_filename,locmagsus_filename,mmom_filename
110 : : character(len=fnlen) :: phon_filename,zeff_filename,zfield_filename
111 : : complex(dpc) :: cplxvar,cplx_weta
112 : : logical :: qeq0
113 : : !arrays
114 : : real(dp) :: qphnrm(3),qphon(3,3)
115 : 3 : real(dp), allocatable :: dint_fsddb(:,:),int_fsddb(:,:,:),int_rsddb(:,:,:)
116 : 3 : real(dp), allocatable :: omega(:),omegacalc(:)
117 : 3 : real(dp), allocatable :: w0hessian(:,:),w0berry(:,:)
118 : 3 : real(dp), allocatable :: displ(:),eigvec(:)
119 : 3 : real(dp), allocatable :: mode_phonspec(:),phonspec(:)
120 : 3 : real(dp), allocatable :: coeffs(:,:,:),pc_in(:),pc_out(:)
121 : 3 : complex(dpc), allocatable :: dummysus(:,:)
122 : 3 : complex(dpc), allocatable :: dummysus1(:,:), dummysus2(:,:)
123 : 3 : complex(dpc), allocatable :: dummymom(:,:),dummymom_tr(:,:),zfield_tr(:,:)
124 : 3 : complex(dpc), allocatable :: phongreen(:,:)
125 : :
126 : : !Pointer-target arrays
127 : 3 : real(dp), allocatable, target :: phfrq(:,:)
128 : 3 : complex(dpc), allocatable, target :: invmagsus(:,:,:), lm_magsus(:,:,:), magsus(:,:,:)
129 : 3 : complex(dpc), allocatable, target :: mmom(:,:,:), mmom_tr(:,:,:)
130 : 3 : complex(dpc), allocatable, target :: zfield(:,:,:)
131 : 3 : complex(dpc), allocatable, target :: ci_alpha(:,:,:),lm_alpha(:,:,:),lm_alpha_nm(:,:,:,:)
132 : 3 : complex(dpc), allocatable, target :: ci_alpha_hc(:,:,:),lm_alpha_hc(:,:,:)
133 : 3 : complex(dpc), allocatable, target :: ci_localpha(:,:,:),lm_localpha(:,:,:)
134 : 3 : complex(dpc), allocatable, target :: ci_locchi(:,:,:),lm_locchi(:,:,:)
135 : 3 : complex(dpc), allocatable, target :: ci_epsilon(:,:,:),lm_epsilon(:,:,:),lm_epsilon_nm(:,:,:,:)
136 : 3 : complex(dpc), allocatable, target :: ci_mchi(:,:,:),lm_mchi(:,:,:),lm_mchi_nm(:,:,:,:)
137 : 3 : complex(dpc), allocatable, target :: modemm(:,:,:),modedisp(:,:,:),modezf(:,:,:)
138 : 3 : complex(dpc), allocatable, target :: modemeff(:,:,:),modezeff(:,:,:)
139 : :
140 : 3 : real(dp), pointer, contiguous :: phfrq_iw(:)
141 : 3 : complex(dpc), pointer, contiguous :: invmagsus_iw(:,:), lm_magsus_iw(:,:), magsus_iw(:,:)
142 : 3 : complex(dpc), pointer, contiguous :: mmom_iw(:,:), mmom_tr_iw(:,:)
143 : 3 : complex(dpc), pointer, contiguous :: zfield_iw(:,:)
144 : 3 : complex(dpc), pointer, contiguous :: ci_alpha_iw(:,:),lm_alpha_iw(:,:),lm_alpha_nm_iw(:,:,:)
145 : 3 : complex(dpc), pointer, contiguous :: ci_alpha_hc_iw(:,:),lm_alpha_hc_iw(:,:)
146 : 3 : complex(dpc), pointer, contiguous :: ci_localpha_iw(:,:),lm_localpha_iw(:,:)
147 : 3 : complex(dpc), pointer, contiguous :: ci_locchi_iw(:,:),lm_locchi_iw(:,:)
148 : 3 : complex(dpc), pointer, contiguous :: ci_epsilon_iw(:,:),lm_epsilon_iw(:,:),lm_epsilon_nm_iw(:,:,:)
149 : 3 : complex(dpc), pointer, contiguous :: ci_mchi_iw(:,:),lm_mchi_iw(:,:),lm_mchi_nm_iw(:,:,:)
150 : 3 : complex(dpc), pointer, contiguous :: modemm_iw(:,:),modedisp_iw(:,:),modezf_iw(:,:)
151 : 3 : complex(dpc), pointer, contiguous :: modemeff_iw(:,:),modezeff_iw(:,:)
152 : :
153 : : ! *********************************************************************
154 : :
155 [ + + ]: 246 : write(msg, '(2a,(80a),4a)' ) ch10,('=',ii=1,80),ch10,ch10,&
156 [ + - ]: 246 : ' Frequency interpolation of constrained DFPT quantities section ',ch10
157 [ + + ]: 9 : call wrtout([std_out, ab_out], msg)
158 : :
159 : : !Identify the calculated omegas
160 : 3 : nwcalc=ddb%nblok
161 [ - + + - : 9 : ABI_MALLOC(omegacalc,(nwcalc))
+ - + - -
+ - + ]
162 [ + + ]: 7 : omegacalc(:)=ddb%omega(1,:)
163 : :
164 : : !Define the omega discretization
165 : 3 : omegastp=(omegamax-omegamin)/(nomega-1)
166 : :
167 : 3 : prtopt=0
168 [ + - ]: 3 : if (magpen<zero) then
169 : : nmat= 1
170 [ + - ]: 3 : else if (magpen>zero) then
171 : 3 : nmat= mpatpol(2) - mpatpol(1) + 1
172 : : end if
173 [ + + ]: 12 : nmdir=sum(mpdir(:))
174 : 3 : ndim=nmat*nmdir
175 : 3 : fs2rs=1
176 : 3 : qphon=zero
177 [ + + ]: 12 : qphon(:,1)=ddb%qpt(1:3,1)
178 [ + + ]: 12 : qeq0=(sqrt(sum(qphon(:,1)**2))<tol8)
179 [ + + ]: 12 : qphnrm(:)=ddb%nrm(1,1)
180 : 3 : optgb=0
181 [ - + + - : 9 : ABI_MALLOC(omega,(nomega))
+ - + - -
+ - + ]
182 [ + - - + : 15 : ABI_MALLOC(phfrq,(3*natom,nomega))
+ - - + +
- + - + -
- + - + ]
183 [ - + + - : 6 : ABI_MALLOC(phonspec,(nomega))
- + - + ]
184 [ - + + - : 12 : ABI_MALLOC(phongreen,(3*natom,3*natom))
- + + - +
- + - - +
- + ]
185 [ - + + - : 9 : ABI_MALLOC(mode_phonspec,(3*natom))
+ - + - -
+ - + ]
186 [ + - - + : 9 : ABI_MALLOC(displ,(2*3*natom*3*natom))
+ - + - +
- - + -
+ ]
187 [ - + + - : 6 : ABI_MALLOC(eigvec,(2*3*natom*3*natom))
- + - + ]
188 [ - + + - : 15 : ABI_MALLOC(modemm,(ndim,3*natom,nomega))
- + + - -
+ + - + -
+ - - + -
+ ]
189 [ - + + - : 15 : ABI_MALLOC(modedisp,(3*natom,3*natom,nomega))
- + + - -
+ + - + -
+ - - + -
+ ]
190 [ - + + - : 12 : ABI_MALLOC(modezf,(ndim,3*natom,nomega))
- + + - -
+ + - - +
- + ]
191 [ - + + - : 9 : ABI_MALLOC(ci_alpha,(3,3,nomega))
+ - + - -
+ - + ]
192 [ - + + - : 6 : ABI_MALLOC(ci_alpha_hc,(3,3,nomega))
- + - + ]
193 [ - + + - : 6 : ABI_MALLOC(lm_alpha,(3,3,nomega))
- + - + ]
194 [ - + + - : 6 : ABI_MALLOC(lm_alpha_hc,(3,3,nomega))
- + - + ]
195 [ - + + - : 12 : ABI_MALLOC(lm_alpha_nm,(3,3,3*natom,nomega))
- + + - +
- + - - +
- + ]
196 [ - + + - : 12 : ABI_MALLOC(ci_localpha,(ndim,3,nomega))
- + + - +
- + - - +
- + ]
197 [ - + + - : 9 : ABI_MALLOC(lm_localpha,(ndim,3,nomega))
- + + - -
+ - + ]
198 [ - + + - : 9 : ABI_MALLOC(ci_locchi,(ndim,3,nomega))
- + + - -
+ - + ]
199 [ - + + - : 9 : ABI_MALLOC(lm_locchi,(ndim,3,nomega))
- + + - -
+ - + ]
200 [ - + + - : 6 : ABI_MALLOC(ci_epsilon,(3,3,nomega))
- + - + ]
201 [ - + + - : 6 : ABI_MALLOC(lm_epsilon,(3,3,nomega))
- + - + ]
202 [ - + + - : 9 : ABI_MALLOC(lm_epsilon_nm,(3,3,3*natom,nomega))
- + + - -
+ - + ]
203 [ - + + - : 6 : ABI_MALLOC(ci_mchi,(3,3,nomega))
- + - + ]
204 [ - + + - : 6 : ABI_MALLOC(lm_mchi,(3,3,nomega))
- + - + ]
205 [ - + + - : 9 : ABI_MALLOC(lm_mchi_nm,(3,3,3*natom,nomega))
- + + - -
+ - + ]
206 [ - + + - : 12 : ABI_MALLOC(modemeff,(3,3*natom,nomega))
- + + - +
- + - - +
- + ]
207 [ - + + - : 9 : ABI_MALLOC(modezeff,(3,3*natom,nomega))
- + + - -
+ - + ]
208 [ - + + - : 12 : ABI_MALLOC(dummysus,(ndim,ndim))
- + + - +
- + - - +
- + ]
209 [ - + + - : 9 : ABI_MALLOC(dummysus1,(ndim,ndim))
- + + - -
+ - + ]
210 [ - + + - : 9 : ABI_MALLOC(dummysus2,(ndim,ndim))
- + + - -
+ - + ]
211 [ - + + - : 15 : ABI_MALLOC(magsus,(ndim,ndim,nomega))
- + + - -
+ + - + -
+ - - + -
+ ]
212 [ - + + - : 12 : ABI_MALLOC(lm_magsus,(ndim,ndim,nomega))
- + + - -
+ + - - +
- + ]
213 [ - + + - : 12 : ABI_MALLOC(invmagsus,(ndim,ndim,nomega))
- + + - -
+ + - - +
- + ]
214 [ - + + - : 12 : ABI_MALLOC(dummymom,(ndim,(natom+5)*3))
+ - - + +
- + - + -
- + - + ]
215 [ - + + - : 9 : ABI_MALLOC(dummymom_tr,(ndim,(natom+5)*3))
- + + - -
+ - + ]
216 [ - + + - : 15 : ABI_MALLOC(mmom,(ndim,(natom+5)*3,nomega))
- + + - -
+ + - + -
+ - - + -
+ ]
217 [ - + + - : 12 : ABI_MALLOC(mmom_tr,((natom+5)*3,ndim,nomega))
- + + - -
+ + - - +
- + ]
218 [ - + + - : 12 : ABI_MALLOC(zfield,(ndim,(natom+5)*3,nomega))
- + + - -
+ + - - +
- + ]
219 [ - + + - : 9 : ABI_MALLOC(zfield_tr,((natom+5)*3,ndim))
- + + - -
+ - + ]
220 [ - + + - : 9 : ABI_MALLOC(dint_fsddb,(2,ddb%msize))
+ - + - -
+ - + ]
221 [ - + + - : 6 : ABI_MALLOC(int_fsddb,(2,ddb%msize,1))
- + - + ]
222 [ - + + - : 6 : ABI_MALLOC(int_rsddb,(2,ddb%msize,1))
- + - + ]
223 : :
224 : : !Compute total cell mass
225 [ + + ]: 21 : mcell= sum(amu(typat(:)))*amu_emass
226 : :
227 : : !For linear interpolation detect the w=0 Hessians and Berry curvatures
228 [ + + ]: 3 : if (omegaflag == 1) then
229 : :
230 [ - + + - : 4 : ABI_MALLOC(w0hessian,(2,ddb%msize))
- + - + ]
231 : 4 : nblok= ddb%nblok
232 : 4 : ifound= 0
233 [ + + ]: 4 : do iblok= 1, nblok
234 [ + - ]: 4 : if (abs(ddb%omega(1,iblok)) < tol12) then
235 [ + + + + ]: 28568 : w0hessian(:,:)= ddb%val_fs(:,:,iblok)
236 : : ifound= 1
237 : : end if
238 : : end do
239 [ - + ]: 2 : if (ifound==0) then
240 : 0 : write(msg, '(3a)' )' No omega=0 block with second-order derivatives', &
241 : 0 : & ' found in the DDB file. This is necessary if omegaflag=1 ',ch10
242 : 0 : ABI_ERROR(msg)
243 : : end if
244 : :
245 [ - + + - : 6 : ABI_MALLOC(w0berry,(2,ddb_lw%msize))
+ - + - -
+ - + ]
246 : 2 : nblok= ddb_lw%nblok
247 : 2 : ifound= 0
248 [ + + ]: 6 : do iblok= 1, nblok
249 [ + - ]: 6 : if (abs(ddb_lw%omega(1,iblok)) < tol12) then
250 [ + + + + ]: 3942112 : w0berry(:,:)= ddb_lw%val_fs(:,:,iblok)
251 : : ifound= 1
252 : : end if
253 : : end do
254 [ - + ]: 2 : if (ifound==0) then
255 : 0 : write(msg, '(3a)' )' No omega=0 block with third-order derivatives', &
256 : 0 : & ' found in the DDB file. This is necessary if omegaflag=1 ',ch10
257 : 0 : ABI_ERROR(msg)
258 : : end if
259 : :
260 : : !Index shift for Berry curvatures
261 : 2 : ipert3= natom + 9
262 : 2 : idir3= 1
263 : 2 : iishift= 3*(mpert*(3*(mpert*((idir3 - 1) + 3*(ipert3 - 1)))))
264 : : end if
265 : :
266 : : !For nwcalc-1 Taylor-expansion interpolation precalculate the coefficients
267 [ + + ]: 3 : if (omegaflag == 3) then
268 [ - + + - : 4 : ABI_MALLOC(coeffs,(2,nwcalc,ddb%msize))
- + + - +
- + - - +
- + ]
269 [ - + + - : 2 : ABI_MALLOC(pc_in,(nwcalc))
- + - + ]
270 [ - + + - : 2 : ABI_MALLOC(pc_out,(nwcalc))
- + - + ]
271 [ + + ]: 4762 : do ii=1,ddb%msize
272 [ + + + + : 8010 : if (all(ddb%flg(ii,:)==1)) then
+ + ]
273 [ + + ]: 4872 : pc_in(:)= ddb%val_fs(1,ii,:)
274 : 1624 : call polcoe(omegacalc,pc_in,nwcalc,pc_out)
275 [ + + ]: 4872 : coeffs(1,:,ii)= pc_out(:)
276 [ + + ]: 4872 : pc_in(:)= ddb%val_fs(2,ii,:)
277 : 1624 : call polcoe(omegacalc,pc_in,nwcalc,pc_out)
278 [ + + ]: 4872 : coeffs(2,:,ii)= pc_out(:)
279 : : end if
280 : : end do
281 [ - + ]: 1 : ABI_FREE(pc_in)
282 [ - + ]: 1 : ABI_FREE(pc_out)
283 : : end if
284 [ - + ]: 3 : if (omegaflag == 2) then
285 [ # # # # : 0 : ABI_MALLOC(pc_in,(nwcalc))
# # # # ]
286 : : endif
287 : :
288 : :
289 : : !Loop over the frequency
290 [ + + ]: 33 : do iw=1,nomega
291 : 30 : omega(iw)=omegamin+omegastp*(iw-1)
292 [ - + ]: 30 : if (nomega==1) omega(iw)=omegamin
293 : :
294 : : !Pointer target associations for iw case
295 : 30 : phfrq_iw => phfrq(:,iw)
296 : 30 : invmagsus_iw => invmagsus(:,:,iw)
297 : 30 : lm_magsus_iw => lm_magsus(:,:,iw)
298 : 30 : magsus_iw => magsus(:,:,iw)
299 : 30 : mmom_iw => mmom(:,:,iw)
300 : 30 : mmom_tr_iw => mmom_tr(:,:,iw)
301 : 30 : zfield_iw => zfield(:,:,iw)
302 : 30 : ci_alpha_iw => ci_alpha(:,:,iw)
303 : 30 : lm_alpha_iw => lm_alpha(:,:,iw)
304 : 30 : lm_alpha_nm_iw => lm_alpha_nm(:,:,:,iw)
305 : 30 : ci_alpha_hc_iw => ci_alpha_hc(:,:,iw)
306 : 30 : lm_alpha_hc_iw => lm_alpha_hc(:,:,iw)
307 : 30 : ci_localpha_iw => ci_localpha(:,:,iw)
308 : 30 : lm_localpha_iw => lm_localpha(:,:,iw)
309 : 30 : ci_locchi_iw => ci_locchi(:,:,iw)
310 : 30 : lm_locchi_iw => lm_locchi(:,:,iw)
311 : 30 : ci_epsilon_iw => ci_epsilon(:,:,iw)
312 : 30 : lm_epsilon_iw => lm_epsilon(:,:,iw)
313 : 30 : lm_epsilon_nm_iw => lm_epsilon_nm(:,:,:,iw)
314 : 30 : ci_mchi_iw => ci_mchi(:,:,iw)
315 : 30 : lm_mchi_iw => lm_mchi(:,:,iw)
316 : 30 : lm_mchi_nm_iw => lm_mchi_nm(:,:,:,iw)
317 : 30 : modemm_iw => modemm(:,:,iw)
318 : 30 : modedisp_iw => modedisp(:,:,iw)
319 : 30 : modezf_iw => modezf(:,:,iw)
320 : 30 : modemeff_iw => modemeff(:,:,iw)
321 : 30 : modezeff_iw => modezeff(:,:,iw)
322 : :
323 : : !Perform the different interpolations
324 : : !Lineal (with analytic Berry curvature) with dissipation if eta/=0
325 [ + + ]: 30 : if (omegaflag==1) then
326 [ + + ]: 95240 : do ii=1,ddb%msize
327 [ + + + + : 125230 : if (all(ddb%flg(ii,:)==1)) then
+ + ]
328 : 29990 : int_fsddb(1,ii,1)= w0hessian(1,ii) + omega(iw)*w0berry(1,ii+iishift) - eta*w0berry(2,ii+iishift)
329 : 29990 : int_fsddb(2,ii,1)= w0hessian(2,ii) + omega(iw)*w0berry(2,ii+iishift) + eta*w0berry(1,ii+iishift)
330 : : end if
331 : : end do
332 : :
333 : : !Polynomial with no dissipation
334 [ - + ]: 10 : else if (omegaflag==2) then
335 [ # # ]: 0 : do ii=1,ddb%msize
336 [ # # # # : 0 : if (all(ddb%flg(ii,:)==1)) then
# # ]
337 [ # # ]: 0 : pc_in(:)= ddb%val_fs(1,ii,:)
338 : 0 : call POLINT(omegacalc,pc_in,nwcalc,omega(iw),int_fsddb(1,ii,1),dint_fsddb(1,ii))
339 [ # # ]: 0 : pc_in(:)= ddb%val_fs(2,ii,:)
340 : 0 : call POLINT(omegacalc,pc_in,nwcalc,omega(iw),int_fsddb(2,ii,1),dint_fsddb(2,ii))
341 [ # # # # : 0 : else if (count(ddb%flg(ii,:)==0)/=nwcalc) then
# # ]
342 : : write(msg,'(a,a,a)')&
343 : 0 : 'ddb_omega_interpol detects differences between the DDB bloks for each frequency.',ch10,&
344 : 0 : & ' The interpolation has been stopped.'
345 : 0 : ABI_ERROR(msg)
346 : : end if
347 : : end do
348 : :
349 : : !Taylor-expansion around w=0 with dissipation if eta/=0
350 [ + - ]: 10 : else if (omegaflag==3) then
351 : 10 : cplx_weta=cmplx(omega(iw),eta,kind=dpc)
352 [ + + ]: 47620 : do ii=1,ddb%msize
353 [ + + + + : 80100 : if (all(ddb%flg(ii,:)==1)) then
+ + ]
354 : 16240 : cplxvar=cmplx(zero,zero,kind=dpc)
355 [ + + ]: 48720 : do jw= 1, nwcalc
356 : 48720 : cplxvar= cplxvar + cplx_weta**(jw-1)*cmplx(coeffs(1,jw,ii),coeffs(2,jw,ii),kind=dpc)
357 : : end do
358 : 16240 : int_fsddb(1,ii,1)=real(cplxvar)
359 : 16240 : int_fsddb(2,ii,1)=aimag(cplxvar)
360 [ + + + - : 94110 : else if (count(ddb%flg(ii,:)==0)/=nwcalc) then
- + ]
361 : : write(msg,'(a,a,a)')&
362 : 0 : 'ddb_omega_interpol detects differences between the DDB bloks for each frequency.',ch10,&
363 : 0 : & ' The interpolation has been stopped.'
364 : 0 : ABI_ERROR(msg)
365 : : end if
366 : : end do
367 : : end if
368 : :
369 : : !Calculate the local spin susceptibilities
370 : : call local_spinsus(dummysus,ddb,1,dummysus1,invmagsus_iw,&
371 : : & dummysus2,magpen,magsus_iw,mpatpol,mpdir,mpert,natom,ndim,nmdir,prtopt,prtvol,qphon,xred,&
372 : 30 : & fs2rs=fs2rs,blkval_fs=int_fsddb)
373 : :
374 : : !Calculate the 1st-order magnetic moments
375 : : call magmom(dummymom,dummymom_tr,ddb,dummysus,dummysus,1,1,1,magpen,&
376 : : & magsus_iw,mmom_iw,mmom_tr_iw,mpatpol,mpdir,mpert,natom,&
377 : : & ndim,nmdir,prtopt,prtvol,qphon,xred,zfield_iw,zfield_tr,&
378 : 30 : & fs2rs=fs2rs,blkval_fs=int_fsddb)
379 : :
380 : : !Now calculate the non-magnetic second-order quantities
381 : 30 : call ddb%to_d2etot(int_fsddb,1,0,qeq0,qphon,qphnrm,ucvol,optgb,omega=omega(iw))
382 : :
383 : : call mp_d2etot(dummysus,ddb,1,magsus_iw,&
384 : : & mpert,mpopt,natom,ndim,zfield_iw,zfield_tr, &
385 : 30 : & fs2rs=fs2rs,blkval_fs=int_fsddb,blkval_rs=int_rsddb)
386 : :
387 : 30 : call ddb%to_d2etot(int_fsddb,1,1,qeq0,qphon,qphnrm,ucvol,optgb,omega=omega(iw))
388 [ + - ]: 30 : if (mpopt==2) call ddb%to_d2etot(int_rsddb,1,1,qeq0,qphon,qphnrm,ucvol,optgb,omega=omega(iw))
389 : :
390 : : !Calculate the phonon propagator (Green's function) and spectral function
391 : : !and the lattice-mediated contributions to the different susceptibilities.
392 [ - + ]: 33 : if (mpopt==1) then
393 : : call phonon_green(amu,displ,eigvec,eta,int_fsddb,&
394 : : & mode_phonspec,mpert,natom,ntypat,omega(iw),&
395 : 0 : & phfrq_iw,phongreen,phonspec(iw),typat)
396 : :
397 : : call ri_d2etot(int_fsddb,ci_alpha_iw,ci_alpha_hc_iw,ci_epsilon_iw,ci_localpha_iw,ci_locchi_iw,ci_mchi_iw,&
398 : : & lm_alpha_iw,lm_alpha_hc_iw,lm_epsilon_iw,lm_localpha_iw,lm_locchi_iw,lm_magsus_iw,lm_mchi_iw,&
399 : 0 : & mpert,mmom_iw,mmom_tr_iw,natom,ndim,phongreen,ucvol)
400 : :
401 : : call lm_normal_modes(amu,int_fsddb,displ,eta,lm_alpha_nm_iw,lm_epsilon_nm_iw,lm_mchi_nm_iw, &
402 : : & mcell,mmom_iw,modemm_iw,modedisp_iw,modemeff_iw,modezeff_iw,modezf_iw,&
403 : 0 : & mpert,natom,ndim,ntypat,omega(iw),phfrq_iw,typat,ucvol,zfield_iw)
404 : :
405 [ + - ]: 30 : else if (mpopt==2) then
406 : : call phonon_green(amu,displ,eigvec,eta,int_rsddb,&
407 : : & mode_phonspec,mpert,natom,ntypat,omega(iw),&
408 : 30 : & phfrq_iw,phongreen,phonspec(iw),typat)
409 : :
410 : : call ri_d2etot(int_rsddb,ci_alpha_iw,ci_alpha_hc_iw,ci_epsilon_iw,ci_localpha_iw,ci_locchi_iw,ci_mchi_iw,&
411 : : & lm_alpha_iw,lm_alpha_hc_iw,lm_epsilon_iw,lm_localpha_iw,lm_locchi_iw,lm_magsus_iw,lm_mchi_iw,&
412 : 30 : & mpert,mmom_iw,mmom_tr_iw,natom,ndim,phongreen,ucvol)
413 : :
414 : : call lm_normal_modes(amu,int_rsddb,displ,eta,lm_alpha_nm_iw,lm_epsilon_nm_iw,lm_mchi_nm_iw, &
415 : : & mcell,mmom_iw,modemm_iw,modedisp_iw,modemeff_iw,modezeff_iw,modezf_iw,&
416 : 30 : & mpert,natom,ndim,ntypat,omega(iw),phfrq_iw,typat,ucvol,zfield_iw)
417 : :
418 : : end if
419 : :
420 : : end do
421 : :
422 : : !!! Print results of interpolation
423 : : !Local magnetic susceptibilities
424 [ + - ]: 3 : locmagsus_filename=trim(outfilename_radix)//"_LOCMAGSUS"
425 [ - + ]: 3 : if (open_file(locmagsus_filename, msg, newunit=locmagsus_unit) /= 0) then
426 : 0 : ABI_ERROR(msg)
427 : : end if
428 : :
429 : 3 : write(locmagsus_unit,*) '#'
430 : 3 : write(locmagsus_unit,*) '# Local magnetic susceptibilities calculated and interpolated by ANADDB'
431 : 3 : write(locmagsus_unit,*) '#'
432 : 3 : write(pfmt, '( "(es15.7, ", I4, "(es15.7))" )' ) ndim**2
433 : :
434 : 3 : write(locmagsus_unit,*) '# Real part of clamped-ion local magnetic susceptibility tensor (at. units)'
435 : 3 : write(msg,'(a,a)') ch10,&
436 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
437 : 3 : call wrtout(locmagsus_unit,msg,'COLL')
438 [ + + ]: 33 : do iw=1,nomega
439 : : write(msg,pfmt) &
440 [ + - + + : 630 : & omega(iw), ((real(magsus(i,j,iw)),j=1,ndim),i=1,ndim)
+ - + + ]
441 : 33 : call wrtout(locmagsus_unit,msg,'COLL')
442 : : end do
443 : :
444 : 3 : write(locmagsus_unit,*) ' '
445 : 3 : write(locmagsus_unit,*) '# Imaginary part of clamped-ion local magnetic susceptibility tensor (at. units)'
446 : 3 : write(msg,'(a,a)') ch10,&
447 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
448 : 3 : call wrtout(locmagsus_unit,msg,'COLL')
449 [ + + ]: 33 : do iw=1,nomega
450 : : write(msg,pfmt) &
451 [ + - + + : 630 : & omega(iw), ((aimag(magsus(i,j,iw)),j=1,ndim),i=1,ndim)
+ - + + ]
452 : 33 : call wrtout(locmagsus_unit,msg,'COLL')
453 : : end do
454 : :
455 : 3 : write(locmagsus_unit,*) ' '
456 : 3 : write(locmagsus_unit,*) '# Real part of relaxed-ion local magnetic susceptibility tensor (at. units)'
457 : 3 : write(msg,'(a,a)') ch10,&
458 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
459 : 3 : call wrtout(locmagsus_unit,msg,'COLL')
460 [ + + ]: 33 : do iw=1,nomega
461 : : write(msg,pfmt) &
462 [ + - + + : 630 : & omega(iw), ((real(magsus(i,j,iw)+lm_magsus(i,j,iw)),j=1,ndim),i=1,ndim)
+ - + + ]
463 : 33 : call wrtout(locmagsus_unit,msg,'COLL')
464 : : end do
465 : :
466 : 3 : write(locmagsus_unit,*) ' '
467 : 3 : write(locmagsus_unit,*) '# Imaginary part of relaxed-ion local magnetic susceptibility tensor (at. units)'
468 : 3 : write(msg,'(a,a)') ch10,&
469 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
470 : 3 : call wrtout(locmagsus_unit,msg,'COLL')
471 [ + + ]: 33 : do iw=1,nomega
472 : : write(msg,pfmt) &
473 [ + - + + : 630 : & omega(iw), ((aimag(magsus(i,j,iw)+lm_magsus(i,j,iw)),j=1,ndim),i=1,ndim)
+ - + + ]
474 : 33 : call wrtout(locmagsus_unit,msg,'COLL')
475 : : end do
476 : :
477 : 3 : write(locmagsus_unit,*) ' '
478 : 3 : write(locmagsus_unit,*) '# Real part of the inverse of the clamped-ion local magnetic susceptibility tensor (at. units)'
479 : 3 : write(msg,'(a,a)') ch10,&
480 : 6 : & ' # At hw X^{-1}_11 X^{-1}_12 ... X^{-1}_21 X^{-1}_22 ...'
481 : 3 : call wrtout(locmagsus_unit,msg,'COLL')
482 [ + + ]: 33 : do iw=1,nomega
483 : : write(msg,pfmt) &
484 [ + - + + : 630 : & omega(iw), ((real(invmagsus(i,j,iw)),j=1,ndim),i=1,ndim)
+ - + + ]
485 : 33 : call wrtout(locmagsus_unit,msg,'COLL')
486 : : end do
487 : :
488 : 3 : write(locmagsus_unit,*) ' '
489 : 3 : write(locmagsus_unit,*) '# Imaginary part of the inverse of the clamped-ion local magnetic susceptibility tensor (at. units)'
490 : 3 : write(msg,'(a,a)') ch10,&
491 : 6 : & ' # At hw X^{-1}_11 X^{-1}_12 ... X^{-1}_21 X^{-1}_22 ...'
492 : 3 : call wrtout(locmagsus_unit,msg,'COLL')
493 [ + + ]: 33 : do iw=1,nomega
494 : : write(msg,pfmt) &
495 [ + - + + : 630 : & omega(iw), ((aimag(invmagsus(i,j,iw)),j=1,ndim),i=1,ndim)
+ - + + ]
496 : 33 : call wrtout(locmagsus_unit,msg,'COLL')
497 : : end do
498 : :
499 : 3 : close (locmagsus_unit)
500 : :
501 : : !Zfields
502 [ + - ]: 3 : zfield_filename=trim(outfilename_radix)//"_LOCZFIELDS"
503 [ - + ]: 3 : if (open_file(zfield_filename, msg, newunit=zfield_unit) /= 0) then
504 : 0 : ABI_ERROR(msg)
505 : : end if
506 : :
507 : 3 : write(pfmt, '( "(es15.7, ", I2, "(es17.7))" )' ) ndim
508 [ + + ]: 57 : do imode= 1, 3*natom
509 : 54 : write(zfield_unit,*) ' '
510 : 54 : write(zfield_unit,'(a,i3)') '# Real part of local Zeeman fields (at. units) induced by phonon mode:', imode
511 : 54 : write(msg,'(a,a,a)') ch10,&
512 : 108 : & ' # At hw Z_{mat_1,1} Z_{mat_1,2} ... Z_{mat_2,1} Z_{mat_2,2}'
513 : 54 : call wrtout(zfield_unit,msg,'COLL')
514 [ + + ]: 594 : do iw=1,nomega
515 : : write(msg,pfmt) &
516 [ + - + + ]: 2700 : & omega(iw), (real(modezf(i,imode,iw)),i=1,ndim)
517 : 594 : call wrtout(zfield_unit,msg,'COLL')
518 : : end do
519 : 54 : write(zfield_unit,*) ' '
520 : 54 : write(zfield_unit,'(a,i3)') '# Imaginary part of local Zeeman fields (at. units) induced by phonon mode:', imode
521 : 54 : write(msg,'(a,a,a)') ch10,&
522 : 108 : & ' # At hw Z_{mat_1,1} Z_{mat_1,2} ... Z_{mat_2,1} Z_{mat_2,2}'
523 : 54 : call wrtout(zfield_unit,msg,'COLL')
524 [ + + ]: 597 : do iw=1,nomega
525 : : write(msg,pfmt) &
526 [ + - + + ]: 2700 : & omega(iw), (aimag(modezf(i,imode,iw)),i=1,ndim)
527 : 594 : call wrtout(zfield_unit,msg,'COLL')
528 : : end do
529 : : end do
530 : :
531 : 3 : write(pfmt, '( "(es15.7, ", I4, "(es15.7))" )' ) ndim*3
532 [ + + ]: 3 : if (qeq0) then
533 : 2 : write(zfield_unit,*) ' '
534 : 2 : write(zfield_unit,*) '# Real part of clamped-ion local Zeeman fields induced by electric field (at. units)'
535 : 2 : write(msg,'(a,a)') ch10,&
536 : 4 : & ' # At hw Z_11 Z_12 Z_13 ... Z_21 Z_22 ...'
537 : 2 : call wrtout(zfield_unit,msg,'COLL')
538 [ + + ]: 22 : do iw=1,nomega
539 [ + - + + : 340 : write(msg,pfmt) omega(iw), ((real(zfield(i,(natom+1)*3+j,iw)),j=1,3),i=1,ndim)
+ - + + ]
540 : 22 : call wrtout(zfield_unit,msg,'COLL')
541 : : end do
542 : :
543 : 2 : write(zfield_unit,*) ' '
544 : 2 : write(zfield_unit,*) '# Imag part of clamped-ion local Zeeman fields induced by electric field (at. units)'
545 : 2 : write(msg,'(a,a)') ch10,&
546 : 4 : & ' # At hw Z_11 Z_12 Z_13 ... Z_21 Z_22 ...'
547 : 2 : call wrtout(zfield_unit,msg,'COLL')
548 [ + + ]: 22 : do iw=1,nomega
549 [ + - + + : 340 : write(msg,pfmt) omega(iw), ((aimag(zfield(i,(natom+1)*3+j,iw)),j=1,3),i=1,ndim)
+ - + + ]
550 : 22 : call wrtout(zfield_unit,msg,'COLL')
551 : : end do
552 : : end if
553 : :
554 : 3 : write(zfield_unit,*) ' '
555 : 3 : write(zfield_unit,*) '# Real part of clamped-ion local Zeeman fields induced by macroscopic Zeeman field (at. units)'
556 : 3 : write(msg,'(a,a)') ch10,&
557 : 6 : & ' # At hw Z_11 Z_12 Z_13 ... Z_21 Z_22 ...'
558 : 3 : call wrtout(zfield_unit,msg,'COLL')
559 [ + + ]: 33 : do iw=1,nomega
560 [ + - + + : 510 : write(msg,pfmt) omega(iw), ((real(zfield(i,(natom+4)*3+j,iw)),j=1,3),i=1,ndim)
+ - + + ]
561 : 33 : call wrtout(zfield_unit,msg,'COLL')
562 : : end do
563 : :
564 : 3 : write(zfield_unit,*) ' '
565 : 3 : write(zfield_unit,*) '# Imag part of clamped-ion local Zeeman fields induced by macroscopic Zeeman field (at. units)'
566 : 3 : write(msg,'(a,a)') ch10,&
567 : 6 : & ' # At hw Z_11 Z_12 Z_13 ... Z_21 Z_22 ...'
568 : 3 : call wrtout(zfield_unit,msg,'COLL')
569 [ + + ]: 33 : do iw=1,nomega
570 [ + - + + : 510 : write(msg,pfmt) omega(iw), ((aimag(zfield(i,(natom+4)*3+j,iw)),j=1,3),i=1,ndim)
+ - + + ]
571 : 33 : call wrtout(zfield_unit,msg,'COLL')
572 : : end do
573 : 3 : close (zfield_unit)
574 : :
575 : : !Magnetic moments
576 [ + - ]: 3 : mmom_filename=trim(outfilename_radix)//"_LOCMAGMOM"
577 [ - + ]: 3 : if (open_file(mmom_filename, msg, newunit=mmom_unit) /= 0) then
578 : 0 : ABI_ERROR(msg)
579 : : end if
580 : :
581 : 3 : write(mmom_unit,*) '#'
582 : 3 : write(mmom_unit,*) '# Local magnetic moments calculated and interpolated by ANADDB'
583 : 3 : write(mmom_unit,*) '#'
584 : :
585 : 3 : write(pfmt, '( "(es15.7, ", I2, "(es17.7))" )' ) ndim
586 [ + + ]: 57 : do imode= 1, 3*natom
587 : 54 : write(mmom_unit,*) ' '
588 : 54 : write(mmom_unit,'(a,i3)') '# Real part of local magnetic moments (at. units) induced by phonon mode:', imode
589 : 54 : write(msg,'(a,a,a)') ch10,&
590 : 108 : & ' # At hw m_{mat_1,1} m_{mat_1,2} ... m_{mat_2,1} m_{mat_2,2}'
591 : 54 : call wrtout(mmom_unit,msg,'COLL')
592 [ + + ]: 594 : do iw=1,nomega
593 : : write(msg,pfmt) &
594 [ + - + + ]: 2700 : & omega(iw), (real(modemm(i,imode,iw)),i=1,ndim)
595 : 594 : call wrtout(mmom_unit,msg,'COLL')
596 : : end do
597 : 54 : write(mmom_unit,*) ' '
598 : 54 : write(mmom_unit,'(a,i3)') '# Imaginary part of local magnetic moments (at. units) induced by phonon mode:', imode
599 : 54 : write(msg,'(a,a,a)') ch10,&
600 : 108 : & ' # At hw m_{mat_1,1} m_{mat_1,2} ... m_{mat_2,1} m_{mat_2,2}'
601 : 54 : call wrtout(mmom_unit,msg,'COLL')
602 [ + + ]: 597 : do iw=1,nomega
603 : : write(msg,pfmt) &
604 [ + - + + ]: 2700 : & omega(iw), (aimag(modemm(i,imode,iw)),i=1,ndim)
605 : 594 : call wrtout(mmom_unit,msg,'COLL')
606 : : end do
607 : : end do
608 : :
609 : 3 : write(pfmt, '( "(es15.7, ", I2, "(es17.7))" )' ) ndim*3
610 [ + + ]: 3 : if (qeq0) then
611 : 2 : write(mmom_unit,*) ' '
612 : 2 : write(mmom_unit,*) '# Real part of the clamped-ion local magnetic moments induced by an electric field (at. units)'
613 : 2 : write(msg,'(a,a,a)') ch10,&
614 : 2 : & ' # At hw m_{mat_1,1}^{Ex} m_{mat_1,1}^{Ey}',&
615 : 4 : & ' ... m_{mat_1,2}^{Ex} ... m_{mat_2,1}^{Ex} ...'
616 : 2 : call wrtout(mmom_unit,msg,'COLL')
617 [ + + ]: 22 : do iw=1,nomega
618 : : write(msg,pfmt) &
619 [ + - + + : 340 : & omega(iw), ((real(ci_localpha(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
620 : 22 : call wrtout(mmom_unit,msg,'COLL')
621 : : end do
622 : :
623 : 2 : write(mmom_unit,*) ' '
624 : 2 : write(mmom_unit,*) '# Imaginary part of the clamped-ion local magnetic moments induced by an electric field (at. units)'
625 : 2 : write(msg,'(a,a,a)') ch10,&
626 : 2 : & ' # At hw m_{mat_1,1}^{Ex} m_{mat_1,1}^{Ey}',&
627 : 4 : & ' ... m_{mat_1,2}^{Ex} ... m_{mat_2,1}^{Ex} ...'
628 : 2 : call wrtout(mmom_unit,msg,'COLL')
629 [ + + ]: 22 : do iw=1,nomega
630 : : write(msg,pfmt) &
631 [ + - + + : 340 : & omega(iw), ((aimag(ci_localpha(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
632 : 22 : call wrtout(mmom_unit,msg,'COLL')
633 : : end do
634 : :
635 : 2 : write(mmom_unit,*) ' '
636 : 2 : write(mmom_unit,*) '# Real part of relaxed-ion local magnetic moments induced by an electric field (at. units)'
637 : 2 : write(msg,'(a,a,a)') ch10,&
638 : 2 : & ' # At hw m_{mat_1,1}^{Ex} m_{mat_1,1}^{Ey}',&
639 : 4 : & ' ... m_{mat_1,2}^{Ex} ... m_{mat_2,1}^{Ex} ...'
640 : 2 : call wrtout(mmom_unit,msg,'COLL')
641 [ + + ]: 22 : do iw=1,nomega
642 : : write(msg,pfmt) &
643 [ + - + + : 340 : & omega(iw), ((real(ci_localpha(i,j,iw)+lm_localpha(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
644 : 22 : call wrtout(mmom_unit,msg,'COLL')
645 : : end do
646 : :
647 : 2 : write(mmom_unit,*) ' '
648 : 2 : write(mmom_unit,*) '# Imaginary part of relaxed-ion local magnetic moments induced by an electric field (at. units)'
649 : 2 : write(msg,'(a,a,a)') ch10,&
650 : 2 : & ' # At hw m_{mat_1,1}^{Ex} m_{mat_1,1}^{Ey}',&
651 : 4 : & ' ... m_{mat_1,2}^{Ex} ... m_{mat_2,1}^{Ex} ...'
652 : 2 : call wrtout(mmom_unit,msg,'COLL')
653 [ + + ]: 22 : do iw=1,nomega
654 : : write(msg,pfmt) &
655 [ + - + + : 340 : & omega(iw), ((aimag(ci_localpha(i,j,iw)+lm_localpha(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
656 : 22 : call wrtout(mmom_unit,msg,'COLL')
657 : : end do
658 : : end if
659 : :
660 : 3 : write(pfmt, '( "(es15.7, ", I2, "(es17.7))" )' ) ndim*3
661 : 3 : write(mmom_unit,*) ' '
662 : 3 : write(mmom_unit,*) '# Real part of the clamped-ion local magnetic moments induced by a Zeeman field (at. units)'
663 : 3 : write(msg,'(a,a,a)') ch10,&
664 : 3 : & ' # At hw m_{mat_1,1}^{Bx} m_{mat_1,1}^{By}',&
665 : 6 : & ' ... m_{mat_1,2}^{Bx} ... m_{mat_2,1}^{Bx} ...'
666 : 3 : call wrtout(mmom_unit,msg,'COLL')
667 [ + + ]: 33 : do iw=1,nomega
668 : : write(msg,pfmt) &
669 [ + - + + : 510 : & omega(iw), ((real(ci_locchi(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
670 : 33 : call wrtout(mmom_unit,msg,'COLL')
671 : : end do
672 : :
673 : 3 : write(mmom_unit,*) ' '
674 : 3 : write(mmom_unit,*) '# Imaginary part of the clamped-ion local magnetic moments induced by a Zeeman field (at. units)'
675 : 3 : write(msg,'(a,a,a)') ch10,&
676 : 3 : & ' # At hw m_{mat_1,1}^{Bx} m_{mat_1,1}^{By}',&
677 : 6 : & ' ... m_{mat_1,2}^{Bx} ... m_{mat_2,1}^{Bx} ...'
678 : 3 : call wrtout(mmom_unit,msg,'COLL')
679 [ + + ]: 33 : do iw=1,nomega
680 : : write(msg,pfmt) &
681 [ + - + + : 510 : & omega(iw), ((aimag(ci_locchi(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
682 : 33 : call wrtout(mmom_unit,msg,'COLL')
683 : : end do
684 : :
685 : 3 : write(mmom_unit,*) ' '
686 : 3 : write(mmom_unit,*) '# Real part of relaxed-ion local magnetic moments induced by a Zeeman field (at. units)'
687 : 3 : write(msg,'(a,a,a)') ch10,&
688 : 3 : & ' # At hw m_{mat_1,1}^{Bx} m_{mat_1,1}^{By}',&
689 : 6 : & ' ... m_{mat_1,2}^{Bx} ... m_{mat_2,1}^{Bx} ...'
690 : 3 : call wrtout(mmom_unit,msg,'COLL')
691 [ + + ]: 33 : do iw=1,nomega
692 : : write(msg,pfmt) &
693 [ + - + + : 510 : & omega(iw), ((real(ci_locchi(i,j,iw)+lm_locchi(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
694 : 33 : call wrtout(mmom_unit,msg,'COLL')
695 : : end do
696 : :
697 : 3 : write(mmom_unit,*) ' '
698 : 3 : write(mmom_unit,*) '# Imaginary part of relaxed-ion local magnetic moments induced by a Zeeman field (at. units)'
699 : 3 : write(msg,'(a,a,a)') ch10,&
700 : 3 : & ' # At hw m_{mat_1,1}^{Bx} m_{mat_1,1}^{By}',&
701 : 6 : & ' ... m_{mat_1,2}^{Bx} ... m_{mat_2,1}^{Bx} ...'
702 : 3 : call wrtout(mmom_unit,msg,'COLL')
703 [ + + ]: 33 : do iw=1,nomega
704 : : write(msg,pfmt) &
705 [ + - + + : 510 : & omega(iw), ((aimag(ci_locchi(i,j,iw)+lm_locchi(i,j,iw)),j=1,3),i=1,ndim)
+ - + + ]
706 : 33 : call wrtout(mmom_unit,msg,'COLL')
707 : : end do
708 : :
709 : 3 : close(mmom_unit)
710 : :
711 : : !Dielectric susceptibility
712 [ + + ]: 3 : if (qeq0) then
713 [ + - ]: 2 : diel_filename=trim(outfilename_radix)//"_DIELTENS"
714 [ - + ]: 2 : if (open_file(diel_filename, msg, newunit=diel_unit) /= 0) then
715 : 0 : ABI_ERROR(msg)
716 : : end if
717 : :
718 [ - + ]: 2 : if (mpopt==1) then
719 [ # # ]: 0 : call wrtout([ab_out,std_out,diel_unit], '# Fixed-spin dielectric tensor calculated and interpolated by ANADDB')
720 [ + - ]: 2 : else if (mpopt==2) then
721 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], '# Relaxed-spin dielectric tensor calculated and interpolated by ANADDB')
722 : : else
723 : 0 : write(msg,'(a)') 'ddb_omega_interpol: variable mpopt just can be 1 or 2'
724 : 0 : ABI_ERROR(msg)
725 : : end if
726 : :
727 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], ' ')
728 : 2 : write(pfmt, '( "(es15.7, ", I4, "(es15.7))" )' ) 9
729 : :
730 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], '# Real part of clamped-ion dielectric tensor')
731 : 2 : write(msg,'(a,a)') ch10,&
732 : 4 : & ' # At hw eps_11 eps_12 ... eps_21 eps_22 ...'
733 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit],msg,'COLL')
734 [ + + ]: 8 : do iw=1,min(nomega,3)
735 : : write(msg,pfmt) &
736 [ + - + + : 78 : & omega(iw), ((real(ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
737 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
738 : : end do
739 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
740 [ + + ]: 22 : do iw=1,nomega
741 : : write(msg,pfmt) &
742 [ + - + + : 260 : & omega(iw), ((real(ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
743 : 22 : call wrtout(diel_unit,msg,'COLL')
744 : : end do
745 : :
746 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], ' ')
747 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], '# Imaginary part of clamped-ion dielectric tensor')
748 : 2 : write(msg,'(a,a)') ch10,&
749 : 4 : & ' # At hw eps_11 eps_12 ... eps_21 eps_22 ...'
750 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit],msg,'COLL')
751 [ + + ]: 8 : do iw=1,min(nomega,3)
752 : : write(msg,pfmt) &
753 [ + - + + : 78 : & omega(iw), ((aimag(ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
754 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
755 : : end do
756 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
757 [ + + ]: 22 : do iw=1,nomega
758 : : write(msg,pfmt) &
759 [ + - + + : 260 : & omega(iw), ((aimag(ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
760 : 22 : call wrtout(diel_unit,msg,'COLL')
761 : : end do
762 : :
763 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], ' ')
764 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], '# Real part of relaxed-ion dielectric tensor ')
765 : 2 : write(msg,'(a,a)') ch10,&
766 : 4 : & ' # At hw eps_11 eps_12 ... eps_21 eps_22 ...'
767 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit],msg,'COLL')
768 [ + + ]: 8 : do iw=1,min(nomega,3)
769 : : write(msg,pfmt) &
770 [ + - + + : 78 : & omega(iw), ((real(lm_epsilon(i,j,iw)+ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
771 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
772 : : end do
773 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
774 [ + + ]: 22 : do iw=1,nomega
775 : : write(msg,pfmt) &
776 [ + - + + : 260 : & omega(iw), ((real(lm_epsilon(i,j,iw)+ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
777 : 22 : call wrtout(diel_unit,msg,'COLL')
778 : : end do
779 : :
780 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], ' ')
781 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit], '# Imaginary part of relaxed-ion dielectric tensor')
782 : 2 : write(msg,'(a,a)') ch10,&
783 : 4 : & ' # At hw eps_11 eps_12 ... eps_21 eps_22 ...'
784 [ + + ]: 8 : call wrtout([ab_out,std_out,diel_unit],msg,'COLL')
785 [ + + ]: 8 : do iw=1,min(nomega,3)
786 : : write(msg,pfmt) &
787 [ + - + + : 78 : & omega(iw), ((aimag(lm_epsilon(i,j,iw)+ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
788 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
789 : : end do
790 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
791 [ + + ]: 22 : do iw=1,nomega
792 : : write(msg,pfmt) &
793 [ + - + + : 260 : & omega(iw), ((aimag(lm_epsilon(i,j,iw)+ci_epsilon(i,j,iw)),j=1,3),i=1,3)
+ - + + ]
794 : 22 : call wrtout(diel_unit,msg,'COLL')
795 : : end do
796 : 2 : write(diel_unit,*) ' '
797 : :
798 : 2 : write(diel_unit,*) '#'
799 : 2 : write(diel_unit,*) '# Phonon modes contribution to dielectric tensor'
800 : 2 : write(diel_unit,*) '#'
801 : :
802 [ + + ]: 38 : do imode= 1, 3*natom
803 : 36 : write(diel_unit,*) ' '
804 : 36 : write(diel_unit,*) '# Real part of dielectric tensor due to phonon mode:', imode
805 : 36 : write(msg,'(a,a)') ch10,&
806 : 72 : & ' # At hw eps_11 eps_12 ... eps_21 eps_22 ...'
807 : 36 : call wrtout(diel_unit,msg,'COLL')
808 [ + + ]: 396 : do iw=1,nomega
809 : : write(msg,pfmt) &
810 [ + - + + : 4680 : & omega(iw), ((real(lm_epsilon_nm(i,j,imode,iw)),j=1,3),i=1,3)
+ - + + ]
811 : 396 : call wrtout(diel_unit,msg,'COLL')
812 : : end do
813 : 36 : write(diel_unit,*) ' '
814 : 36 : write(diel_unit,*) '# Imaginary part of dielectric tensor due to phonon mode:', imode
815 : 36 : write(msg,'(a,a)') ch10,&
816 : 72 : & ' # At hw eps_11 eps_12 ... eps_21 eps_22 ...'
817 : 36 : call wrtout(diel_unit,msg,'COLL')
818 [ + + ]: 398 : do iw=1,nomega
819 : : write(msg,pfmt) &
820 [ + - + + : 4680 : & omega(iw), ((aimag(lm_epsilon_nm(i,j,imode,iw)),j=1,3),i=1,3)
+ - + + ]
821 : 396 : call wrtout(diel_unit,msg,'COLL')
822 : : end do
823 : : end do
824 : :
825 : 2 : close(diel_unit)
826 : : end if
827 : :
828 : : !Magnetoelectric tensor mu_0 dM/dE
829 [ + + ]: 3 : if (qeq0) then
830 [ + - ]: 2 : alpha_filename=trim(outfilename_radix)//"_MAGNETOELTENS_dMdE"
831 [ - + ]: 2 : if (open_file(alpha_filename, msg, newunit=alpha_unit) /= 0) then
832 : 0 : ABI_ERROR(msg)
833 : : end if
834 : :
835 [ + + ]: 6 : call wrtout([ab_out,std_out], ' ')
836 [ - + ]: 2 : if (mpopt==1) then
837 [ # # ]: 0 : call wrtout([ab_out,std_out,alpha_unit], '# Fixed-spin magnetoelectric tensor calculated and interpolated by ANADDB')
838 [ + - ]: 2 : else if (mpopt==2) then
839 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],'# Relaxed-spin magnetoelectric tensor calculated and interpolated by ANADDB')
840 : : else
841 : 0 : write(msg,'(a)') 'ddb_omega_interpol: variable mpopt just can be 1 or 2'
842 : 0 : ABI_ERROR(msg)
843 : : end if
844 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],'# (Magnetization induced by electric field: mu_0 dM_a/dE_b)')
845 : :
846 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],' ')
847 : 2 : write(pfmt, '( "(es15.7, ", I4, "(es15.7))" )' ) 9
848 : :
849 : 2 : convfac= mu_0_SI * e_Cb * BField_Tesla / Bohr_meter**2 * 1.d12
850 : :
851 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],'# Real part of clamped-ion magnetoelectric tensor (ps/m)')
852 : 2 : write(msg,'(a,a)') ch10,&
853 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
854 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],msg,'COLL')
855 [ + + ]: 8 : do iw=1,min(nomega,3)
856 : : write(msg,pfmt) &
857 [ + - + + : 78 : & omega(iw), ((real(ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
858 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
859 : : end do
860 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
861 [ + + ]: 22 : do iw=1,nomega
862 : : write(msg,pfmt) &
863 [ + - + + : 260 : & omega(iw), ((real(ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
864 : 22 : call wrtout(alpha_unit,msg,'COLL')
865 : : end do
866 : :
867 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit], ' ')
868 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit], '# Imaginary part of clamped-ion magnetoelectric tensor (ps/m)')
869 : 2 : write(msg,'(a,a)') ch10,&
870 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
871 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],msg,'COLL')
872 [ + + ]: 8 : do iw=1,min(nomega,3)
873 : : write(msg,pfmt) &
874 [ + - + + : 78 : & omega(iw), ((aimag(ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
875 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
876 : : end do
877 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
878 [ + + ]: 22 : do iw=1,nomega
879 : : write(msg,pfmt) &
880 [ + - + + : 260 : & omega(iw), ((aimag(ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
881 : 22 : call wrtout(alpha_unit,msg,'COLL')
882 : : end do
883 : :
884 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],' ')
885 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],'# Real part of relaxed-ion magnetoelectric tensor (ps/m)')
886 : 2 : write(msg,'(a,a)') ch10,&
887 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
888 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],msg,'COLL')
889 [ + + ]: 8 : do iw=1,min(nomega,3)
890 : : write(msg,pfmt) &
891 [ + - + + : 78 : & omega(iw), ((real(lm_alpha(i,j,iw)+ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
892 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
893 : : end do
894 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
895 [ + + ]: 22 : do iw=1,nomega
896 : : write(msg,pfmt) &
897 [ + - + + : 260 : & omega(iw), ((real(lm_alpha(i,j,iw)+ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
898 : 22 : call wrtout(alpha_unit,msg,'COLL')
899 : : end do
900 : :
901 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],' ')
902 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],'# Imaginary part of relaxed-ion magnetoelectric tensor (ps/m)')
903 : 2 : write(msg,'(a,a)') ch10,&
904 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
905 [ + + ]: 8 : call wrtout([ab_out,std_out,alpha_unit],msg,'COLL')
906 [ + + ]: 8 : do iw=1,min(nomega,3)
907 : : write(msg,pfmt) &
908 [ + - + + : 78 : & omega(iw), ((aimag(lm_alpha(i,j,iw)+ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
909 [ + + ]: 20 : call wrtout([ab_out,std_out],msg,'COLL')
910 : : end do
911 [ + + ]: 6 : call wrtout([ab_out,std_out],'...')
912 [ + + ]: 22 : do iw=1,nomega
913 : : write(msg,pfmt) &
914 [ + - + + : 260 : & omega(iw), ((aimag(lm_alpha(i,j,iw)+ci_alpha(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
915 : 22 : call wrtout(alpha_unit,msg,'COLL')
916 : : end do
917 : 2 : write(alpha_unit,*) ''
918 : :
919 : 2 : write(alpha_unit,*) '#'
920 : 2 : write(alpha_unit,*) '# Phonon modes contribution to magnetoelectric tensor'
921 : 2 : write(alpha_unit,*) '#'
922 : :
923 [ + + ]: 38 : do imode= 1, 3*natom
924 : 36 : write(alpha_unit,*) ' '
925 : 36 : write(alpha_unit,*) '# Real part of magnetoelectric tensor (ps/m) due to phonon mode:', imode
926 : 36 : write(msg,'(a,a)') ch10,&
927 : 72 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
928 : 36 : call wrtout(alpha_unit,msg,'COLL')
929 [ + + ]: 396 : do iw=1,nomega
930 : : write(msg,pfmt) &
931 [ + - + + : 4680 : & omega(iw), ((real(lm_alpha_nm(i,j,imode,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
932 : 396 : call wrtout(alpha_unit,msg,'COLL')
933 : : end do
934 : 36 : write(alpha_unit,*) ' '
935 : 36 : write(alpha_unit,*) '# Imaginary part of magnetoelectric tensor (ps/m) due to phonon mode:', imode
936 : 36 : write(msg,'(a,a)') ch10,&
937 : 72 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
938 : 36 : call wrtout(alpha_unit,msg,'COLL')
939 [ + + ]: 398 : do iw=1,nomega
940 : : write(msg,pfmt) &
941 [ + - + + : 4680 : & omega(iw), ((aimag(lm_alpha_nm(i,j,imode,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
942 : 396 : call wrtout(alpha_unit,msg,'COLL')
943 : : end do
944 : : end do
945 : :
946 : 2 : close(alpha_unit)
947 : : end if
948 : :
949 : : !Magnetoelectric tensor dP/dH
950 [ + + ]: 3 : if (qeq0) then
951 [ + - ]: 2 : alpha_hc_filename=trim(outfilename_radix)//"_MAGNETOELTENS_dPdH"
952 [ - + ]: 2 : if (open_file(alpha_hc_filename, msg, newunit=alpha_hc_unit) /= 0) then
953 : 0 : ABI_ERROR(msg)
954 : : end if
955 : :
956 [ - + ]: 2 : if (mpopt==1) then
957 : 0 : write(alpha_hc_unit,*) '# Fixed-spin magnetoelectric tensor calculated and interpolated by ANADDB'
958 [ + - ]: 2 : else if (mpopt==2) then
959 : 2 : write(alpha_hc_unit,*) '# Relaxed-spin magnetoelectric tensor calculated and interpolated by ANADDB'
960 : : else
961 : 0 : write(msg,'(a)') 'ddb_omega_interpol: variable mpopt just can be 1 or 2'
962 : 0 : ABI_ERROR(msg)
963 : : end if
964 : 2 : write(alpha_hc_unit,*) '# (Polarization induced by Zeeman field strength: dP_a/dH_b)'
965 : :
966 : 2 : write(alpha_hc_unit,*) '#'
967 : 2 : write(pfmt, '( "(es15.7, ", I4, "(es15.7))" )' ) 9
968 : :
969 : 2 : convfac= mu_0_SI * e_Cb * BField_Tesla / Bohr_meter**2 * 1.d12
970 : :
971 : 2 : write(alpha_hc_unit,*) '# Real part of clamped-ion magnetoelectric tensor (ps/m)'
972 : 2 : write(msg,'(a,a)') ch10,&
973 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
974 : 2 : call wrtout(alpha_hc_unit,msg,'COLL')
975 [ + + ]: 22 : do iw=1,nomega
976 : : write(msg,pfmt) &
977 [ + - + + : 260 : & omega(iw), ((real(ci_alpha_hc(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
978 : 22 : call wrtout(diel_unit,msg,'COLL')
979 : : end do
980 : :
981 : 2 : write(alpha_hc_unit,*) ' '
982 : 2 : write(alpha_hc_unit,*) '# Imaginary part of clamped-ion magnetoelectric tensor (ps/m)'
983 : 2 : write(msg,'(a,a)') ch10,&
984 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
985 : 2 : call wrtout(alpha_hc_unit,msg,'COLL')
986 [ + + ]: 22 : do iw=1,nomega
987 : : write(msg,pfmt) &
988 [ + - + + : 260 : & omega(iw), ((aimag(ci_alpha_hc(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
989 : 22 : call wrtout(alpha_hc_unit,msg,'COLL')
990 : : end do
991 : :
992 : 2 : write(alpha_hc_unit,*) ' '
993 : 2 : write(alpha_hc_unit,*) '# Real part of relaxed-ion magnetoelectric tensor (ps/m)'
994 : 2 : write(msg,'(a,a)') ch10,&
995 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
996 : 2 : call wrtout(alpha_hc_unit,msg,'COLL')
997 [ + + ]: 22 : do iw=1,nomega
998 : : write(msg,pfmt) &
999 [ + - + + : 260 : & omega(iw), ((real(lm_alpha_hc(i,j,iw)+ci_alpha_hc(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1000 : 22 : call wrtout(alpha_hc_unit,msg,'COLL')
1001 : : end do
1002 : :
1003 : 2 : write(alpha_hc_unit,*) ' '
1004 : 2 : write(alpha_hc_unit,*) '# Imaginary part of relaxed-ion magnetoelectric tensor (ps/m)'
1005 : 2 : write(msg,'(a,a)') ch10,&
1006 : 4 : & ' # At hw alpha_11 alpha_12 ... alpha_21 alpha_22 ...'
1007 : 2 : call wrtout(alpha_hc_unit,msg,'COLL')
1008 [ + + ]: 22 : do iw=1,nomega
1009 : : write(msg,pfmt) &
1010 [ + - + + : 260 : & omega(iw), ((aimag(lm_alpha_hc(i,j,iw)+ci_alpha_hc(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1011 : 22 : call wrtout(alpha_hc_unit,msg,'COLL')
1012 : : end do
1013 : :
1014 : 2 : close(alpha_hc_unit)
1015 : : end if
1016 : :
1017 : : !Magnetic susceptibility
1018 [ + - ]: 3 : mmag_filename=trim(outfilename_radix)//"_MAGSUS"
1019 [ - + ]: 3 : if (open_file(mmag_filename, msg, newunit=mmag_unit) /= 0) then
1020 : 0 : ABI_ERROR(msg)
1021 : : end if
1022 : :
1023 [ + + ]: 9 : call wrtout([ab_out,std_out], ' ')
1024 [ - + ]: 3 : if (mpopt==1) then
1025 [ # # ]: 0 : call wrtout([ab_out,std_out,mmag_unit],'# Fixed-spin magnetic susceptibility calculated and interpolated by ANADDB')
1026 [ + - ]: 3 : else if (mpopt==2) then
1027 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],'# Relaxed-spin magnetic susceptibility calculated and interpolated by ANADDB')
1028 : : else
1029 : 0 : write(msg,'(a)') 'ddb_omega_interpol: variable mpopt just can be 1 or 2'
1030 : 0 : ABI_ERROR(msg)
1031 : : end if
1032 : :
1033 : 3 : write(pfmt, '( "(es15.7, ", I4, "(es15.7))" )' ) 9
1034 : :
1035 : 3 : convfac= two * mu_0_SI * mu_B_SI * BField_Tesla / Bohr_meter**3
1036 : :
1037 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],' ')
1038 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],'# Real part of clamped-ion magnetic susceptibility (SI adim)')
1039 : 3 : write(msg,'(a,a)') ch10,&
1040 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
1041 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],msg,'COLL')
1042 [ + + ]: 12 : do iw=1,min(nomega,3)
1043 : : write(msg,pfmt) &
1044 [ + - + + : 117 : & omega(iw), ((real(ci_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1045 [ + + ]: 30 : call wrtout([ab_out,std_out],msg,'COLL')
1046 : : end do
1047 [ + + ]: 9 : call wrtout([ab_out,std_out],'...')
1048 [ + + ]: 33 : do iw=1,nomega
1049 : : write(msg,pfmt) &
1050 [ + - + + : 390 : & omega(iw), ((real(ci_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1051 : 33 : call wrtout(mmag_unit,msg,'COLL')
1052 : : end do
1053 : :
1054 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],' ')
1055 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],'# Imaginary part of clamped-ion magnetic susceptibility (SI adim)')
1056 : 3 : write(msg,'(a,a)') ch10,&
1057 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
1058 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],msg,'COLL')
1059 [ + + ]: 12 : do iw=1,min(nomega,3)
1060 : : write(msg,pfmt) &
1061 [ + - + + : 117 : & omega(iw), ((aimag(ci_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1062 [ + + ]: 30 : call wrtout([ab_out,std_out],msg,'COLL')
1063 : : end do
1064 [ + + ]: 9 : call wrtout([ab_out,std_out],'...')
1065 [ + + ]: 33 : do iw=1,nomega
1066 : : write(msg,pfmt) &
1067 [ + - + + : 390 : & omega(iw), ((aimag(ci_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1068 : 33 : call wrtout(mmag_unit,msg,'COLL')
1069 : : end do
1070 : :
1071 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],' ')
1072 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],'# Real part of relaxed-ion magnetic susceptibility (SI adim)')
1073 : 3 : write(msg,'(a,a)') ch10,&
1074 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
1075 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],msg,'COLL')
1076 [ + + ]: 12 : do iw=1,min(nomega,3)
1077 : : write(msg,pfmt) &
1078 [ + - + + : 117 : & omega(iw), ((real(ci_mchi(i,j,iw)+lm_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1079 [ + + ]: 30 : call wrtout([ab_out,std_out],msg,'COLL')
1080 : : end do
1081 [ + + ]: 9 : call wrtout([ab_out,std_out],'...')
1082 [ + + ]: 33 : do iw=1,nomega
1083 : : write(msg,pfmt) &
1084 [ + - + + : 390 : & omega(iw), ((real(ci_mchi(i,j,iw)+lm_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1085 : 33 : call wrtout(mmag_unit,msg,'COLL')
1086 : : end do
1087 : :
1088 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],' ')
1089 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],'# Imaginary part of relaxed-ion magnetic susceptibility (SI adim)')
1090 : 3 : write(msg,'(a,a)') ch10,&
1091 : 6 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
1092 [ + + ]: 12 : call wrtout([ab_out,std_out,mmag_unit],msg,'COLL')
1093 [ + + ]: 12 : do iw=1,min(nomega,3)
1094 : : write(msg,pfmt) &
1095 [ + - + + : 117 : & omega(iw), ((aimag(ci_mchi(i,j,iw)+lm_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1096 [ + + ]: 30 : call wrtout([ab_out,std_out],msg,'COLL')
1097 : : end do
1098 [ + + ]: 9 : call wrtout([ab_out,std_out],'...')
1099 [ + + ]: 33 : do iw=1,nomega
1100 : : write(msg,pfmt) &
1101 [ + - + + : 390 : & omega(iw), ((aimag(ci_mchi(i,j,iw)+lm_mchi(i,j,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1102 : 33 : call wrtout(mmag_unit,msg,'COLL')
1103 : : end do
1104 : :
1105 : 3 : write(mmag_unit,*) ''
1106 : :
1107 : 3 : write(mmag_unit,*) '#'
1108 : 3 : write(mmag_unit,*) '# Phonon modes contribution to magnetic susceptibility tensor'
1109 : 3 : write(mmag_unit,*) '#'
1110 : :
1111 [ + + ]: 57 : do imode= 1, 3*natom
1112 : 54 : write(mmag_unit,*) ' '
1113 : 54 : write(mmag_unit,*) '# Real part of magnetic susceptibility (SI adim) due to phonon mode:', imode
1114 : 54 : write(msg,'(a,a)') ch10,&
1115 : 108 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
1116 : 54 : call wrtout(mmag_unit,msg,'COLL')
1117 [ + + ]: 594 : do iw=1,nomega
1118 : : write(msg,pfmt) &
1119 [ + - + + : 7020 : & omega(iw), ((real(lm_mchi_nm(i,j,imode,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1120 : 594 : call wrtout(mmag_unit,msg,'COLL')
1121 : : end do
1122 : 54 : write(mmag_unit,*) ' '
1123 : 54 : write(mmag_unit,*) '# Imaginary part of magnetic susceptibility (SI adim) due to phonon mode:', imode
1124 : 54 : write(msg,'(a,a)') ch10,&
1125 : 108 : & ' # At hw X_11 X_12 ... X_21 X_22 ...'
1126 : 54 : call wrtout(mmag_unit,msg,'COLL')
1127 [ + + ]: 597 : do iw=1,nomega
1128 : : write(msg,pfmt) &
1129 [ + - + + : 7020 : & omega(iw), ((aimag(lm_mchi_nm(i,j,imode,iw))*convfac,j=1,3),i=1,3)
+ - + + ]
1130 : 594 : call wrtout(mmag_unit,msg,'COLL')
1131 : : end do
1132 : : end do
1133 : :
1134 : 3 : close(mmag_unit)
1135 : :
1136 : : !Phonon spectral function
1137 [ + - ]: 3 : phon_filename=trim(outfilename_radix)//"_SPECTRAL_PHONON"
1138 [ - + ]: 3 : if (open_file(phon_filename, msg, newunit=phon_unit) /= 0) then
1139 : 0 : ABI_ERROR(msg)
1140 : : end if
1141 : :
1142 : 3 : write(phon_unit,*) '#'
1143 [ - + ]: 3 : if (mpopt==1) then
1144 : 0 : write(phon_unit,*) '# Fixed-spin phonon spectral function calculated and interpolated by ANADDB'
1145 [ + - ]: 3 : else if (mpopt==2) then
1146 : 3 : write(phon_unit,*) '# Relaxed-spin phonon spectral function calculated and interpolated by ANADDB'
1147 : : else
1148 : 0 : write(msg,'(a)') 'ddb_omega_interpol: variable mpopt just can be 1 or 2'
1149 : 0 : ABI_ERROR(msg)
1150 : : end if
1151 : :
1152 : 3 : write(msg,'(a,a)') ch10, ' # At hw Phonon SF '
1153 : 3 : call wrtout(phon_unit,msg,'COLL')
1154 [ + + ]: 33 : do iw=1,nomega
1155 : 30 : write(msg,*) omega(iw), phonspec(iw)
1156 : 33 : call wrtout(phon_unit,msg,'COLL')
1157 : : end do
1158 : :
1159 : 3 : close(phon_unit)
1160 : :
1161 : : !Phonon frequencies
1162 [ + - ]: 3 : phon_filename=trim(outfilename_radix)//"_PHFRW"
1163 [ - + ]: 3 : if (open_file(phon_filename, msg, newunit=phon_unit) /= 0) then
1164 : 0 : ABI_ERROR(msg)
1165 : : end if
1166 : :
1167 : 3 : write(phon_unit,*) '#'
1168 [ - + ]: 3 : if (mpopt==1) then
1169 : 0 : write(phon_unit,*) '# Frozen-spin phonon frequencies calculated and interpolated by ANADDB'
1170 [ + - ]: 3 : else if (mpopt==2) then
1171 : 3 : write(phon_unit,*) '# Relaxed-Spin phonon frequencies calculated and interpolated by ANADDB'
1172 : : else
1173 : 0 : write(msg,'(a)') 'ddb_omega_interpol: variable mpopt just can be 1 or 2'
1174 : 0 : ABI_ERROR(msg)
1175 : : end if
1176 : :
1177 : 3 : write(pfmt, '( "(es15.7, ", I4, "(es15.7))" )' ) natom*3
1178 : 3 : write(msg,'(a,a)') ch10,&
1179 : 6 : & ' # At hw eval(1) eval(2) ...'
1180 : 3 : call wrtout(phon_unit,msg,'COLL')
1181 [ + + ]: 33 : do iw=1,nomega
1182 : 30 : write(msg,pfmt) omega(iw), phfrq(:,iw)
1183 : 33 : call wrtout(phon_unit,msg,'COLL')
1184 : : end do
1185 : :
1186 : 3 : close(phon_unit)
1187 : :
1188 : : !Phonon eigendisplacements
1189 [ + - ]: 3 : phon_filename=trim(outfilename_radix)//"_PHDISP"
1190 [ - + ]: 3 : if (open_file(phon_filename, msg, newunit=phon_unit) /= 0) then
1191 : 0 : ABI_ERROR(msg)
1192 : : end if
1193 : :
1194 : 3 : write(phon_unit,*) '#'
1195 [ - + ]: 3 : if (mpopt==1) then
1196 : 0 : write(phon_unit,*) '# Frozen-spin phonon eigendisplacements calculated and interpolated by ANADDB'
1197 [ + - ]: 3 : else if (mpopt==2) then
1198 : 3 : write(phon_unit,*) '# Relaxed-spin phonon eigendisplacements calculated and interpolated by ANADDB'
1199 : : else
1200 : 0 : write(msg,'(a)') 'ddb_omega_interpol: variable mpopt just can be 1 or 2'
1201 : 0 : ABI_ERROR(msg)
1202 : : end if
1203 : :
1204 : 3 : write(pfmt, '( "(es15.7, ", I2, "(es17.7))" )' ) 3*natom
1205 [ + + ]: 57 : do imode= 1, 3*natom
1206 : 54 : write(phon_unit,*) ' '
1207 : 54 : write(phon_unit,'(a,i3)') '# Real part of phonon eigendisplacement:', imode
1208 [ + + ]: 594 : do iw=1,nomega
1209 : : write(msg,pfmt) &
1210 [ + - + + ]: 10260 : & omega(iw), (real(modedisp(i,imode,iw)),i=1,3*natom)
1211 : 594 : call wrtout(phon_unit,msg,'COLL')
1212 : : end do
1213 : :
1214 : 54 : write(phon_unit,*) ' '
1215 : 54 : write(phon_unit,'(a,i3)') '# Imaginary part of phonon eigendisplacement:', imode
1216 [ + + ]: 597 : do iw=1,nomega
1217 : : write(msg,pfmt) &
1218 [ + - + + ]: 10260 : & omega(iw), (aimag(modedisp(i,imode,iw)),i=1,3*natom)
1219 : 594 : call wrtout(phon_unit,msg,'COLL')
1220 : : end do
1221 : : end do
1222 : :
1223 : 3 : close(phon_unit)
1224 : : !Born effective charges
1225 [ + - ]: 3 : zeff_filename=trim(outfilename_radix)//"_ZEFF"
1226 [ - + ]: 3 : if (open_file(zeff_filename, msg, newunit=zeff_unit) /= 0) then
1227 : 0 : ABI_ERROR(msg)
1228 : : end if
1229 : :
1230 [ + + ]: 3 : if (qeq0) then
1231 : 2 : write(zeff_unit,*) '#'
1232 : 2 : write(zeff_unit,*) '# Electric and magnetic Born effective charges calculated and interpolated by ANADDB'
1233 : 2 : write(zeff_unit,*) '#'
1234 : :
1235 : 2 : write(pfmt, '( "(es15.7, ", I2, "(es17.7))" )' ) 3
1236 [ + + ]: 38 : do imode= 1, 3*natom
1237 : 36 : write(zeff_unit,*) ' '
1238 : 36 : write(zeff_unit,'(a,i3)') '# Real part of electric Born charge (at. units) induced by phonon mode:', imode
1239 : 36 : write(msg,'(a,a)') ch10,&
1240 : 72 : & ' # At hw Z^x_{n} Z^y_{n} Z^z_{n}'
1241 : 36 : call wrtout(zeff_unit,msg,'COLL')
1242 [ + + ]: 396 : do iw=1,nomega
1243 : : write(msg,pfmt) &
1244 [ + - + + ]: 1440 : & omega(iw), (real(modezeff(i,imode,iw)),i=1,3)
1245 : 396 : call wrtout(zeff_unit,msg,'COLL')
1246 : : end do
1247 : 36 : write(zeff_unit,*) ' '
1248 : 36 : write(zeff_unit,'(a,i3)') '# Imaginary part of electric Born charge (at. units) induced by phonon mode:', imode
1249 : 36 : write(msg,'(a,a)') ch10,&
1250 : 72 : & ' # At hw Z^x_{n} Z^y_{n} Z^z_{n}'
1251 : 36 : call wrtout(zeff_unit,msg,'COLL')
1252 [ + + ]: 398 : do iw=1,nomega
1253 : : write(msg,pfmt) &
1254 [ + - + + ]: 1440 : & omega(iw), (aimag(modezeff(i,imode,iw)),i=1,3)
1255 : 396 : call wrtout(zeff_unit,msg,'COLL')
1256 : : end do
1257 : : end do
1258 : : end if
1259 : :
1260 [ + + ]: 57 : do imode= 1, 3*natom
1261 : 54 : write(zeff_unit,*) ' '
1262 : 54 : write(zeff_unit,'(a,i3)') '# Real part of magnetic Born charge (at. units) induced by phonon mode:', imode
1263 : 54 : write(msg,'(a,a)') ch10,&
1264 : 108 : & ' # At hw M^x_{n} M^y_{n} M^z_{n}'
1265 : 54 : call wrtout(zeff_unit,msg,'COLL')
1266 [ + + ]: 594 : do iw=1,nomega
1267 : : write(msg,pfmt) &
1268 [ + - + + ]: 2160 : & omega(iw), (real(modemeff(i,imode,iw)),i=1,3)
1269 : 594 : call wrtout(zeff_unit,msg,'COLL')
1270 : : end do
1271 : 54 : write(zeff_unit,*) ' '
1272 : 54 : write(zeff_unit,'(a,i3)') '# Imaginary part of magnetic Born charge (at. units) induced by phonon mode:', imode
1273 : 54 : write(msg,'(a,a)') ch10,&
1274 : 108 : & ' # At hw M^x_{n} M^y_{n} M^z_{n}'
1275 : 54 : call wrtout(zeff_unit,msg,'COLL')
1276 [ + + ]: 597 : do iw=1,nomega
1277 : : write(msg,pfmt) &
1278 [ + - + + ]: 2160 : & omega(iw), (aimag(modemeff(i,imode,iw)),i=1,3)
1279 : 594 : call wrtout(zeff_unit,msg,'COLL')
1280 : : end do
1281 : : end do
1282 : :
1283 : 3 : close(zeff_unit)
1284 : :
1285 [ - + ]: 3 : ABI_FREE(dint_fsddb)
1286 [ - + ]: 3 : ABI_FREE(int_fsddb)
1287 [ - + ]: 3 : ABI_FREE(int_rsddb)
1288 [ - + ]: 3 : ABI_FREE(dummysus)
1289 [ - + ]: 3 : ABI_FREE(dummysus1)
1290 [ - + ]: 3 : ABI_FREE(dummysus2)
1291 [ - + ]: 3 : ABI_FREE(dummymom_tr)
1292 [ - + ]: 3 : ABI_FREE(dummymom)
1293 [ - + ]: 3 : ABI_FREE(omegacalc)
1294 [ - + ]: 3 : ABI_FREE(magsus)
1295 [ - + ]: 3 : ABI_FREE(lm_magsus)
1296 [ - + ]: 3 : ABI_FREE(invmagsus)
1297 [ - + ]: 3 : ABI_FREE(mmom)
1298 [ - + ]: 3 : ABI_FREE(mmom_tr)
1299 [ - + ]: 3 : ABI_FREE(zfield)
1300 [ - + ]: 3 : ABI_FREE(zfield_tr)
1301 [ - + ]: 3 : ABI_FREE(ci_epsilon)
1302 [ - + ]: 3 : ABI_FREE(lm_epsilon)
1303 [ - + ]: 3 : ABI_FREE(lm_epsilon_nm)
1304 [ - + ]: 3 : ABI_FREE(ci_mchi)
1305 [ - + ]: 3 : ABI_FREE(lm_mchi)
1306 [ - + ]: 3 : ABI_FREE(lm_mchi_nm)
1307 [ - + ]: 3 : ABI_FREE(ci_alpha)
1308 [ - + ]: 3 : ABI_FREE(ci_alpha_hc)
1309 [ - + ]: 3 : ABI_FREE(lm_alpha)
1310 [ - + ]: 3 : ABI_FREE(lm_alpha_hc)
1311 [ - + ]: 3 : ABI_FREE(lm_alpha_nm)
1312 [ - + ]: 3 : ABI_FREE(ci_localpha)
1313 [ - + ]: 3 : ABI_FREE(lm_localpha)
1314 [ - + ]: 3 : ABI_FREE(ci_locchi)
1315 [ - + ]: 3 : ABI_FREE(lm_locchi)
1316 [ - + ]: 3 : ABI_FREE(omega)
1317 [ - + ]: 3 : ABI_FREE(phfrq)
1318 [ - + ]: 3 : ABI_FREE(phongreen)
1319 [ - + ]: 3 : ABI_FREE(phonspec)
1320 [ - + ]: 3 : ABI_FREE(mode_phonspec)
1321 [ - + ]: 3 : ABI_FREE(displ)
1322 [ - + ]: 3 : ABI_FREE(eigvec)
1323 [ - + ]: 3 : ABI_FREE(modemm)
1324 [ - + ]: 3 : ABI_FREE(modedisp)
1325 [ - + ]: 3 : ABI_FREE(modezf)
1326 [ - + ]: 3 : ABI_FREE(modezeff)
1327 [ - + ]: 3 : ABI_FREE(modemeff)
1328 [ + + ]: 3 : ABI_SFREE(w0hessian)
1329 [ + + ]: 3 : ABI_SFREE(w0berry)
1330 [ + + ]: 3 : ABI_SFREE(coeffs)
1331 [ - + ]: 3 : ABI_SFREE(pc_in)
1332 : :
1333 : 3 : end subroutine ddb_omega_interpol
1334 : : !!***
1335 : :
1336 : : !!****f* ABINIT/phonon_green
1337 : : !! NAME
1338 : : !! phonon_green
1339 : : !!
1340 : : !! FUNCTION
1341 : : !! Computes the phonon Green's function and spectral function at
1342 : : !! a given value of frequency and imaginary damping
1343 : : !!
1344 : : !! COPYRIGHT
1345 : : !! Copyright (C) 2024 ABINIT group (FIXME: add author)
1346 : : !! This file is distributed under the terms of the
1347 : : !! GNU General Public License, see ~abinit/COPYING
1348 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1349 : : !!
1350 : : !! INPUTS
1351 : : !! amu(ntypat)= atomic masses
1352 : : !! ifc(3*natom,3*natom)= Interatomic-force constants calculated at a given omega
1353 : : !! natom= number of atoms in the cell
1354 : : !! ntypat= number of atom types in the cell
1355 : : !! omega= frequency at which the IFCs have been calculated
1356 : : !! typat(natom)= array with the type of atoms in the cell
1357 : : !!
1358 : : !! OUTPUT
1359 : : !! phonspec= phonon spectral function at the input omega
1360 : : !! phfrq= phonon frequencies calculated with the IFCs at the input omega
1361 : : !! phfrq= phonon eigenvectors calculated with the IFCs at the input omega
1362 : : !!
1363 : : !!
1364 : : !! SIDE EFFECTS
1365 : : !!
1366 : : !! NOTES
1367 : : !!
1368 : : !! PARENTS
1369 : : !!
1370 : : !! CHILDREN
1371 : : !!
1372 : : !! SOURCE
1373 : :
1374 : : #if defined HAVE_CONFIG_H
1375 : : #include "config.h"
1376 : : #endif
1377 : :
1378 : : #include "abi_common.h"
1379 : :
1380 : :
1381 : 30 : subroutine phonon_green(amu,displ,eigvec,eta,blkval,&
1382 : 30 : & mode_phonspec,mpert,natom,ntypat,omega,&
1383 : 30 : & phfrq,phongreen,phonspec,typat)
1384 : :
1385 : : use defs_basis
1386 : : use m_errors
1387 : : use m_profiling_abi
1388 : :
1389 : : implicit none
1390 : :
1391 : : !Arguments ------------------------------------
1392 : : !scalars
1393 : : integer, intent(in) :: mpert,natom,ntypat
1394 : : real(dp), intent(in) :: eta,omega
1395 : : real(dp), intent(out) :: phonspec
1396 : : !arrays
1397 : : integer, intent(in) :: typat(natom)
1398 : : real(dp), intent(in) :: amu(ntypat)
1399 : : real(dp), intent(in) :: blkval(2,3,mpert,3,mpert,1)
1400 : : real(dp), intent(out) :: displ(2*3*natom*3*natom)
1401 : : real(dp), intent(out) :: eigvec(2*3*natom*3*natom)
1402 : : real(dp), intent(out) :: phfrq(3*natom)
1403 : : real(dp), intent(out) :: mode_phonspec(3*natom)
1404 : : complex(dpc), intent(out) :: phongreen(3*natom,3*natom)
1405 : :
1406 : : !Local variables-------------------------------
1407 : : !scalars
1408 : : integer :: iat1,iat2,idir1,idir2,ipert1,ipert2
1409 : : integer :: icol,ier,imode,info,irow,lwork,pdim
1410 : : real(dp) :: mfac1, mfac2
1411 : : complex(dpc) :: cplx_eta
1412 : : !arrays
1413 : 30 : integer, allocatable :: ipiv(:)
1414 : : real(dp) :: dum(2,0)
1415 : 30 : real(dp), allocatable :: invmassfac(:,:)
1416 : 30 : real(dp), allocatable :: matrx(:,:),zhpev1(:,:),zhpev2(:)
1417 : 30 : real(dp), allocatable :: eigval(:)
1418 : : ! real(dp), allocatable, save :: delta_asrw0(:,:)
1419 : 30 : complex(dpc), allocatable :: ifc(:,:)
1420 : 30 : complex(dpc), allocatable :: dynmat(:,:),ifc_w2mass(:,:)
1421 : 30 : complex(dpc),allocatable :: work(:),work1(:,:)
1422 : 30 : complex(dpc),allocatable :: mass_phongreen(:,:)
1423 : : !character(len=500) :: msg
1424 : :
1425 : : ! *************************************************************************
1426 : :
1427 : : DBG_ENTER("COLL")
1428 : :
1429 : : !Extract the IFCs.
1430 : 30 : pdim=3*natom
1431 [ - + + - : 120 : ABI_MALLOC(ifc,(pdim,pdim))
- + + - +
- + - - +
- + ]
1432 [ + + ]: 210 : do ipert2= 1, natom
1433 [ + + ]: 750 : do idir2= 1, 3
1434 : 540 : icol= (ipert2-1)*3 + idir2
1435 [ + + ]: 3960 : do ipert1= 1, natom
1436 [ + + ]: 13500 : do idir1= 1, 3
1437 : 9720 : irow= (ipert1-1)*3 + idir1
1438 : : ifc(irow,icol)= &
1439 : : & cmplx(blkval(1,idir1,ipert1,idir2,ipert2,1), &
1440 : 12960 : & blkval(2,idir1,ipert1,idir2,ipert2,1),kind=dpc)
1441 : : end do
1442 : : end do
1443 : : end do
1444 : : end do
1445 : :
1446 : : !Apply ASR: it has weird consequences on the intensities of the spectral function
1447 : : !better not applied.
1448 : : ! ABI_MALLOC_IFNOT(delta_asrw0,(3*natom,3))
1449 : : ! if (omega < tol14) then
1450 : : ! call asrw0(delta_asrw0,ifc,natom,0)
1451 : : ! else
1452 : : ! call asrw0(delta_asrw0,ifc,natom,1)
1453 : : ! end if
1454 : :
1455 : : !Build an array with the inverse mass factors
1456 [ - + + - : 120 : ABI_MALLOC(invmassfac,(natom,natom))
- + + - +
- + - - +
- + ]
1457 [ + + ]: 210 : do iat2= 1, natom
1458 [ + + ]: 1290 : do iat1= 1, natom
1459 : 1260 : invmassfac(iat1,iat2)=one/sqrt(amu(typat(iat1))*amu(typat(iat2)))/amu_emass
1460 : : end do
1461 : : end do
1462 : :
1463 : : !Build the dynamical and (Phi-M(w+eta)**2) matrices
1464 [ - + + - : 90 : ABI_MALLOC(dynmat,(pdim,pdim))
- + + - -
+ - + ]
1465 [ - + + - : 90 : ABI_MALLOC(ifc_w2mass,(pdim,pdim))
- + + - -
+ - + ]
1466 : :
1467 : 30 : cplx_eta=cmplx(0.0_dp,eta)
1468 [ + + ]: 210 : do iat2= 1, natom
1469 [ + + ]: 750 : do idir2= 1, 3
1470 : 540 : icol= (iat2-1)*3 + idir2
1471 [ + + ]: 3960 : do iat1= 1, natom
1472 [ + + ]: 13500 : do idir1= 1, 3
1473 : 9720 : irow= (iat1-1)*3 + idir1
1474 : 9720 : dynmat(irow,icol)= invmassfac(iat1,iat2)*ifc(irow,icol)
1475 : 9720 : ifc_w2mass(irow,icol)= ifc(irow,icol)
1476 [ + + ]: 12960 : if (irow==icol) then
1477 : : ifc_w2mass(irow,icol)= ifc_w2mass(irow,icol) - &
1478 : 540 : & amu(typat(iat1))*amu_emass*(omega+cplx_eta)**2
1479 : : end if
1480 : : end do
1481 : : end do
1482 : : end do
1483 : : end do
1484 : :
1485 : : !Invert to obtain the phonon Green's function
1486 [ - + + - : 90 : ABI_MALLOC(work1,(pdim,pdim))
- + + - -
+ - + ]
1487 [ - + + + : 10320 : work1=ifc_w2mass
+ + ]
1488 : :
1489 [ - + + - : 90 : ABI_MALLOC(ipiv,(pdim))
+ - + - -
+ - + ]
1490 : 30 : call zgetrf( pdim, pdim, work1, pdim, ipiv, info )
1491 [ - + ]: 30 : ABI_CHECK(info == 0, sjoin('zgetrf returned:', itoa(info)))
1492 : :
1493 [ - + ]: 30 : ABI_MALLOC(work,(2))
1494 : 30 : call zgetri( pdim, work1, pdim, ipiv, work, -1, info )
1495 [ - + ]: 30 : ABI_CHECK(info == 0, sjoin('zgetri returned:', itoa(info)))
1496 : 30 : lwork=int(work(1))
1497 : :
1498 [ - + + - : 90 : ABI_REMALLOC(work,(lwork))
+ - + - -
+ - + -
+ ]
1499 : 30 : call zgetri( pdim, work1, pdim, ipiv, work, lwork, info )
1500 [ - + ]: 30 : ABI_CHECK(info == 0, sjoin('zgetri returned:', itoa(info)))
1501 [ - + ]: 30 : ABI_FREE(work)
1502 : :
1503 [ + + + + ]: 10290 : phongreen=-work1
1504 : :
1505 : : !Now apply the mass factors
1506 [ - + + - : 120 : ABI_MALLOC(mass_phongreen,(pdim,pdim))
- + + - +
- + - - +
- + ]
1507 [ + + ]: 570 : do icol= 1, pdim
1508 [ + + ]: 540 : iat2= ceiling(icol/three)
1509 : 540 : mfac2= sqrt(amu(typat(iat2))*amu_emass)
1510 [ + + ]: 10290 : do irow= 1, pdim
1511 [ + + ]: 9720 : iat1= ceiling(irow/three)
1512 : 9720 : mfac1= sqrt(amu(typat(iat1))*amu_emass)
1513 : 10260 : mass_phongreen(irow,icol)= -mfac1*work1(irow,icol)*mfac2
1514 : : end do
1515 : : end do
1516 : :
1517 : : !Finally extract the spectral function from the trace
1518 [ + + ]: 570 : do irow= 1, pdim
1519 : 570 : mode_phonspec(irow)= -one/pi * aimag(two*cmplx(omega,eta,kind=dpc)*mass_phongreen(irow,irow))
1520 : : end do
1521 [ + + ]: 570 : phonspec= sum(mode_phonspec(:))
1522 [ - + ]: 30 : ABI_FREE(mass_phongreen)
1523 : :
1524 : : !Diagonalize the Dynamical matrix
1525 [ + - - + : 90 : ABI_MALLOC(matrx,(2,(3*natom*(3*natom+1))/2))
+ - + - +
- - + -
+ ]
1526 [ - + + - : 90 : ABI_MALLOC(eigval,(pdim))
+ - + - -
+ - + ]
1527 [ + + ]: 570 : do icol= 1, pdim
1528 [ + + ]: 5700 : do irow= 1, icol
1529 : 5130 : matrx(1,irow + (icol-1)*icol/2)=real(dynmat(irow,icol))
1530 : 5670 : matrx(2,irow + (icol-1)*icol/2)=aimag(dynmat(irow,icol))
1531 : : end do
1532 : : end do
1533 : :
1534 [ + - - + : 90 : ABI_MALLOC(zhpev1,(2,2*3*natom-1))
+ - + - +
- - + -
+ ]
1535 [ + - - + : 90 : ABI_MALLOC(zhpev2,(3*3*natom-2))
+ - + - +
- - + -
+ ]
1536 : :
1537 : 30 : call ZHPEV ('V','U',3*natom,matrx,eigval,eigvec,3*natom,zhpev1,zhpev2,ier)
1538 [ - + ]: 30 : ABI_CHECK(ier == 0, sjoin('zhpev returned:', itoa(ier)))
1539 : :
1540 [ - + ]: 30 : ABI_FREE(matrx)
1541 [ - + ]: 30 : ABI_FREE(dynmat)
1542 [ - + ]: 30 : ABI_FREE(zhpev1)
1543 [ - + ]: 30 : ABI_FREE(zhpev2)
1544 : :
1545 : : !Get the phonon frequencies (negative by convention, if the eigenvalue of the dynamical matrix is negative)
1546 [ + + ]: 570 : do imode=1,3*natom
1547 [ + + ]: 570 : if(eigval(imode)>=1.0d-16)then
1548 : 520 : phfrq(imode)=sqrt(eigval(imode))
1549 [ - + ]: 20 : else if(eigval(imode)>=-1.0d-16)then
1550 : 0 : phfrq(imode)=zero
1551 : : else
1552 : 20 : phfrq(imode)=-sqrt(-eigval(imode))
1553 : : end if
1554 : : end do
1555 : :
1556 : : !Fix the phase of the eigenvectors
1557 : 30 : call fxphas_seq(eigvec,dum, 0, 0, 1, 3*natom*3*natom, 0, 3*natom, 3*natom, 0)
1558 : :
1559 : : !Normalise the eigenvectors
1560 : 30 : call pheigvec_normalize(natom, eigvec)
1561 : :
1562 : : ! Get the phonon displacements
1563 : 30 : call phdispl_from_eigvec(natom, ntypat, typat, amu, eigvec, displ)
1564 : :
1565 [ - + ]: 30 : ABI_FREE(ifc)
1566 [ - + ]: 30 : ABI_FREE(ifc_w2mass)
1567 [ - + ]: 30 : ABI_FREE(ipiv)
1568 [ - + ]: 30 : ABI_FREE(work1)
1569 [ - + ]: 30 : ABI_FREE(invmassfac)
1570 [ - + ]: 30 : ABI_FREE(eigval)
1571 : : ! ABI_FREE(delta_asrw0)
1572 : :
1573 : : DBG_EXIT("COLL")
1574 : :
1575 : 30 : end subroutine phonon_green
1576 : : !!***
1577 : :
1578 : : !!****f* ABINIT/lm_normal_modes
1579 : : !! NAME
1580 : : !! lm_normal_modes
1581 : : !!
1582 : : !! FUNCTION
1583 : : !! Calculates the lattice-mediated contributions of the different type of
1584 : : !! susceptibilities by projecting the calculation on the contributions of
1585 : : !! the phonon modes.
1586 : : !!
1587 : : !! COPYRIGHT
1588 : : !! Copyright (C) 2024 ABINIT group (FIXME: add author)
1589 : : !! This file is distributed under the terms of the
1590 : : !! GNU General Public License, see ~abinit/COPYING
1591 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1592 : : !!
1593 : : !! INPUTS
1594 : : !! amu(ntypat)= atomic masses
1595 : : !! blkval(2,3,mpert,3,mpert,1)= array with second-order derivatives
1596 : : !! displ(2,3,natom,3,natom)= phonon eigendisplacements
1597 : : !! mmom(ndim,(natom+2)*3)= first-order magnetic moments
1598 : : !! natom= number of atoms in the cell
1599 : : !! ndim= dimension of the penalized degrees of freedom
1600 : : !! ntypat= number of atom types in the cell
1601 : : !! typat(natom)= array with the type of atoms in the cell
1602 : : !!
1603 : : !! OUTPUT
1604 : : !! modemm(ndim,3*natom)= mode-resolved local magnetic moments
1605 : : !! modezf(ndim,3*natom)= mode-resolved local Zeeman fields
1606 : : !! modemeff(3,3*natom)= mode-resolved magnetic Born charges
1607 : : !! modezeff(3,3*natom)= mode-resolved electric Born charges
1608 : : !!
1609 : : !! SIDE EFFECTS
1610 : : !!
1611 : : !! NOTES
1612 : : !!
1613 : : !! PARENTS
1614 : : !!
1615 : : !! CHILDREN
1616 : : !!
1617 : : !! SOURCE
1618 : :
1619 : : #if defined HAVE_CONFIG_H
1620 : : #include "config.h"
1621 : : #endif
1622 : :
1623 : : #include "abi_common.h"
1624 : :
1625 : 30 : subroutine lm_normal_modes(amu,blkval,displ,eta,lm_alpha_nm,lm_epsilon_nm,lm_mchi_nm, &
1626 : 30 : & mcell,mmom,modemm,modedisp,modemeff,modezeff,modezf,mpert,natom,ndim,ntypat,omega,phfrq,typat,ucvol,zfield)
1627 : :
1628 : : use defs_basis
1629 : : use m_errors
1630 : : use m_profiling_abi
1631 : :
1632 : : implicit none
1633 : :
1634 : : !Arguments ------------------------------------
1635 : : !scalars
1636 : : integer, intent(in) :: mpert,natom,ndim,ntypat
1637 : : real(dp), intent(in) :: eta,mcell,omega,ucvol
1638 : : !arrays
1639 : : integer, intent(in) :: typat(natom)
1640 : : real(dp), intent(in) :: amu(ntypat)
1641 : : real(dp), intent(in) :: blkval(2,3,mpert,3,mpert,1)
1642 : : real(dp), intent(inout) :: displ(2*3*natom*3*natom)
1643 : : real(dp), intent(in) :: phfrq(3*natom)
1644 : : complex(dpc), intent(in) :: mmom(ndim,(natom+5)*3)
1645 : : complex(dpc), intent(in) :: zfield(ndim,(natom+5)*3)
1646 : : complex(dpc), intent(out) :: lm_alpha_nm(3,3,3*natom)
1647 : : complex(dpc), intent(out) :: lm_epsilon_nm(3,3,3*natom)
1648 : : complex(dpc), intent(out) :: lm_mchi_nm(3,3,3*natom)
1649 : : complex(dpc), intent(out) :: modemm(ndim,3*natom)
1650 : : complex(dpc), intent(out) :: modedisp(3*natom,3*natom)
1651 : : complex(dpc), intent(out) :: modemeff(3,3*natom)
1652 : : complex(dpc), intent(out) :: modezeff(3,3*natom)
1653 : : complex(dpc), intent(out) :: modezf(ndim,3*natom)
1654 : :
1655 : : !Local variables-------------------------------
1656 : : !scalars
1657 : : integer :: i1,iat1,iat2,idir1,idir2,im,imode,index,ipert,ipert1,irow,jpert
1658 : : integer :: pdim
1659 : : real(dp) :: fac
1660 : : complex(dpc) :: cplx_eta,cplx_w2
1661 : : !arrays
1662 : 30 : complex(dpc),allocatable :: c_blkval(:,:,:,:) !,norm(:)
1663 : 30 : complex(dpc),allocatable :: zeff(:,:), zeff_tr(:,:), modezeff_tr(:,:)
1664 : 30 : complex(dpc),allocatable :: meff(:,:), meff_tr(:,:), modemeff_tr(:,:)
1665 : : !character(len=500) :: msg
1666 : :
1667 : : !TMP:
1668 : : ! complex(dpc) :: basein(2,2),baseout(2,2),hmat(2,2),rmat(2,2)
1669 : : ! complex(dpc) :: totbasein(3*natom,2), totbaseout(3*natom,2)
1670 : : complex*16,parameter :: ure=(1.d0,0.d0),uim=(0.d0,1.d0)
1671 : :
1672 : :
1673 : : ! *************************************************************************
1674 : :
1675 : : DBG_ENTER("COLL")
1676 : :
1677 : : !Rotate doubly degenerated modes
1678 : 30 : call alignph(amu,displ,blkval,mpert,natom,ntypat,phfrq,typat,silent=1)
1679 : :
1680 : : !Define the complex eigendisplacementes array
1681 [ + + ]: 570 : do imode=1,3*natom
1682 [ + + ]: 2190 : do idir1=1,3
1683 [ + + ]: 11880 : do ipert1=1,natom
1684 : 9720 : i1=idir1+(ipert1-1)*3
1685 : 9720 : index=i1+3*natom*(imode-1)
1686 : 11340 : modedisp(i1,imode)= cmplx(displ(2*index-1),displ(2*index),kind=dpc)
1687 : : end do
1688 : : end do
1689 : : end do
1690 : :
1691 : : ! ABI_MALLOC(norm,(3*natom))
1692 : : ! do imode= 1, 3*natom
1693 : : ! norm(imode)= sqrt(dot_product(modedisp(:,imode),modedisp(:,imode)))
1694 : : ! end do
1695 : :
1696 : : !Compute the mode-resolved macroscopic quantities
1697 : : !(Born and magnetic charges)
1698 : 30 : pdim= 3*natom
1699 [ - + + - : 150 : ABI_MALLOC(c_blkval,(3,mpert,3,mpert))
+ - - + +
- + - + -
- + - + ]
1700 [ - + + + : 193260 : c_blkval= cmplx(blkval(1,:,:,:,:,1),blkval(2,:,:,:,:,1),kind=dpc)
+ + + + +
+ ]
1701 : :
1702 : : !Born charges
1703 [ - + + - : 90 : ABI_MALLOC(zeff,(3,pdim))
+ - + - -
+ - + ]
1704 [ - + + - : 60 : ABI_MALLOC(zeff_tr,(pdim,3))
- + - + ]
1705 [ - + + - : 60 : ABI_MALLOC(modezeff_tr,(pdim,3))
- + - + ]
1706 [ - + + - : 60 : ABI_MALLOC(meff,(3,pdim))
- + - + ]
1707 [ - + + - : 60 : ABI_MALLOC(meff_tr,(pdim,3))
- + - + ]
1708 [ - + + - : 60 : ABI_MALLOC(modemeff_tr,(pdim,3))
- + - + ]
1709 [ + + + + ]: 2190 : modezeff(:,:)=(zero,zero)
1710 [ + + + + ]: 2190 : modemeff(:,:)=(zero,zero)
1711 [ + + + + ]: 1740 : modezeff_tr(:,:)=(zero,zero)
1712 [ + + + + ]: 1740 : modemeff_tr(:,:)=(zero,zero)
1713 : 30 : ipert= natom + 2
1714 : 30 : jpert= natom + 5
1715 [ + + ]: 120 : do im= 1, 3
1716 [ + + ]: 660 : do iat2= 1, natom
1717 [ + + ]: 2250 : do idir2= 1, 3
1718 : 1620 : imode= (iat2-1)*3 + idir2
1719 [ + + ]: 11880 : do iat1= 1, natom
1720 [ + + ]: 40500 : do idir1= 1, 3
1721 : 29160 : irow= (iat1-1)*3 + idir1
1722 : :
1723 : : !Electric Born charges
1724 : 29160 : zeff(im,irow)= c_blkval(im,ipert,idir1,iat1)
1725 : 29160 : zeff_tr(irow,im)= c_blkval(idir1,iat1,im,ipert)
1726 : : modezeff(im,imode)= modezeff(im,imode) + zeff(im,irow)* &
1727 : 29160 : & modedisp(irow,imode)
1728 : : modezeff_tr(imode,im)= modezeff_tr(imode,im) + zeff_tr(irow,im)* &
1729 : 29160 : & conjg(modedisp(irow,imode))
1730 : :
1731 : : !Magnetic Born charges
1732 : 29160 : meff(im,irow)= c_blkval(im,jpert,idir1,iat1)
1733 : 29160 : meff_tr(irow,im)= c_blkval(idir1,iat1,im,jpert)
1734 : : modemeff(im,imode)= modemeff(im,imode) + meff(im,irow)* &
1735 : 29160 : & modedisp(irow,imode)
1736 : : modemeff_tr(imode,im)= modemeff_tr(imode,im) + meff_tr(irow,im)* &
1737 : 38880 : & conjg(modedisp(irow,imode))
1738 : :
1739 : : end do
1740 : : end do
1741 : : end do
1742 : : end do
1743 : : end do
1744 [ - + ]: 30 : ABI_FREE(zeff)
1745 [ - + ]: 30 : ABI_FREE(zeff_tr)
1746 [ - + ]: 30 : ABI_FREE(meff)
1747 [ - + ]: 30 : ABI_FREE(meff_tr)
1748 : :
1749 : : !Compute the mode-resolved local magnetic moments and fields
1750 [ + + + + ]: 2730 : modemm(:,:)=(zero,zero)
1751 [ + + + + ]: 2730 : modezf(:,:)=(zero,zero)
1752 [ + + ]: 150 : do im= 1, ndim
1753 [ + + ]: 870 : do iat2= 1, natom
1754 [ + + ]: 3000 : do idir2= 1, 3
1755 : 2160 : imode= (iat2-1)*3 + idir2
1756 [ + + ]: 15840 : do iat1= 1, natom
1757 [ + + ]: 54000 : do idir1= 1, 3
1758 : 38880 : irow= (iat1-1)*3 + idir1
1759 : 38880 : modemm(im,imode)= modemm(im,imode) + mmom(im,irow)*modedisp(irow,imode)
1760 : 51840 : modezf(im,imode)= modezf(im,imode) + zfield(im,irow)*modedisp(irow,imode)
1761 : : end do
1762 : : end do
1763 : : end do
1764 : : end do
1765 : : end do
1766 : :
1767 : : !Compute the normal modes contribution to the susceptibilities
1768 : 30 : cplx_eta= cmplx(0.0_dp,eta)
1769 : 30 : cplx_w2= (omega+cplx_eta)**2
1770 : :
1771 : : !Dielectric tensor
1772 : 30 : fac= -four_pi/ucvol
1773 [ + + ]: 120 : do idir1= 1, 3
1774 [ + + ]: 390 : do idir2= 1, 3
1775 [ + + ]: 5220 : do imode= 1, pdim
1776 : : lm_epsilon_nm(idir1,idir2,imode)= fac*modezeff(idir1,imode)*modezeff_tr(imode,idir2)/ &
1777 : 5130 : & (cplx_w2 - phfrq(imode)**2)
1778 : : end do
1779 : : end do
1780 : : end do
1781 : :
1782 : : !Magnetoelectric susceptibility
1783 : 30 : fac= -one/ucvol
1784 [ + + ]: 120 : do idir1= 1, 3
1785 [ + + ]: 390 : do idir2= 1, 3
1786 [ + + ]: 5220 : do imode= 1, pdim
1787 : : lm_alpha_nm(idir1,idir2,imode)= fac*modemeff(idir1,imode)*modezeff_tr(imode,idir2)/ &
1788 : 5130 : & (cplx_w2 - phfrq(imode)**2)
1789 : : end do
1790 : : end do
1791 : : end do
1792 : :
1793 : : !Magnetic susceptibility
1794 : 120 : fac= -one/ucvol
1795 [ + + ]: 120 : do idir1= 1, 3
1796 [ + + ]: 390 : do idir2= 1, 3
1797 [ + + ]: 5220 : do imode= 1, pdim
1798 : : lm_mchi_nm(idir1,idir2,imode)= fac*modemeff(idir1,imode)*modemeff_tr(imode,idir2)/ &
1799 : 5130 : & (cplx_w2 - phfrq(imode)**2)
1800 : : end do
1801 : : end do
1802 : : end do
1803 : :
1804 : : !Normalize mode-projected quantities
1805 : : ! do imode= 1, natom*3
1806 : : ! modezeff(:,imode)= modezeff(:,imode) / norm(imode)
1807 : : ! modemeff(:,imode)= modemeff(:,imode) / norm(imode)
1808 : : ! modemm(:,imode)= modemm(:,imode) / norm(imode)
1809 : : ! modezf(:,imode)= modezf(:,imode) / norm(imode)
1810 : : ! end do
1811 : :
1812 : : !Apply a cell mass factor to the mode-resolved Born charges
1813 : 30 : fac= sqrt(mcell)
1814 [ + + + + ]: 2190 : modezeff= modezeff*fac
1815 [ + + + + ]: 2190 : modemeff= modemeff*fac
1816 [ + + + + ]: 2730 : modemm= modemm*fac
1817 [ + + + + ]: 2730 : modezf= modezf*fac
1818 : :
1819 : :
1820 : : ! ABI_FREE(norm)
1821 [ - + ]: 30 : ABI_FREE(c_blkval)
1822 [ - + ]: 30 : ABI_FREE(modezeff_tr)
1823 [ - + ]: 30 : ABI_FREE(modemeff_tr)
1824 : :
1825 : : DBG_EXIT("COLL")
1826 : :
1827 : 30 : end subroutine lm_normal_modes
1828 : : !!***
1829 : :
1830 : : !!****f* ABINIT/ri_d2etot
1831 : : !! NAME
1832 : : !! ri_d2etot
1833 : : !!
1834 : : !! FUNCTION
1835 : : !! Extracts the CI 2nd-order susceptibilities and calculates the corresponding
1836 : : !! lattice-mediated contributions
1837 : : !!
1838 : : !! COPYRIGHT
1839 : : !! Copyright (C) 2024 ABINIT group (FIXME: add author)
1840 : : !! This file is distributed under the terms of the
1841 : : !! GNU General Public License, see ~abinit/COPYING
1842 : : !! or http://www.gnu.org/copyleft/gpl.txt .
1843 : : !!
1844 : : !! INPUTS
1845 : : !! blkval(2,3,mpert,3,mpert,1)= array with second-order derivatives
1846 : : !! magsus(ndim,ndim)= local magnetic susceptibility (RS) or its inverse (FS)
1847 : : !! mpert= maximum number of perturbations
1848 : : !! mcoup(ndim,natom+5)= magnetic Zeman fields (FS) or moments (RS)
1849 : : !! mcoup_tr(natom+5,ndim)= hermitian conjugate of mmom
1850 : : !! natom= number of atoms in the cell
1851 : : !! ndim= dimension of the local magnetic degrees of freedom
1852 : : !! phongreen(3*natom,3*natom)= Phonon Green's function
1853 : : !!
1854 : : !! OUTPUT
1855 : : !! blkval_lm(2,3,mpert,3,mpert,1)= array with the lattice-mediated
1856 : : !! second-order derivatives
1857 : : !!
1858 : : !! SIDE EFFECTS
1859 : : !!
1860 : : !! NOTES
1861 : : !!
1862 : : !! PARENTS
1863 : : !!
1864 : : !! CHILDREN
1865 : : !!
1866 : : !! SOURCE
1867 : :
1868 : : #if defined HAVE_CONFIG_H
1869 : : #include "config.h"
1870 : : #endif
1871 : :
1872 : : #include "abi_common.h"
1873 : :
1874 : :
1875 : 30 : subroutine ri_d2etot(blkval,ci_alpha,ci_alpha_hc,ci_epsilon,ci_localpha,ci_locchi,ci_mchi,&
1876 : 30 : & lm_alpha,lm_alpha_hc,lm_epsilon,lm_localpha,lm_locchi,lm_magsus,lm_mchi,mpert,mcoup, &
1877 : 30 : & mcoup_tr,natom,ndim,phongreen,ucvol)
1878 : :
1879 : : !Arguments ------------------------------------
1880 : : !scalars
1881 : : integer, intent(in) :: mpert,natom,ndim
1882 : : real(dp), intent(in) :: ucvol
1883 : : !arrays
1884 : : real(dp), intent(in) :: blkval(2,3,mpert,3,mpert,1)
1885 : : complex(dpc), intent(out) :: ci_alpha(3,3)
1886 : : complex(dpc), intent(out) :: ci_alpha_hc(3,3)
1887 : : complex(dpc), intent(out) :: lm_alpha(3,3)
1888 : : complex(dpc), intent(out) :: lm_alpha_hc(3,3)
1889 : : complex(dpc), intent(out) :: ci_localpha(ndim,3)
1890 : : complex(dpc), intent(out) :: lm_localpha(ndim,3)
1891 : : complex(dpc), intent(out) :: ci_locchi(ndim,3)
1892 : : complex(dpc), intent(out) :: lm_locchi(ndim,3)
1893 : : complex(dpc), intent(out) :: ci_epsilon(3,3)
1894 : : complex(dpc), intent(out) :: lm_epsilon(3,3)
1895 : : complex(dpc), intent(out) :: ci_mchi(3,3)
1896 : : complex(dpc), intent(out) :: lm_mchi(3,3)
1897 : : complex(dpc), intent(out) :: lm_magsus(ndim,ndim)
1898 : : complex(dpc), intent(in) :: mcoup(ndim,(natom+5)*3)
1899 : : complex(dpc), intent(in) :: mcoup_tr((natom+5)*3,ndim)
1900 : : complex(dpc), intent(in) :: phongreen(3*natom,3*natom)
1901 : :
1902 : : !Local variables-------------------------------
1903 : : !scalars
1904 : : integer :: icol,idir1,idir2,ipert1,ipert2,irow
1905 : : integer :: pdim
1906 : : real(dp) :: fac
1907 : : !arrays
1908 : 30 : complex(dpc),allocatable :: c_blkval(:,:,:,:)
1909 : 30 : complex(dpc),allocatable :: coup(:,:),coup_tr(:,:)
1910 : :
1911 : : ! *************************************************************************
1912 : :
1913 : : DBG_ENTER("COLL")
1914 : :
1915 : 30 : pdim= 3*natom
1916 [ - + + - : 150 : ABI_MALLOC(c_blkval,(3,mpert,3,mpert))
+ - - + +
- + - + -
- + - + ]
1917 [ - + + + : 193260 : c_blkval= cmplx(blkval(1,:,:,:,:,1),blkval(2,:,:,:,:,1),kind=dpc)
+ + + + +
+ ]
1918 : :
1919 : : !Dielectric tensor
1920 : 30 : fac= -four_pi/ucvol
1921 [ - + + - : 90 : ABI_MALLOC(coup,(3,pdim))
+ - + - -
+ - + ]
1922 [ - + + - : 60 : ABI_MALLOC(coup_tr,(pdim,3))
- + - + ]
1923 : 30 : ipert1= natom + 2
1924 [ + + ]: 120 : do idir1= 1, 3
1925 : : irow= idir1
1926 [ + + ]: 660 : do ipert2= 1, natom
1927 [ + + ]: 2250 : do idir2= 1, 3
1928 : 1620 : icol= (ipert2-1)*3 + idir2
1929 : 1620 : coup(irow,icol)= c_blkval(idir1,ipert1,idir2,ipert2)
1930 : 2160 : coup_tr(icol,irow)= c_blkval(idir2,ipert2,idir1,ipert1)
1931 : : end do
1932 : : end do
1933 : : end do
1934 : 30 : ipert2= natom + 2
1935 [ - + + - : 44520 : lm_epsilon= fac*matmul(coup(:,:),matmul(phongreen,coup_tr(:,:)))
+ - + - +
- - + - +
+ + + + +
+ + + + +
+ - + + +
+ + + + -
+ - + - +
+ + + ]
1936 : :
1937 [ + + ]: 120 : do idir1= 1, 3
1938 [ + + ]: 390 : do idir2= 1, 3
1939 : 360 : ci_epsilon(idir1,idir2)= c_blkval(idir1,ipert1,idir2,ipert2)
1940 : : end do
1941 : : end do
1942 : :
1943 : : !Magnetoelectric susceptibility
1944 : : !i) mu_0 dM/dE
1945 : 30 : fac= -one/ucvol
1946 : 30 : ipert1= natom + 5
1947 [ + + ]: 120 : do idir1= 1, 3
1948 : : irow= idir1
1949 [ + + ]: 660 : do ipert2= 1, natom
1950 [ + + ]: 2250 : do idir2= 1, 3
1951 : 1620 : icol= (ipert2-1)*3 + idir2
1952 : 2160 : coup(irow,icol)= c_blkval(idir1,ipert1,idir2,ipert2)
1953 : : end do
1954 : : end do
1955 : : end do
1956 : :
1957 : 120 : ipert1= natom + 2
1958 [ + + ]: 120 : do idir1= 1, 3
1959 : : irow= idir1
1960 [ + + ]: 660 : do ipert2= 1, natom
1961 [ + + ]: 2250 : do idir2= 1, 3
1962 : 1620 : icol= (ipert2-1)*3 + idir2
1963 : 2160 : coup_tr(icol,irow)= c_blkval(idir2,ipert2,idir1,ipert1)
1964 : : end do
1965 : : end do
1966 : : end do
1967 : :
1968 [ - + + - : 44520 : lm_alpha= fac*matmul(coup(:,:),matmul(phongreen,coup_tr(:,:)))
+ - + - +
- - + - +
+ + + + +
+ + + + +
+ - + + +
+ + + + -
+ - + - +
+ + + ]
1969 : :
1970 : 120 : ipert1= natom + 5
1971 : 120 : ipert2= natom + 2
1972 [ + + ]: 120 : do idir1= 1, 3
1973 [ + + ]: 390 : do idir2= 1, 3
1974 : 360 : ci_alpha(idir1,idir2)= c_blkval(idir1,ipert1,idir2,ipert2)
1975 : : end do
1976 : : end do
1977 : :
1978 : : !ii) dP/dH
1979 : 120 : fac= -one/ucvol
1980 : 120 : ipert1= natom + 2
1981 [ + + ]: 120 : do idir1= 1, 3
1982 : : irow= idir1
1983 [ + + ]: 660 : do ipert2= 1, natom
1984 [ + + ]: 2250 : do idir2= 1, 3
1985 : 1620 : icol= (ipert2-1)*3 + idir2
1986 : 2160 : coup(irow,icol)= c_blkval(idir1,ipert1,idir2,ipert2)
1987 : : end do
1988 : : end do
1989 : : end do
1990 : :
1991 : 120 : ipert1= natom + 5
1992 [ + + ]: 120 : do idir1= 1, 3
1993 : : irow= idir1
1994 [ + + ]: 660 : do ipert2= 1, natom
1995 [ + + ]: 2250 : do idir2= 1, 3
1996 : 1620 : icol= (ipert2-1)*3 + idir2
1997 : 2160 : coup_tr(icol,irow)= c_blkval(idir2,ipert2,idir1,ipert1)
1998 : : end do
1999 : : end do
2000 : : end do
2001 : :
2002 [ - + + - : 44550 : lm_alpha_hc= fac*matmul(coup(:,:),matmul(phongreen,coup_tr(:,:)))
+ - + - +
- - + - +
+ + + + +
+ + + + +
+ - + + +
+ + + + -
+ - + - +
+ + + ]
2003 : :
2004 : 120 : ipert1= natom + 2
2005 : 120 : ipert2= natom + 5
2006 [ + + ]: 120 : do idir1= 1, 3
2007 [ + + ]: 390 : do idir2= 1, 3
2008 : 360 : ci_alpha_hc(idir1,idir2)= c_blkval(idir1,ipert1,idir2,ipert2)
2009 : : end do
2010 : : end do
2011 : :
2012 : : !Local magnetoelectric susceptibilty
2013 [ + + + + ]: 480 : ci_localpha(:,:)= mcoup(:,(natom+1)*3+1:(natom+2)*3)
2014 [ - + + - : 48450 : lm_localpha(:,:)= -matmul(mcoup(:,1:natom*3),matmul(phongreen,coup_tr(:,:)))
+ - + - +
- - + - +
+ + + + +
+ + + + +
+ - - + +
- - + + -
+ - + - -
+ - + + +
+ + + + +
+ + + + -
+ - + - +
+ + + +
- ]
2015 : :
2016 : : !Magnetic susceptibility
2017 : 120 : fac= -one/ucvol
2018 : 120 : ipert1= natom + 5
2019 [ + + ]: 120 : do idir1= 1, 3
2020 : : irow= idir1
2021 [ + + ]: 660 : do ipert2= 1, natom
2022 [ + + ]: 2250 : do idir2= 1, 3
2023 : 1620 : icol= (ipert2-1)*3 + idir2
2024 : 1620 : coup(irow,icol)= c_blkval(idir1,ipert1,idir2,ipert2)
2025 : 2160 : coup_tr(icol,irow)= c_blkval(idir2,ipert2,idir1,ipert1)
2026 : : end do
2027 : : end do
2028 : : end do
2029 : 30 : ipert2= natom + 5
2030 [ - + + - : 44550 : lm_mchi= fac*matmul(coup(:,:),matmul(phongreen,coup_tr(:,:)))
+ - + - +
- - + - +
+ + + + +
+ + + + +
+ - + + +
+ + + + -
+ - + - +
+ + + ]
2031 : :
2032 [ + + ]: 120 : do idir1= 1, 3
2033 [ + + ]: 390 : do idir2= 1, 3
2034 : 360 : ci_mchi(idir1,idir2)= c_blkval(idir1,ipert1,idir2,ipert2)
2035 : : end do
2036 : : end do
2037 : :
2038 : : !Local magnetic moments induced by macroscopic Zeeman
2039 [ + + + + ]: 480 : ci_locchi(:,:)= mcoup(:,(natom+4)*3+1:(natom+5)*3)
2040 [ - + + - : 48450 : lm_locchi(:,:)= -matmul(mcoup(:,1:natom*3),matmul(phongreen,coup_tr(:,:)))
+ - + - +
- - + - +
+ + + + +
+ + + + +
+ - - + +
- - + + -
+ - + - -
+ - + + +
+ + + + +
+ + + + -
+ - + - +
+ + + +
- ]
2041 : :
2042 : : !Local magnetic susceptibility
2043 [ - + + - : 64530 : lm_magsus(:,:)= -matmul(mcoup(:,1:natom*3),matmul(phongreen,mcoup_tr(1:natom*3,:)))
- + + - +
- + - - +
- + + + +
+ + + + +
+ + + - -
+ + - - +
+ - + - +
- - + - +
+ + + + +
+ + + + +
+ - + - +
- + + + +
+ - ]
2044 : :
2045 [ - + ]: 30 : ABI_FREE(c_blkval)
2046 [ - + ]: 30 : ABI_FREE(coup)
2047 [ - + ]: 30 : ABI_FREE(coup_tr)
2048 : :
2049 : : DBG_EXIT("COLL")
2050 : :
2051 : 30 : end subroutine ri_d2etot
2052 : : !!***
2053 : :
2054 [ + - + - : 210 : end module m_ddb_omega_interpol
+ - + - +
- + - +
- ]
2055 : : !!***
|
