Miedema et al.’s Enthalpy code — 25 years after..
July 28, 2008 Posted by 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!
Sub-groups of space groups
July 3, 2008 Posted by Emre S. Tasci
ISOTROPY is a highly effective program written by the guru of applications of Group Theory on Phase Transitions, Harold T. Stokes (For those who are interested in the topic, I recommend his Introduction to Isotropy Subgroups and Displacive Phase Transitions titled paper(co-authored by Dorian M. Hatch) and Structures and phase transitions in perovskites – a group-theoretical approach paper (co-authored by Christopher J. Howard).
I’m pretty new and noobie in Group Theory – it was one of the subjects I found myself always running away. Anyway, I wanted to limit my search for some specific phase transition search in the database. Say, if an A-B binary has been reported in S1 and S2 structures, I didn’t want to go all the way looking for possible transition mechanisms if S1->S2 isn’t allowed by the group theory at all!
In the meantime, I’ve begun experimenting with the ISOTROPY program, yet I’m still at the very beginning. You can find detailed information and an online version in the related website, but to give an example, suppose we’d like to know about the subgroups of spacegroup 123 (P4/mmm). We enter the following commands:
VALUE PARENT 123
SHOW PARENT
SHOW SUBGROUP
SHOW LATTICE
DISPLAY ISOTROPY
and it displays a list of which a portion is presented below:
Parent Lattice Subgroup Lattice
123 P4/mmm q 123 P4/mmm q
123 P4/mmm q 47 Pmmm o
123 P4/mmm q 83 P4/m q
123 P4/mmm q 65 Cmmm o-b
123 P4/mmm q 10 P2/m m
123 P4/mmm q 12 C2/m m-b
123 P4/mmm q 2 P-1 t
123 P4/mmm q 89 P422 q
123 P4/mmm q 111 P-42m q
123 P4/mmm q 99 P4mm q
123 P4/mmm q 115 P-4m2 q
123 P4/mmm q 25 Pmm2 o
123 P4/mmm q 38 Amm2 o-b
123 P4/mmm q 6 Pm m
123 P4/mmm q 123 P4/mmm q
123 P4/mmm q 127 P4/mbm q
123 P4/mmm q 127 P4/mbm q
You may have noticed that, there are some recurring lines – this is due to some other properties that we have not specifically selected for display, i.e., it would display lines only containing "123 P4/mmm" if we hadn’t specifically asked for subgroup and lattice information (in analogy with the degeneracies in QM).
I needed a table containing all the subgroups of each space-group but there were two problems:
1. As far as I know, ISOTROPY does not accept an input file where you can pass the commands you’ve prepared via a script/macro. So, I was afraid that I would type "VALUE PARENT 1 – DISPLAY ISOTROPY – VALUE PARENT 2 – DISPLAY ISOTROPY – …" (by the way, ISOTROPY has a clever interpreter so that you can just type the first letters of commands and parameters as long as they are not dubious).
2. Again, as far as I know, you can’t direct the output to a file, meaning you have to copy the output manually and paste into a file (x230 in my case).
Luckily, it turned out that:
1. You could not supply an input file, but you could instead, paste the contents of such an input file and ISOTROPY would process it line by line.
2. Stokes had included a very useful and helpful property that all the sessions was recorded automatically in the "iso.log" file.
So, all I had to do was
1. Write a script that would produce the query to retrieve the subgroups of all 230 space-groups.
2. Paste the outcome of this script to ISOTROPY.
3. Write another script that would parse the "iso.log" file and split these reported subgroups per spacegroup (and while I’m at it, remove the duplicate lines and order the subgroups).
Maybe you won’t believe me but I am also not happy to paste lines and lines of code to this blog’s entries. I will deal it soon but in the mean time, please bear with them a little longer.. Ath the moment, I will describe the process and will include them at the end.
* build the query via query_builder.php
* feed this query to ISOTROPY
* split the space groups from the "iso.log" via split_spacegroups.php
* run sortall.php which :
* runs the sortkillduplines.php on each of these files
* sorts again via the sort.pl wrt to the 4th col (space number of the subgroup)
* removes the blank lines
* And that’s it..
Here comes the codes (php+perl)
query_builder.php:
<?
// Generates the query that will be "pasted" into the ISOTROPY prompt
// The query will be stored in isotropy_query.txt
$fp = fopen("isotropy_query.txt","w");
fwrite($fp,"PAGE 1000\nSHOW SUBGROUP\nSHOW LATTICE\nSHOW PARENT\n");
for($i=1;$i<231;$i++)
fwrite($fp,"VALUE PARENT $i\nDISPLAY ISOTROPY\n");
fwrite($fp,"\n\n");
fclose($fp);
?>
sortkillduplines:
<?
// PERL sorts the mentioned file (via $_GET["file"]) and deletes duplicate lines
require("commandline.inc.php");
if($_GET["file"]) $file = $_GET["file"];
else
{
$in = fopen(âphp://stdinâ, ârâ);
$fw = fopen("sortkilldup_tmp.tmp","w");
while(!feof($in))
{
fwrite($fw, rtrim(fgets($in, 4096))."\n");
}
fclose($fw);
$file = "sortkilldup_tmp.tmp";
}
if($_GET["seperator"])$f_seperator = TRUE;//Means that an extra seperator has been placed for this order purpose and is seperated from the text by the parameter defined with seperator=xxx.
if($_GET["n"]||$_GET["numeric"]||$_GET["num"])$f_numeric = TRUE; // If it is TRUE than it is assumed that the first word of each line is a number and will be sort accordingly.
// Reference : http://www.stonehenge.com/merlyn/UnixReview/col06.html
if($f_numeric) $exec = exec("perl -e âprint sort numerically <>;\nsub numerically { \$a <=> \$b; }â ".$file." > sorttemp.txt");
else $exec = exec("perl -e âprint sort <>â ".$file." > sorttemp.txt");
$fi = file("sorttemp.txt");
unlink("sorttemp.txt");
$curlinecont = "wrwerwer";
if(!$f_seperator)
{
for($i=0; $i<sizeof($fi); $i++)
{
$line = rtrim($fi[$i]);
if($curlinecont!=$line)
{
echo $line."\n";
$curlinecont = $line;
}
}
}
else
{
$frs = fopen("sorttemp2.txt","w"); // will capture the dup killed output in this file and resort it later.
for($i=0; $i<sizeof($fi); $i++)
{
$lineorg = rtrim($fi[$i]);
$linearr = explode($_GET["seperator"],$lineorg);
if($curlinecont!=$linearr[0])
{
fwrite($frs, $linearr[1]."\n");
$curlinecont = $linearr[0];
}
}
fclose($frs);
// resort the file
passthru("perl -e âprint sort <>â sorttemp2.txt");
unlink("sorttemp2.txt");
}
// unlink("sortkilldup_tmp.tmp");
?>
split_spacegroups.php:
<?
// splits the isotropy.log file generated by the query :
// PAGE 1000
// SHOW SUBGROUP
// SHOW LATTICE
// SHOW PARENT
// VALUE PARENT 1
// DISPLAY ISOTROPY
// â¦
// VALUE PARENT 230
// DISPLAY ISOTROPY
// with respect to the parent groups
$file = fopen("iso.log", "r") or exit("Unable to open file!");
//Output a line of the file until the end is reached
$cur = 0;
$fp = fopen("dummy.txt","w"); // we need to open a dummy file in order to include the fclose() in the iteration
while(!feof($file))
{
// Read the line
$line = fgets($file);
// Check if it is a spacegroup information entry by looking at the first word
// - if this word is a number, than process that line
if(ereg("^([0-9]+)",$line,$num))
{
$num = $num[0];
if($cur!=$num)
{
// We have a new space group, so create its file:
fclose($fp); // First close the previous space groupâs file
$filename = sprintf("%03d",$num)."_sg.txt";
$fp = fopen($filename,"w");
$cur = $num; // Set the current check data to the sg number
}
fwrite($fp,$line);
}
}
fclose($file);
unlink("dummy.txt");
?>
sortall.php:
<?
// sorts & kills dupes in all the spacegroup files
for($num=1;$num<231;$num++)
{
$filename = sprintf("%03d",$num)."_sg.txt";
$exec = exec("php ../../toolbox/sortkillduplines.php file=".$filename." > temp");
$exec = exec("./sort.pl temp > temp2");
$exec = exec("perl -pi -e \"s/^\\n//\" < temp2 > ".$filename);
}
unlink("temp");
unlink("temp2");
?>
or you can directly skip all these steps and download the generated spacegroup files from here. 8)
Lattice values are given in Schoenflies Notation where it corresponds to (with Pearson notation):
T : triclinic (AP)
M : prmitive monoclinic (MP)
M-B : base-centered monoclinic (MC)
O : primitive orthorhombic (OP)
O-B : base-centered orthorhombic (OC)
O-V : body-centered orthorhombic (OI)
O-F : face-centered orthorhombic (OF)
Q : primitive tetragonal (TP)
Q-V : body-centered tetragonal (TI)
RH : trigonal (HR)
H : hexagonal (HP)
C : primitive cubic (CP)
C-F : face-centered cubic (CF)
C-V : body-centered cubic (CI)
Quick reference on Bravais Lattices and info on how/why they are called after Bravais but not Frankenheim..
