Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_phi2
2 : : !! NAME
3 : : !! m_tdep_phi2
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the data types holding the second-order IFC
7 : : !! (Phi2_type) and the phonon frequencies and eigenvectors (Eigen_type),
8 : : !! and the routines writing the first- and second-order IFC, building the
9 : : !! dynamical matrix and writing the phonon properties.
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_tdep_phi2
26 : :
27 : : use defs_basis
28 : : use m_errors
29 : : use m_abicore
30 : : use m_xmpi
31 : : use m_io_tools
32 : : use m_sort, only : sort_dp
33 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
34 : : use m_tdep_shell, only : Shell_type
35 : : use m_tdep_sym, only : Symmetries_type
36 : : use m_tdep_qpt, only : Qpoints_type
37 : : use m_tdep_sampling, only : tdep_Sampling_type
38 : : use m_tdep_supercell, only : tdep_Supercell_type
39 : :
40 : : implicit none
41 : : !!***
42 : :
43 : : !----------------------------------------------------------------------
44 : : !!****t* m_tdep_phi2/Phi2_type
45 : : !! NAME
46 : : !! Phi2_type
47 : : !!
48 : : !! FUNCTION
49 : : !! Second-order IFC
50 : : !!
51 : : !! SOURCE
52 : :
53 : : type Phi2_type
54 : :
55 : : double precision, allocatable :: SR(:,:)
56 : : ! SR(3*natom, 3*natom)
57 : : ! Short-range part of the second-order IFC
58 : :
59 : : double precision, allocatable :: LR(:,:)
60 : : ! LR(3*natom, 3*natom)
61 : : ! Long-range part of the second-order IFC
62 : :
63 : : double precision, allocatable :: Tot(:,:)
64 : : ! Tot(3*natom, 3*natom)
65 : : ! Second-order IFC (SR + LR)
66 : :
67 : : contains
68 : :
69 : : procedure :: init => tdep_init_phi2
70 : : procedure :: free => tdep_destroy_phi2
71 : :
72 : : end type Phi2_type
73 : : !!***
74 : :
75 : : !----------------------------------------------------------------------
76 : :
77 : : !!****t* m_tdep_phi2/Eigen_type
78 : : !! NAME
79 : : !! Eigen_type
80 : : !!
81 : : !! FUNCTION
82 : : !! Holds the phonon frequencies, eigenvectors and dynamical matrices
83 : : !! on a set of q-points.
84 : : !!
85 : : !! SOURCE
86 : :
87 : : type Eigen_type
88 : :
89 : : double precision, allocatable :: eigenval(:,:)
90 : : ! eigenval(3*natom_unitcell,nqpt)
91 : : ! Phonon frequencies (Ha)
92 : :
93 : : double precision, allocatable :: eigenvec(:,:,:,:,:,:)
94 : : ! eigenvec(2,3,natom_unitcell,3,natom_unitcell,nqpt)
95 : : ! Phonon eigenvectors (real and imaginary parts)
96 : :
97 : : double precision, allocatable :: dynmat(:,:,:,:,:,:)
98 : : ! dynmat(2,3,natom_unitcell,3,natom_unitcell,nqpt)
99 : : ! Dynamical matrices (real and imaginary parts)
100 : :
101 : : contains
102 : :
103 : : procedure :: init => tdep_init_eigen2nd
104 : : procedure :: free => tdep_destroy_eigen2nd
105 : :
106 : : end type Eigen_type
107 : : !!***
108 : :
109 : : public :: tdep_write_phi1
110 : : public :: tdep_write_phi2
111 : : public :: tdep_init_phi2
112 : : public :: tdep_destroy_phi2
113 : : public :: tdep_build_phi2_33
114 : : public :: tdep_calc_dij
115 : : public :: tdep_write_dij
116 : : public :: tdep_init_eigen2nd
117 : : public :: tdep_destroy_eigen2nd
118 : : public :: tdep_write_yaml
119 : :
120 : : contains
121 : :
122 : : !=====================================================================================================
123 : :
124 : : !!****f* ABINIT/m_tdep_phi2/tdep_write_phi1
125 : : !! NAME
126 : : !! tdep_write_phi1
127 : : !!
128 : : !! FUNCTION
129 : : !! Write the first-order IFC of each atom of the unitcell in the output file.
130 : : !!
131 : : !! INPUTS
132 : : !! Phi1(3*natom) = First-order IFC.
133 : : !! natom = Number of atoms in the supercell.
134 : : !! natom_unitcell = Number of atoms in the unitcell.
135 : : !!
136 : : !! OUTPUT
137 : : !! Only writing.
138 : : !!
139 : : !! SOURCE
140 : :
141 : 46 : subroutine tdep_write_phi1(Phi1, natom, natom_unitcell)
142 : :
143 : : integer,intent(in) :: natom
144 : : integer,intent(in) :: natom_unitcell
145 : : real(dp),intent(in) :: Phi1(3*natom)
146 : :
147 : : integer :: iatcell,ii
148 : :
149 : 46 : write(ab_out,*) ' '
150 : 46 : write(ab_out,*) '#############################################################################'
151 : 46 : write(ab_out,*) '#### For each shell, list of coefficients (IFC), number of neighbours... ####'
152 : 46 : write(ab_out,*) '#############################################################################'
153 : :
154 : : ! Write the IFCs in the data.out file (with others specifications:
155 : : ! number of atoms in a shell, distance, Trace...)
156 [ + + ]: 156 : do iatcell=1,natom_unitcell
157 : 110 : write(ab_out,'(a,i4)') ' ############# List of (first order) IFC for the reference atom=',iatcell
158 [ + - + + ]: 440 : write(ab_out,'(2x,3(f9.6,1x))') (Phi1((iatcell-1)*3+ii),ii=1,3)
159 : 156 : write(ab_out,*) ' '
160 : : end do !iatcell
161 : :
162 : 46 : end subroutine tdep_write_phi1
163 : : !!***
164 : :
165 : : !=====================================================================================================
166 : :
167 : : !!****f* ABINIT/m_tdep_phi2/tdep_write_phi2
168 : : !! NAME
169 : : !! tdep_write_phi2
170 : : !!
171 : : !! FUNCTION
172 : : !! Write the second-order IFC in the output file: for each atom of the
173 : : !! unitcell, and for each shell (sorted by increasing distance), the 3x3 IFC
174 : : !! of every atom of the shell, the interatomic vector and the trace of the IFC.
175 : : !! If prtvol > 1, also write the IFC in the _Phi2_unitcell.dat and _Phi2.dat files.
176 : : !!
177 : : !! INPUTS
178 : : !! Invar = Input object containing the input variables.
179 : : !! MPIdata = Info on MPI parallelism.
180 : : !! SC = Supercell object, used for the interatomic distances (SC%dij).
181 : : !! Phi2(3*natom,3*natom) = Second-order IFC, in the supercell.
182 : : !! Shell2at = Shell object for second-order interactions.
183 : : !!
184 : : !! OUTPUT
185 : : !! Only writing.
186 : : !!
187 : : !! SOURCE
188 : :
189 : 46 : subroutine tdep_write_phi2(Invar,MPIdata,SC,Phi2,Shell2at)
190 : :
191 : : type(atdep_dataset_type),intent(in) :: Invar
192 : : type(MPI_enreg_type), intent(in) :: MPIdata
193 : : type(tdep_Supercell_type),intent(in) :: SC
194 : : type(Shell_type),intent(in) :: Shell2at
195 : : double precision,intent(in) :: Phi2(3*Invar%natom,3*Invar%natom)
196 : :
197 : : logical :: debug
198 : : integer :: iatcell,ishell,jatom
199 : : integer :: nshell,ii,this_shell,iatshell
200 : : double precision :: dist_shell,tmp1,tmp2,tmp3
201 : 46 : integer,allocatable :: iperm(:)
202 : 46 : double precision,allocatable :: dist(:)
203 : :
204 : 46 : debug = (Invar%prtvol > 1)
205 : :
206 : 46 : write(ab_out,*) ' '
207 : 46 : write(ab_out,*) '#############################################################################'
208 : 46 : write(ab_out,*) '#### For each shell, list of coefficients (IFC), number of neighbours... ####'
209 : 46 : write(ab_out,*) '#############################################################################'
210 : :
211 : 46 : nshell=Shell2at%nshell
212 : : ! Write the IFCs in the data.out file (with others specifications:
213 : : ! number of atoms in a shell, distance, Trace...)
214 : : ! Sort the shells with distance in increasing order.
215 : : ! Shells at the same distance (within tol6) are kept in the order of their index.
216 [ - + + - : 138 : ABI_MALLOC(dist,(nshell))
+ - + - -
+ - + ]
217 [ - + + - : 138 : ABI_MALLOC(iperm,(nshell))
+ - + - -
+ - + ]
218 [ + + ]: 568 : do ishell=1,nshell
219 : 522 : dist(ishell)=SC%dij(Shell2at%iatref(ishell),Shell2at%jatref(ishell),1)
220 : 568 : iperm(ishell)=ishell
221 : : end do
222 : 46 : call sort_dp(nshell,dist,iperm,tol6)
223 [ + + ]: 156 : do iatcell=1,Invar%natom_unitcell
224 : 110 : write(ab_out,'(a,i4)') ' ############# List of (second order) IFC for the reference atom=',iatcell
225 [ + + ]: 2486 : do ishell=1,nshell
226 : 2330 : this_shell=iperm(ishell)
227 : 2330 : dist_shell=dist(ishell)
228 : :
229 : : ! Write the IFC properly
230 [ + + ]: 2440 : if (Shell2at%neighbours(iatcell,this_shell)%n_interactions.ne.0) then
231 : 1143 : write(ab_out,'(a,i4,a,i4,a,f9.6)') ' ======== NEW SHELL (ishell=',this_shell,&
232 : 2286 : & '): There are',Shell2at%neighbours(iatcell,this_shell)%n_interactions,' atoms on this shell at distance=',dist_shell
233 [ + + ]: 6819 : do iatshell=1,Shell2at%neighbours(iatcell,this_shell)%n_interactions
234 : 5676 : jatom=Shell2at%neighbours(iatcell,this_shell)%atomj_in_shell(iatshell)
235 : 5676 : write(ab_out,'(a,i4,a,i4)') ' For jatom=',jatom,' ,with type=',mod(jatom-1,Invar%natom_unitcell)+1
236 [ + + ]: 22704 : do ii=1,3
237 [ + + ]: 17028 : if (abs(Phi2((iatcell-1)*3+ii,(jatom-1)*3+1)).lt.5.d-7) then
238 : 2772 : tmp1=0.d0
239 : : else
240 : 14256 : tmp1=Phi2((iatcell-1)*3+ii,(jatom-1)*3+1)
241 : : end if
242 [ + + ]: 17028 : if (abs(Phi2((iatcell-1)*3+ii,(jatom-1)*3+2)).lt.5.d-7) then
243 : 2718 : tmp2=0.d0
244 : : else
245 : 14310 : tmp2=Phi2((iatcell-1)*3+ii,(jatom-1)*3+2)
246 : : end if
247 [ + + ]: 17028 : if (abs(Phi2((iatcell-1)*3+ii,(jatom-1)*3+3)).lt.5.d-7) then
248 : 2992 : tmp3=0.d0
249 : : else
250 : 14036 : tmp3=Phi2((iatcell-1)*3+ii,(jatom-1)*3+3)
251 : : end if
252 : 22704 : write(ab_out,'(2x,3(f9.6,1x))') tmp1,tmp2,tmp3
253 : : end do
254 : 5676 : write(ab_out,'(a,3(1x,f11.6))') ' The components of the vector are:', SC%dij(iatcell,jatom,2:4)
255 : 5676 : write(ab_out,'(a,(1x,f9.6))') ' Trace=',Phi2((iatcell-1)*3+1,(jatom-1)*3+1)+Phi2((iatcell-1)*3+2,&
256 : 11352 : & (jatom-1)*3+2)+Phi2((iatcell-1)*3+3,(jatom-1)*3+3)
257 : 6819 : write(ab_out,*) ' '
258 : : end do
259 : : end if
260 : : end do !ishell
261 : : end do !iatcell
262 [ - + ]: 46 : ABI_FREE(dist)
263 [ - + ]: 46 : ABI_FREE(iperm)
264 : :
265 : : ! Write the Phi2_unitcell.dat and Phi2.dat files
266 [ - + - - ]: 46 : if (debug.and.MPIdata%iam_master) then
267 : 0 : write(ab_out,'(a)') ' See the Phi2*.dat file'
268 : 0 : open(unit=52,file=trim(Invar%output_prefix)//'_Phi2_unitcell.dat')
269 : 0 : open(unit=55,file=trim(Invar%output_prefix)//'_Phi2.dat')
270 [ # # ]: 0 : do jatom=1,3*Invar%natom
271 [ # # ]: 0 : if (jatom.le.3*Invar%natom_unitcell) then
272 : 0 : write(52,'(10000(f10.6,1x))') Phi2(jatom,:)
273 : : end if
274 : 0 : write(55,'(10000(f10.6,1x))') Phi2(jatom,:)
275 : : end do
276 : 0 : close(52)
277 : 0 : close(55)
278 : : end if
279 : :
280 : 46 : end subroutine tdep_write_phi2
281 : : !!***
282 : :
283 : : !=====================================================================================================
284 : :
285 : : !!****f* ABINIT/m_tdep_phi2/tdep_init_phi2
286 : : !! NAME
287 : : !! tdep_init_phi2
288 : : !!
289 : : !! FUNCTION
290 : : !! Initialize the Phi2 object: allocate (and set to zero) the short-range
291 : : !! part of the second-order IFC and, if dipdip is true, its long-range part
292 : : !! and the total IFC.
293 : : !!
294 : : !! INPUTS
295 : : !! dipdip = Whether the dipole-dipole interaction is included.
296 : : !! natom = Number of atoms in the supercell.
297 : : !!
298 : : !! OUTPUT
299 : : !!
300 : : !! SIDE EFFECTS
301 : : !! Phi2 = Second-order IFC object, initialized.
302 : : !!
303 : : !! SOURCE
304 : :
305 : 48 : subroutine tdep_init_phi2(Phi2,dipdip,natom)
306 : :
307 : : class(Phi2_type),intent(inout) :: Phi2
308 : : logical, intent(in) :: dipdip
309 : : integer, intent(in) :: natom
310 : :
311 [ + - - + : 14237442 : ABI_CALLOC(Phi2%SR ,(3*natom,3*natom))
+ - - + +
- + - + -
- + - + -
+ + + +
+ ]
312 [ + + ]: 48 : if (dipdip) then
313 [ - + + - : 923727 : ABI_CALLOC(Phi2%Tot,(3*natom,3*natom))
- + + - -
+ - + - +
+ + + + ]
314 [ - + + - : 923727 : ABI_CALLOC(Phi2%LR ,(3*natom,3*natom))
- + + - -
+ - + - +
+ + + + ]
315 : : end if
316 : :
317 : 48 : end subroutine tdep_init_phi2
318 : : !!***
319 : :
320 : : !=====================================================================================================
321 : :
322 : : !!****f* ABINIT/m_tdep_phi2/tdep_destroy_phi2
323 : : !! NAME
324 : : !! tdep_destroy_phi2
325 : : !!
326 : : !! FUNCTION
327 : : !! Free the memory allocated in the Phi2 object.
328 : : !!
329 : : !! SIDE EFFECTS
330 : : !! Phi2 = Second-order IFC object. All its allocated arrays are deallocated.
331 : : !!
332 : : !! SOURCE
333 : :
334 : 48 : subroutine tdep_destroy_phi2(Phi2)
335 : :
336 : : class(Phi2_type),intent(inout) :: Phi2
337 : :
338 [ - + ]: 48 : ABI_FREE(Phi2%SR)
339 [ + + ]: 48 : ABI_SFREE(Phi2%Tot)
340 [ + + ]: 48 : ABI_SFREE(Phi2%LR)
341 : :
342 : 48 : end subroutine tdep_destroy_phi2
343 : : !!***
344 : :
345 : : !=====================================================================================================
346 : :
347 : : !!****f* ABINIT/m_tdep_phi2/tdep_calc_dij
348 : : !! NAME
349 : : !! tdep_calc_dij
350 : : !!
351 : : !! FUNCTION
352 : : !! Compute the dynamical matrix at a given q-point from the second-order IFC:
353 : : !! D_{ij}(q) = sum_{j'} Phi2_{ij'} exp(2 i pi q.R_{ij'})
354 : : !! where j' runs over the periodic images of the unitcell atom j in the
355 : : !! supercell. Check that it is hermitian, then diagonalize D_{ij}/sqrt(M_i M_j)
356 : : !! to obtain the phonon frequencies (negative for imaginary frequencies)
357 : : !! and the normalized eigenvectors.
358 : : !!
359 : : !! INPUTS
360 : : !! iqpt = Index of the q-point (the hermiticity warning is only printed for iqpt=1).
361 : : !! Invar = Input object containing the input variables.
362 : : !! Phi2(3*natom,3*natom) = Second-order IFC, in the supercell.
363 : : !! qpt_cart(3) = Cartesian coordinates of the q-point.
364 : : !! Rlatt_cart(3,natom_unitcell,natom) = Cartesian lattice vector between
365 : : !! a unitcell reference atom and the
366 : : !! ideal site of each atom.
367 : : !!
368 : : !! OUTPUT
369 : : !! dij(3*natom_unitcell,3*natom_unitcell) = Dynamical matrix (not divided by
370 : : !! the masses). Its contribution is
371 : : !! added to the input value, which
372 : : !! must therefore be set to zero by
373 : : !! the caller.
374 : : !! eigenV(3*natom_unitcell,3*natom_unitcell) = Phonon eigenvectors.
375 : : !! omega(3*natom_unitcell) = Phonon frequencies (Ha).
376 : : !!
377 : : !! SOURCE
378 : :
379 : 0 : subroutine tdep_calc_dij(dij,eigenV,iqpt,Invar,omega,Phi2,qpt_cart,Rlatt_cart)
380 : :
381 : : type(atdep_dataset_type),intent(in) :: Invar
382 : : integer,intent(in) :: iqpt
383 : : double precision,intent(in) :: Phi2(3*Invar%natom,3*Invar%natom)
384 : : double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
385 : : double precision,intent(in) :: qpt_cart(3)
386 : : double precision,intent(out) :: omega (3*Invar%natom_unitcell)
387 : : double complex ,intent(out) :: dij (3*Invar%natom_unitcell,3*Invar%natom_unitcell)
388 : : double complex ,intent(out) :: eigenV(3*Invar%natom_unitcell,3*Invar%natom_unitcell)
389 : :
390 : : logical :: debug
391 : : integer :: LWORK,ii,jj,kk,iatom,jatom,iatcell,jatcell,itypat,jtypat,iat_mod,INFO,itemp,imode,nmode
392 : : double precision :: phase
393 : : double complex :: norm
394 : 0 : double precision, allocatable :: RWORK(:)
395 : 0 : double complex, allocatable :: WORKC(:)
396 : : ! double complex, allocatable :: mass_mat(:,:)
397 : :
398 : 0 : debug = (Invar%prtvol > 1)
399 : :
400 : : ! Calculation of the dynamical matrix (Dij)
401 [ # # ]: 0 : do iatcell=1,Invar%natom_unitcell
402 [ # # ]: 0 : do jatom=1,Invar%natom
403 : 0 : iat_mod=mod(jatom+Invar%natom_unitcell-1,Invar%natom_unitcell)+1
404 : 0 : phase=0.d0
405 [ # # ]: 0 : do kk=1,3
406 : 0 : phase=phase+2*pi*Rlatt_cart(kk,iatcell,jatom)*qpt_cart(kk)
407 : : end do
408 [ # # ]: 0 : do ii=1+(iatcell-1)*3,3+(iatcell-1)*3
409 [ # # ]: 0 : do jj=1,3
410 : 0 : dij(ii,3*(iat_mod-1)+jj)=dij(ii,3*(iat_mod-1)+jj)+dcmplx(Phi2(ii,3*(jatom-1)+jj),0.d0)*exp(dcmplx(0.d0,phase))
411 : : end do !jj
412 : : end do !ii
413 : : end do !jatom
414 : : end do !iatcell
415 : :
416 : : ! The Dij has to be an hermitian matrix
417 : 0 : itemp=0
418 [ # # ]: 0 : do ii=1,3*Invar%natom_unitcell
419 [ # # ]: 0 : do jj=ii,3*Invar%natom_unitcell
420 [ # # # # ]: 0 : if ((abs(real(dij(ii,jj))-real(dij(jj,ii))).gt.tol10).or.(abs(aimag(dij(ii,jj))+aimag(dij(jj,ii))).gt.tol10)) then
421 [ # # ]: 0 : if (debug) then
422 : 0 : write (ab_out,'(a,1x,2(i4,1x))') 'for ii,jj=',ii,jj
423 : 0 : write (ab_out,'(a,1x,1(f12.8,1x))') 'abs(realij-realji)=',abs(real(dij(ii,jj))-real(dij(jj,ii)))
424 : 0 : write (ab_out,'(a,1x,1(f12.8,1x))') 'abs(imagij+imagji)=',abs(aimag(dij(ii,jj))+aimag(dij(jj,ii)))
425 : : end if
426 : 0 : itemp=itemp+1
427 : : end if
428 : : end do
429 : : end do
430 [ # # # # ]: 0 : if (itemp.ne.0.and.iqpt.eq.1) then
431 : 0 : write(ab_out,*) 'WARNING: The Dij matrix is not hermitian'
432 : 0 : write(ab_out,*) ' Probably: one shell may not have the whole number of atoms'
433 : 0 : write(ab_out,*) ' The Dij matrix is symetrized'
434 : : end if
435 : :
436 : : ! Diagonalization of dynamical matrix Dij/sqrt(Mi*Mj)
437 : 0 : LWORK=2*3*Invar%natom_unitcell-1
438 [ # # # # : 0 : ABI_MALLOC(WORKC,(LWORK)); WORKC(:)=czero
# # # # #
# # # #
# ]
439 [ # # # # : 0 : ABI_MALLOC(RWORK,(3*3*Invar%natom_unitcell-2)); RWORK(:)=zero
# # # # #
# # # # #
# # ]
440 [ # # ]: 0 : do iatcell=1,Invar%natom_unitcell
441 : 0 : itypat=Invar%typat_unitcell(iatcell)
442 [ # # ]: 0 : do jatcell=1,Invar%natom_unitcell
443 : 0 : jtypat=Invar%typat_unitcell(jatcell)
444 [ # # ]: 0 : do ii=1,3
445 [ # # ]: 0 : do jj=1,3
446 : : eigenV(ii+(iatcell-1)*3,jj+(jatcell-1)*3)=dij(ii+(iatcell-1)*3,jj+(jatcell-1)*3)/&
447 : 0 : & dcmplx(dsqrt(Invar%amu(itypat)*Invar%amu(jtypat))*amu_emass,0.d0)
448 : : end do !jj
449 : : end do !ii
450 : : end do !jatcell
451 : : end do !iatcell
452 : 0 : call ZHEEV('V','U',3*Invar%natom_unitcell,eigenV(:,:),3*Invar%natom_unitcell,omega(:),WORKC,LWORK,RWORK,INFO)
453 : :
454 : : ! Normalization of the eigenvectors
455 : 0 : nmode=3*Invar%natom_unitcell
456 [ # # ]: 0 : do imode=1,nmode
457 : : norm=zero
458 [ # # ]: 0 : do iatom=1,Invar%natom_unitcell
459 [ # # ]: 0 : do ii=1,3
460 : 0 : norm=norm+eigenV(3*(iatom-1)+ii,imode)*conjg(eigenV(3*(iatom-1)+ii,imode))
461 : : end do
462 : : end do
463 [ # # ]: 0 : eigenV(:,imode)=eigenV(:,imode)/dsqrt(real(norm))
464 : : end do
465 : :
466 : : ! Remove the squared-negative frequencies
467 [ # # ]: 0 : do ii=1,Invar%natom_unitcell
468 [ # # ]: 0 : do jj=1,3
469 [ # # ]: 0 : if (omega((ii-1)*3+jj).lt.0.d0) then
470 : 0 : omega((ii-1)*3+jj)=-dsqrt(-omega((ii-1)*3+jj))
471 : : else
472 : 0 : omega((ii-1)*3+jj)=dsqrt(omega((ii-1)*3+jj))
473 : : end if
474 : : end do
475 : : end do
476 [ # # ]: 0 : ABI_FREE(WORKC)
477 [ # # ]: 0 : ABI_FREE(RWORK)
478 : :
479 : 0 : end subroutine tdep_calc_dij
480 : : !!***
481 : :
482 : : !=====================================================================================================
483 : : !FB subroutine tdep_write_dij(Eigen2nd,iqpt,Invar,qpt_cart)
484 : :
485 : : !!****f* ABINIT/m_tdep_phi2/tdep_write_dij
486 : : !! NAME
487 : : !! tdep_write_dij
488 : : !!
489 : : !! FUNCTION
490 : : !! Write, for a given q-point, the dynamical matrix (on unit 52), the phonon
491 : : !! frequencies (on unit 53, in the unit selected by Invar%enunit) and the
492 : : !! normalized phonon eigenvectors (on unit 51). These units must be opened by
493 : : !! the caller (see tdep_calc_phdos).
494 : : !!
495 : : !! INPUTS
496 : : !! Eigen2nd = Eigen object containing the phonon frequencies, eigenvectors
497 : : !! and dynamical matrices.
498 : : !! iqpt = Index of the q-point in Eigen2nd.
499 : : !! Invar = Input object containing the input variables.
500 : : !! qpt(3) = Reduced coordinates of the q-point.
501 : : !!
502 : : !! OUTPUT
503 : : !! Only writing.
504 : : !!
505 : : !! SOURCE
506 : :
507 : 9290 : subroutine tdep_write_dij(Eigen2nd,iqpt,Invar,qpt)
508 : :
509 : : type(Eigen_type),intent(in) :: Eigen2nd
510 : : type(atdep_dataset_type),intent(in) :: Invar
511 : : integer,intent(in) :: iqpt
512 : : !FB double precision,intent(in) :: qpt_cart(3)
513 : : double precision,intent(in) :: qpt(3)
514 : :
515 : 9290 : double precision, allocatable :: omega (:)
516 : 9290 : double complex, allocatable :: dij (:,:)
517 : 9290 : double complex, allocatable :: eigenV(:,:)
518 : : double precision :: norm_eigenV
519 : : integer :: ii,jj,iatcell,jatcell
520 : :
521 [ + - - + : 85890 : ABI_MALLOC(omega ,(3*Invar%natom_unitcell)) ; omega(:) = zero
+ - + - +
- - + - +
+ + ]
522 [ - + + - : 776642 : ABI_MALLOC(dij ,(3*Invar%natom_unitcell,3*Invar%natom_unitcell)); dij(:,:) =czero
- + + - +
- + - - +
- + + + +
+ ]
523 [ - + + - : 767352 : ABI_MALLOC(eigenV,(3*Invar%natom_unitcell,3*Invar%natom_unitcell)); eigenV(:,:)=czero
- + + - -
+ - + + +
+ + ]
524 [ + + ]: 67310 : omega(:)=Eigen2nd%eigenval(:,iqpt)
525 [ + + ]: 28630 : do iatcell=1,Invar%natom_unitcell
526 [ + + ]: 104348 : do jatcell=1,Invar%natom_unitcell
527 [ + + ]: 322212 : do jj=1,3
528 : 227154 : norm_eigenV = 0.0d0
529 [ + + ]: 908616 : do ii=1,3
530 : : dij ((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%dynmat (1,ii,iatcell,jj,jatcell,iqpt),&
531 : 681462 : & Eigen2nd%dynmat (2,ii,iatcell,jj,jatcell,iqpt))
532 : : eigenV((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%eigenvec(1,ii,iatcell,jj,jatcell,iqpt),&
533 : 681462 : & Eigen2nd%eigenvec(2,ii,iatcell,jj,jatcell,iqpt))
534 : : norm_eigenV = norm_eigenV + real(eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj))**2 + &
535 : 908616 : aimag(eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj))**2
536 : : end do !ii
537 : 227154 : norm_eigenV = dsqrt(norm_eigenV)
538 : : ! if (norm_eigenV.lt.tol8) then
539 : : ! write(6,*) "iqpt, norm=",iqpt,norm_eigenV
540 : : ! ABI_ERROR(' STOP: THE NORM OF EIGENVECTOR IS ZERO')
541 : : ! end if
542 [ + + ]: 302872 : if (norm_eigenV.gt.tol8) then
543 [ + + ]: 901016 : do ii=1,3
544 : 901016 : eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj) = eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj) / dcmplx(norm_eigenV,0.d0)
545 : : end do !ii
546 : : end if
547 : : end do !jj
548 : : end do !jatcell
549 : : end do !iatcell
550 : :
551 : : ! Print the dynamical matrix (Dij)
552 : 9290 : write(52,'(a,1x,3(f10.6,1x))') 'For qpt=',qpt(:)
553 : 9290 : write(52,'(a,i4,a)') ' Dij(',iqpt,'real)='
554 [ + + ]: 28630 : do iatcell=1,Invar%natom_unitcell
555 [ + + ]: 246494 : write(52,'(100(f10.6,1x))') real(dij(1+(iatcell-1)*3,:))
556 [ + + ]: 246494 : write(52,'(100(f10.6,1x))') real(dij(2+(iatcell-1)*3,:))
557 [ + + ]: 255784 : write(52,'(100(f10.6,1x))') real(dij(3+(iatcell-1)*3,:))
558 : : end do
559 : 9290 : write(52,'(a,i4,a)') ' Dij(',iqpt,'imag)='
560 [ + + ]: 28630 : do iatcell=1,Invar%natom_unitcell
561 [ + + ]: 246494 : write(52,'(100(f10.6,1x))') aimag(dij(1+(iatcell-1)*3,:))
562 [ + + ]: 246494 : write(52,'(100(f10.6,1x))') aimag(dij(2+(iatcell-1)*3,:))
563 [ + + ]: 255784 : write(52,'(100(f10.6,1x))') aimag(dij(3+(iatcell-1)*3,:))
564 : : end do
565 : 9290 : write(52,*)' '
566 : :
567 : : ! Print the frequencies (omega)
568 [ + + + - : 18830 : if (Invar%enunit.eq.0) write(53,'(i5,1x,100(f15.3,1x))') iqpt,(omega(ii)*Ha_eV*1000,ii=1,3*Invar%natom_unitcell)
+ + ]
569 [ + + + - : 52205 : if (Invar%enunit.eq.1) write(53,'(i5,1x,100(f15.3,1x))') iqpt,(omega(ii)*Ha_cmm1 ,ii=1,3*Invar%natom_unitcell)
+ + ]
570 [ + + + - : 9866 : if (Invar%enunit.eq.2) write(53,'(i5,1x,100(f15.3,1x))') iqpt,(omega(ii)*1000 ,ii=1,3*Invar%natom_unitcell)
+ + ]
571 [ + + + - : 14279 : if (Invar%enunit.eq.3) write(53,'(i5,1x,100(f15.3,1x))') iqpt,(omega(ii)*Ha_THz ,ii=1,3*Invar%natom_unitcell)
+ + ]
572 : :
573 : : ! Print the eigenvectors (eigenV)
574 : : ! write(51,*) 'For iqpt=',iqpt
575 : : ! do ii=1,3*Invar%natom_unitcell
576 : : ! write(51,*) 'Mode number',ii,' energy',omega(ii)
577 : : ! write(51,*) ' Real:'
578 : : ! write(51,*) real(eigenV(:,ii))
579 : : ! write(51,*) ' Imag:'
580 : : ! write(51,*) aimag(eigenV(:,ii))
581 : : ! end do
582 : : ! write(51,*) ' '
583 : :
584 : 9290 : write(51,'(a,4x,i3,4(f15.6,1x))') 'q-pt=',iqpt,(qpt(ii),ii=1,3),qpt(1)*qpt(2)*qpt(3)
585 : :
586 [ + + ]: 67310 : do ii=1,3*Invar%natom_unitcell
587 : 67310 : write(51,'(i5,5x,f15.6)') ii,omega(ii)*Ha_cmm1
588 : : enddo
589 : 9290 : write(51,*) 'Phonon Eigenvectors'
590 : 9290 : write(51,'(a,17x,a,33x,a,34x,a)') 'Mode Ion','X','Y','Z'
591 [ + + ]: 67310 : do ii=1,3*Invar%natom_unitcell
592 [ + + ]: 294464 : do jj=1,Invar%natom_unitcell
593 : 227154 : write(51,'(i2,3x,i2,3x,6(f15.12,1x))') ii,jj,real(eigenV(3*(jj-1)+1,ii)),&
594 : 227154 : & aimag(eigenV(3*(jj-1)+1,ii)),&
595 : 227154 : & real(eigenV(3*(jj-1)+2,ii)),&
596 : 227154 : & aimag(eigenV(3*(jj-1)+2,ii)),&
597 : 227154 : & real(eigenV(3*(jj-1)+3,ii)),&
598 : 512328 : & aimag(eigenV(3*(jj-1)+3,ii))
599 : : enddo
600 : : enddo
601 : :
602 : :
603 [ - + ]: 9290 : ABI_FREE(omega)
604 [ - + ]: 9290 : ABI_FREE(dij)
605 [ - + ]: 9290 : ABI_FREE(eigenV)
606 : :
607 : 9290 : end subroutine tdep_write_dij
608 : : !!***
609 : :
610 : : !=====================================================================================================
611 : :
612 : : !!****f* ABINIT/m_tdep_phi2/tdep_build_phi2_33
613 : : !! NAME
614 : : !! tdep_build_phi2_33
615 : : !!
616 : : !! FUNCTION
617 : : !! Build the 3x3 second-order IFC of an interaction from the IFC of the
618 : : !! reference interaction of its shell, by applying the symmetry operation
619 : : !! isym: Phi2_33 = S Phi2_ref S^-1, and transposing the result if the
620 : : !! interaction is reversed with respect to the reference one (itrans=2).
621 : : !!
622 : : !! INPUTS
623 : : !! isym = Index of the symmetry operation.
624 : : !! Phi2_ref(3,3) = Second-order IFC of the reference interaction.
625 : : !! Sym = Symmetries object describing all the symmetry operations of the crystal.
626 : : !! itrans = 1 for a direct transformation, 2 if the pair of atoms is reversed.
627 : : !!
628 : : !! OUTPUT
629 : : !! Phi2_33(3,3) = Second-order IFC of the interaction.
630 : : !!
631 : : !! SOURCE
632 : :
633 : 320932 : subroutine tdep_build_phi2_33(isym,Phi2_ref,Phi2_33,Sym,itrans)
634 : :
635 : : type(Symmetries_type),intent(in) :: Sym
636 : : ! type(atdep_dataset_type),intent(in) :: Invar
637 : : double precision, intent(in) :: Phi2_ref(3,3)
638 : : double precision, intent(out) :: Phi2_33(3,3)
639 : : integer,intent(in) :: isym,itrans
640 : :
641 : : double precision :: Phi2_tmp(3,3),tmp1(3,3)
642 : :
643 : : ! Transform in the new basis wrt S_ref
644 : 320932 : call DGEMM('N','N',3,3,3,1.d0,Sym%S_ref(:,:,isym,1),3,Phi2_ref,3,0.d0,Phi2_tmp,3)
645 : 320932 : call DGEMM('N','N',3,3,3,1.d0,Phi2_tmp,3,Sym%S_inv(:,:,isym,1),3,0.d0,Phi2_33,3)
646 : :
647 [ - + ]: 320932 : if ((itrans.lt.1).or.(itrans.gt.2)) then
648 : 0 : ABI_BUG('This value of the symmetry index is not permitted')
649 : : end if
650 : : ! Transpose the 3x3 matrix if required
651 [ + + ]: 320932 : if (itrans.eq.2) then
652 : 48632 : tmp1(:,:)=Phi2_33(:,:)
653 : 48632 : Phi2_33(1,2)=tmp1(2,1)
654 : 48632 : Phi2_33(1,3)=tmp1(3,1)
655 : 48632 : Phi2_33(2,3)=tmp1(3,2)
656 : 48632 : Phi2_33(2,1)=tmp1(1,2)
657 : 48632 : Phi2_33(3,1)=tmp1(1,3)
658 : 48632 : Phi2_33(3,2)=tmp1(2,3)
659 : : end if
660 : :
661 : 320932 : end subroutine tdep_build_phi2_33
662 : : !!***
663 : :
664 : : !=====================================================================================================
665 : :
666 : : !!****f* ABINIT/m_tdep_phi2/tdep_init_eigen2nd
667 : : !! NAME
668 : : !! tdep_init_eigen2nd
669 : : !!
670 : : !! FUNCTION
671 : : !! Initialize the Eigen object: allocate (and set to zero) the arrays holding
672 : : !! the phonon frequencies, eigenvectors and dynamical matrices.
673 : : !!
674 : : !! INPUTS
675 : : !! natom_unitcell = Number of atoms in the unitcell.
676 : : !! nqpt = Number of q-points.
677 : : !!
678 : : !! OUTPUT
679 : : !!
680 : : !! SIDE EFFECTS
681 : : !! Eigen2nd = Eigen object, initialized.
682 : : !!
683 : : !! SOURCE
684 : :
685 : 92 : subroutine tdep_init_eigen2nd(Eigen2nd,natom_unitcell,nqpt)
686 : :
687 : : class(Eigen_type),intent(inout) :: Eigen2nd
688 : : integer, intent(in) :: natom_unitcell,nqpt
689 : :
690 [ + - - + : 203146 : ABI_MALLOC(Eigen2nd%eigenval,(3*natom_unitcell,nqpt)); Eigen2nd%eigenval(:,:) =zero
+ - - + +
- + - + -
- + - + -
+ + + +
+ ]
691 [ - + + - : 9573306 : ABI_MALLOC(Eigen2nd%eigenvec,(2,3,natom_unitcell,3,natom_unitcell,nqpt)); Eigen2nd%eigenvec(:,:,:,:,:,:)=zero
+ - - + +
- - + + -
+ - + - -
+ - + - +
+ + + + +
+ + + + +
+ + ]
692 [ - + + - : 9573122 : ABI_MALLOC(Eigen2nd%dynmat ,(2,3,natom_unitcell,3,natom_unitcell,nqpt)); Eigen2nd%dynmat (:,:,:,:,:,:)=zero
- + + - -
+ + - - +
- + - + +
+ + + + +
+ + + + +
+ ]
693 : :
694 : 92 : end subroutine tdep_init_eigen2nd
695 : : !!***
696 : :
697 : : !=====================================================================================================
698 : :
699 : : !!****f* ABINIT/m_tdep_phi2/tdep_destroy_eigen2nd
700 : : !! NAME
701 : : !! tdep_destroy_eigen2nd
702 : : !!
703 : : !! FUNCTION
704 : : !! Free the memory allocated in the Eigen object.
705 : : !!
706 : : !! SIDE EFFECTS
707 : : !! Eigen2nd = Eigen object. All its allocated arrays are deallocated.
708 : : !!
709 : : !! SOURCE
710 : :
711 : 92 : subroutine tdep_destroy_eigen2nd(Eigen2nd)
712 : :
713 : : class(Eigen_type),intent(inout) :: Eigen2nd
714 : :
715 [ - + ]: 92 : ABI_FREE(Eigen2nd%eigenval)
716 [ - + ]: 92 : ABI_FREE(Eigen2nd%eigenvec)
717 [ - + ]: 92 : ABI_FREE(Eigen2nd%dynmat)
718 : :
719 : 92 : end subroutine tdep_destroy_eigen2nd
720 : : !!***
721 : :
722 : : !=====================================================================================================
723 : :
724 : : !!****f* ABINIT/m_tdep_phi2/tdep_write_yaml
725 : : !! NAME
726 : : !! tdep_write_yaml
727 : : !!
728 : : !! FUNCTION
729 : : !! Write the phonon band structure along the q-point path (q-points,
730 : : !! distances, labels of the special q-points, frequencies in THz and
731 : : !! eigenvectors) in the _phonon-bands.yaml file, in the phonopy format.
732 : : !!
733 : : !! INPUTS
734 : : !! Eigen2nd = Eigen object containing the phonon frequencies and eigenvectors
735 : : !! along the q-point path.
736 : : !! Qpt = Qpoints object containing the q-point path.
737 : : !! Prefix = Root name of the output file.
738 : : !!
739 : : !! OUTPUT
740 : : !! Only writing.
741 : : !!
742 : : !! SOURCE
743 : :
744 : 46 : subroutine tdep_write_yaml(Eigen2nd,Qpt,Prefix)
745 : :
746 : : type(Eigen_type),intent(in) :: Eigen2nd
747 : : type(Qpoints_type),intent(in) :: Qpt
748 : : character(len=*) :: Prefix
749 : :
750 : : integer :: ii,jj,iatcell,jatcell,iqpt,imode,nmode
751 : : double precision :: distance
752 : 46 : double complex, allocatable :: eigenV(:,:)
753 : :
754 : 46 : nmode=size(Eigen2nd%eigenval,dim=1)
755 : 46 : open(unit=52,file=trim(Prefix)//'_phonon-bands.yaml')
756 : 46 : write(52,'(a,i4)') 'nqpoint:',Qpt%nqpt
757 : 46 : write(52,'(a,i4)') 'npath:',Qpt%qpt_tot-1
758 : 46 : write(52,'(a)') 'segment_nqpoint:'
759 [ + + ]: 285 : do ii=1,Qpt%qpt_tot-1
760 : 285 : write(52,'(a,i4)') '- ',Qpt%lgth_segments(ii)
761 : : end do
762 : 46 : write(52,'(a,i4)') 'natom:',nmode/3
763 : 46 : write(52,'(a)') 'phonon:'
764 : 46 : distance=0.d0
765 [ + + ]: 9336 : do iqpt=1,Qpt%nqpt
766 : 9290 : write(52,'(a,3(f15.6,1x,a))') '- q-position: [',Qpt%qpt_red(1,iqpt),',',Qpt%qpt_red(2,iqpt),',',Qpt%qpt_red(3,iqpt),']'
767 [ + + + + ]: 37022 : if (iqpt.gt.1) distance=distance+sqrt(sum((Qpt%qpt_cart(:,iqpt)-Qpt%qpt_cart(:,iqpt-1))**2))
768 : 9290 : write(52,'(a,f15.6)') ' distance:',distance
769 [ + + ]: 71942 : do ii=1,Qpt%qpt_tot
770 [ + + + + ]: 260753 : if (sum(abs(Qpt%qpt_red(:,iqpt)-Qpt%special_red(ii,:))).lt.tol8) then
771 [ + - ]: 285 : write(52,'(3a)') " label: '",trim(Qpt%special_qpt(ii)),"'"
772 : 285 : exit
773 : : end if
774 : : end do !ii
775 [ - + + - : 776642 : ABI_MALLOC(eigenV,(nmode,nmode)); eigenV(:,:)=czero
- + + - +
- + - - +
- + + + +
+ ]
776 [ + + ]: 28630 : do iatcell=1,nmode/3
777 [ + + ]: 104348 : do jatcell=1,nmode/3
778 [ + + ]: 322212 : do ii=1,3
779 [ + + ]: 984334 : do jj=1,3
780 : : eigenV((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%eigenvec(1,ii,iatcell,jj,jatcell,iqpt),&
781 : 908616 : & Eigen2nd%eigenvec(2,ii,iatcell,jj,jatcell,iqpt))
782 : : end do !ii
783 : : end do !jj
784 : : end do !jatcell
785 : : end do !iatcell
786 : 9290 : write(52,'(a)') ' band:'
787 [ + + ]: 67310 : do imode=1,nmode
788 : 58020 : write(52,'(a,i4)') ' - #',imode
789 : 58020 : write(52,'(a,f15.6)') ' frequency:',Eigen2nd%eigenval(imode,iqpt)*Ha_THz
790 : 58020 : write(52,'(a)') ' eigenvector:'
791 [ + + ]: 294464 : do iatcell=1,nmode/3
792 : 227154 : write(52,'(a,i4)') " - # atom ", iatcell
793 [ + + ]: 966636 : do ii=1,3
794 : 681462 : write(52,'(a,f18.9,a,f18.9,a)') " - [",real(eigenV((iatcell-1)*3+ii,imode)),','&
795 : 1590078 : & ,aimag(eigenV((iatcell-1)*3+ii,imode)),']'
796 : : end do
797 : : end do
798 : : end do !nmode
799 [ - + ]: 9290 : ABI_FREE(eigenV)
800 : 9336 : write(52,'(a)') ''
801 : : end do
802 : :
803 : 46 : close(52)
804 : :
805 : :
806 : 46 : end subroutine tdep_write_yaml
807 : : !!***
808 : :
809 : : !=====================================================================================================
810 : :
811 [ # # # # : 0 : end module m_tdep_phi2
# # # # #
# # # # #
# # # # ]
812 : : !!***
|
