Line data Source code
1 :
2 : #if defined HAVE_CONFIG_H
3 : #include "config.h"
4 : #endif
5 :
6 : #include "abi_common.h"
7 :
8 : module m_tdep_phi2
9 :
10 : use defs_basis
11 : use m_errors
12 : use m_abicore
13 : use m_xmpi
14 : use m_io_tools
15 : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
16 : use m_tdep_shell, only : Shell_type
17 : use m_tdep_sym, only : Symmetries_type
18 : use m_tdep_qpt, only : Qpoints_type
19 : use m_tdep_sampling, only : tdep_Sampling_type
20 : use m_tdep_solver, only : tdep_Solver_type
21 : use m_tdep_model, only : tdep_Model_type
22 :
23 : implicit none
24 :
25 : type Eigen_type
26 :
27 : double precision, allocatable :: eigenval(:,:)
28 : double precision, allocatable :: eigenvec(:,:,:,:,:,:)
29 : double precision, allocatable :: dynmat(:,:,:,:,:,:)
30 :
31 : end type Eigen_type
32 :
33 : public :: tdep_calc_ftot2
34 : public :: tdep_calc_phi1
35 : public :: tdep_write_phi1
36 : public :: tdep_calc_phi2
37 : public :: tdep_write_phi2
38 : public :: tdep_build_phi2_33
39 : public :: tdep_calc_dij
40 : public :: tdep_write_dij
41 : public :: tdep_init_eigen2nd
42 : public :: tdep_destroy_eigen2nd
43 : public :: tdep_write_yaml
44 :
45 : contains
46 :
47 : !====================================================================================================
48 49 : subroutine tdep_calc_ftot2(Model,Invar,Phi2,ucart)
49 :
50 : type(tdep_Model_type),intent(inout) :: Model
51 : type(atdep_dataset_type),intent(in) :: Invar
52 : double precision, intent(in) :: Phi2(3*Invar%natom,3*Invar%natom)
53 : double precision, intent(in) :: ucart(3,Invar%natom,Invar%my_nstep)
54 :
55 : integer :: jj,istep,jatom
56 49 : double precision, allocatable :: ucart_blas(:)
57 49 : double precision, allocatable :: ftot2(:)
58 :
59 : ! Compute Forces of the model TDEP
60 23003 : ABI_CALLOC(ucart_blas,(3*Invar%natom))
61 22905 : ABI_CALLOC(ftot2 ,(3*Invar%natom))
62 1101 : do istep=1,Invar%my_nstep
63 460748 : ucart_blas(:)=0.d0
64 460748 : ftot2 (:)=0.d0
65 154284 : do jatom=1,Invar%natom
66 613980 : do jj=1,3
67 612928 : ucart_blas(3*(jatom-1)+jj)=ucart(jj,jatom,istep)
68 : end do
69 : end do
70 460748 : Model%Phi1Ui(istep)=sum(Model%Phi1(:)*ucart_blas(:))
71 1052 : call DGEMM('N','N',3*Invar%natom,1,3*Invar%natom,1.d0,Phi2,3*Invar%natom,ucart_blas,3*Invar%natom,0.d0,ftot2(:),3*Invar%natom)
72 1052 : call DGEMM('T','N',1,1,3*Invar%natom,1./2.d0,ftot2,3*Invar%natom,ucart_blas,3*Invar%natom,0.d0,Model%Phi2UiUj(istep),3*Invar%natom)
73 460797 : Model%Forces(3*Invar%natom*(istep-1)+1:3*Invar%natom*istep)=-Model%Phi1(:)-ftot2(:)
74 : end do !istep
75 49 : ABI_FREE(ucart_blas)
76 49 : ABI_FREE(ftot2)
77 :
78 49 : end subroutine tdep_calc_ftot2
79 :
80 : !=====================================================================================================
81 43 : subroutine tdep_calc_phi1(Solver,Shell1at,Sym,Phi1)
82 :
83 : type(tdep_Solver_type),intent(in) :: Solver
84 : type(Shell_type),intent(in) :: Shell1at
85 : type(Symmetries_type),intent(in) :: Sym
86 : double precision,intent(out) :: Phi1(3*Solver%natom)
87 :
88 : integer :: ishell,isym,iatom,ncoeff,ncoeff_prev
89 : integer :: nshell,ii,iatshell,iat_mod
90 43 : double precision,allocatable :: Phi1_coeff(:)
91 43 : double precision,allocatable :: Phi1_3(:),Phi1_ref(:,:)
92 :
93 43 : nshell=Shell1at%nshell
94 129 : ABI_MALLOC(Phi1_coeff,(Solver%ncoeff1st))
95 43 : Phi1_coeff(:) = Solver%theta(1:Solver%ncoeff1st)
96 377 : ABI_MALLOC(Phi1_ref,(3,nshell)); Phi1_ref(:,:)=zero
97 215 : ABI_MALLOC(Phi1_3,(3)) ; Phi1_3(:)=0.d0
98 105 : do ishell=1,nshell
99 : ! Build the 3x3 IFC per shell
100 62 : ncoeff =Shell1at%ncoeff(ishell)
101 62 : ncoeff_prev=Shell1at%ncoeff_prev(ishell)
102 248 : do ii=1,3
103 248 : Phi1_ref(ii,ishell)=sum(Shell1at%proj(ii,1:ncoeff,ishell)*Phi1_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
104 : end do
105 : ! Build the vector-IFC of an atom in this shell
106 62 : if (Shell1at%neighbours(1,ishell)%n_interactions.eq.0) cycle
107 43 : do iatshell=1,Shell1at%neighbours(1,ishell)%n_interactions
108 0 : iatom=Shell1at%neighbours(1,ishell)%atomj_in_shell(iatshell)
109 0 : isym =Shell1at%neighbours(1,ishell)%sym_in_shell(iatshell)
110 0 : do ii=1,3
111 0 : Phi1_3(ii)=sum(Sym%S_ref(ii,:,isym,1)*Phi1_ref(:,ishell))
112 : end do
113 62 : Phi1((iatom-1)*3+1:(iatom-1)*3+3)=Phi1_3(:)
114 : end do !iatshell
115 : end do !ishell
116 : ! Acoustic sum rule
117 172 : do ii=1,3
118 19282 : do iatom=1,Solver%natom
119 19110 : iat_mod=mod(iatom+Solver%natom_unitcell-1,Solver%natom_unitcell)+1
120 19110 : if (iat_mod==1) cycle
121 19239 : Phi1((iatom-iat_mod+1)*3+ii)=Phi1((iatom-iat_mod+1)*3+ii)-Phi1((iatom-1)*3+ii)
122 : end do
123 : end do
124 43 : ABI_FREE(Phi1_coeff)
125 43 : ABI_FREE(Phi1_3)
126 43 : ABI_FREE(Phi1_ref)
127 :
128 : ! Remove the rounding errors before writing (for non regression testing purposes)
129 19153 : do ii=1,3*Solver%natom
130 19153 : if (abs(Phi1(ii)).lt.tol8) Phi1(ii)=zero
131 : end do
132 :
133 43 : end subroutine tdep_calc_phi1
134 :
135 : !=====================================================================================================
136 44 : subroutine tdep_write_phi1(Invar,Phi1)
137 :
138 : type(atdep_dataset_type),intent(in) :: Invar
139 : double precision,intent(in) :: Phi1(3*Invar%natom)
140 :
141 : integer :: iatcell,ii
142 :
143 44 : write(Invar%stdout,*) ' '
144 44 : write(Invar%stdout,*) '#############################################################################'
145 44 : write(Invar%stdout,*) '#### For each shell, list of coefficients (IFC), number of neighbours... ####'
146 44 : write(Invar%stdout,*) '#############################################################################'
147 :
148 : ! Write the IFCs in the data.out file (with others specifications:
149 : ! number of atoms in a shell, distance, Trace...)
150 150 : do iatcell=1,Invar%natom_unitcell
151 106 : write(Invar%stdout,'(a,i4)') ' ############# List of (first order) IFC for the reference atom=',iatcell
152 424 : write(Invar%stdout,'(2x,3(f9.6,1x))') (Phi1((iatcell-1)*3+ii),ii=1,3)
153 150 : write(Invar%stdout,*) ' '
154 : end do !iatcell
155 :
156 44 : end subroutine tdep_write_phi1
157 :
158 : !=====================================================================================================
159 43 : subroutine tdep_calc_phi2(Solver,Shell2at,Sym,Phi2)
160 :
161 : type(tdep_Solver_type),intent(in) :: Solver
162 : type(Shell_type),intent(in) :: Shell2at
163 : type(Symmetries_type),intent(in) :: Sym
164 : double precision,intent(out) :: Phi2(3*Solver%natom,3*Solver%natom)
165 :
166 : integer :: ishell,isym,eatom,fatom,ncoeff,ncoeff_prev
167 : integer :: nshell,ii,jj,kk,kappa,iatshell,itrans
168 43 : double precision,allocatable :: Phi2_33(:,:),Phi2_ref(:,:,:)
169 43 : double precision,allocatable :: Phi2_coeff(:)
170 :
171 43 : nshell=Shell2at%nshell
172 129 : ABI_MALLOC(Phi2_coeff,(Solver%ncoeff2nd))
173 1968 : Phi2_coeff(:) = Solver%theta(Solver%ncoeff1st+1:Solver%ncoeff1st+Solver%ncoeff2nd)
174 6356 : ABI_MALLOC(Phi2_ref,(3,3,nshell)); Phi2_ref(:,:,:)=zero
175 602 : ABI_MALLOC(Phi2_33,(3,3)) ; Phi2_33(:,:)=0.d0
176 522 : do ishell=1,nshell
177 : ! Build the 3x3 IFC per shell
178 479 : ncoeff =Shell2at%ncoeff(ishell)
179 479 : ncoeff_prev=Shell2at%ncoeff_prev(ishell)
180 479 : kappa=0
181 1916 : do ii=1,3
182 6227 : do jj=1,3
183 4311 : kappa=kappa+1
184 23073 : Phi2_ref(ii,jj,ishell)=sum(Shell2at%proj(kappa,1:ncoeff,ishell)*Phi2_coeff(ncoeff_prev+1:ncoeff_prev+ncoeff))
185 : end do
186 : end do
187 68344 : do eatom=1,Solver%natom
188 : ! Build the 3x3 IFC of an atom in this shell
189 67822 : if (Shell2at%neighbours(eatom,ishell)%n_interactions.eq.0) cycle
190 462411 : do iatshell=1,Shell2at%neighbours(eatom,ishell)%n_interactions
191 405630 : fatom=Shell2at%neighbours(eatom,ishell)%atomj_in_shell(iatshell)
192 405630 : isym =Shell2at%neighbours(eatom,ishell)%sym_in_shell(iatshell)
193 405630 : itrans=Shell2at%neighbours(eatom,ishell)%transpose_in_shell(iatshell)
194 405630 : if (fatom.lt.eatom) cycle
195 206000 : call tdep_build_phi2_33(isym,Phi2_ref(:,:,ishell),Phi2_33,Sym,itrans)
196 : ! Symetrization of the Phi2 matrix
197 2678000 : Phi2((eatom-1)*3+1:(eatom-1)*3+3,3*(fatom-1)+1:3*(fatom-1)+3)=Phi2_33(:,:)
198 891822 : do ii=1,3
199 2877630 : do jj=1,3
200 2472000 : Phi2((fatom -1)*3+ii,3*(eatom-1)+jj)=Phi2_33(jj,ii)
201 : end do
202 : end do
203 : end do !iatshell
204 : end do !eatom
205 : end do !ishell
206 : ! Acoustic sum rule
207 6413 : do eatom=1,Solver%natom
208 25523 : do jj=1,3
209 82810 : do kk=1,3
210 10967916 : do fatom=1,Solver%natom
211 10891476 : if (fatom==eatom) cycle
212 : Phi2((eatom-1)*3+jj,(eatom-1)*3+kk)=Phi2((eatom-1)*3+jj,3*(eatom-1)+kk)&
213 10948806 : & -Phi2((eatom-1)*3+jj,3*(fatom-1)+kk)
214 : enddo
215 : enddo
216 : enddo
217 : enddo
218 43 : ABI_FREE(Phi2_coeff)
219 43 : ABI_FREE(Phi2_33)
220 43 : ABI_FREE(Phi2_ref)
221 :
222 : ! Remove the rounding errors before writing (for non regression testing purposes)
223 19153 : do ii=1,3*Solver%natom
224 10910629 : do jj=1,3*Solver%natom
225 10910586 : if (abs(Phi2(ii,jj)).lt.tol8) Phi2(ii,jj)=zero
226 : end do
227 : end do
228 :
229 43 : end subroutine tdep_calc_phi2
230 :
231 : !=====================================================================================================
232 44 : subroutine tdep_write_phi2(distance,Invar,MPIdata,Phi2,Shell2at)
233 :
234 : type(atdep_dataset_type),intent(in) :: Invar
235 : type(Shell_type),intent(in) :: Shell2at
236 : type(MPI_enreg_type), intent(in) :: MPIdata
237 : double precision,intent(in) :: distance(Invar%natom,Invar%natom,4)
238 : double precision,intent(in) :: Phi2(3*Invar%natom,3*Invar%natom)
239 :
240 : integer :: iatcell,ishell,jshell,jatom
241 : integer :: nshell,ii,this_shell,iatshell
242 : double precision :: max_bound,min_bound,dist_shell,tmp1,tmp2,tmp3
243 44 : integer,allocatable :: tab_shell(:)
244 :
245 44 : write(Invar%stdout,*) ' '
246 44 : write(Invar%stdout,*) '#############################################################################'
247 44 : write(Invar%stdout,*) '#### For each shell, list of coefficients (IFC), number of neighbours... ####'
248 44 : write(Invar%stdout,*) '#############################################################################'
249 :
250 44 : nshell=Shell2at%nshell
251 : ! Write the IFCs in the data.out file (with others specifications:
252 : ! number of atoms in a shell, distance, Trace...)
253 632 : ABI_MALLOC(tab_shell,(nshell)); tab_shell(:)=0
254 150 : do iatcell=1,Invar%natom_unitcell
255 2392 : tab_shell(:)=0
256 106 : write(Invar%stdout,'(a,i4)') ' ############# List of (second order) IFC for the reference atom=',iatcell
257 : ! Sort the IFC with distance in increasing order
258 106 : min_bound=-1.d0
259 2436 : do ishell=1,nshell
260 113058 : do jshell=1,nshell
261 113058 : if ((distance(Shell2at%iatref(jshell),Shell2at%jatref(jshell),1).ge.min_bound).and.(tab_shell(jshell).eq.0)) then
262 56529 : max_bound=distance(Shell2at%iatref(jshell),Shell2at%jatref(jshell),1)
263 56529 : this_shell=jshell
264 : end if
265 : end do
266 :
267 113058 : do jshell=1,nshell
268 : if ((distance(Shell2at%iatref(jshell),Shell2at%jatref(jshell),1).lt.max_bound).and.&
269 : & (distance(Shell2at%iatref(jshell),Shell2at%jatref(jshell),1).ge.min_bound).and.&
270 : & (distance(Shell2at%iatref(jshell),Shell2at%jatref(jshell),1).ne.&
271 110772 : & distance(Shell2at%iatref(this_shell),Shell2at%jatref(this_shell),1)).and.&
272 2286 : & (tab_shell(jshell).eq.0)) then
273 4843 : max_bound=distance(Shell2at%iatref(jshell),Shell2at%jatref(jshell),1)
274 4843 : this_shell=jshell
275 : end if
276 : end do
277 2286 : tab_shell(this_shell)=1
278 2286 : min_bound=max_bound
279 2286 : dist_shell=distance(Shell2at%iatref(this_shell),Shell2at%jatref(this_shell),1)
280 :
281 : ! Write the IFC properly
282 2392 : if (Shell2at%neighbours(iatcell,this_shell)%n_interactions.ne.0) then
283 1099 : write(Invar%stdout,'(a,i4,a,i4,a,f9.6)') ' ======== NEW SHELL (ishell=',this_shell,&
284 2198 : & '): There are',Shell2at%neighbours(iatcell,this_shell)%n_interactions,' atoms on this shell at distance=',dist_shell
285 6367 : do iatshell=1,Shell2at%neighbours(iatcell,this_shell)%n_interactions
286 5268 : jatom=Shell2at%neighbours(iatcell,this_shell)%atomj_in_shell(iatshell)
287 5268 : write(Invar%stdout,'(a,i4,a,i4)') ' For jatom=',jatom,' ,with type=',mod(jatom-1,Invar%natom_unitcell)+1
288 21072 : do ii=1,3
289 15804 : if (abs(Phi2((iatcell-1)*3+ii,(jatom-1)*3+1)).lt.5.d-7) then
290 2660 : tmp1=0.d0
291 : else
292 13144 : tmp1=Phi2((iatcell-1)*3+ii,(jatom-1)*3+1)
293 : end if
294 15804 : if (abs(Phi2((iatcell-1)*3+ii,(jatom-1)*3+2)).lt.5.d-7) then
295 2622 : tmp2=0.d0
296 : else
297 13182 : tmp2=Phi2((iatcell-1)*3+ii,(jatom-1)*3+2)
298 : end if
299 15804 : if (abs(Phi2((iatcell-1)*3+ii,(jatom-1)*3+3)).lt.5.d-7) then
300 2864 : tmp3=0.d0
301 : else
302 12940 : tmp3=Phi2((iatcell-1)*3+ii,(jatom-1)*3+3)
303 : end if
304 21072 : write(Invar%stdout,'(2x,3(f9.6,1x))') tmp1,tmp2,tmp3
305 : end do
306 5268 : write(Invar%stdout,'(a,3(1x,f11.6))') ' The components of the vector are:', distance(iatcell,jatom,2:4)
307 5268 : write(Invar%stdout,'(a,(1x,f9.6))') ' Trace=',Phi2((iatcell-1)*3+1,(jatom-1)*3+1)+Phi2((iatcell-1)*3+2,&
308 10536 : & (jatom-1)*3+2)+Phi2((iatcell-1)*3+3,(jatom-1)*3+3)
309 6367 : write(Invar%stdout,*) ' '
310 : end do
311 : end if
312 : end do !ishell
313 : end do !iatcell
314 44 : ABI_FREE(tab_shell)
315 :
316 : ! Write the Phi2_unitcell.dat and Phi2.dat files
317 44 : if (Invar%debug.and.MPIdata%iam_master) then
318 2 : write(Invar%stdout,'(a)') ' See the Phi2*.dat file'
319 2 : open(unit=52,file=trim(Invar%output_prefix)//'_Phi2_unitcell.dat')
320 2 : open(unit=55,file=trim(Invar%output_prefix)//'_Phi2.dat')
321 674 : do jatom=1,3*Invar%natom
322 672 : if (jatom.le.3*Invar%natom_unitcell) then
323 30 : write(52,'(10000(f10.6,1x))') Phi2(jatom,:)
324 : end if
325 674 : write(55,'(10000(f10.6,1x))') Phi2(jatom,:)
326 : end do
327 2 : close(52)
328 2 : close(55)
329 : end if
330 :
331 44 : end subroutine tdep_write_phi2
332 :
333 : !=====================================================================================================
334 0 : subroutine tdep_calc_dij(dij,eigenV,iqpt,Invar,omega,Phi2,qpt_cart,Rlatt_cart)
335 :
336 : type(atdep_dataset_type),intent(in) :: Invar
337 : integer,intent(in) :: iqpt
338 : double precision,intent(in) :: Phi2(3*Invar%natom,3*Invar%natom)
339 : double precision,intent(in) :: Rlatt_cart(3,Invar%natom_unitcell,Invar%natom)
340 : double precision,intent(in) :: qpt_cart(3)
341 : double precision,intent(out) :: omega (3*Invar%natom_unitcell)
342 : double complex ,intent(out) :: dij (3*Invar%natom_unitcell,3*Invar%natom_unitcell)
343 : double complex ,intent(out) :: eigenV(3*Invar%natom_unitcell,3*Invar%natom_unitcell)
344 :
345 : integer :: LWORK,ii,jj,kk,iatom,jatom,iatcell,jatcell,itypat,jtypat,iat_mod,INFO,itemp,imode,nmode
346 : double precision :: phase
347 : double complex :: norm
348 0 : double precision, allocatable :: RWORK(:)
349 0 : double complex, allocatable :: WORKC(:)
350 : ! double complex, allocatable :: mass_mat(:,:)
351 :
352 : ! Calculation of the dynamical matrix (Dij)
353 0 : do iatcell=1,Invar%natom_unitcell
354 0 : do jatom=1,Invar%natom
355 0 : iat_mod=mod(jatom+Invar%natom_unitcell-1,Invar%natom_unitcell)+1
356 0 : phase=0.d0
357 0 : do kk=1,3
358 0 : phase=phase+2*pi*Rlatt_cart(kk,iatcell,jatom)*qpt_cart(kk)
359 : end do
360 0 : do ii=1+(iatcell-1)*3,3+(iatcell-1)*3
361 0 : do jj=1,3
362 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))
363 : end do !jj
364 : end do !ii
365 : end do !jatom
366 : end do !iatcell
367 :
368 : ! The Dij has to be an hermitian matrix
369 0 : itemp=0
370 0 : do ii=1,3*Invar%natom_unitcell
371 0 : do jj=ii,3*Invar%natom_unitcell
372 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
373 0 : if (Invar%debug) then
374 0 : write (Invar%stdout,'(a,1x,2(i4,1x))') 'for ii,jj=',ii,jj
375 0 : write (Invar%stdout,'(a,1x,1(f12.8,1x))') 'abs(realij-realji)=',abs(real(dij(ii,jj))-real(dij(jj,ii)))
376 0 : write (Invar%stdout,'(a,1x,1(f12.8,1x))') 'abs(imagij+imagji)=',abs(aimag(dij(ii,jj))+aimag(dij(jj,ii)))
377 : end if
378 0 : itemp=itemp+1
379 : end if
380 : end do
381 : end do
382 0 : if (itemp.ne.0.and.iqpt.eq.1) then
383 0 : write(Invar%stdout,*) 'WARNING: The Dij matrix is not hermitian'
384 0 : write(Invar%stdout,*) ' Probably: one shell may not have the whole number of atoms'
385 0 : write(Invar%stdout,*) ' The Dij matrix is symetrized'
386 : end if
387 :
388 : ! Diagonalization of dynamical matrix Dij/sqrt(Mi*Mj)
389 0 : LWORK=2*3*Invar%natom_unitcell-1
390 0 : ABI_MALLOC(WORKC,(LWORK)); WORKC(:)=czero
391 0 : ABI_MALLOC(RWORK,(3*3*Invar%natom_unitcell-2)); RWORK(:)=zero
392 0 : do iatcell=1,Invar%natom_unitcell
393 0 : itypat=Invar%typat_unitcell(iatcell)
394 0 : do jatcell=1,Invar%natom_unitcell
395 0 : jtypat=Invar%typat_unitcell(jatcell)
396 0 : do ii=1,3
397 0 : do jj=1,3
398 : eigenV(ii+(iatcell-1)*3,jj+(jatcell-1)*3)=dij(ii+(iatcell-1)*3,jj+(jatcell-1)*3)/&
399 0 : & dcmplx(dsqrt(Invar%amu(itypat)*Invar%amu(jtypat))*amu_emass,0.d0)
400 : end do !jj
401 : end do !ii
402 : end do !jatcell
403 : end do !iatcell
404 0 : call ZHEEV('V','U',3*Invar%natom_unitcell,eigenV(:,:),3*Invar%natom_unitcell,omega(:),WORKC,LWORK,RWORK,INFO)
405 :
406 : ! Normalization of the eigenvectors
407 0 : nmode=3*Invar%natom_unitcell
408 0 : do imode=1,nmode
409 : norm=zero
410 0 : do iatom=1,Invar%natom_unitcell
411 0 : do ii=1,3
412 0 : norm=norm+eigenV(3*(iatom-1)+ii,imode)*conjg(eigenV(3*(iatom-1)+ii,imode))
413 : end do
414 : end do
415 0 : eigenV(:,imode)=eigenV(:,imode)/dsqrt(real(norm))
416 : end do
417 :
418 : ! Remove the squared-negative frequencies
419 0 : do ii=1,Invar%natom_unitcell
420 0 : do jj=1,3
421 0 : if (omega((ii-1)*3+jj).lt.0.d0) then
422 0 : omega((ii-1)*3+jj)=-dsqrt(-omega((ii-1)*3+jj))
423 : else
424 0 : omega((ii-1)*3+jj)=dsqrt(omega((ii-1)*3+jj))
425 : end if
426 : end do
427 : end do
428 0 : ABI_FREE(WORKC)
429 0 : ABI_FREE(RWORK)
430 :
431 0 : end subroutine tdep_calc_dij
432 :
433 : !=====================================================================================================
434 : !FB subroutine tdep_write_dij(Eigen2nd,iqpt,Invar,qpt_cart)
435 22588 : subroutine tdep_write_dij(Eigen2nd,iqpt,Invar,qpt)
436 :
437 : type(atdep_dataset_type),intent(in) :: Invar
438 : integer,intent(in) :: iqpt
439 : !FB double precision,intent(in) :: qpt_cart(3)
440 : double precision,intent(in) :: qpt(3)
441 : type(Eigen_type),intent(in) :: Eigen2nd
442 :
443 22588 : double precision, allocatable :: omega (:)
444 22588 : double complex, allocatable :: dij (:,:)
445 22588 : double complex, allocatable :: eigenV(:,:)
446 : double precision :: norm_eigenV
447 : integer :: ii,jj,iatcell,jatcell
448 :
449 221145 : ABI_MALLOC(omega ,(3*Invar%natom_unitcell)) ; omega(:) = zero
450 2053372 : ABI_MALLOC(dij ,(3*Invar%natom_unitcell,3*Invar%natom_unitcell)); dij(:,:) =czero
451 2030784 : ABI_MALLOC(eigenV,(3*Invar%natom_unitcell,3*Invar%natom_unitcell)); eigenV(:,:)=czero
452 175969 : omega(:)=Eigen2nd%eigenval(:,iqpt)
453 73715 : do iatcell=1,Invar%natom_unitcell
454 274786 : do jatcell=1,Invar%natom_unitcell
455 855411 : do jj=1,3
456 603213 : norm_eigenV = 0.0d0
457 2412852 : do ii=1,3
458 : dij ((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%dynmat (1,ii,iatcell,jj,jatcell,iqpt),&
459 1809639 : & Eigen2nd%dynmat (2,ii,iatcell,jj,jatcell,iqpt))
460 : eigenV((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%eigenvec(1,ii,iatcell,jj,jatcell,iqpt),&
461 1809639 : & Eigen2nd%eigenvec(2,ii,iatcell,jj,jatcell,iqpt))
462 : norm_eigenV = norm_eigenV + real(eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj))**2 + &
463 2412852 : aimag(eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj))**2
464 : end do !ii
465 603213 : norm_eigenV = dsqrt(norm_eigenV)
466 : ! if (norm_eigenV.lt.tol8) then
467 : ! write(6,*) "iqpt, norm=",iqpt,norm_eigenV
468 : ! ABI_ERROR(' STOP: THE NORM OF EIGENVECTOR IS ZERO')
469 : ! end if
470 804284 : if (norm_eigenV.gt.tol8) then
471 2390116 : do ii=1,3
472 2390116 : eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj) = eigenV((iatcell-1)*3+ii, (jatcell-1)*3+jj) / dcmplx(norm_eigenV,0.d0)
473 : end do !ii
474 : end if
475 : end do !jj
476 : end do !jatcell
477 : end do !iatcell
478 :
479 : ! Print the dynamical matrix (Dij)
480 22588 : write(52,'(a,1x,3(f10.6,1x))') 'For qpt=',qpt(:)
481 22588 : write(52,'(a,i4,a)') ' Dij(',iqpt,'real)='
482 73715 : do iatcell=1,Invar%natom_unitcell
483 654340 : write(52,'(100(f10.6,1x))') real(dij(1+(iatcell-1)*3,:))
484 654340 : write(52,'(100(f10.6,1x))') real(dij(2+(iatcell-1)*3,:))
485 676928 : write(52,'(100(f10.6,1x))') real(dij(3+(iatcell-1)*3,:))
486 : end do
487 22588 : write(52,'(a,i4,a)') ' Dij(',iqpt,'imag)='
488 73715 : do iatcell=1,Invar%natom_unitcell
489 654340 : write(52,'(100(f10.6,1x))') aimag(dij(1+(iatcell-1)*3,:))
490 654340 : write(52,'(100(f10.6,1x))') aimag(dij(2+(iatcell-1)*3,:))
491 676928 : write(52,'(100(f10.6,1x))') aimag(dij(3+(iatcell-1)*3,:))
492 : end do
493 22588 : write(52,*)' '
494 :
495 : ! Print the frequencies (omega)
496 58192 : 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)
497 124216 : if (Invar%enunit.eq.1) write(53,'(i5,1x,100(f15.3,1x))') iqpt,(omega(ii)*Ha_cmm1 ,ii=1,3*Invar%natom_unitcell)
498 24022 : if (Invar%enunit.eq.2) write(53,'(i5,1x,100(f15.3,1x))') iqpt,(omega(ii)*1000 ,ii=1,3*Invar%natom_unitcell)
499 37303 : if (Invar%enunit.eq.3) write(53,'(i5,1x,100(f15.3,1x))') iqpt,(omega(ii)*Ha_THz ,ii=1,3*Invar%natom_unitcell)
500 :
501 : ! Print the eigenvectors (eigenV)
502 : ! write(51,*) 'For iqpt=',iqpt
503 : ! do ii=1,3*Invar%natom_unitcell
504 : ! write(51,*) 'Mode number',ii,' energy',omega(ii)
505 : ! write(51,*) ' Real:'
506 : ! write(51,*) real(eigenV(:,ii))
507 : ! write(51,*) ' Imag:'
508 : ! write(51,*) aimag(eigenV(:,ii))
509 : ! end do
510 : ! write(51,*) ' '
511 :
512 22588 : write(51,'(a,4x,i3,4(f15.6,1x))') 'q-pt=',iqpt,(qpt(ii),ii=1,3),qpt(1)*qpt(2)*qpt(3)
513 :
514 175969 : do ii=1,3*Invar%natom_unitcell
515 175969 : write(51,'(i5,5x,f15.6)') ii,omega(ii)*Ha_cmm1
516 : enddo
517 22588 : write(51,*) 'Phonon Eigenvectors'
518 22588 : write(51,'(a,17x,a,33x,a,34x,a)') 'Mode Ion','X','Y','Z'
519 175969 : do ii=1,3*Invar%natom_unitcell
520 779182 : do jj=1,Invar%natom_unitcell
521 603213 : write(51,'(i2,3x,i2,3x,6(f15.12,1x))') ii,jj,real(eigenV(3*(jj-1)+1,ii)),&
522 603213 : & aimag(eigenV(3*(jj-1)+1,ii)),&
523 603213 : & real(eigenV(3*(jj-1)+2,ii)),&
524 603213 : & aimag(eigenV(3*(jj-1)+2,ii)),&
525 603213 : & real(eigenV(3*(jj-1)+3,ii)),&
526 1359807 : & aimag(eigenV(3*(jj-1)+3,ii))
527 : enddo
528 : enddo
529 :
530 :
531 22588 : ABI_FREE(omega)
532 22588 : ABI_FREE(dij)
533 22588 : ABI_FREE(eigenV)
534 :
535 22588 : end subroutine tdep_write_dij
536 :
537 : !=====================================================================================================
538 300160 : subroutine tdep_build_phi2_33(isym,Phi2_ref,Phi2_33,Sym,itrans)
539 :
540 : type(Symmetries_type),intent(in) :: Sym
541 : ! type(atdep_dataset_type),intent(in) :: Invar
542 : double precision, intent(in) :: Phi2_ref(3,3)
543 : double precision, intent(out) :: Phi2_33(3,3)
544 : integer,intent(in) :: isym,itrans
545 :
546 : double precision :: Phi2_tmp(3,3),tmp1(3,3)
547 :
548 : ! Transform in the new basis wrt S_ref
549 300160 : call DGEMM('N','N',3,3,3,1.d0,Sym%S_ref(:,:,isym,1),3,Phi2_ref,3,0.d0,Phi2_tmp,3)
550 300160 : call DGEMM('N','N',3,3,3,1.d0,Phi2_tmp,3,Sym%S_inv(:,:,isym,1),3,0.d0,Phi2_33,3)
551 :
552 300160 : if ((itrans.lt.1).or.(itrans.gt.2)) then
553 0 : ABI_BUG('This value of the symmetry index is not permitted')
554 : end if
555 : ! Transpose the 3x3 matrix if required
556 300160 : if (itrans.eq.2) then
557 46327 : tmp1(:,:)=Phi2_33(:,:)
558 46327 : Phi2_33(1,2)=tmp1(2,1)
559 46327 : Phi2_33(1,3)=tmp1(3,1)
560 46327 : Phi2_33(2,3)=tmp1(3,2)
561 46327 : Phi2_33(2,1)=tmp1(1,2)
562 46327 : Phi2_33(3,1)=tmp1(1,3)
563 46327 : Phi2_33(3,2)=tmp1(2,3)
564 : end if
565 :
566 300160 : end subroutine tdep_build_phi2_33
567 :
568 : !=====================================================================================================
569 88 : subroutine tdep_init_eigen2nd(Eigen2nd,natom_unitcell,nqpt)
570 :
571 : integer, intent(in) :: natom_unitcell,nqpt
572 : type(Eigen_type),intent(out) :: Eigen2nd
573 :
574 309521 : ABI_MALLOC(Eigen2nd%eigenval,(3*natom_unitcell,nqpt)); Eigen2nd%eigenval(:,:) =zero
575 13432522 : ABI_MALLOC(Eigen2nd%eigenvec,(2,3,natom_unitcell,3,natom_unitcell,nqpt)); Eigen2nd%eigenvec(:,:,:,:,:,:)=zero
576 13432346 : ABI_MALLOC(Eigen2nd%dynmat ,(2,3,natom_unitcell,3,natom_unitcell,nqpt)); Eigen2nd%dynmat (:,:,:,:,:,:)=zero
577 :
578 88 : end subroutine tdep_init_eigen2nd
579 :
580 : !=====================================================================================================
581 88 : subroutine tdep_destroy_eigen2nd(Eigen2nd)
582 :
583 : type(Eigen_type),intent(inout) :: Eigen2nd
584 :
585 88 : ABI_FREE(Eigen2nd%eigenval)
586 88 : ABI_FREE(Eigen2nd%eigenvec)
587 88 : ABI_FREE(Eigen2nd%dynmat)
588 :
589 88 : end subroutine tdep_destroy_eigen2nd
590 :
591 : !=====================================================================================================
592 44 : subroutine tdep_write_yaml(Eigen2nd,Qpt,Prefix)
593 :
594 : type(Eigen_type),intent(in) :: Eigen2nd
595 : type(Qpoints_type),intent(in) :: Qpt
596 : character(len=*) :: Prefix
597 :
598 : integer :: ii,jj,iatcell,jatcell,iqpt,imode,nmode
599 : double precision :: distance
600 44 : double complex, allocatable :: eigenV(:,:)
601 :
602 44 : nmode=size(Eigen2nd%eigenval,dim=1)
603 44 : open(unit=52,file=trim(Prefix)//'_phonon-bands.yaml')
604 44 : write(52,'(a,i4)') 'nqpoint:',Qpt%nqpt
605 44 : write(52,'(a,i4)') 'npath:',Qpt%qpt_tot-1
606 44 : write(52,'(a)') 'segment_nqpoint:'
607 270 : do ii=1,Qpt%qpt_tot-1
608 270 : write(52,'(a,i4)') '- ',Qpt%lgth_segments(ii)
609 : end do
610 44 : write(52,'(a,i4)') 'natom:',nmode/3
611 44 : write(52,'(a)') 'phonon:'
612 44 : distance=0.d0
613 22632 : do iqpt=1,Qpt%nqpt
614 22588 : 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),']'
615 90220 : if (iqpt.gt.1) distance=distance+sqrt(sum((Qpt%qpt_cart(:,iqpt)-Qpt%qpt_cart(:,iqpt-1))**2))
616 22588 : write(52,'(a,f15.6)') ' distance:',distance
617 178068 : do ii=1,Qpt%qpt_tot
618 645318 : if (sum(abs(Qpt%qpt_red(:,iqpt)-Qpt%special_red(ii,:))).lt.tol8) then
619 270 : write(52,'(3a)') " label: '",trim(Qpt%special_qpt(ii)),"'"
620 270 : exit
621 : end if
622 : end do !ii
623 2053372 : ABI_MALLOC(eigenV,(nmode,nmode)); eigenV(:,:)=czero
624 73715 : do iatcell=1,nmode/3
625 274786 : do jatcell=1,nmode/3
626 855411 : do ii=1,3
627 2613923 : do jj=1,3
628 : eigenV((iatcell-1)*3+ii,(jatcell-1)*3+jj)=dcmplx(Eigen2nd%eigenvec(1,ii,iatcell,jj,jatcell,iqpt),&
629 2412852 : & Eigen2nd%eigenvec(2,ii,iatcell,jj,jatcell,iqpt))
630 : end do !ii
631 : end do !jj
632 : end do !jatcell
633 : end do !iatcell
634 22588 : write(52,'(a)') ' band:'
635 175969 : do imode=1,nmode
636 153381 : write(52,'(a,i4)') ' - #',imode
637 153381 : write(52,'(a,f15.6)') ' frequency:',Eigen2nd%eigenval(imode,iqpt)*Ha_THz
638 153381 : write(52,'(a)') ' eigenvector:'
639 779182 : do iatcell=1,nmode/3
640 603213 : write(52,'(a,i4)') " - # atom ", iatcell
641 2566233 : do ii=1,3
642 1809639 : write(52,'(a,f18.9,a,f18.9,a)') " - [",real(eigenV((iatcell-1)*3+ii,imode)),','&
643 4222491 : & ,aimag(eigenV((iatcell-1)*3+ii,imode)),']'
644 : end do
645 : end do
646 : end do !nmode
647 22588 : ABI_FREE(eigenV)
648 22632 : write(52,'(a)') ''
649 : end do
650 :
651 44 : close(52)
652 :
653 :
654 44 : end subroutine tdep_write_yaml
655 :
656 : !=====================================================================================================
657 0 : end module m_tdep_phi2
|
