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 [100], [110] and [111] 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 [100] [110] and [111] 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 '
          READ *, N

    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          [100]    FREQUENCIES'

          DO 10 KX = 0.0D0, PIOA, STEP

    CCCC Calculated wave frequencies for [100] propagation direction

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

          PRINT *,'      KX          [110]    FREQUENCIES'

          DO 20 KX = 0.0D0, PIOA, STEP

    CCCC Calculated wave frequencies for [110] propagation direction

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

          PRINT *,'      KX          [111]    FREQUENCIES'

          DO 30 KX = 0.0D0, PIOA, STEP

    CCCC Calculated wave frequencies for [111] 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
    \\everypar{\\advance\\lineno by1 \\llap{\\fiverm\\ \\ }}}

    \\def\\setupverbatim{\\tt \\lineno=0
      \\def\\par{\\leavevmode\\endgraf}  \\catcode`\\`=\\active
        \\obeylines \\uncatcodespecials \\obeyspaces
        \\everypar{\\advance\\lineno by1 \\llap{\\fiverm\\the\\lineno\\ \\ }}}
{\\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
Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Create a free website or blog at WordPress.com.

%d bloggers like this: