# Rcjp's Weblog

## April 30, 1993

### BSc Fortran Numerics Course – Phonon Dispersion Curves

Filed under: physics — rcjp @ 11:39 am

Amongst the programs were those to calculate phonon dispersion curves and had a note…

“The frequency of oscillation of an atom in a crystal can be calculated from an equation involving a dynamical matrix which incorportes the mass of the atom and the forces on it.

The frequency of oscillation in three crystal directions ,  and  are calculated for iron. The program includes data such as the atomic weights and elastic bulk modulus and requests input of the number of atoms in the unit cell which is 2 for BCC and 4 for FCC.”


CCCC*******************************************************************
CCCC   Calculation of phonon frequencies travelling through a
CCCC   iron crystal lattice, for   and  directions
CCCC*******************************************************************
C  NOTES: User input determines number of atoms per unit cell
C         BCC=2 FCC=4
C
C  VARIABLES:
C  RHO      - Density
C  K        - Bulk modulus
C  N        - Number of atoms per unit cell
C  KX,KY,KZ - Wave direction components
C  A        - Unit cell size
C  PIOA     - Pi over A
C  STEP     - Stepsize for increment to Kx
C  AT       - Atomic weight
C  GAMMAM   - Force constant divided by mass
C
C  SUBROUTINES:
C  F02ABF   -  NAG Routine calculates the eigenvalues and eigenvectors
C              of a symmetric matrix.
C  PH       -  Programs the dynamical matrix and calls F02ABF.
C
CCCC*******************************************************************
REAL*8 RHO,K,KX,KY,KZ,A,PIOA,STEP,AT,GAMMAM,R

CCCC Data for Iron

DATA AT,RHO,K,R/55.847,7860.0D0,170.0D9,0.8D0/

CCCC Read number of atoms in the unit cell

PRINT *,'ENTER THE NUMBER OF ATOMS IN UNIT CELL '

CCCC Calculate unit cell size

A      = (N*AT/(1000.0*6.023D23*RHO))**(1.0/3.0)
GAMMAM = K/(3.0D0*(4.0D0+4.0D0*R)*RHO*A*A)
PIOA   = 4.0d0*DATAN(1.0d0)/A
STEP   = PIOA/12.0D0

PRINT *,'      KX              FREQUENCIES'

DO 10 KX = 0.0D0, PIOA, STEP

CCCC Calculated wave frequencies for  propagation direction

KY = 0.0D0
KZ = 0.0D0
CALL PH(KX,KY,KZ,AT,GAMMAM,R,A)
10 CONTINUE

PRINT *,'      KX              FREQUENCIES'

DO 20 KX = 0.0D0, PIOA, STEP

CCCC Calculated wave frequencies for  propagation direction

KY = KX
KZ = 0.0D0
CALL PH(KX,KY,KZ,AT,GAMMAM,R,A)
20  CONTINUE

PRINT *,'      KX              FREQUENCIES'

DO 30 KX = 0.0D0, PIOA, STEP

CCCC Calculated wave frequencies for  propagation direction

KY = KX
KZ = KX
CALL PH(KX,KY,KZ,AT,GAMMAM,R,A)
30  CONTINUE
END

SUBROUTINE PH(KX,KY,KZ,AT,GAMMAM,R,A)
CCCC*******************************************************************
C    Finds the eigenvalues - the wave frequencies for a specified
C    wave direction
CCCC*******************************************************************
C VARIABLES:
C  A        - Unit cell size
C  AT       - Atomic Weight
C  GAMMAM   - Force constant divided by mass
C  R        - Force constant ratio
C  COEFF    - Coefficient
C  DELTA    - Constant
C  D(,)     - Dynamical Matrix
C  KX,KY,KZ - Wave direction components
C
C  ADDITIONAL VARIABLES FOR NAG ROUTINE:
C  IA        - Specifies the first dimension of D(,)
C  ORDER     - Specifies the order of the matrix
C  IFAIL     - ERROR INDICATOR
C  EIGVAL()  - Contains eigenvalues
C  EIGVEC(,) - Contains eigenvectors
C  IV        - First dimension of eigenvector array
C  E         - Work space
C
C  SUBROUTINES:
C  F02ABF   -  NAG Routine calculates the eigenvalues and eigenvectors
C              of a symmetric matrix.
CCCC*******************************************************************

EXTERNAL F02ABF
REAL*8 KX,KY,KZ,A,GAMMAM,D(3,3),EIGVAL(3),COEFF,
*                 EIGVEC(3,3),DELTA,R,AT,E(3)
INTEGER ORDER

CCCC  Set values for NAG routine

IA    = 3
ORDER = 3
IFAIL = 0
IV    = 3

CCCC Calculate Constants

COEFF  = 8.0D0/3.0D0
DELTA  = COEFF+2.0D0*R - COEFF*DCOS(KX*A)*DCOS(KY*A)*DCOS(KZ*A)

CCCC Calculated dynamical matrix elements

