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