Line data Source code
1 : module m_dfpt_dynmat_eig
2 :
3 : use m_types
4 :
5 : implicit none
6 :
7 : contains
8 :
9 0 : subroutine DiagonalizeDynMat(atoms, qvec, calcEv, dynMat, w, a, iqpt, l_scalemass, add_tag,l_sumrule)
10 :
11 : USE m_juDFT_stop
12 : implicit none
13 :
14 : ! Type parameters
15 : type(t_atoms), intent(in) :: atoms
16 : real, intent(in) :: qvec(3)
17 :
18 : ! Scalar parameters
19 : logical, intent(in) :: calcEv
20 :
21 : ! Array parameters
22 : complex, intent(in) :: dynMat(:, :)
23 : real, allocatable, intent(out) :: w(:)
24 : complex, allocatable, intent(out) :: a(:, :)
25 : logical, intent(in) :: l_scalemass
26 : character(len=*), intent(in) :: add_tag
27 : logical, intent(in) :: l_sumrule
28 :
29 : ! Array parameters
30 :
31 : ! Scalar variables
32 : ! jobVl : (LAPACK) switch for calculation of left eigenvector
33 : ! jobVr : (LAPACK) switch for calculation of right eigenvector
34 : ! n : (LAPACK) order of the matrix a to diagonalize
35 : ! ldA : (LAPACK) leading dimension of the array a
36 : ! ldVl : (LAPACK) leading dimension of the array vL
37 : ! ldVr : (LAPACK) leading dimension of the array vR
38 : ! lWork : (LAPACK) dimension of the array work or if equals -1 switch for LAPACK to determine optimal size of work array
39 : ! info : (LAPACK) status of calculation
40 : character :: jobz
41 : character :: uplo
42 : integer :: n
43 : integer :: ldA
44 : integer :: lWork
45 : integer :: info
46 : integer :: ii
47 : integer :: jj
48 : integer :: itype
49 : integer :: ieqat
50 : integer :: iatom
51 : integer, intent(in) :: iqpt
52 0 : complex, allocatable :: dynMat0(:, :)
53 : character(len=100) :: trash
54 : integer :: iread, iDir
55 0 : real :: numbers(3*atoms%nat,6*atoms%nat)
56 :
57 : ! Array variables
58 : ! a : (LAPACK) matrix to diagonalize
59 : ! w : (LAPACK) computed eigenvalues
60 : ! work : (LAPACK) swap array for LAPACK diagonalization routine
61 : ! rwork : (LAPACK) swap array for LAPACK diagonalization routine
62 0 : real, allocatable :: rwork(:)
63 0 : complex, allocatable :: work(:)
64 : character(len=:), allocatable :: filename
65 : character(len=50) :: filenameTemp
66 : REAL :: atomic_mass_array(118)
67 :
68 : atomic_mass_array = [1.01, 4.00, 6.94, 9.01, 10.81, 12.01, 14.01, 16.00, 19.00, 20.18, & ! up to neon
69 : & 22.99, 24.31, 26.98, 28.09, 30.97, 32.06, 35.45, 39.95, & ! up to argon
70 : & 39.10, 40.08, 44.96, 47.87, 50.94, 52.00, 54.94, 55.85, 58.93, & ! up to cobalt
71 : & 58.69, 63.55, 65.38, 69.72, 72.63, 74.92, 78.97, 79.90, 83.80, & ! up to krypton
72 : & 85.47, 87.62, 88.91, 91.22, 92.91, 95.95, 97.40, 101.07, 102.91, & ! up to ruthenium
73 : & 106.42, 107.87, 112.41, 114.82, 118.71, 121.76, 127.60, 126.90, 131.29, & ! up to xenon
74 : & 132.91, 137.33, 138.91, 140.12, 140.91, 144.24, 146.00, 150.36, 151.96, & ! up to europium
75 : & 157.25, 158.93, 162.50, 164.93, 167.26, 168.93, 173.05, 174.97, 178.49, & ! up to hafnium
76 : & 180.95, 183.84, 186.21, 190.23, 192.22, 195.08, 196.97, 200.59, 204.38, & ! up to thallium
77 : & 207.20, 208.98, 209.98, 210.00, 222.00, 223.00, 226.00, 227.00, 232.04, & ! up to thorium
78 : & 231.04, 238.03, 237.00, 244.00, 243.00, 247.00, 247.00, 251.00, 252.00, & ! up to einsteinium
79 : & 257.00, 258.00, 259.00, 262.00, 267.00, 269.00, 270.00, 272.00, 273.00, & ! up to hassium
80 0 : & 277.00, 281.00, 281.00, 285.00, 286.00, 289.00, 288.00, 293.00, 294.00, 294.00]
81 :
82 : ! TODO: This is ridiculous. Remove asap.
83 0 : if (iqpt < 10) then
84 0 : write(filenameTemp, '("dynMatq=000",i1)') iqpt
85 0 : else if (iqpt > 9 .and. iqpt < 100) then
86 0 : write(filenameTemp, '("dynMatq=00",i2)') iqpt
87 0 : else if (iqpt > 99 .and. iqpt < 1000) then
88 0 : write(filenameTemp, '("dynMatq=0",i3)') iqpt
89 0 : else if (iqpt > 999 .and. iqpt < 10000) then
90 0 : write(filenameTemp, '("dynMatq=",i4)') iqpt
91 : end if
92 0 : if ((TRIM(add_tag)).EQ."full".OR.(TRIM(add_tag).EQ."band")) filenameTemp = TRIM(add_tag)//"_"//filenameTemp
93 : ! filename = trim(filenameTemp)
94 : ! write(*, *) filename
95 0 : open( 109, file=filenameTemp, status='replace', action='write', form='formatted')
96 :
97 : ! Set parameter for LAPACK diagonalization routine
98 0 : if (calcEv ) then
99 0 : jobz = 'V'
100 : else
101 0 : jobz = 'N'
102 : end if
103 :
104 : ! Evaluate the lower triangle of the Dynamical matrix
105 0 : uplo = 'U'
106 :
107 : ! Dimensions of the Dynamical matrix
108 0 : n = 3 * atoms%nat
109 0 : lda = 3 * atoms%nat
110 :
111 : ! Copy Dynamical matrix into work array as it is destroyed if not calculating eigenvectors
112 0 : allocate( a( lda, n ) )
113 0 : do jj = 1, n
114 0 : do ii = 1, lda
115 0 : a(ii, jj) = (dynMat(ii, jj) + conjg(dynMat(jj, ii)))/2.0
116 0 : IF (l_scalemass) a(ii, jj) = a(ii, jj)/SQRT(atomic_mass_array(atoms%nz(CEILING(jj/3.0)))*atomic_mass_array(atoms%nz(CEILING(ii/3.0))))
117 : end do
118 : end do
119 :
120 0 : IF (l_sumrule) THEN
121 0 : IF (iqpt/=1) THEN
122 0 : ALLOCATE(dynmat0,mold=dynmat)
123 0 : OPEN( 110, file="dynMatq=gamma", status="old")
124 0 : DO iread = 1, 3 + 3*atoms%nat ! Loop over dynmat rows
125 0 : IF (iread<4) THEN
126 0 : READ( 110,*) trash
127 0 : write(*,*) iread, trash
128 : ELSE
129 0 : READ( 110,*) numbers(iread-3,:)
130 0 : write(*,*) iread, numbers(iread-3,:)
131 0 : dynmat0(iread-3,:) = CMPLX(numbers(iread-3,::2),numbers(iread-3,2::2))
132 : END IF
133 : END DO ! iread
134 0 : CLOSE(110)
135 : ELSE
136 0 : ALLOCATE(dynmat0,mold=dynmat)
137 0 : dynmat0 = a
138 : ! Write non-corrected dynMatq for sumrule at q!=gamma
139 0 : OPEN( 111, file="dynMatq=gamma",status="replace",action="write",form="formatted")
140 0 : write(111, '(a,3f9.3)') 'q =', qvec
141 0 : write(111, '(a)') '==================================='
142 0 : write(111, '(a)')
143 0 : write(111, '(a)') 'Original Dynamical Matrix [mass corrected]'
144 0 : DO ii = 1, lda
145 0 : write(111, '(3(2(es16.8,1x),3x))') dynmat0(ii, :)
146 : END DO
147 0 : CLOSE(111)
148 : END IF
149 :
150 0 : allocate( w(n))
151 0 : w = 0.
152 :
153 0 : allocate( rwork(3 * n))
154 0 : rwork = 0.
155 :
156 0 : lwork = 2 * n
157 0 : allocate(work(lwork))
158 0 : work = 0.
159 :
160 : ! Diagonalize Gamma dynamical matrix
161 0 : call zheev( jobz, uplo, n, dynmat0, lda, w, work, lwork, rwork, info )
162 :
163 : ! w: eigenvalues q=0, a: eigenvectors q=0
164 0 : DO jj = 1, lda
165 0 : DO ii = 1, lda
166 0 : DO iDir = 1, 3
167 0 : a(jj, ii) = a(jj, ii) - w(iDir)*dynmat0(jj,iDir)*conjg(dynmat0(ii,iDir))
168 : END DO
169 : END DO
170 : END DO
171 0 : DEALLOCATE(w, rwork, work)
172 0 : DEALLOCATE(dynmat0)
173 : END IF
174 :
175 0 : write(*, '(a,3f9.3)') 'q =', qvec
176 0 : write(*, '(a)') '==================================='
177 0 : write(*, '(a)')
178 0 : write(*, '(a)') 'Original Dynamical Matrix [mass corrected]'
179 0 : DO ii = 1, lda
180 0 : write(*, '(3(2(es16.8,1x),3x))') a(ii, :)
181 : END DO
182 :
183 0 : write(109, '(a,3f9.3)') 'q =', qvec
184 0 : write(109, '(a)') '==================================='
185 0 : write(109, '(a)')
186 0 : write(109, '(a)') 'Original Dynamical Matrix [mass corrected]'
187 0 : DO ii = 1, lda
188 0 : write(109, '(3(2(es16.8,1x),3x))') a(ii, :)
189 : END DO
190 :
191 0 : allocate( w(n))
192 0 : w = 0.
193 :
194 0 : allocate( rwork(3 * n))
195 0 : rwork = 0.
196 :
197 0 : lwork = 2 * n
198 0 : allocate(work(lwork))
199 0 : work = 0.
200 :
201 : ! Diagonalize dynamical matrix
202 0 : call zheev( jobz, uplo, n, a, lda, w, work, lwork, rwork, info )
203 :
204 0 : if ( info < 0 ) then
205 0 : write(*, *) 'The ', info, '(st/nd/)th argument had an illegal value'
206 0 : CALL juDFT_error("Illegal argument value.",calledby="DiagonalizeDynMat")
207 0 : else if ( info > 0 ) then
208 0 : write(*, *) 'The diagonalization algorithm failed to converge; ', info, ' off-diagonal elements of an intermediate tridiaonal form&
209 0 : & did not converge to zero.'
210 0 : CALL juDFT_error("Diagonalization failed.",calledby="DiagonalizeDynMat")
211 : end if
212 :
213 0 : write(*, '(a)') 'The eigenvalues of the Dynamical matrix are:'
214 0 : iatom = 0
215 0 : do itype = 1, atoms%ntype
216 0 : do ieqat = 1, atoms%neq(itype)
217 0 : iatom = iatom + 1
218 0 : write(*, "(a,i2,a,1x,3(es16.8,1x),',',5x)") 'Atom', iatom, ':', w((iatom - 1)* 3 + 1:(iatom - 1)* 3 + 3)
219 0 : write(109, "(a,i2,a,1x,3(es16.8,1x),',',5x)") 'Atom', iatom, ':', w((iatom - 1)* 3 + 1:(iatom - 1)* 3 + 3)
220 : end do ! ieqat
221 : end do ! itype
222 :
223 0 : if ( calcEv ) then
224 : !todo make the output nicer
225 0 : DO ii = 1, lda
226 0 : write(*, '(3(2(es16.8,1x),3x))') a(ii, :)
227 0 : write(*, *)
228 0 : write(109, '(3(2(es16.8,1x),3x))') a(ii, :)
229 0 : write(109, *)
230 : END DO
231 : else
232 0 : a(:, :) = cmplx(0., 0.)
233 : end if
234 :
235 0 : close( 109 )
236 :
237 :
238 :
239 0 : end subroutine diagonalizeDynMat
240 :
241 0 : subroutine CalculateFrequencies( atoms, iqpt, eigenVals, eigenFreqs, add_tag )
242 :
243 : implicit none
244 :
245 : ! Type parameter
246 : type(t_atoms), intent(in) :: atoms
247 :
248 : ! Scalar parameter
249 : integer, intent(in) :: iqpt
250 :
251 : ! Array parameter
252 : real, intent(in) :: eigenVals(:)
253 : complex, allocatable, intent(out) :: eigenFreqs(:)
254 : character(len=*), intent(in) :: add_tag
255 :
256 : ! Scalar variables
257 : integer :: itype
258 : integer :: ieqat
259 : integer :: iatom
260 : integer :: idir
261 : real :: massInElectronMasses
262 : real :: convFact
263 :
264 : ! Array variables
265 : character(len=50) :: filenameTemp
266 : REAL :: atomic_mass_array(118)
267 :
268 : ! TODO: This is ridiculous. Remove asap.
269 0 : if (iqpt < 10) then
270 0 : write(filenameTemp, '("dynMatq=000",i1)') iqpt
271 0 : else if (iqpt > 9 .and. iqpt < 100) then
272 0 : write(filenameTemp, '("dynMatq=00",i2)') iqpt
273 0 : else if (iqpt > 99 .and. iqpt < 1000) then
274 0 : write(filenameTemp, '("dynMatq=0",i3)') iqpt
275 0 : else if (iqpt > 999 .and. iqpt < 10000) then
276 0 : write(filenameTemp, '("dynMatq=",i4)') iqpt
277 : end if
278 0 : if ((TRIM(add_tag)).EQ."full".OR.(TRIM(add_tag).EQ."band")) filenameTemp = TRIM(add_tag)//"_"//filenameTemp
279 0 : open( 109, file=filenameTemp, status='old', action='write', form='formatted', position='append')
280 :
281 0 : allocate(eigenFreqs(3*atoms%nat))
282 0 : eigenFreqs = 0.
283 0 : itype = 1
284 :
285 : atomic_mass_array = [1.01, 4.00, 6.94, 9.01, 10.81, 12.01, 14.01, 16.00, 19.00, 20.18, & ! up to neon
286 : & 22.99, 24.31, 26.98, 28.09, 30.97, 32.06, 35.45, 39.95, & ! up to argon
287 : & 39.10, 40.08, 44.96, 47.87, 50.94, 52.00, 54.94, 55.85, 58.93, & ! up to cobalt
288 : & 58.69, 63.55, 65.38, 69.72, 72.63, 74.92, 78.97, 79.90, 83.80, & ! up to krypton
289 : & 85.47, 87.62, 88.91, 91.22, 92.91, 95.95, 97.40, 101.07, 102.91, & ! up to ruthenium
290 : & 106.42, 107.87, 112.41, 114.82, 118.71, 121.76, 127.60, 126.90, 131.29, & ! up to xenon
291 : & 132.91, 137.33, 138.91, 140.12, 140.91, 144.24, 146.00, 150.36, 151.96, & ! up to europium
292 : & 157.25, 158.93, 162.50, 164.93, 167.26, 168.93, 173.05, 174.97, 178.49, & ! up to hafnium
293 : & 180.95, 183.84, 186.21, 190.23, 192.22, 195.08, 196.97, 200.59, 204.38, & ! up to thallium
294 : & 207.20, 208.98, 209.98, 210.00, 222.00, 223.00, 226.00, 227.00, 232.04, & ! up to thorium
295 : & 231.04, 238.03, 237.00, 244.00, 243.00, 247.00, 247.00, 251.00, 252.00, & ! up to einsteinium
296 : & 257.00, 258.00, 259.00, 262.00, 267.00, 269.00, 270.00, 272.00, 273.00, & ! up to hassium
297 : & 277.00, 281.00, 281.00, 285.00, 286.00, 289.00, 288.00, 293.00, 294.00, 294.00]
298 :
299 0 : massInElectronMasses = 1836.15
300 :
301 0 : convFact = 4.35974472220e-18 / (5.2918e-11)**2 / 9.1093837015e-31 / massInElectronMasses
302 :
303 0 : write(*, *)
304 0 : write(*, '(a)') 'Eigenfrequencies in THz'
305 0 : write(109, *)
306 0 : write(109, '(a)') 'Eigenfrequencies in THz'
307 0 : iatom = 0
308 0 : do itype = 1, atoms%ntype
309 0 : do ieqat = 1, atoms%neq(itype)
310 0 : iatom = iatom + 1
311 0 : do idir = 1, 3
312 : !if (abs(eigenVals((iatom - 1) * 3 + idir)) < 1e-5) then
313 : ! eigenFreqs((iatom - 1) * 3 + idir) = cmplx(0., 0.)
314 0 : if (eigenVals((iatom - 1) * 3 + idir) < 0 ) then
315 0 : eigenFreqs((iatom - 1) * 3 + idir) = cmplx(0., -sqrt(abs(eigenVals((iatom - 1) * 3 + idir)) * convFact) / tpi_const * 1e-12)
316 : else
317 0 : eigenFreqs((iatom - 1) * 3 + idir) = sqrt(eigenVals((iatom - 1) * 3 + idir) * convFact) / tpi_const * 1e-12
318 : end if
319 : end do ! idir
320 0 : write(*, "(a,i2,a,1x,3(2es16.8,1x),',',5x)") 'Atom', iatom, ':', eigenFreqs((iatom - 1)* 3 + 1:(iatom - 1)* 3 + 3)
321 0 : write(109, "(a,i2,a,1x,3(2es16.8,1x),',',5x)") 'Atom', iatom, ':', eigenFreqs((iatom - 1)* 3 + 1:(iatom - 1)* 3 + 3)
322 : end do ! ieqat
323 : end do ! itype
324 :
325 :
326 0 : write(*, *)
327 0 : write(*, '(a)') 'Eigenfrequencies in 1/cm'
328 0 : write(109, *)
329 0 : write(109, '(a)') 'Eigenfrequencies in 1/cm'
330 0 : iatom = 0
331 0 : do itype = 1, atoms%ntype
332 0 : do ieqat = 1, atoms%neq(itype)
333 0 : iatom = iatom + 1
334 0 : do idir = 1, 3
335 0 : eigenFreqs((iatom - 1) * 3 + idir) = eigenFreqs((iatom - 1) * 3 + idir) * 33
336 : end do ! idir
337 0 : write(*, "(a,i2,a,1x,3(2es16.8,1x),',',5x)") 'Atom', iatom, ':', eigenFreqs((iatom - 1)* 3 + 1:(iatom - 1)* 3 + 3)
338 0 : write(109, "(a,i2,a,1x,3(2es16.8,1x),',',5x)") 'Atom', iatom, ':', eigenFreqs((iatom - 1)* 3 + 1:(iatom - 1)* 3 + 3)
339 : end do ! ieqat
340 : end do ! itype
341 0 : close( 109 )
342 :
343 0 : end subroutine CalculateFrequencies
344 :
345 : end module m_dfpt_dynmat_eig
|
