Hex, Bugs and More Physics | Emre S. Tasci

Yes, it’s been 25 years since A.K. Niessen, A.R. Miedema et al.’s "Model Predictions for the Enthalpy of Formation of Transition Metal Alloys II" titled article (CALPHAD 7(1) 51-70 1983) was published. It can be thought as a sequel to Miedema et al.’s 1979 (Calphad 1 341-359 1979) and 1980 dated (Physica 100 1-28 1980) papers with some update on the data as well as a computer code to calculate the enthalpies of formation written in ALGOL.

I have been currently working on these papers and by the way ported the code presented in CALPHAD 1983 to FORTRAN and here it is:

C      THE IMPLEMENTATION OF THE ALGOL CODE FROM
C      A.K. Niessen, F.R. de Boer, R. Boom, P.F. de Chatel
C      W.C.M. Mattens and A.R. Miedema
C      CALPHAD Vol.7 No.1 pp. 51-70, 1983
C      Model Predictions for the Enthalpy of Formation of Transition 
C        Metal Alloys II"
C
C      by Emre S. Tasci, TUDelft (2008)

      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
C      DOUBLE PRECISION va,aa,p,r,x,dn,dg,ng,pq,xa,xm,csa,csm
C      DOUBLE PRECISION cas,cms,fam,dph,vm,am,fma,var,vmr
C      DOUBLE PRECISION vas, vms, xas, xms, xva,xvm, dgl,pqrs, pqrl
      INTEGER arrout, z1,z2,a,m,n,w,dh,ga,gm,gal,gml,nr,dHtrans
      LOGICAL detp, detr, bool
      CHARACTER(LEN=5) Symbol
      DIMENSION arrout(15)
      DIMENSION Phi(75), acf(75), dHtrans(75), Nws113(75), Rcf(75),
     + V213(75), xxx(9), Symbol(75)
      COMMON /CDH/ x, dh
      COMMON /CDPH/ dph, Phi, Nws113, ng, r, Rcf, p, pq, dn, pqrs, pqrl

      OPEN(UNIT=3,FILE='inpdata.txt',STATUS='OLD')
      DO I=1,75
        READ(3,*)nr,Symbol(nr),Phi(nr),Nws113(nr),V213(nr),
     +  acf(nr),Rcf(nr),dHtrans(nr)
C       WRITE(*,*)nr,Symbol(nr),Phi(nr),Nws113(nr),V213(nr),
C     + acf(nr),Rcf(nr),dHtrans(nr)
      END DO
      CLOSE(UNIT=3)

      xxx(1) = 0.001
      xxx(2) = 1.0/6.0
      xxx(3) = 1.0/4.0
      xxx(4) = 1.0/3.0
      xxx(5) = 1.0/2.0
      xxx(6) = 2.0/3.0
      xxx(7) = 3.0/4.0
      xxx(8) = 5.0/6.0
      xxx(9) = 0.999

      DO z1 = 1,46 
          WRITE(*,200)z1,Symbol(z1),Phi(z1),Nws113(z1)**3,
     +      V213(z1)**(3.0/2),dHtrans(z1)
200       FORMAT(I2," ",A5," Phi: ",F5.2,"V  Nws: ",F5.2,
     +     "d.u.  Vmole: "
     +     ,F5.2,"cm3  DeltaHtrans:",I3,"kJ/mole",/)
          WRITE(*,250)"M","AM5","AM3","AM2","AM",
     +     "MA2","MA3","MA5","AinM","AM","MinA"
250       FORMAT(A,"      ",9(A4,"  "),A4,/)
          DO z2 = 1, 75
             m = z2
             a = z1
             CALL PQRSL(a,m)

             va = V213(a)
             aa = acf(a)
             ga = dHtrans(a)
             gal = ga

             vm = V213(m)
             am = acf(m)
             gm = dHtrans(m)

             IF (m.EQ.67.OR.m.EQ.68) THEN
                 gml = 0
             ELSE
                 gml = gm
             END IF

             n = 1

             DO I = 1,9
              xa = xxx(I)
              IF(xa.EQ.0.001.OR.xa.EQ.0.999) THEN
                  bool = .TRUE.
              ELSE
                  bool = .FALSE.
              END IF

              xm = 1 - xa
              xva = xa * va
              xvm = xm * vm

              dgl = xa * gal + xm * gml
              dg = xa * ga + xm * gm
              csa = xva / (xva + xvm)

              csm = 1 - csa
              fam = csm * (1 + 8 * csa * csa * csm * csm)
              fma = fam * csa / csm

              var = va * (1 + aa * csm * dph)
              vmr = vm * (1 - am * csa * dph)

              vas = va * (1 + aa * fam * dph)
              vms = vm * (1 - am * fma * dph)

              xar = xa * var
              xas = xa * vas
              xmr = xm * vmr
              xms = xm * vms

              cas = xas / (xas+xms)
              cms = 1 - cas
              csa = xar / (xar + xmr)
              csm = 1 - csa

              fam = cms * (1 + 8 * cas * cas * cms * cms)
              fma = fam * cas / cms

              IF(bool) THEN
                  IF(n.EQ.1) THEN 
                      GOTO 292
                  ELSE
                      GOTO 392
                  END IF
              END IF
                  IF(xa.EQ.0.5) THEN
                      x = csm * xar * p * pqrl + dgl
                      CALL MAXI
                      arrout(n+5) = dh
                  END IF
                  x = fam * xas * p * pqrs + dg
                  CALL MAXI
                  arrout(n) = dh
                  arrout(11) = gml
                  GOTO 492
