File: phdos.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (170 lines) | stat: -rw-r--r-- 4,176 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170

! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine phdos
use modmain
use modphonon
use modtest
implicit none
! local variables
integer iq,i,iw
integer i1,i2,i3
real(8) wmin,wmax,wd,dw
real(8) tmax,temp(ntemp),s(ntemp)
real(8) v(3),t1,t2
! allocatable arrays
real(8), allocatable :: wp(:),w(:),gw(:),f(:)
complex(8), allocatable :: dynq(:,:,:),dynr(:,:,:)
complex(8), allocatable :: dynp(:,:),ev(:,:)
! external functions
real(8) splint
external splint
! initialise universal variables
call init0
call init2
allocate(wp(nbph),w(nwplot),gw(nwplot),f(nwplot))
allocate(dynq(nbph,nbph,nqpt),dynr(nbph,nbph,nqptnr))
allocate(dynp(nbph,nbph),ev(nbph,nbph))
! read in the dynamical matrices
call readdyn(dynq)
! apply the acoustic sum rule
call sumrule(dynq)
! Fourier transform the dynamical matrices to real-space
call dynqtor(dynq,dynr)
! find the minimum and maximum frequencies
wmin=0.d0
wmax=0.d0
do iq=1,nqpt
  call dynev(dynq(:,:,iq),wp,ev)
  wmin=min(wmin,wp(1))
  wmax=max(wmax,wp(nbph))
end do
wmax=wmax+(wmax-wmin)*0.1d0
wmin=wmin-(wmax-wmin)*0.1d0
wd=wmax-wmin
if (wd.lt.1.d-8) wd=1.d0
dw=wd/dble(nwplot)
! generate energy grid
do iw=1,nwplot
  w(iw)=dw*dble(iw-1)+wmin
end do
gw(:)=0.d0
do i1=0,ngrkf-1
  v(1)=dble(i1)/dble(ngrkf)
  do i2=0,ngrkf-1
    v(2)=dble(i2)/dble(ngrkf)
    do i3=0,ngrkf-1
      v(3)=dble(i3)/dble(ngrkf)
! compute the dynamical matrix at this particular q-point
      call dynrtoq(v,dynr,dynp)
! find the phonon frequencies
      call dynev(dynp,wp,ev)
      do i=1,nbph
        t1=(wp(i)-wmin)/dw+1.d0
        iw=nint(t1)
        if ((iw.ge.1).and.(iw.le.nwplot)) then
          gw(iw)=gw(iw)+1.d0
        end if
      end do
    end do
  end do
end do
t1=1.d0/(dw*dble(ngrkf)**3)
gw(:)=t1*gw(:)
! smooth phonon DOS if required
if (nswplot.gt.0) call fsmooth(nswplot,nwplot,gw)
! write phonon DOS to file
open(50,file='PHDOS.OUT',form='FORMATTED')
do iw=1,nwplot
  write(50,'(2G18.10)') w(iw),gw(iw)
end do
close(50)
write(*,*)
write(*,'("Info(phdos):")')
write(*,'(" phonon density of states written to PHDOS.OUT")')
!-------------------------------------------!
!     thermodynamic properties from DOS     !
!-------------------------------------------!
! maximum temperature
tmax=wmax/kboltz
! temperature grid
do i=1,ntemp
  temp(i)=tmax*dble(i)/dble(ntemp)
end do
open(50,file='THERMO.OUT',form='FORMATTED')
! zero point energy
do iw=1,nwplot
  f(iw)=gw(iw)*w(iw)
end do
t1=splint(nwplot,w,f)
t1=0.5d0*dble(natmtot)*t1
write(50,*)
write(50,'("Zero-point energy : ",G18.10)') t1
! vibrational energy
write(50,*)
write(50,'("Vibrational energy vs. temperature :")')
do i=1,ntemp
  do iw=1,nwplot
    t1=w(iw)/(2.d0*kboltz*temp(i))
    if (t1.gt.0.d0) then
      f(iw)=gw(iw)*w(iw)*cosh(t1)/sinh(t1)
    else
      f(iw)=0.d0
    end if
  end do
  t1=splint(nwplot,w,f)
  t1=0.5d0*dble(natmtot)*t1
  write(50,'(2G18.10)') temp(i),t1
  s(i)=t1
end do
! free energy
write(50,*)
write(50,'("Free energy vs. temperature :")')
do i=1,ntemp
  do iw=1,nwplot
    t1=2.d0*sinh(w(iw)/(2.d0*kboltz*temp(i)))
    if (t1.gt.0.d0) then
      f(iw)=gw(iw)*log(t1)
    else
      f(iw)=0.d0
    end if
  end do
  t1=splint(nwplot,w,f)
  t1=dble(natmtot)*kboltz*temp(i)*t1
  write(50,'(2G18.10)') temp(i),t1
! compute entropy from S = (F-E)/T
  s(i)=(s(i)-t1)/temp(i)
end do
! entropy
write(50,*)
write(50,'("Entropy vs. temperature :")')
do i=1,ntemp
  write(50,'(2G18.10)') temp(i),s(i)
end do
! heat capacity
write(50,*)
write(50,'("Heat capacity vs. temperature :")')
do i=1,ntemp
  do iw=1,nwplot
    t1=w(iw)/(kboltz*temp(i))
    t2=exp(t1)-1.d0
    if (abs(t2).gt.1.d-14) then
      f(iw)=gw(iw)*(t1**2)*(t2+1.d0)/t2**2
    else
      f(iw)=0.d0
    end if
  end do
  t1=splint(nwplot,w,f)
  t1=dble(natmtot)*kboltz*t1
  write(50,'(2G18.10)') temp(i),t1
end do
close(50)
write(*,'(" thermodynamic properties written to THERMO.OUT")')
! write phonon DOS to test file
call writetest(210,'phonon DOS',nv=nwplot,tol=1.d-2,rva=gw)
deallocate(wp,w,gw,f,dynq,dynr,dynp,ev)
return
end subroutine