Branch data 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
# # # # #
# # # # #
# # # # ]
|
