mirror of https://gitlab.com/QEF/q-e.git
*** empty log message ***
git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@633 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz
parent
e1739befad
commit
710179c08f
|
|
@ -1,354 +0,0 @@
|
|||
!
|
||||
! Copyright (C) 2001 PWSCF group
|
||||
! This file is distributed under the terms of the
|
||||
! GNU General Public License. See the file `License'
|
||||
! in the root directory of the present distribution,
|
||||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
!
|
||||
!--------------------------------------------------------------------
|
||||
program dynmat
|
||||
!--------------------------------------------------------------------
|
||||
!
|
||||
! read a dynamical matrix at q=0 , diagonalise it, calculate IR
|
||||
! Input data (namelist "input")
|
||||
!
|
||||
! flmat character input file containing the dynamical matrix
|
||||
! (default: flmat='dynmat')
|
||||
! q(3) real calculate LO modes (add nonanalytic terms) along
|
||||
! the direction q (default: q=(0,0,0) )
|
||||
! amass(nt) real mass for atom type nt, a.u.
|
||||
! (default: amass is read from file flmat)
|
||||
! asr logical impose Acoustic Sum Rule (default:asr=.true.)
|
||||
! flout character output file containing frequencies and modes
|
||||
! (default: flle='dynout')
|
||||
! flmol character as above, in a format suitable for 'molden'
|
||||
! (default: flmol='moldout')
|
||||
!
|
||||
implicit none
|
||||
integer nax
|
||||
parameter (nax=30)
|
||||
character(len=50) flmat, flout, flmol
|
||||
character(len=3) atm(nax)
|
||||
logical asr, lread, gamma
|
||||
complex(kind=8) dyn(3,3,nax,nax), z(3*nax,3*nax)
|
||||
real(kind=8) tau(3,nax), amass(nax), amass_(nax), zstar(3,3,nax),&
|
||||
eps0(3,3), a0, omega, amconv, q(3), q_(3), w2(3*nax)
|
||||
real(kind=8) deps_dtau(3,3,3,nax)
|
||||
integer ityp(nax), itau(nax), nat, na, nt, ntyp, nu
|
||||
external latgen, recips, error, dyndiag, nonanal, readmat
|
||||
namelist /input/ amass,asr,flmat,flout,flmol,q
|
||||
!
|
||||
!
|
||||
asr =.true.
|
||||
flout=' '
|
||||
flmat='dynmat'
|
||||
flout='dynout'
|
||||
flmol='moldout'
|
||||
do nt=1, ntyp
|
||||
amass(nt)=0.0
|
||||
end do
|
||||
q(1)=0.0
|
||||
q(2)=0.0
|
||||
q(3)=0.0
|
||||
!
|
||||
read (5,input)
|
||||
!
|
||||
inquire(file=flmat,exist=lread)
|
||||
if (lread) then
|
||||
write(6,'(/5x,a,a)') 'Reading Dynamical Matrix from file ',&
|
||||
& flmat
|
||||
else
|
||||
write(6,'(/5x,a)') 'file not found', flmat
|
||||
stop
|
||||
end if
|
||||
!
|
||||
call readmat (flmat,asr,nax,nat,ntyp,ityp,atm,a0,omega, amass_&
|
||||
&,tau,zstar,eps0,dyn,deps_dtau,q_)
|
||||
!
|
||||
gamma = abs(q_(1)**2+q_(2)**2+q_(3)**2).lt.1.0e-8
|
||||
amconv = 1.66042e-24/9.1095e-28*0.5
|
||||
do nt=1, ntyp
|
||||
if (amass(nt).gt.0.0) then
|
||||
amass(nt)=amass(nt)*amconv
|
||||
else
|
||||
amass(nt)=amass_(nt)
|
||||
end if
|
||||
end do
|
||||
!
|
||||
if (gamma) then
|
||||
do na=1,nat
|
||||
itau(na)=na
|
||||
end do
|
||||
call nonanal (nax,nat,dyn,q,itau,nax,eps0,zstar,omega)
|
||||
end if
|
||||
!
|
||||
call dyndiag(nax,nat,amass,ityp,dyn,w2,z)
|
||||
!
|
||||
call writemodes(nax,nat,q_,w2,z,flout)
|
||||
!
|
||||
call writemolden(nax,nat,atm,a0,tau,ityp,w2,z,flmol)
|
||||
!
|
||||
if (gamma) call RamanIR &
|
||||
(nax, nat, omega, w2, z, zstar, eps0, deps_dtau)
|
||||
!
|
||||
stop
|
||||
end
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
subroutine readmat (flmat,asr,nax,nat,ntyp,ityp,atm,a0, &
|
||||
omega,amass,tau,zstar,eps0,dynr,deps_dtau,q)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
implicit none
|
||||
character(len=*) flmat
|
||||
integer nax, nat, ntyp, ityp(nax)
|
||||
character(len=3) atm(ntyp)
|
||||
real(kind=8) amass(ntyp), tau(3,nax), a0, omega
|
||||
real(kind=8) dynr(2,3,3,nax,nax), eps0(3,3), zstar(3,3,nax), &
|
||||
deps_dtau(3,3,3,nax), q(3)
|
||||
logical asr
|
||||
!
|
||||
character(len=80) line
|
||||
real(kind=8) at(3,3), celldm(6), sum
|
||||
integer ibrav, nt, na, nb, naa, nbb, i, j, k
|
||||
logical qfinito, noraman
|
||||
!
|
||||
!
|
||||
noraman=.true.
|
||||
open (unit=1,file=flmat,status='old',form='formatted')
|
||||
read(1,'(a)') line
|
||||
read(1,'(a)') line
|
||||
read(1,*) ntyp,nat,ibrav,celldm
|
||||
if (nat.gt.nax) stop ' too many atoms '
|
||||
a0=celldm(1)
|
||||
call latgen(ibrav,celldm,at(1,1),at(1,2),at(1,3))
|
||||
call volume(a0,at(1,1),at(1,2),at(1,3),omega)
|
||||
do nt=1,ntyp
|
||||
read(1,*) i,atm(nt),amass(nt)
|
||||
end do
|
||||
do na=1,nat
|
||||
read(1,*) i,ityp(na), (tau(j,na),j=1,3)
|
||||
end do
|
||||
read(1,'(a)') line
|
||||
read(1,'(a)') line
|
||||
read(1,'(a)') line
|
||||
read(1,'(a)') line
|
||||
read(line(11:80),*) (q(i),i=1,3)
|
||||
qfinito=q(1).ne.0.0 .or. q(2).ne.0.0 .or. q(3).ne.0.0
|
||||
read(1,'(a)') line
|
||||
do na = 1,nat
|
||||
do nb = 1,nat
|
||||
read (1,*) naa, nbb
|
||||
if (na.ne.naa .or. nb.ne.nbb) then
|
||||
print *, 'na, nb, naa, nbb =', na,nb,naa,nbb
|
||||
stop ' sgrunt!'
|
||||
end if
|
||||
read (1,*) ((dynr(1,i,j,na,nb), &
|
||||
& dynr(2,i,j,na,nb), j=1,3), i=1,3)
|
||||
end do
|
||||
end do
|
||||
!
|
||||
if (.not.qfinito) then
|
||||
read(1,*)
|
||||
read(1,'(a)') line
|
||||
if (line(1:23).ne.' Dielectric Tensor:') then
|
||||
do na=1,nat
|
||||
do j=1,3
|
||||
do i=1,3
|
||||
zstar(i,j,na)=0.0
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
do j=1,3
|
||||
do i=1,3
|
||||
eps0(i,j)=0.0
|
||||
end do
|
||||
eps0(j,j)=1.0
|
||||
end do
|
||||
else
|
||||
read(1,*)
|
||||
read(1,*) ((eps0(i,j), j=1,3), i=1,3)
|
||||
read(1,*)
|
||||
read(1,*)
|
||||
read(1,*)
|
||||
do na = 1,nat
|
||||
read(1,*)
|
||||
read(1,*) ((zstar(i,j,na), j=1,3),i=1,3)
|
||||
end do
|
||||
read(1,*,end=10,err=10)
|
||||
read(1,*,end=10,err=10)
|
||||
read(1,*,end=10,err=10)
|
||||
do na = 1,nat
|
||||
do i = 1, 3
|
||||
read(1,*,end=10,err=10)
|
||||
read(1,*,end=10,err=10) &
|
||||
((deps_dtau(k,j,i,na), j=1,3), k=1,3)
|
||||
end do
|
||||
end do
|
||||
write(6,'(/5x,a)') 'Raman cross sections read'
|
||||
noraman=.false.
|
||||
10 continue
|
||||
end if
|
||||
end if
|
||||
if (noraman) deps_dtau=0.d0
|
||||
if (asr) then
|
||||
!
|
||||
! ASR sulle cariche efficaci
|
||||
!
|
||||
do i=1,3
|
||||
do j=1,3
|
||||
sum=0.0
|
||||
do na=1,nat
|
||||
sum = sum + zstar(i,j,na)
|
||||
end do
|
||||
do na=1,nat
|
||||
zstar(i,j,na) = zstar(i,j,na) - sum/nat
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!
|
||||
! ASR sulla matrice dinamica
|
||||
!
|
||||
do i=1,3
|
||||
do j=1,3
|
||||
do na=1,nat
|
||||
sum=0.0
|
||||
do nb=1,nat
|
||||
if (na.ne.nb) &
|
||||
& sum=sum+dynr(1,i,j,na,nb)
|
||||
end do
|
||||
dynr(1,i,j,na,na) = -sum
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
!
|
||||
close(unit=1)
|
||||
!
|
||||
return
|
||||
end
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
subroutine RamanIR (nax, nat, omega, w2, z, zstar, eps0, deps_dtau)
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
! write IR cross sections
|
||||
! on input: z = eigendisplacements
|
||||
! zstar = effective charges
|
||||
! deps_dtau = derivatives of epsilon wrt atomic displacement
|
||||
use allocate
|
||||
implicit none
|
||||
! input
|
||||
integer nax, nat
|
||||
real(kind=8) omega, w2(3*nat), zstar(3,3,nat), eps0(3,3), &
|
||||
deps_dtau(3,3,3,nat)
|
||||
complex(kind=8) z(3*nax,3*nat)
|
||||
! local
|
||||
integer na, nu, ipol, jpol, lpol
|
||||
logical noraman
|
||||
real(kind=8), pointer :: infrared(:), raman(:,:,:)
|
||||
real(kind=8):: polar(3), rydcm1, cm1thz, freq, u2, r1fac, r2fac, irfac
|
||||
real(kind=8):: alpha, beta2
|
||||
!
|
||||
! conversion factors RYD=>THZ, RYD=>1/CM e 1/CM=>THZ
|
||||
!
|
||||
rydcm1 = 13.6058*8065.5
|
||||
cm1thz = 241.796/8065.5
|
||||
!
|
||||
! conversion factor for IR cross sections from
|
||||
! (Ry atomic units * e^2) to (Debye/A)^2/amu
|
||||
! 1 Ry mass unit = 2 * mass of one electron = 2 amu
|
||||
! 1 e = 4.80324x10^(-10) esu = 4.80324 Debye/A
|
||||
! (1 Debye = 10^(-18) esu*cm = 0.2081928 e*A)
|
||||
!
|
||||
irfac = 4.80324**2/2.d0
|
||||
!
|
||||
! conversion factor for Raman cross sections from Ry atomic units
|
||||
! to A^4/amu
|
||||
!
|
||||
r1fac = 0.529177**4/2.d0
|
||||
!
|
||||
! derivatives of epsilon are translated into derivatives of molecular
|
||||
! polarizabilities by assuming a Clausius-Mossotti behavior
|
||||
! (for anisotropic system epsilon is replaced by its trace)
|
||||
!
|
||||
r2fac = 3.d0*omega/(4.d0*3.1415926) * &
|
||||
3.d0 / ( 2.d0 + (eps0(1,1) + eps0(2,2) + eps0(3,3))/3.d0 )
|
||||
!
|
||||
call mallocate(infrared,3*nat)
|
||||
call mallocate(raman,3,3,3*nat)
|
||||
!
|
||||
noraman=.true.
|
||||
do nu = 1,3*nat
|
||||
do ipol=1,3
|
||||
polar(ipol)=0.0
|
||||
end do
|
||||
do na=1,nat
|
||||
do ipol=1,3
|
||||
do jpol=1,3
|
||||
polar(ipol) = polar(ipol) + &
|
||||
zstar(ipol,jpol,na)*z((na-1)*3+jpol,nu)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
!
|
||||
! infrared is in units of e^2/(Ry mass unit)
|
||||
!
|
||||
infrared(nu) = polar(1)**2+polar(2)**2+polar(3)**2
|
||||
!
|
||||
do ipol=1,3
|
||||
do jpol=1,3
|
||||
raman(ipol,jpol,nu)=0.0
|
||||
do na=1,nat
|
||||
do lpol=1,3
|
||||
raman(ipol,jpol,nu) = raman(ipol,jpol,nu) + r2fac * &
|
||||
deps_dtau(ipol,jpol,lpol,na) * z((na-1)*3+lpol,nu)
|
||||
end do
|
||||
end do
|
||||
noraman=noraman .and. abs(raman(ipol,jpol,nu)).lt.1.d-12
|
||||
end do
|
||||
end do
|
||||
! Raman cross sections are in units of bohr^4/(Ry mass unit)
|
||||
end do
|
||||
!
|
||||
write (6,'(/5x,''IR cross sections are in (D/A)^2/amu units'')')
|
||||
if (noraman) then
|
||||
write (6,'(/''# mode [cm-1] [THz] IR'')')
|
||||
else
|
||||
write (6,'(5x,''Raman cross sections are in A^4/amu units'')')
|
||||
write (6,'(/''# mode [cm-1] [THz] IR*10^-3 Raman'')')
|
||||
end if
|
||||
!
|
||||
do nu = 1,3*nat
|
||||
!
|
||||
u2=0.d0
|
||||
do na=1,nat
|
||||
do ipol=1,3
|
||||
u2 = u2 + abs(z((na-1)*3+ipol,nu))**2
|
||||
end do
|
||||
end do
|
||||
!
|
||||
freq = sqrt(abs(w2(nu)))*rydcm1
|
||||
if (w2(nu).lt.0.0) freq = -freq
|
||||
!
|
||||
! alpha, beta2: see PRB 54, 7830 (1996) and refs quoted therein
|
||||
!
|
||||
if (noraman) then
|
||||
write (6,'(i5,f10.2,f12.4,2f10.4)') &
|
||||
nu, freq, freq*cm1thz, infrared(nu)*irfac/u2
|
||||
else
|
||||
alpha = (raman(1,1,nu) + raman(2,2,nu) + raman(3,3,nu))/3.d0
|
||||
beta2 = ( (raman(1,1,nu) - raman(2,2,nu))**2 + &
|
||||
(raman(1,1,nu) - raman(3,3,nu))**2 + &
|
||||
(raman(2,2,nu) - raman(3,3,nu))**2 + 6.d0 * &
|
||||
(raman(1,2,nu)**2 + raman(1,3,nu)**2 + raman(2,3,nu)**2) )/2.d0
|
||||
write (6,'(i5,f10.2,f12.4,2f10.4)') &
|
||||
nu, freq, freq*cm1thz, infrared(nu)*irfac/u2, &
|
||||
(45.d0*alpha**2 + 7.0d0*beta2)*r1fac
|
||||
end if
|
||||
end do
|
||||
!
|
||||
return
|
||||
!
|
||||
end subroutine RamanIR
|
||||
Loading…
Reference in New Issue