292               x = (csm * xar * p * pqrl + dgl) / xa
                  CALL MAXI
                  arrout(n) = dh
                  arrout(11) = gml
                  GOTO 492
392               x = (csm * xar * p * pqrl + dgl) / xm
                  CALL MAXI
                  arrout(n) = dh
                  arrout(12) = gal
492               n = n + 1
             END DO
C            DO J=1,15
C              WRITE(*,500) J, arrout(J)
C            END DO
             WRITE(*,600)Symbol(m),arrout(2),arrout(3),arrout(4),
     +        arrout(5),
     +        arrout(6),arrout(7),arrout(8),arrout(1),arrout(10),
     +        arrout(9)
             IF(Symbol(m).EQ."Ni".OR.Symbol(m).EQ."Pd".OR.Symbol(m).EQ.
     +       "Lu".OR.Symbol(m).EQ."Pt".OR.Symbol(m).EQ."Pu".OR.
     +       Symbol(m).EQ."Au".OR.Symbol(m).EQ."H".OR.Symbol(m).EQ.
     +       "Cs".OR.Symbol(m).EQ."Mg".OR.Symbol(m).EQ."Ba".OR.
     +       Symbol(m).EQ."Hg".OR.Symbol(m).EQ."Tl".OR.Symbol(m).EQ.
     +       "Pb".OR.Symbol(m).EQ."Bi")WRITE(*,*)""
          END DO
      WRITE(*,"(2A,/)")"-----------------------------------",
     + "----------------------------"
      END DO

C500   FORMAT("O[",I3,"] = ",I50)
600   FORMAT(A5,"  ",9(I4,"  "),I4)

      STOP
      END

      SUBROUTINE MAXI
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      INTEGER dh
      COMMON /CDH/ x, dh
      IF(x.LT.-999.OR.x.GT.999) THEN
        dh = 999
      ELSE
        dh = NINT(x)
      END IF
C      WRITE(*,*)"x:",dh
      RETURN
      END

      SUBROUTINE PQRSL(a,b)
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      INTEGER a,b
      LOGICAL detp,detr
      DIMENSION Phi(75), acf(75), dHtrans(75), Nws113(75), Rcf(75),
     + V213(75)
      COMMON /CDPH/ dph, Phi, Nws113, ng, r, Rcf, p, pq, dn, pqrs, pqrl

          dph = Phi(a) - Phi(b)
          dn = Nws113(a)-Nws113(b)
          ng = (1/Nws113(a) + 1/Nws113(b))/2

          IF (detr(a).EQV.detr(b)) THEN
              r = 0
          ELSE
              r = Rcf(a)*Rcf(b)
          END IF

          IF (detp(a).AND.detp(b)) THEN
              p = 1.15 * 12.35
C             BOTH TRUE
          ELSE IF (detp(a).EQV.detp(b)) THEN
              p = 12.35 / 1.15
C             BOTH FALSE
          ELSE
              p = 12.35
          END IF
          p = p / ng
          pq = -dph * dph + 9.4 * dn * dn
          pqrs = pq - r
          pqrl = pq - 0.73 * r

      RETURN
      END

      LOGICAL FUNCTION  detp(z)
      INTEGER z
      IF (z.LT.47.AND.(z.NE.23.AND.z.NE.31)) THEN
          detp = .true.
      ELSE
          detp = .false.
      END IF
      RETURN
      END

      LOGICAL FUNCTION  detr(z)
      INTEGER z
      IF (z.LT.47.OR.(z.GT.54.AND.z.LT.58)) THEN
          detr = .true.
      ELSE
          detr = .false.
      END IF
      RETURN
      END

Even though ALGOL doesn’t exist anymore, being a highly semantic language, it is still used in forms of pseudo-code to present algorithms and because of this reason, it didn’t take much to implement the code in FORTRAN. You can download the datafile from here. And also if you don’t want to run the code, here is the result file. And, as you can see, I haven’t trimmed or commented the code since I’ll feed it the updated values and move on forward, sorry for that!