Hex, Bugs and More Physics | Emre S. Tasci

Converts POSCAR files to CIF format. Uses Octave (latvec2par) to convert the unit cell vectors to lattice cell parameters. It is assumed (compulsory) that the atom coordinates in the POSCAR file are in fractional (direct) coordinates.

Parameters

f         : input file  (Default = POSCAR)
o         : output file (Default = <inputfilename>.cif)
t         : if set (t=1), then the output is written to the screen
speclist  : Define the name of the species, seperated by 'xx' or ',' to be used
            with neighbour listing (Default = AxxBxxCxx...)
            (Labels can also be seperated via [space] as long as speclist is given
            in quotation -- Example: ... speclist=\"Au Si\" )

Example

sururi@husniya OUTCARs_final_structures $ cat POSCAR_Pd3S
Pd3S
   1.00000000000000
     4.7454403619558345    0.0098182468538853    0.0000000000000000
    -1.4895902658503473    4.5055989020479856    0.0000000000000000
     0.0000000000000000    0.0000000000000000    7.2191545483192190
   6   2
Direct
  0.1330548697855782  0.5102695954022698  0.2500000000000000
  0.4897304045977232  0.8669451452144267  0.7500000000000000
  0.5128136657304309  0.1302873334247993  0.2500000000000000
  0.8697126665752007  0.4871863042695738  0.7500000000000000
  0.0013210250693640  0.9986789749306360  0.0000000000000000
  0.0013210250693640  0.9986789749306360  0.5000000000000000
  0.6856021298170862  0.6825558526447357  0.2500000000000000
  0.3174442073552622  0.3143978111829124  0.7500000000000000

sururi@husniya OUTCARs_final_structures $ php POSCAR2CIF.php f=POSCAR_Pd3S t=1 speclist=Pd,S
data_
loop_
_symmetry_equiv_pos_as_xyz
x,y,z
_cell_length_a  4.745451
_cell_length_b  4.745451
_cell_length_c  7.219155
_cell_angle_alpha       90.000000
_cell_angle_beta        90.000000
_cell_angle_gamma       108.175792
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
Pd01    Pd        0.1330548697855782  0.5102695954022698  0.2500000000000000
Pd02    Pd        0.4897304045977232  0.8669451452144267  0.7500000000000000
Pd03    Pd        0.5128136657304309  0.1302873334247993  0.2500000000000000
Pd04    Pd        0.8697126665752007  0.4871863042695738  0.7500000000000000
Pd05    Pd        0.0013210250693640  0.9986789749306360  0.0000000000000000
Pd06    Pd        0.0013210250693640  0.9986789749306360  0.5000000000000000
S01     S         0.6856021298170862  0.6825558526447357  0.2500000000000000
S02     S         0.3174442073552622  0.3143978111829124  0.7500000000000000

Code

#!/usr/bin/php
<?PHP
/* Emre S. Tasci <e.tasci@tudelft.nl>                    *
 * Script name: POSCAR2CIF.php                           *
 * Converts the specified POSCAR file to CIF format      *
 * Uses octave to convert the unit cell vectors          *
 * to lattice parameters.                                *
 *                                              23/05/09 */
 
 
//require("/code/toolbox/commandline.inc.php"); # Equivalent handling of $_GET / $_POST
# ============/code/toolbox/commandline.inc.php========================================
// Enables the same treat for command line parameters
// as those of GET parameters.
// Also sets the $first, $last ranges including the gall parameter
 
//Support for command line variable passing:
for($i=1;$i<sizeof($_SERVER["argv"]);$i++)
{
  list($var0,$val0) = explode("=",$_SERVER["argv"][$i]);
  $_GET[$var0] = $val0;
}
 
$first = -2; //disabled by default
$last  = -3;
 
if($_GET["gfirst"]) $first = $_GET["gfirst"];
if($_GET["glast"] ) $last  = $_GET["glast"];
if($_GET["gall"] ) {$first = $_GET["gall"]; $last = $first;}
# ============/code/toolbox/commandline.inc.php========================================
 
 
if($_GET["h"]||$_GET["help"])
{
$help = <<<lol
 * Script name: POSCAR2CIF.php                           *
 * Converts the specified POSCAR file to CIF format      *
 * Uses octave to convert the unit cell vectors          *
 * to lattice parameters.                                *
 
 * Emre S. Tasci <e.tasci@tudelft.nl>                    *
 *                                              23/05/09 *
 
 Usage:
 f         : input file  (Default = POSCAR)
 o         : output file (Default = <inputfilename>.cif)
 t         : if set (t=1), then the output is written to the screen
 speclist  : Define the name of the species, seperated by 'xx' or ',' to be used
             with neighbour listing (Default = AxxBxxCxx...)
             (Labels can also be seperated via [space] as long as speclist is given
             in quotation -- Example: ... speclist=\"Au Si\" )
 
 Example : php POSCAR2CIF.php f=POSCAR_RhBi4 speclist=Bi,Rh
lol;
 echo $help."\n";
 exit;
}
 
 
if($_GET["file"])$inputfile = $_GET["file"];
else if($_GET["f"])$inputfile = $_GET["f"];
else $inputfile="POSCAR";
 
$outfile = $inputfile.".cif";
if($_GET["o"]) $outfile=$_GET["o"];
if($_GET["t"]) $outfile = "php://stdout";
 
$species_type_template = Array("A","B","C","D","E","F","G","H","I","J","K","L","M","N","O","P","Q","R","S","T","U","V","W","X","Y","Z");
if($_GET["speclist"])$species_type_template = split("xx| |,",$_GET["speclist"]);
 
//echo sizeof($_GET["speclist"])."\n";
exec('sed -n "3,6p" '.$inputfile.'|awk \'{print $1"\t"$2"\t"$3}\'',$out);
$octfp = fopen("octrungetvolnpar.m","w");
 
fwrite($octfp,"kini=[\n");
for ($i=0;$i<3;$i++)
    fwrite($octfp,$out[$i]."\n");
fwrite($octfp,"];\n\n");
 
fwrite($octfp,"numatom = [".$out[3]."];\n\n");
$numatoms = split("[ \t]+",trim($out[3]));
 
// Write the octave operations
fwrite($octfp,"k1 = latvec2par (kini);totatom=sum(numatom );\n");
fwrite($octfp,"for i=1:6;printf(\"%f\\n\",k1(i));endfor;\nprintf(\"%d\\n\",totatom)\n");
 
fclose($octfp);
 
// Execute the octave file to find the lattice parameters and volume
exec('octave -q '.$thisdir.'octrungetvolnpar.m',$out3);
unlink("octrungetvolnpar.m");
//echo $out3[0]."\n";
 
/*
for ($i=0;$i<sizeof($out3);$i++)
    echo $i."\t".$out3[$i]."\n";
 */
 
exec('grep "Direct" '.$inputfile.' -n|sed "s:^\([0-9]\+\).*:\1:"',$start);# Line number of the Direct prompt
$start = $start[0];
$start++;
$finish = $start+$out3[sizeof($out3)-1]-1;
exec('sed -n "'.$start.','.$finish.'p" '.$inputfile,$outr);
$fpo = fopen($outfile, "w");
fwrite($fpo, "data_\nloop_\n_symmetry_equiv_pos_as_xyz\nx,y,z\n");
fwrite($fpo, "_cell_length_a\t".$out3[0]."\n");
//echo $out3[0];
fwrite($fpo, "_cell_length_b\t".$out3[1]."\n");
fwrite($fpo, "_cell_length_c\t".$out3[2]."\n");
fwrite($fpo, "_cell_angle_alpha\t".$out3[3]."\n");
fwrite($fpo, "_cell_angle_beta\t".$out3[4]."\n");
fwrite($fpo, "_cell_angle_gamma\t".$out3[5]."\n");
fwrite($fpo, "loop_\n_atom_site_label\n_atom_site_type_symbol\n_atom_site_fract_x\n_atom_site_fract_y\n_atom_site_fract_z\n");
$k = 0;
for($specie=0;$specie<sizeof($numatoms);$specie++)
    for($i=1;$i<=$numatoms[$specie];$i++)
    {
        fwrite($fpo, $species_type_template[$specie].sprintf("%02d",$i)."\t".$species_type_template[$specie]."\t".$outr[$k]."\n");
        $k++;
    }
 
 
fclose($fpo);
?>

A very simple code I wrote while studying the Python language that fetches the information from a VASP POSCAR file and prepares (and feeds) the input for Harold Stokes’ findsym (you should actually download the executable included in the ISOTROPY package instead of the referred web version) program to find the unit cell symmetry.

#!/usr/bin/python
# -*- coding: utf-8 -*-

"""
POSCAR2findsym.py
/* Emre S. Tasci <e.tasci@tudelft.nl>                    *
   A simple python code that parses the POSCAR file,
   prepares an input file for Harold Stokes' findsym
   (http://stokes.byu.edu/isotropy.html)
   executable.

   Then, checks the current system to see if findsym is
   accessible: if it is, then feeds the prepared file 
   and directly outputs the result; if findsym is not 
   executable then outputs the generated input file.

   If you are interested only in the results section of 
   the findsymm output, you should append "| grep "\-\-\-\-\-\-" -A1000"
   to the execution of this code.

   Accepts alternative POSCAR filename for argument
   (default = POSCAR)

   If a second argument is proposed then acts as if
   findsym is not accessible (ie, prints out the input
   file instead of feeding it to findsymm)

                                             10/03/09 */
"""

import sys
import re
from numpy import *
from os import popen
import commands

force_output = False

if len(sys.argv)<2:
    filename="POSCAR"
elif len(sys.argv)<3:
    filename=sys.argv[1]
else:
    filename=sys.argv[1]
    force_output = True

f = open(filename,'r')

title = f.readline().strip()
tolerance = 0.00001
latt_vec_type = 1 # We will be specifying in vectors

f.readline() # Read Dummy
lat_vec = array([])
for i in range(0,3):
    lat_vec = append(lat_vec,f.readline().strip())

# Read atom species
species_count = f.readline()
species_count = re.split('[\ ]+', species_count.strip())
species_count = map(int,species_count)
number_of_species = len(species_count)
total_atoms = sum(species_count)

initial_guess = "P" # For the centering

species_designation = array([])
for i in range(1,number_of_species+1):
    for j in range(0,species_count[i-1]):
        species_designation = append(species_designation,i)

species_designation = map(str,map(int,species_designation))
# print species_designation
sep = " "
# print sep.join(species_designation)

# Read Coordinates
f.readline() # Read Dummy
pos_vec = array([])
for i in range(0,total_atoms):
    pos_vec = append(pos_vec,f.readline().strip())

# print pos_vec

findsym_input = array([])
findsym_input = append(findsym_input, title)
findsym_input = append(findsym_input, tolerance)
findsym_input = append(findsym_input, latt_vec_type)
findsym_input = append(findsym_input, lat_vec)
findsym_input = append(findsym_input, 2)
findsym_input = append(findsym_input, initial_guess)
findsym_input = append(findsym_input, total_atoms)
findsym_input = append(findsym_input, species_designation)
findsym_input = append(findsym_input, pos_vec)
# print findsym_input
sep = "\n"
findsym_input_txt = sep.join(findsym_input)
# print findsym_input_txt

# Check if findsym is accessible:
status,output =  commands.getstatusoutput("echo \"\n\n\"|findsym ")
if(output.find("command not found")!=-1 or force_output):
    # if it is not there, then just output the input
    print findsym_input_txt
elif(status==6144):
    # feed it to findsym
    pipe = popen("echo \""+findsym_input_txt+"\" | findsym").readlines()
    for line in pipe:
        print line.strip()
quit()

Example Usage:

sururi@husniya POSCAR2findsym $ cat POSCAR
Si(S)-Au(Pd) -- Pd3S
   1.00000000000000
     4.7454403619558345    0.0098182468538853    0.0000000000000000
    -1.4895902658503473    4.5055989020479856    0.0000000000000000
     0.0000000000000000    0.0000000000000000    7.2191545483192190
   6   2
Direct
  0.1330548697855782  0.5102695954022698  0.2500000000000000
  0.4897304045977232  0.8669451452144267  0.7500000000000000
  0.5128136657304309  0.1302873334247993  0.2500000000000000
  0.8697126665752007  0.4871863042695738  0.7500000000000000
  0.0013210250693640  0.9986789749306360  0.0000000000000000
  0.0013210250693640  0.9986789749306360  0.5000000000000000
  0.6856021298170862  0.6825558526447357  0.2500000000000000
  0.3174442073552622  0.3143978111829124  0.7500000000000000

sururi@husniya POSCAR2findsym $ python POSCAR2findsym.py
FINDSYM, Version 3.2.3, August 2007
Written by Harold T. Stokes and Dorian M. Hatch
Brigham Young University
 
Si(S)-Au(Pd) -- Pd3S
Tolerance:    0.00001
Lattice vectors in cartesian coordinates:
4.74544   0.00982   0.00000
-1.48959   4.50560   0.00000
0.00000   0.00000   7.21915
Lattice parameters, a,b,c,alpha,beta,gamma:
4.74545   4.74545   7.21915  90.00000  90.00000 108.17579
Centering: P
Number of atoms in unit cell:
8
Type of each atom:
1  1  1  1  1  1  2  2
Position of each atom (dimensionless coordinates)
1   0.13305   0.51027   0.25000
2   0.48973   0.86695   0.75000
3   0.51281   0.13029   0.25000
4   0.86971   0.48719   0.75000
5   0.00132   0.99868   0.00000
6   0.00132   0.99868   0.50000
7   0.68560   0.68256   0.25000
8   0.31744   0.31440   0.75000
------------------------------------------
Space Group  40  C2v-16    Ama2
Origin at    0.00000   0.00000   0.50000
Vectors a,b,c:
0.00000   0.00000   1.00000
-1.00000  -1.00000   0.00000
1.00000  -1.00000   0.00000
Values of a,b,c,alpha,beta,gamma:
7.21915   5.56683   7.68685  90.00000  90.00000  90.00000
Atomic positions in terms of a,b,c:
Wyckoff position b, y = -0.17834, z =  0.31139
1   0.75000   0.17834   0.31139
2   0.25000   0.82166   0.31139
Wyckoff position b, y =  0.32155, z =  0.19126
3   0.75000   0.17845   0.69126
4   0.25000   0.82155   0.69126
Wyckoff position a, z =  0.00132
5   0.50000   0.00000   0.00132
6   0.00000   0.00000   0.00132
Wyckoff position b, y = -0.31592, z =  0.00152
7   0.75000   0.81592   0.50152
8   0.25000   0.18408   0.50152

After migrating this blog to a new server, I had some difficulties in writing up the equations. It was due to the fact that on this server, I wasn’t able to execute the necessary commands to parse the latex code and convert it to PNG format.

Last week, some need arouse for this and I finally coded the necessary way around to be able to convert LaTeX definitions to their PNG counterparts. Before I start to quote the necessary code and explanations, let me tell that, this method only will work if you can successfully manage to produce this conversation at your -local- computer since what the code does is to generate the image on your local computer and then copy it to the server where your blog (/or whatever you have) resides. So, make sure:

• You can convert LaTeX code to PNG on your -own- computer.
• You are running a httpd server so that the code can copy the output PNG file from your local computer

What the code does:
It has 3 seperate processes. Let’s designate them A, B & C.

• Part A is called on the remote server. It presents a form with textbox in it where you can type the LaTeX code. The form’s target is the same file but this time on your local server, so when you submit the data —
• Part B is called on your local server. Using the texvc code, an open source software that also comes packed with mediawiki, it converts the LaTeX code into a PNG image and places it under the relative math/ directory of your local server. After doing this, via the META tag ‘refresh’ it causes the page to refresh itself loading its correspondence in the remote server, hence initiating —
• Part C where the generated image identified by the MD5 hashname is fetched from the local server and copied unto the relative math/ directory on the remote server. The LaTeX code that was used is also passed in the GET method and it is included in the generated code as the title of the image, so if necessary, the equation can be modified in the future.

The code is as follows:

<?PHP if($_GET["com"]||(!$_POST["formula"] && !$_GET["file"])){?>
<body onload="document.postform.formula.focus();">
<form action="http://yourlocalserver.dns.name/tex/index.php" method="post" name=postform id=postform>
<textarea name="formula" cols=40 rows=5 id=formula>
</textarea><br>
<input type="submit" value="Submit" />
<?PHP
}
if($_POST["formula"])
{
    // NL
    $lol = './texvc  ./mathtemp ./math "'.$_POST["formula"].'" iso-8859-1';
    exec($lol,$out);
    $file = substr($out[0],1,32);
    $formula = base64_encode($_POST["formula"]);
    $html = '<HTML>
<HEAD>
<META http-equiv="refresh" content="0;URL=http://remoteserver.dns.name/tex/index.php?formula='.$formula.'&file='.$file.'">
</HEAD>';
    echo $html;
}
else if($_GET["file"])
{
    // COM
    $file = "math/".$_GET["file"].".png";
    if(!copy("yourlocalserver.dns.name/tex/".$file, $file)) exit( "copy failed<br>");
    $formula = base64_decode($_GET["formula"]);
?>
<form action="http://yourlocalserver.dns.name/tex/index.php" method="post">
<textarea name="formula" cols=40 rows=5><?PHP echo $formula; ?>
</textarea><br>
<input type="submit" value="Submit" />
<?
 
    echo "<br><br><img src=\"".$file."\" title=\"".$formula."\"><br>";
    echo "&lt;img src=\"/tex/".$file."\" title=\"".$formula."\"&gt;<br>";
}
?>

1. Create a directory called ‘tex’ on both servers.
2. Create subdir ‘math’ on both servers
3. Create subdir ‘mathtemp’ on the local server
4. Modify the permissions of these subdirs accordingly
5. Place the executable ‘texvc’ in the ‘tex’ directory
6. Place the code quoted above in the ‘tex’ directories of both servers in a file called ‘index.php’

If everything well up to this point, let’s give it a try and you shall have some output similar to this:

For an input LaTeX code of:
\sum_{n=0}^\infty\frac{x^n}{2n!}

where you can immediately see the resulting image, modify the code if there’s something wrong or copy the IMG tag and put it to your raw HTML input box in your text…

It’s of great importance to have sufficient number of k-points when doing DFT calculations. The most practical (but CPU time costly) way is to check energies for different number of k-points until some kind of "convergence" is reached.

Here are the two scripts I’ve recently written for VASP to check for suitable k-point values.

The first code runs the calculation with increasing k-points (I have a cubic cell, so using same number of k-points for the 3 lattice axes).

Before running this code, make sure that the essential files (KPOINTS POSCAR POTCAR INCAR) are in the directory, also copy your KPOINTS file to a file named "KP" which has KPOINTS mesh defined as 2x2x2:

sururi@husniya Al3Li $ cat KP
Automatic
0
Gamma
2 2 2
0 0 0

The following script will then run the calculations, incrementing the number of K-Points by 1 per each axis. It will store the output (and the 4 input files) in archive files in the format k{%03d}.tar.bz2 where the value in the curly brackets is actually the number of k-points per axis. As you can easily see, it runs from k= 1x1x1 (silly but since this is tutorial) to k=13x13x13.

#!/bin/bash

# A script file to run the system for various k-points in order to find the optimal setting
# Emre S. Tasci, 20/10/2008

for (( i = 1; i <= 13; i++ )) 
do  
    date
    echo $i

    # change the KPOINTS (KP is a copy of the KPOINTS with k=2:
    cat KP|sed s/"2 2 2"/"$i $i $i"/ > KPOINTS

    # run the vasp
    vasp > vasp.out

    # archive the existing files (except the archive files):
    tar -cvjf $(printf "k%02d.tar.bz2" $i) $(ls |grep -v "\(.bz2\)\|\(runrunrun.sh\)\|\(KP$\)\|\(work\)")

    # remove the output files:
    rm $(ls *|grep -v "\(.bz2\)\|\(KPOINTS\)\|\(POSCAR\)\|\(POTCAR\)\|\(INCAR\)\|\(runrunrun.sh\)\|\(KP$\)") -f
    date
    echo "=================================================="
done

After it ends, you should see something like

-rw-r----- 1 sururi users 299K 2008-10-20 10:14 POTCAR
-rw-r--r-- 1 sururi users  283 2008-10-20 10:14 POSCAR
-rw-r--r-- 1 sururi users  604 2008-10-20 10:14 INCAR
-rw-r--r-- 1 sururi users   31 2008-10-20 11:06 KP
-rw-r--r-- 1 sururi users   34 2008-10-20 14:33 KPOINTS
-rw-r--r-- 1 sururi users 195K 2008-10-20 11:08 k01.tar.bz2
-rw-r--r-- 1 sururi users 196K 2008-10-20 11:08 k02.tar.bz2
-rw-r--r-- 1 sururi users 193K 2008-10-20 11:08 k03.tar.bz2
-rw-r--r-- 1 sururi users 195K 2008-10-20 11:08 k04.tar.bz2
-rw-r--r-- 1 sururi users 195K 2008-10-20 11:08 k05.tar.bz2
-rw-r--r-- 1 sururi users 197K 2008-10-20 11:09 k06.tar.bz2
-rw-r--r-- 1 sururi users 197K 2008-10-20 11:09 k07.tar.bz2
-rw-r--r-- 1 sururi users 200K 2008-10-20 11:10 k08.tar.bz2
-rw-r--r-- 1 sururi users 201K 2008-10-20 11:10 k09.tar.bz2
-rw-r--r-- 1 sururi users 205K 2008-10-20 11:11 k10.tar.bz2
-rw-r--r-- 1 sururi users 205K 2008-10-20 11:12 k11.tar.bz2
-rw-r--r-- 1 sururi users 211K 2008-10-20 11:13 k12.tar.bz2
-rw-r--r-- 1 sururi users 211K 2008-10-20 11:14 k13.tar.bz2
-rwxr--r-- 1 sururi users  732 2008-10-20 14:02 runrunrun.sh

Now, it’s time for the second script. Create a new directory (say “work”) and then type the following into a script file (don’t forget to set it as executable (or you can also always “sh ” 😉 ) )

#!/bin/bash
 
# Script name: energy_extract_from_OUTCAR.sh
# Emre S. Tasci, 20/10/2008
 
# Reads the OUTCAR files from the archive files in the parent directory
# Stores the energies in corresponding ENE files.
 
# (Assuming filenames are given like "k05.tar.bz2")
 
k=00
 
for i in ../k*.bz2
do
    kpre=$k
    enepre="0.00000000"
 
    # Extract OUTCAR
    tar -xjf $i OUTCAR
 
    # Parse file counter from filename 
    k=$(echo $(basename $i)|sed "s:k\([0-9]*\).*:\1:")
 
    # write the energies calculated in this run 
    cat OUTCAR | grep "energy without"|awk '{print $8}' > $(printf "%s_ENE" $k)
 
    #calculate the energy difference between this and the last run
    if [ -e "$(printf "%s_ENE" $kpre)" ] 
    then
        enepre=$(tail -n1 $(printf "%s_ENE" $kpre));
    fi
 
    enethis=$(tail -n1 $(printf "%s_ENE" $k));
 
    # Using : awk '{ print ($1 >= 0) ? $1 : 0 - $1}' : for absolute value
 
    echo -e $(printf "%s_ENE" $kpre) " :\t" $enepre "\t" $(printf "%s_ENE" $k) ":\t" $enethis "\t" $(echo $(awk "BEGIN { print $enepre - $enethis}") | awk '{ print ($1 >= 0) ? $1 : 0 - $1}')
 
    rm -f OUTCAR
done;

when runned, this script will produce an output similar to the following:

sururi@husniya work $ ./energy_extract_from_OUTCAR.sh
00_ENE  :        0.00000000      01_ENE :        -6.63108952     6.63109
01_ENE  :        -6.63108952     02_ENE :        -11.59096452    4.95988
02_ENE  :        -11.59096452    03_ENE :        -12.96519853    1.37423
03_ENE  :        -12.96519853    04_ENE :        -13.20466179    0.239463
04_ENE  :        -13.20466179    05_ENE :        -13.26411934    0.0594576
05_ENE  :        -13.26411934    06_ENE :        -13.26528991    0.00117057
06_ENE  :        -13.26528991    07_ENE :        -13.40540825    0.140118
07_ENE  :        -13.40540825    08_ENE :        -13.35505746    0.0503508
08_ENE  :        -13.35505746    09_ENE :        -13.38130280    0.0262453
09_ENE  :        -13.38130280    10_ENE :        -13.36356457    0.0177382
10_ENE  :        -13.36356457    11_ENE :        -13.37065368    0.00708911
11_ENE  :        -13.37065368    12_ENE :        -13.37249683    0.00184315
12_ENE  :        -13.37249683    13_ENE :        -13.38342842    0.0109316

Which shows the final energies of the sequential runs as well as the difference of them. You can easily plot the energies in the gnuplot via as an example “plot “12_ENE”“. You can also plot the evolution of the energy difference with help from awk. To do this, make sure you first pipe the output of the script to a file (say “energies.txt”):

./energy_extract_from_OUTCAR.sh > energies.txt

And then, obtaining the last column via awk

cat energies.txt |awk '{print $7}' > enediff.txt

Now you can also easily plot the difference file.

Hope this scripts will be of any help.

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!