Hex, Bugs and More Physics | Emre S. Tasci

a blog about physics, computation, computational physics and materials…

Some unit cell functions for Octave

May 23, 2009 Posted by Emre S. Tasci

I’ve written some functions to aid in checking the transformations as well as switching between the parameter ↔ vector representations.

latpar2vec([a,b,c,alpha,beta,gamma])
Given the a,b,c,α,β and γ, the function calculates the lattice vectors

function f=latpar2vec(X)
% function f=latpar2vec([a,b,c,alpha,beta,gamma])
% Calculates the lattice vectors from the given lattice paramaters.
  
a     = X(1);
b     = X(2);
c     = X(3);
alpha = X(4);
beta  = X(5);
gamma = X(6);
 
% Converting the angles to radians:
alpha = d2r(alpha);
beta  = d2r(beta);
gamma = d2r(gamma);
 
%Aligning a with the x axis:
ax = a;
ay = 0;
az = 0;
 
%Orienting b to lie on the x-y plane:
bz = 0;
 
% Now we have 6 unknowns for 6 equations..
bx = b*cos(gamma);
by = (b**2-bx**2)**.5;
cx = c*cos(beta);
cy = (b*c*cos(alpha)-bx*cx)/by;
cz = (c**2-cx**2-cy**2)**.5;
f=[ax ay az
   bx by bz
   cx cy cz];
endfunction

latvec2par([ax ay az; bx by bz; cx cy cz])
Calculates the lattice parameters a,b,c,α,β and γ from the lattice vectors

function f=latvec2par(x)
% function f=latvec2par(x)
%     ax ay az
% x = bx by bz
%     cx cy cz
%
% takes the unit cell vectors and calculates the lattice parameters
% Emre S. Tasci <e.tasci@tudelft.nl>   03/10/08
 
av = x(1,:);
bv = x(2,:);
cv = x(3,:);
 
a = norm(av);
b = norm(bv);
c = norm(cv);
 
alpha = acos((bv*cv')/(b*c))*180/pi;
beta  = acos((av*cv')/(a*c))*180/pi;
gamma = acos((av*bv')/(a*b))*180/pi;
 
f = [a b c alpha beta gamma];
endfunction

volcell([a,b,c,A,B,G])
Calculates the volume of the cell from the given lattice parameters, which is the determinant of the matrice build from the lattice vectors.

function f=volcell(X)
% function f=volcell([a,b,c,A,B,G])
% Calculate the cell volume from lattice parameters
%
% Emre S. Tasci, 09/2008
a=X(1);
b=X(2);
c=X(3);
A=X(4); % alpha
B=X(5); % beta
G=X(6); % gamma
 f=a*b*c*(1-cos(d2r(A))^2-cos(d2r(B))^2-cos(d2r(G))^2+2*cos(d2r(A))*cos(d2r(B))*cos(d2r(G)))^.5;
endfunction

Why’s there no “volcell” function for the unit cell vectors? You’re joking, right? (det(vector)) ! 🙂

Example

octave:13> % Define the unit cell for PtSn4 :
octave:13> A = latpar2vec([ 6.41900 11.35700  6.38800  90.0000  90.0000  90.0000 ])
A =
 
    6.41900    0.00000    0.00000
    0.00000   11.35700    0.00000
    0.00000    0.00000    6.38800
 
octave:14> % Cell volume :
octave:14> Apar = [ 6.41900 11.35700  6.38800  90.0000  90.0000  90.0000 ]
Apar =
 
    6.4190   11.3570    6.3880   90.0000   90.0000   90.0000
 
octave:15> % Define the unit cell for PtSn4 :
octave:15> Apar=[ 6.41900 11.35700  6.38800  90.0000  90.0000  90.0000 ]
Apar =
 
    6.4190   11.3570    6.3880   90.0000   90.0000   90.0000
 
octave:16> % Cell volume :
octave:16> Avol = volcell (Apar)
Avol =  465.69
 
octave:17> % Calculate the lattice vectors :
octave:17> A = latpar2vec (Apar)
A =
 
    6.41900    0.00000    0.00000
    0.00000   11.35700    0.00000
    0.00000    0.00000    6.38800
 
octave:18> % Verify the volume :
octave:18> det(A)
ans =  465.69
 
octave:19> % Define the transformation matrix :
octave:19> R = [ 0 0 -1 ; -1 0 0 ; .5 .5 0]
R =
 
   0.00000   0.00000  -1.00000
  -1.00000   0.00000   0.00000
   0.50000   0.50000   0.00000
 
octave:21> % The reduced unit cell volume will be half of the original one as is evident from :
octave:21> det(R)
ans =  0.50000
 
octave:22> % Do the transformation :
octave:22> N = R*A
N =
 
  -0.00000  -0.00000  -6.38800
  -6.41900   0.00000   0.00000
   3.20950   5.67850   0.00000
 
octave:23> % The reduced cell parameters :
octave:23> Npar = latvec2par (N)
Npar =
 
     6.3880     6.4190     6.5227   119.4752    90.0000    90.0000
 
octave:24> % And the volume :
octave:24> det(N), volcell (Npar)
ans =  232.84
ans =  232.84

POSCAR2Cif

 Posted by 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 = &lt;inputfilename&gt;.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);
?>