D(1,1) = DELTA-2.0D0*R*DCOS(2.0*KX*A)
D(2,1) = COEFF*DSIN(KX*A)*DSIN(KY*A)*DCOS(KZ*A)
D(3,1) = COEFF*DSIN(KX*A)*DCOS(KY*A)*DSIN(KZ*A)
D(1,2) = D(2,1)
D(2,2) = DELTA-2.0D0*R*DCOS(2.0*KY*A)
D(3,2) = COEFF*DCOS(KX*A)*DSIN(KY*A)*DSIN(KZ*A)
D(1,3) = D(3,1)
D(2,3) = D(3,2)
D(3,3) = DELTA-2.0D0*R*DCOS(2.0*KZ*A)

CALL F02ABF(D,IA,ORDER,EIGVAL,EIGVEC,IV,E,IFAIL)

IF (IFAIL.NE.0) THEN
PRINT *,'ERROR IN NAG ROUTINE IFAIL = ',IFAIL
STOP
ENDIF

CCCC Multiply eigenvalues by force constant divided by mass

PRINT  '(1X,4(D12.5,1X))',KX,(EIGVAL(I)*GAMMAM,I=1,3)
END


I think the programs were part of some coursework I did during my BSc in Computational Physics at Salford University.

It must have been after my stint at BNFL because I was using raw TeX as I thought Latex looked too babyish. I’d put together some basic macros…


\\font\\ninerm=cmr9
\\font\\twelverm=cmr12
\\font\\seventeenrm=cmr17

\\font\\smallcaps=cmcsc10
\\font\\bigtitle=cmss12
\\font\\part=cmbx12
\\parindent 0mm
\\parskip 1em
\\baselineskip=1.2em

\\def\\E{{\\bf E}}
\\def\\B{{\\bf B}}
\\def\\H{{\\bf H}}
\\def\\P{{\\bf P}}
\\def\\Pl{{\\bf P}_L}
\\def\\Pnl{{\\bf P}_{NL}}
\\def\\D{{\\bf D}}
\\def\\r{{\\bf r}}
\\def\\del{\\nabla}
\\def\\intf{\\intop^{+\\infty}_{-\\infty}}
\\def\\chione{\\chi^{(1)}}
\\def\\chithree{\\chi^{(3)}}
\\def\\muo{\\mu_0}
\\def\\ko{k_0}
\\def\\wo{w_0}
\\def\\betao{\\beta_0}
\\def\\epsilono{\\epsilon_0}
\\def\\Eh{\\tilde E}
\\def\\Ah{\\tilde A}
\\def\\Eb{{\\overline E}}
\\def\\nb{{\\overline n}}
\\def\\betab{{\\overline\\beta}}
\\def\\pone#1#2{{\\partial#1\\over\\partial#2}}
\\def\\ptwo#1#2{{\\partial^2 #1\\over\\partial #2^2}}
\\def\\Tos{T_0^2}
\\def\\sqp{\\sqrt{P_0}}
\\def\\fig#1{{\\bf Figure #1. }}

%\\def\\section#1{{\\vskip 5mm\\seventeenrm #1\\vskip 1em}}
\\def\\section#1{\\vfill\\eject{\\seventeenrm #1\\vskip 1em}}
\\def\\subsection#1{\\vskip 5mm{\\bigtitle #1\\vskip 1em}}
\\def\\ve{\\vfill\\eject}
\\def\\fortran{{\\tt FORTRAN }}
\\def\\pagenumbers{\\footline={\\hss\\tenrm\\folio\\hss}}

\\def\\boxit#1{\\vbox{\\hrule\\hbox{\\vrule\\kern3pt
\\vbox{\\kern3pt#1\\kern3pt}\\kern3pt\\vrule}\\hrule}}

\\def\\uncatcodespecials{\\def\\do##1{\\catcode##1=12 }\\dospecials}

\\newcount\\lineno % the number of file lines

{\\catcode\\=\\active \\gdef{\\relax\\lq}}

\\def\\myupverbatim{\\tt \\lineno=0
\\def\\par{\\leavevmode\\endgraf}  \\catcode\\=\\active
\\obeylines \\uncatcodespecials \\obeyspaces

\\def\\setupverbatim{\\tt \\lineno=0
\\def\\par{\\leavevmode\\endgraf}  \\catcode\\=\\active
\\obeylines \\uncatcodespecials \\obeyspaces
{\\obeyspaces\\global\\let =\\ } % let active space = control space

\\def\\listing#1{\\par\\parskip -0.3em\\begingroup\\setupverbatim\\input#1 \\endgroup
\\parskip 1em}
\\def\\results#1{\\par\\parskip -0.3em\\begingroup\\myupverbatim\\input#1 \\endgroup
\\parskip 1em}


and I was using the fabulous GLE (Graphics Language Editor) by Chris Pugmire for my graphs (sigh, I still dont think anything looks as good, it had a cool trick to use TeX to do the axes captions). A typical gle input file was


size 21 29
amove 3 14
begin scale 0.6 0.6
begin graph
nobox
size 15 15
title "Current in a Diode Circuit" font texcmr
xtitle "Voltage /V" font texcmtt
ytitle "Current /Amps" font texcmtt
data vidata.dat
d1 lstyle 1 smooth
end graph
rmove 0 -14
begin graph
nobox
size 15 15
title "Solution Times for NAG and Bisection" font texcmr
xtitle "Voltage /V" font texcmtt
ytitle "Time for Solution /Secs" font texcmtt
data tdata.dat
yaxis max 1e-2
d1 lstyle 1 smooth key NAG
d2 lstyle 2 smooth key Bisection
end graph
end scale