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 = <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);
?>
