Consider the following 16×16 matrix:
0.06483 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.65809 0.65809 0.60365 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.10993 0.00000 0.00000 0.37862 0.37862 0.00000 0.00000 0.00000 0.00000 0.86568 0.00000 0.00000 0.42424 0.00000 0.00000 0.00000 0.00000 0.10993 0.00000 0.00000 0.12949 0.00000 0.00000 0.00000 0.00000 0.00000 -0.86568 0.00000 0.00000 0.42424 0.00000 0.00000 0.00000 0.00000 0.10259 0.00000 0.00000 0.84469 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.73695 0.00000 0.37862 0.00000 0.00000 0.01584 0.00000 0.00000 0.00000 0.00000 0.00000 0.84812 0.00000 0.00000 0.87983 0.00000 0.00000 0.00000 0.37862 0.12949 0.00000 0.00000 0.01584 0.00000 0.00000 0.00000 0.00000 0.00000 -0.84812 0.00000 0.00000 0.87983 0.00000 0.00000 0.00000 0.00000 0.84469 0.00000 0.00000 0.01168 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.87050 0.65809 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.88440 0.28986 0.45572 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.65809 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.28986 0.88440 0.45572 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.60365 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.45572 0.45572 0.89801 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.86568 0.00000 0.00000 0.84812 0.00000 0.00000 0.00000 0.00000 0.00000 0.35993 0.00000 0.00000 0.19336 0.00000 0.00000 0.00000 0.00000 -0.86568 0.00000 0.00000 -0.84812 0.00000 0.00000 0.00000 0.00000 0.00000 0.35993 0.00000 0.62239 -0.62239 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.18908 0.00000 0.00000 0.00000 0.00000 0.42424 0.00000 0.00000 0.87983 0.00000 0.00000 0.00000 0.00000 0.00000 0.19336 0.62239 0.00000 0.54930 0.00000 0.00000 0.00000 0.00000 0.42424 0.00000 0.00000 0.87983 0.00000 0.00000 0.00000 0.00000 0.00000 -0.62239 0.00000 0.00000 0.54930 0.00000 0.00000 0.00000 0.00000 0.73695 0.00000 0.00000 0.87050 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.55094
Equality relations:
(1,9) = (1,8)
(3,3) = (2,2)
(2,6) = (2,5)
(3,12) = -(2,11)
(3,15) = (2,14)
(6,6) = (5,5)
(6,12) = -(5,11)
(6,15) = (5,14)
(9,9) = (8,8)
(9,10) = (8,10)
(12,12) = (11,11)
(12,15) = -(12,14)
(15,15) = (14,14)
These kind of relations are sought upon tensor matrices (at least, in group theory 8), kind of postulating the imposed limits and… (have I said “relations”?). So, the sort-of-important-like thingy is to be able to spot them after evaluating some calculations in the computers.
Have you verified yet? Good! Now let’s take a look at this:
So, this is what we’ll be discussing in this post – preparation of such graphs. We’ll be using the PHP module Imagick for ImageMagick functions and objects. On Ubuntu, you can install the package from the command line as:
sudo apt-get install php5-imagick
Before going further into the graph production, here is a description of producing a random matrix compatible with the given conditions in Octave:
clear; # Declare the matrix T=[]; # == >> === non-zero coefficients ===== 0 ===== # I used the shorcuts and scripting of Vim to easily format the like lines.. non_zero_coeff = []; non_zero_coeff = [non_zero_coeff; 1 1]; non_zero_coeff = [non_zero_coeff; 1 8]; non_zero_coeff = [non_zero_coeff; 1 9]; non_zero_coeff = [non_zero_coeff; 1 10]; non_zero_coeff = [non_zero_coeff; 2 2]; non_zero_coeff = [non_zero_coeff; 2 5]; non_zero_coeff = [non_zero_coeff; 2 11]; non_zero_coeff = [non_zero_coeff; 2 14]; non_zero_coeff = [non_zero_coeff; 3 3]; non_zero_coeff = [non_zero_coeff; 3 6]; non_zero_coeff = [non_zero_coeff; 3 12]; non_zero_coeff = [non_zero_coeff; 3 15]; non_zero_coeff = [non_zero_coeff; 4 4]; non_zero_coeff = [non_zero_coeff; 4 7]; non_zero_coeff = [non_zero_coeff; 3 12]; non_zero_coeff = [non_zero_coeff; 4 16]; non_zero_coeff = [non_zero_coeff; 5 5]; non_zero_coeff = [non_zero_coeff; 5 11]; non_zero_coeff = [non_zero_coeff; 1 9]; non_zero_coeff = [non_zero_coeff; 5 14]; non_zero_coeff = [non_zero_coeff; 6 6]; non_zero_coeff = [non_zero_coeff; 6 12]; non_zero_coeff = [non_zero_coeff; 1 9]; non_zero_coeff = [non_zero_coeff; 6 15]; non_zero_coeff = [non_zero_coeff; 7 7]; non_zero_coeff = [non_zero_coeff; 7 16]; non_zero_coeff = [non_zero_coeff; 8 8]; non_zero_coeff = [non_zero_coeff; 8 9]; non_zero_coeff = [non_zero_coeff; 8 10]; non_zero_coeff = [non_zero_coeff; 6 15]; non_zero_coeff = [non_zero_coeff; 9 9]; non_zero_coeff = [non_zero_coeff; 9 10]; non_zero_coeff = [non_zero_coeff; 10 10]; non_zero_coeff = [non_zero_coeff; 11 11]; non_zero_coeff = [non_zero_coeff; 11 14]; non_zero_coeff = [non_zero_coeff; 12 12]; non_zero_coeff = [non_zero_coeff; 12 14]; non_zero_coeff = [non_zero_coeff; 13 13]; non_zero_coeff = [non_zero_coeff; 14 14]; non_zero_coeff = [non_zero_coeff; 15 15]; non_zero_coeff = [non_zero_coeff; 16 16]; # Assign random values to the non-zero coefficients: for i=1:rows(non_zero_coeff) T(non_zero_coeff(i,1),non_zero_coeff(i,2)) = rand; endfor # == << === non-zero coefficients ===== 1 ===== # == >> === Relations ===== 0 ===== T(1,9)=T(1,8); T(3,3)=T(2,2); T(2,6)=T(2,5); T(3,12)=-T(2,11); T(3,15)=T(2,14); T(6,6)=T(5,5); T(6,12)=-T(5,11); T(6,15)=T(5,14); T(9,9)=T(8,8); T(9,10)=T(8,10); T(12,12)=T(11,11); T(12,15)=-T(12,14); T(15,15)=T(14,14); # == << === Relations ===== 1 ===== # == >> === Symmetry ===== 0 ===== for i=1:rows(T) for j=i:columns(T) T(j,i) = T(i,j); endfor endfor # == << === Symmetry ===== 1 ===== # Export the matrix to a file: save("input.data.txt","T");
Doing a double loop over i=1:rows and j=i:cols, for each of the coefficients in and above the diagonal, taking one as reference (checking if abs(T(i,j))>0 and recording to an array, say “non-zero-coeff-array” if so) and then comparing it with the rest of the elements (that comes after it) with another double loop i2=i:rows and j2=j:cols. If T(i2,j2) is equal to +-T(i,j), we also note this (three values for each i2,j2 equivalency identification: i2,j2 and +/-1) in another array (say “equal-array”).
For example, according to our case the (9,9) element of the non-zero-coeff-array would be TRUE and the same element of the equal-array would be a sub-array holding the values (8,8,1).
Once you have access to the ImageMagick library through Imagick, we can start “painting” a canvas by first defining the canvas (width, height and background color):
# Declare the object (/holder)
$image = new Imagick();
# Define the canvas (size 640x640 pixels with white background)
$image->newImage( 640, 640, new ImagickPixel( 'white' ) );
Let’s put a circle (filled, black) in the center with a radius of 30 pixels:
# Define the circle:
$circle = new ImagickDraw();
$circle->setFillColor("#000000");
$circle->circle(320,320,350,320);
# Add (/register) it to the canvas:
$image->drawImage($circle);
As for the production of the image file, we need to designate it as a PNG image and then export it to a file:
<code># Designate the image type as PNG:
$image->setImageFormat( "png" );
# Export (/write) it to a file "out.png":
$fp = fopen("out.png","w");
fwrite($fp,$image);
fclose($fp);
This is what you -should- get:
not bad for beginning, right?
The circle function has 4 parameters: the (x,y) coordinates of the center and the coordinates of a point that is on the edge of the circle. Let’s draw another circle, this time an unfilled one below this filled one:
# Declare the object (/holder)
$image = new Imagick();
# Define the canvas (size 640x640 pixels with white background)
$image --->newImage( 640, 640, new ImagickPixel( 'white' ) );
# Define the 1st, filled circle:
$circle = new ImagickDraw();
$circle->setFillColor("#000000");
$circle->circle(320,320,350,320);
# Add (/register) it to the canvas:
$image->drawImage($circle);
# Define the 2nd, unfilled circle:
$circle = new ImagickDraw();
$circle->setStrokeColor("#000000");
$circle->setStrokeWidth(2);
$circle->setFillColor("none");
$circle->circle(320,350,350,350);
# Add (/register) it to the canvas:
$image->drawImage($circle);
# Designate the image type as PNG:
$image->setImageFormat( "png" );
# Export (/write) it to a file "out.png":
$fp = fopen("out.png","w");
fwrite($fp,$image);
fclose($fp);
The “stroke” methods are responsible for the perimeter and the fillcolor setting to “none” ensures that we have an otherwise invisible circle.
So far, so good.. Let’s add a line that passes between the circles from (120,50) to (520,50):
# Define the line:
$line = new ImagickDraw();
$line->setStrokeColor("#000000");
$line->setStrokeWidth(2);
$line->line(120,350,520,350);
$image->drawImage($line);
Suppose that we have a 16×16 matrix to be mapped onto a 640×640 canvas. We could map (i=0,j=0) to (x=0,y=0) and (i=16,j=16) to (x=640,y=640) but it wouldn’t be right:
doesn’t look right, right? And here’s the code that produced it:
# Declare the object (/holder)
$image = new Imagick();
# Define the canvas (size 640x640 pixels with white background)
$image--->newImage( 640, 640, new ImagickPixel( 'white' ) );
# Define the 1st, filled circle:
$dxy = 640 / 16;
for($i=1;$i<=16;$i++)
{
for($j=1;$j<=16;$j++)
{
$circle = new ImagickDraw();
$circle->setFillColor("#000000");
$circle->circle($dxy*($j-1),$dxy*($i-1),$dxy*($j-1)+$dxy/6,$dxy*($i-1));
# Add (/register) it to the canvas:
$image->drawImage($circle);
}
}
# Designate the image type as PNG:
$image->setImageFormat( "png" );
# Export (/write) it to a file "out.png":
$fp = fopen("out.png","w");
fwrite($fp,$image);
fclose($fp);
What we need to do is to shift it half of the periodicity in the x & y directions:
That’s more like it — this shift was introduced to the code by changing the coordinate parameters of the circle(s):
$circle->circle($dxy*($j-1/2),$dxy*($i-1/2),$dxy*($j-1/2)+$dxy/6,$dxy*($i-1/2));
So, depending on the value (zero/non-zero) of an element of our tensor, we can put it as a big circle (r=$dxy/6) or a small one (r=$dxy/12). In addition, if it’s negative, we can have it drawn as an unfilled one (
$circle->setStrokeColor("#000000"); $circle->setStrokeWidth(2); $circle->setFillColor("none");
Suppose that we’d like to connect the circles corresponding to the (i1,j1) and the (i2,j2) elements of the tensor. On the canvas, their centers will be given as ($dxy*($j1-1/2),$dxy*($i1-1/2)) and ($dxy*($j2-1/2),$dxy*($i2-1/2)) where $dxy is the conversion unit of 1 step in the row-col matrix to that of the canvas (for simplicity I assumed a square canvas which is not necessary — then we’d have a $dx=$width/$num_cols and a $dy=$height/$num_rows. Also note that the rows are associated with the height ($i <-> vertical) and the columns are with the width ($j <-> horizontal) (that’s why $j’s come before $i’s in the center equations).
So, let’s join (11,5) to (9,3): in a 640×640 canvas, (11,5) is centered at ((640/16)*(5-1/2) = 180, (640/16)*(11-1/2) = 420) and similarly (9,3) is centered at (100,340). Let’s put them on the map and connect their centers with a line:
# Declare the object (/holder)
$image = new Imagick();
# Define the canvas (size 640x640 pixels with white background)
$image->newImage( 640, 640, new ImagickPixel( 'white' ) );
$dxy = 640/16; # 40px
$radius = $dxy/6; #6.6667px
# Define the unfilled (11,5) circle:
$circle = new ImagickDraw();
$circle->setStrokeColor("#000000");
$circle->setStrokeWidth(2);
$circle->setFillColor("none");
$circle->circle(180,420,180+$radius,420);
# Add (/register) it to the canvas:
$image->drawImage($circle);
# Define the unfilled (9,3) circle:
$circle = new ImagickDraw();
$circle->setStrokeColor("#000000");
$circle->setStrokeWidth(2);
$circle->setFillColor("none");
$circle->circle(100,340,100+$radius,340);
# Add (/register) it to the canvas:
$image->drawImage($circle);
# Connect them with a line:
$line = new ImagickDraw();
$line->setStrokeColor("#000000");
$line->setStrokeWidth(2);
$line->line(180,420,100,340);
$image->drawImage($line);
# Designate the image type as PNG:
$image->setImageFormat( "png" );
# Export (/write) it to a file "out.png":
$fp = fopen("out.png","w");
fwrite($fp,$image);
fclose($fp);
By now, hopefully you’ve got the idea of how to swing things in your way. In my case, if I spent the 40% of the time doing the things I described, the remaining 60% went in adjusting the line such that it starts on the perimeter of the circle, not from the center of it, i.e.,
and this requires the determination of the angle of the line connecting our nodes. The angle of a line passing through two points (x1,y1) & (x2,y2) is given by: arctan[(y2-y1)/(x2-x1)], where the value is nothing but the A of the y=Ax+B line equation, and we need to translate the beginning and the final points of our line in this direction by an amount of a radius:
# Find the tangent of the line:
$tg = (340-420)/(100-180);
$alpha_radian = atan($tg);
$alpha_deg = $alpha_radian * 180 / M_PI;
# Connect them with a line:
$line = new ImagickDraw();
$line->setStrokeColor("#000000");
$line->setStrokeWidth(2);
$line->line(180-$radius*cos($alpha_radian),420-$radius*sin($alpha_radian),100+ $radius*cos($alpha_radian),340+$radius*sin($alpha_radian));
$image->drawImage($line);
and that’s more or less all about it!
You can have a sneak peak of the working code via here: http://144.122.31.125/cgi-bin/cryst/programs/tensor_graph/
]]>Abstract
This how-to is intended to guide the reader during the preparation of the input file for an intended surface calculation study. Mainly VESTA software is used to manipulate and transform initial structure file along with some home brewed scripts to do simple operations such as translations.
[1] Mois Ilia Aroyo, Juan Manuel Perez-Mato, Cesar Capillas, Eli Kroumova, Svetoslav Ivantchev, Gotzon Madariaga, Asen Kirov, and Hans Wondratschek. Bilbao crystallographic server: I. databases and crystallographic computing programs. Zeitschrift für Kristallographie, 221(1):1527, 01 2006. doi: 10.1524/zkri.2006.221.1.15.
[2] Mois I. Aroyo, Asen Kirov, Cesar Capillas, J. M. Perez-Mato, and Hans Wondratschek. Bilbao crystallographic server. ii. representations of crystallographic point groups and space groups. Acta Crystallographica Section A, 62(2):115128, Mar 2006. http://dx.doi.org/10.1107/S0108767305040286
[3] M I Aroyo, J M Perez-Mato, D Orobengoa, E Tasci, G De La Flor, and A Kirov. Crystallography online: Bilbao crystallographic server. Bulg Chem Commun, 43(2):183197, 2011.
[4] S. R. Hall, F. H. Allen, and I. D. Brown. The crystallographic information le (cif): a new standard archive file for crystallography. Acta Crystallographica Section A, 47(6):655685, Nov 1991.
[5] G. Kresse and J. Furthmüller. Ecient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B, 54:1116911186, Oct 1996.
[6] Koichi Momma and Fujio Izumi. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. Journal of Applied Crystallography, 44(6):12721276, Dec 2011.
[7] A Zaoui and M Ferhat. Dynamical stability and high pressure phases of platinum carbide. Solid State Communications, 151(12):867869, 6 2011.
[8] Saulius Grazulis, Adriana Daskevic, Andrius Merkys, Daniel Chateigner, Luca Lutterotti, Miguel Quirós, Nadezhda R. Serebryanaya, Peter Moeck, Robert T. Downs, and Armel Le Bail. Crystallography open database (cod): an open-access collection of crystal structures and platform for world-wide collaboration. Nucleic Acids Research, 40(D1):D420D427, 2012.
[9] Alec Belsky, Mariette Hellenbrandt, Vicky Lynn Karen, and Peter Luksch. New developments in the inorganic crystal structure database (icsd): accessibility in support of materials research and design. Acta Crystallographica Section B, 58(3 Part 1):364369, Jun 2002.
[10] Springer, editor. The Landolt-Boernstein Database (http://www.springermaterials.com/navigation/). Springer Materials.
[11] Yongsheng Zhang. First-principles statistical mechanics approach to step decoration at solid surfaces . PhD thesis, Frei Universitat Berlin, 2008.
Suppose that we want to calculate the energies corresponding to various values of the lattice parameter a, ranging from 3.2 to 4, incrementing by 2: 3.2, 3.4, … , 4.0
Incrementation:
In bash, you can define a for loop to handle this:
for (( i = 1; i <= 9; i++ ))
do
echo $i;
done
1
2
3
4
5
6
7
8
9
for (( i = 1; i <= 5; i++ ))
do
val=$(awk "BEGIN {print 3+$i/5 }")
echo $val;
done
3.2
3.4
3.6
3.8
4
octave:1> [3.2:0.2:4] ans = 3.2000 3.4000 3.6000 3.8000 4.0000
sururi@husniya:~/shared/qe/diamond_tutorial$ echo "[3.2:0.2:4]"|octave -q ans = 3.2000 3.4000 3.6000 3.8000 4.0000
octave:4> [3:0.1:4] ans = Columns 1 through 9: 3.0000 3.1000 3.2000 3.3000 3.4000 3.5000 3.6000 3.7000 3.8000 Columns 10 and 11: 3.9000 4.0000
sururi@vala:/tmp$ echo "[3.2:0.1:4]"|octave -q|sed -n "2,1000p"| grep -v "Column" 3.2000 3.3000 3.4000 3.5000 3.6000 3.7000 3.8000 3.9000 4.0000 sururi@vala:/tmp$ echo "[3.2:0.1:4]"|octave -q|sed -n "2,1000p"| grep -v "Column"|grep -v -e "^$" 3.2000 3.3000 3.4000 3.5000 3.6000 3.7000 3.8000 3.9000 4.0000
cat diamond.scf.in | sed "s:#REPLACE#:$a:g">tmp_$a.in
Parsing the Energy:
We submit this job via mpirun (in our case via mpi-pw.x of the previous post) and parse the output for the final energy:
grep -e ! tmp_$a.out|tail -n1|awk '{print $(NF-1)}'
The Script:
Putting all these together, here is the script in its full glory:
#!/bin/bash
# Takes in the template Quantum Espresso input file, runs it iteratively
# changing the #REPLACE# parameters iteratively with the range elements,
# grepping the energy and putting it along with the current value of the
# #REPLACE# parameter.
#
# Emre S. Tasci, 03/2012
# Example: Suppose that the file "diamond.scf.in" contains the following line:
# A = #REPLACE#
# Calling "qe-opti.sh diamond.scf.in a 3.2 0.2 4.2" would cause that line to be modified as
# A = 3.2
# and submitted, followed by
# A = 3.4, 3.6, .., 4.2 in the sequential runs.
#
# at the end, having constructed the E_vs_a.txt file.
if [ $# -ne 5 ]
then
echo "Usage: qe-opti.sh infile_template.in variable_name initial_value increment final_value"
echo "Example: qe-opti.sh diamond.scf.in a 3.2 0.2 4.2"
echo -e "\nNeeds octave to be accessible via \"octave\" command for the range to work"
echo -e "\nEmre S. Tasci, 03/2012"
exit
fi
logext=".txt"
logfilename=E_vs_$2$logext
rm -f $logfilename
range=$(echo "[$3:$4:$5]"|octave -q|sed -n "2,1000p"|grep -vr "^$"|grep -v "Column")
for a in $range
do
# echo $a
cat $1|sed "s:#REPLACE#:$a:g">tmp_$a.in
/home/sururi/bin/mpi-pw.x 4 tmp_$a.in
energ=$(grep -e ! tmp_$a.out|tail -n1|awk '{print $(NF-1)}')
# echo $energ
echo -e $a " " $energ >> $logfilename
echo -e $a " " $energ
done
echo "Results stored in: "$logfilename
sururi@bebop:~/shared/qe/diamond_tutorial$ grep "REPLACE" -B2 -A2 diamond.scf.in # a,b,c in Angstrom ======= #A = 3.56712 A = #REPLACE# B = 2.49 C = 2.49
sururi@bebop:~/shared/qe/diamond_tutorial$ qe-opti.sh diamond.scf.in a 3.2 0.2 4.2 Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.2000.in > tmp_3.2000.out Start : Tue Mar 13 23:21:07 CET 2012 Finish: Tue Mar 13 23:21:08 CET 2012 Elapsed time:0:00:01 3.2000 -22.55859979 Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.4000.in > tmp_3.4000.out Start : Tue Mar 13 23:21:08 CET 2012 Finish: Tue Mar 13 23:21:09 CET 2012 Elapsed time:0:00:01 3.4000 -22.67717405 Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.6000.in > tmp_3.6000.out Start : Tue Mar 13 23:21:09 CET 2012 Finish: Tue Mar 13 23:21:10 CET 2012 Elapsed time:0:00:01 3.6000 -22.69824702 Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.8000.in > tmp_3.8000.out Start : Tue Mar 13 23:21:10 CET 2012 Finish: Tue Mar 13 23:21:11 CET 2012 Elapsed time:0:00:01 3.8000 -22.65805473 Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_4.0000.in > tmp_4.0000.out Start : Tue Mar 13 23:21:11 CET 2012 Finish: Tue Mar 13 23:21:13 CET 2012 Elapsed time:0:00:02 4.0000 -22.57776354 Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_4.2000.in > tmp_4.2000.out Start : Tue Mar 13 23:21:13 CET 2012 Finish: Tue Mar 13 23:21:14 CET 2012 Elapsed time:0:00:01 4.2000 -22.47538891 Results stored in: E_vs_a.txt
sururi@bebop:~/shared/qe/diamond_tutorial$ cat E_vs_a.txt 3.2000 -22.55859979 3.4000 -22.67717405 3.6000 -22.69824702 3.8000 -22.65805473 4.0000 -22.57776354 4.2000 -22.47538891
sururi@husniya:~/shared/qe/diamond_tutorial$ plot E_vs_a.txt
Other than the ‘plot’, I’ve also written another one, fit for data sets just like these where you have a parabol-like data and you’d like to try fitting it and here is what it does (automatically):
sururi@husniya:~/shared/qe/diamond_tutorial$ plotfitparabol.sh E_vs_a.txt func f(x)=0.674197*x*x-4.881275*x-13.855234 min, f(min) 3.620066 -22.690502
Here are the source codes of the ‘plot’ and ‘plotfitparabol.sh’ scripts. Doei!
sururi@mois:/vala/bin$ cat plot
#!/bin/bash
gnuplot -persist <<END
set term postscript enhanced color
set output "$1.ps"
plot "$1" $2 $3
set term wxt
replot
END
sururi@mois:/vala/bin$ cat plotfitparabol.sh
#!/bin/bash
# Usage : plotfitparabol.sh <datafile>
# Reads the data points from the given file
# then uses Octave's "polyfit" function to calculate the coefficients
# for the quadratic fit and also the minimum (maximum) value.
# Afterwards, prepares a GNUPlot script to display:
# the data set (scatter)
# fit function
# the extremum point's value
#
# Written with Quantum Espresso optimizing procedures in mind.
if [ $# -ne 1 ]
then
echo "Usage: plotfitparabol.sh <datafile>"
echo "Example: plotfitparabol.sh E_vs_A.txt"
echo -e "\nNeeds octave to be accessible via \"octave\" command to work"
echo -e "\nEmre S. Tasci <emre.tasci@ehu.es>, 03/2012"
exit
fi
octave_script="data=load(\"$1\");m=polyfit(data(:,1),data(:,2),2);min=-m(2)*0.5/m(1);fmin=min*min*m(1)+min*m(2)+m(3);printf(\"f(x)=%f*x*x%+f*x%+f|%f\t%f\",m(1),m(2),m(3),min,fmin)"
oct=$(echo $octave_script|octave -q)
#echo $oct
func=${oct/|*/ }
minfmin=${oct/*|/ }
echo -e "func \t" $func
echo -e "min, f(min) \t" $minfmin
echo '$func
set term postscript enhanced color
set label "$1" at graph 0.03, graph 0.94
set output "$1.ps"
plot "$1" title "data", f(x) title "$func", "<echo '$minfmin'" title "$minfmin"
set term wxt
replot'
gnuplot -persist <<END
$func
set term postscript enhanced color
set label "$1" at graph 0.03, graph 0.94
set output "$1.ps"
plot "$1" title "data", f(x) title "$func", "<echo '$minfmin'" title "$minfmin"
set term wxt
replot
END
]]>As a change, I selected and installed the Quantum Espresso (QE) package instead of VASP that I was using in the past, main reason QE being open source.
As a newbie to QE and long-gone returner to DFT, I started to practice by mostly ye olde trial/error approach.
The first thing that surprised me the lack of a commenting option in regard to QE input files! (Actually, I’m almost sure that it surely is supported but I have yet to find!). Here’s my sample input file:
&control title='Diamond Relaxing', # Calculation Type ========== calculation='scf' #calculation='relax', #calculation='vc-relax' #calculation='cp' # Restart =================== restart_mode='from_scratch', #restart_mode='restart' # Read/Write units ========== #ndr = 51, #ndw = 51, #nstep = 100, #nstep = 1000, nstep = 2000, # max_seconds = 600 # jobs stops after max_seconds CPU time. iprint = 10, # band energies are written every iprint iterations #tstress = .TRUE., # calculate stress. It is set to .TRUE. automatically if # calculation='vc-md' or 'vc-relax' #tprnfor = .TRUE., # print forces. Set to .TRUE. if calculation='relax','md','vc-md' #dt = 5.0d0, # time step for molecular dynamics, in Rydberg atomic units # (1 a.u.=4.8378 * 10^-17 s : beware, the CP code use # Hartree atomic units, half that much!!!) prefix='diamond_rel', # prepended to input/output filenames pseudo_dir = '/usr/share/espresso-4.3.2/pseudo/', outdir='/home/sururi/shared/tmp/' etot_conv_thr = 1.d-6 # convergence threshold on total energy (a.u) for ionic minimization forc_conv_thr = 1.d-3 # convergence threshold on forces (a.u) for ionic minimization / &system ## ibrav ## 0 "free", see above not used ## 1 cubic P (sc) not used ## 2 cubic F (fcc) not used ## 3 cubic I (bcc) not used ## 4 Hexagonal and Trigonal P celldm(3)=c/a ## 5 Trigonal R, 3fold axis c celldm(4)=cos(alpha) ## -5 Trigonal R, 3fold axis <111> celldm(4)=cos(alpha) ## 6 Tetragonal P (st) celldm(3)=c/a ## 7 Tetragonal I (bct) celldm(3)=c/a ## 8 Orthorhombic P celldm(2)=b/a,celldm(3)=c/a ## 9 Orthorhombic base-centered(bco) celldm(2)=b/a,celldm(3)=c/a ## 10 Orthorhombic face-centered celldm(2)=b/a,celldm(3)=c/a ## 11 Orthorhombic body-centered celldm(2)=b/a,celldm(3)=c/a ## 12 Monoclinic P, unique axis c celldm(2)=b/a,celldm(3)=c/a, ## celldm(4)=cos(ab) ## -12 Monoclinic P, unique axis b celldm(2)=b/a,celldm(3)=c/a, ## celldm(5)=cos(ac) ## 13 Monoclinic base-centered celldm(2)=b/a,celldm(3)=c/a, ## celldm(4)=cos(ab) ## 14 Triclinic celldm(2)= b/a, ## celldm(3)= c/a, ## celldm(4)= cos(bc), ## celldm(5)= cos(ac), ## celldm(6)= cos(ab ibrav= 2, # Cell dimensions ======================== # a,b,c in Angstrom ======= A = 3.56712 B = 2.49 C = 2.49 cosAB=0.5 cosAC=0.5 cosBC=0.5 ## OR ## ## celldm(1:6) in Bohr ===== ## 1 Bohr = 0.529177249 A ## 1 A = 1.889725989 Bohr # celldm(1) = 6.7409 # alat # celldm(2) = 4.7054 # b/a # celldm(3) = 4.7054 # c/a # celldm(4) = 0.5 # cos(ab) # celldm(5) = 0.5 # cos(ac) # celldm(6) = 0.5 # cos(ab) nat=2, ntyp= 1, ecutwfc = 40, # kinetic energy cutoff (Ry) for wavefunctions # ecutrho = 160 # kinetic energy cutoff (Ry) for charge density and potential. # For norm-conserving pseudopotential you should stick to the # default value # nbnd = 6 # number of electronic states (bands) to be calculated. # Default: # for an insulator, nbnd = number of valence bands (nbnd = # of electrons /2); # for a metal, 20% more (minimum 4 more) / &electrons emass = 400.d0 emass_cutoff = 2.5d0 #electron_dynamics = 'sd' # CP specific #electron_dynamics = 'bfgs' # CP specific conv_thr = 1D-6 # Convergence threshold for selfconsistency / &ions #ion_dynamics = 'none' #ion_dynamics = 'bfgs' #ion_dynamics = 'sd' # === damping === 0 == # ion_dynamics = 'damp' # ion_damping = 0.2 # ion_velocities = 'zero' # CP specific -- initial ionic velocities # === damping === 1 == #ion_nstepe = 10 # CP specific -- number of electronic steps per ionic step. (Def: 1) / &cell cell_dynamics = 'none', #cell_dynamics = 'bfgs', # ion_dynamics must be 'bfgs' too #cell_dofree = 'xyz' #cell_factor = 3 press = 0.0d, # press_conv_thr = 0.05 # Convergence threshold on the pressure for variable cell relaxation / ATOMIC_SPECIES C 12.0107 C.pz-vbc.UPF ATOMIC_POSITIONS C 0.0 0.0 0.0 C 0.25 0.25 0.25 K_POINTS automatic 4 4 4 1 1 1
And then the problem of (yet undiscovered) comments in the QE input files hits us. But with bash, it’s pretty easy. So, if you save the above contents to a file called, say, “diamond.scf.in” and filter it via grep and sed, you’re done:
sururi@husniya:~/shared/qe/diamond_tutorial$ cat diamond.scf.in | grep -v -E "^[\t ]*#" | sed "s:\([^#]*\)#.*:\1:" &control title='Diamond Relaxing', calculation='scf' restart_mode='from_scratch', nstep = 2000, iprint = 10, prefix='diamond_rel', pseudo_dir = '/usr/share/espresso-4.3.2/pseudo/', outdir='/home/sururi/shared/tmp/' etot_conv_thr = 1.d-6 forc_conv_thr = 1.d-3 / &system ibrav= 2, A = 3.56712 B = 2.49 C = 2.49 cosAB=0.5 cosAC=0.5 cosBC=0.5 nat=2, ntyp= 1, ecutwfc = 40, / &electrons emass = 400.d0 emass_cutoff = 2.5d0 conv_thr = 1D-6 / &ions / &cell cell_dynamics = 'none', press = 0.0d, / ATOMIC_SPECIES C 12.0107 C.pz-vbc.UPF ATOMIC_POSITIONS C 0.0 0.0 0.0 C 0.25 0.25 0.25 K_POINTS automatic 4 4 4 1 1 1
mpirun -n #_of_processes /usr/share/espresso-4.3.2/bin/pw.x < input_file > output_file
which I guess is more or less similar to what you are doing to run your jobs.
So, combining this “comment filtering” and job submission, I’ve written the following bash script (“mpi-pw.x“) to automatize (and keep log of) the job:
sururi@cowboy:~$ cat /home/sururi/bin/mpi-pw.x
#!/bin/bash
# Script to submit Quantum Espresso jobs via mpi
# EST, Fri Mar 9 17:02:10 CET 2012
function timer()
{
if [[ $# -eq 0 ]]; then
echo $(date '+%s')
else
local stime=$1
etime=$(date '+%s')
if [[ -z "$stime" ]]; then stime=$etime; fi
dt=$((etime - stime))
ds=$((dt % 60))
dm=$(((dt / 60) % 60))
dh=$((dt / 3600))
printf '%d:%02d:%02d' $dh $dm $ds
fi
}
if [ $# -ne 2 ]
then
echo "Usage: mpi-pw.x num_of_processes infile.in"
exit
fi
job=$2
job_wo_ext=${job%.*}
outext=".out"
logext=".log"
echo -e "Running: mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < $2 > $job_wo_ext$outext\n"
echo -e "Running: mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < $2 > $job_wo_ext$outext\n">$job_wo_ext$logext
echo "Start : "$(date)
echo "Start : "$(date)>>$job_wo_ext$logext
# Take out the comments from the input file
#grep -v -E "^[\t ]*#" $2 > in_wo_comments.in.tmp
grep -v -E "^[\t ]*#" $2 | sed "s:\([^#]*\)#.*:\1:" > in_wo_comments.in.tmp
t=$(timer)
mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < in_wo_comments.in.tmp > $job_wo_ext$outext
echo "Finish: "$(date)
echo "Finish: "$(date)>>$job_wo_ext$logext
printf 'Elapsed time:%s\n\n' $(timer $t)
printf 'Elapsed time:%s\n\n' $(timer $t) >>$job_wo_ext$logext
rm -f in_wo_comments.in.tmp
(I’ve picked the “timer” function from Mitch Frazier’s entry in LinuxJournal)
What it does is, really simple: it strips out the extension of the input file ($2 – the second argument in calling); filters out the comments, constructing the temporary input file “in_wo_comments.in.tmp”; starts the timer; submits the job using the number of processes passed in the calling ($1 – first argument) and directs the output to a file whose name is the extension stripped input filename + “.out”; when the job is (this way or that) done, stops the timer and calculates the elapsed time via the ‘timer’ function. During all this time, it prints these information both on screen and to a file called extension stripped input filename + “.log”.
I have a very similar another script called “mpi-cp.x” that I’ve forked for Carr-Parinelli jobs (just search & replace all the occurrences of “mpi-pw.x” with “mpi-cp.x” and you’re done 8).
And here is the script in action:
sururi@cowboy:~/shared/qe/diamond_tutorial$ ll total 16 drwxr-xr-x 2 sururi sururi 4096 2012-03-13 15:59 ./ drwxr-xr-x 7 sururi sururi 4096 2012-03-13 15:10 ../ -rw-r--r-- 1 sururi sururi 4954 2012-03-13 15:59 diamond.scf.in sururi@cowboy:~/shared/qe/diamond_tutorial$ mpi-pw.x 4 diamond.scf.in Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < diamond.scf.in > diamond.scf.out Start : Tue Mar 13 16:18:20 CET 2012 Finish: Tue Mar 13 16:18:21 CET 2012 Elapsed time:0:00:01 sururi@cowboy:~/shared/qe/diamond_tutorial$ grep -e ! diamond.scf.out ! total energy = -22.70063050 Ry sururi@cowboy:~/shared/qe/diamond_tutorial$ cat diamond.scf.log Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < diamond.scf.in > diamond.scf.out Start : Tue Mar 13 16:18:20 CET 2012 Finish: Tue Mar 13 16:18:21 CET 2012 Elapsed time:0:00:01
Now that we have the mpi-pw.x and mpi-cp.x, our next entry in the series will be the automated search for optimized values.
]]>The red arrows indicate the mercurial transfers. You can see that, the changes are only applied to the servers, i.e., nobody codes directly there (no push, only pulls and ups). My two computers can also be thought of two different developers but that part actually isn’t very relevant to this entry’s message. The thing I want to happen is: whenever I push a change into the central repository, I would like to see that it is applied on the servers. In other words I want the following procedures to be automated after the PC@wherever: hg push operation:
Even though I’ll be using Mercurial in this entry, SVN implementation is pretty much similar. After all, we’ll be using the so-called hooks, i.e. the procedures that takes place upon a triggered event such as push.
Since we will be automating (and I’ll be using SSH for the communication between the nodes), it is essential that the nodes can communicate freely via the help of the ssh-keys (all the relevant information including the usage of “hg-ssh” can be found in a previous entry titled ‘Accessing Mercurial with limited SSH access using key and hg-ssh’).
The trigger event is the incoming data to the central repository, so we edit the .hg/hgrc file of the project on the central repository and type:
[hooks]
incoming = .hg/incominghook.sh > /dev/null
#!/bin/sh
# Go to the project's dir no matter where you are
cd /path/to/the/project/repository/
# Update it (since this is the central repo, there is no need to 'pull',
# everybody is pushing to here)
hg up
# Connect to the servers using the id files specifically generated for this purpose
# so that, upon connection, filtered by the command option
ssh -i ~/.ssh/id_rsa_specific_id_for_the_project__central_repo server1
ssh -i ~/.ssh/id_rsa_specific_id_for_the_project__central_repo server2
The ~/.ssh/authorized_keys file on server1 contains the following related entry corresponding to the supplied id being used in the connection made from central repo:
command="/home/sururi/bin/hgpull_project_x.sh",
no-port-forwarding,no-X11-forwarding,no-agent-forwarding,no-pty
ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAnq.....==sururi@centralRepo (project X)
cd /path/to/the/project/on/server1
hg pull --ssh "ssh -i ~/.ssh/id_rsa_specific_id_for_the_project__server1" ssh://sururi@centralrepo//path/to/project
hg up
command="hg-ssh /path/to/project/",no-port-forwarding, no-X11-forwarding,no-agent-forwarding,no-pty ssh-rsa AAAAB...
To limit an SSH connection, you use ssh-keys: you create a pair of keys, private and public, using the ‘ssh-keygen’ command and then adding the public one to the ~/.ssh/authorized_keys. As an example, consider two computers ‘local’ and ‘remote’ and we want to connect to remote from local without having to enter password every time. So, from the console of local, first I create the key pair:
sururi@local:/tmp/tmp$ ssh-keygen
Generating public/private rsa key pair.
Enter file in which to save the key (/home/sururi/.ssh/id_rsa): ./id_rsa_example
Enter passphrase (empty for no passphrase):
Enter same passphrase again:
Your identification has been saved in ./id_rsa_example.
Your public key has been saved in ./id_rsa_example.pub.
The key fingerprint is:
84:59:8d:8c:43:4a:bc:31:ff:5a:12:23:34:45:56:67 sururi@remote
The key's randomart image is:
+--[ RSA 2048]----+
|  oAE*=A+o      |
| .o.=..+=..     |
|  .o...       |
|Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â |
|Â Â Â Â Â Â Â Â XXXÂ Â |
+-----------------+
Now, I should append the contents of the ‘id_rsa_example.pub’ file to the ./ssh/authorized_keys at the remote computer. To do this, of course, I should be able to connect it via SSH by normal means. So I execute (from the local computer):
ssh your_username@remote 'cat >> ~/.ssh/authorized_keys' < ./id_rsa_example.pub
We can limit the things an SSH-connected user can do to a single thing via the authorized keys option. Checking this file on the remote computer’s ~/.ssh/authorized_keys file, you should see something like:
ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAnq5rMLnoab+2F28g/nb58RBENWtX395TuyDFsYkalGaZxrziwDoau/wglkU19DbcAVKgw0p6lMEIuh2iALOppRzxrTgFFhJkL1dxzkugbbPEoSWyfrj9FivzpnxHWgRHQApQeWUBOZhroDTURwfqcyC9SW020CR57jLWfgw+idqwtCu+ZBYmEyHSJcZIH2mWXLrUQ8OalxCFVaLKL50Lpc7V8XJPs+Pg6MPVgfDUqMdjrGkAF7j4viOHTjDWP1h4Ngim70dOeyxWtuqbCbxM4APTShaqET42sj1jHxL2m1dJzXX8s/gEdN0O09hZPhI6rlC+ANWIdJ1vJfODMXWaQ== sururi@local
command="date" ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAnq5rMLnoab+2F28g/nb58RBENWtX395TuyDFsYkalGaZxrziwDoau/wglkU19DbcAVKgw0p6lMEIuh2iALOppRzxrTgFFhJkL1dxzkugbbPEoSWyfrj9FivzpnxHWgRHQApQeWUBOZhroDTURwfqcyC9SW020CR57jLWfgw+idqwtCu+ZBYmEyHSJcZIH2mWXLrUQ8OalxCFVaLKL50Lpc7V8XJPs+Pg6MPVgfDUqMdjrGkAF7j4viOHTjDWP1h4Ngim70dOeyxWtuqbCbxM4APTShaqET42sj1jHxL2m1dJzXX8s/gEdN0O09hZPhI6rlC+ANWIdJ1vJfODMXWaQ== sururi@local== sururi@vala
sururi@local:~$ ssh -i /tmp/tmp/id_rsa_example sururi@remote
Sun Aug 28 22:05:54 CEST 2011
Connection to remote closed.
command="hg-ssh /path/to/repository" ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAnq5rMLnoab+2F28g/nb58RBENWtX395TuyDFsYkalGaZxrziwDoau/wglkU19DbcAVKgw0p6lMEIuh2iALOppRzxrTgFFhJkL1dxzkugbbPEoSWyfrj9FivzpnxHWgRHQApQeWUBOZhroDTURwfqcyC9SW020CR57jLWfgw+idqwtCu+ZBYmEyHSJcZIH2mWXLrUQ8OalxCFVaLKL50Lpc7V8XJPs+Pg6MPVgfDUqMdjrGkAF7j4viOHTjDWP1h4Ngim70dOeyxWtuqbCbxM4APTShaqET42sj1jHxL2m1dJzXX8s/gEdN0O09hZPhI6rlC+ANWIdJ1vJfODMXWaQ== sururi@local
command="hg-ssh /path/to/repository",no-port-forwarding,no-X11-forwarding,no-agent-forwarding,no-pty ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAnq5rMLnoab+2F28g/nb58RBENWtX395TuyDFsYkalGaZxrziwDoau/wglkU19DbcAVKgw0p6lMEIuh2iALOppRzxrTgFFhJkL1dxzkugbbPEoSWyfrj9FivzpnxHWgRHQApQeWUBOZhroDTURwfqcyC9SW020CR57jLWfgw+idqwtCu+ZBYmEyHSJcZIH2mWXLrUQ8OalxCFVaLKL50Lpc7V8XJPs+Pg6MPVgfDUqMdjrGkAF7j4viOHTjDWP1h4Ngim70dOeyxWtuqbCbxM4APTShaqET42sj1jHxL2m1dJzXX8s/gEdN0O09hZPhI6rlC+ANWIdJ1vJfODMXWaQ== sururi@local
sururi@local:~/project$ hg pull --ssh "ssh -i /tmp/tmp/id_rsa_example" ssh://sururi@remote//path/to/repo
pulling from ssh://sururi@remote//path/to/repo
searching for changes
no changes found
]]>I had heard of git and bazaar but -for some reasons I don’t know why- never had taken the step to give them a try. Very recently, in a scientific gathering, one of the participants suggested me to use mercurial. The interesting thing about this is, he was suggesting it as a means to enable an “undo” option for a -scientifing- refinement program which has the tendency to update your files in according to the changes you applied. I know this sounds natural but when you tried some method and it didn’t worked, you couldn’t go just one step back and use a different method – you have to start from scratch with your initial input file and apply the methods so far but the last consecutively over again. So, my friend was actually recommending me to use a revision control system for the work files (which can be done perfectly well via SVN, by the way) and he happened to be using mercurial. He suggested to check this website by Joel Spolsky (http://hginit.com/) to get it going and the website is very intriguing, indeed. So I have been running some tests to get my feet wet and it was going so well until I found that the log command didn’t reveal the list of the changes made upon the files per changeset (‘revision’ in terms of SVN). I searched and searched and at the end I found the solution on hg tip (http://hgtip.com/tips/advanced/2010-01-15-styling-mercurials-cli/) website by Steve Losh. His ‘Nice Log’ script was exactly what I was looking for plus some additional effective logging templates.
I don’t think for our big project I can make the switch from SVN to Mercurial, but from now on, I’ll be ‘hg init’ing instead of ‘svnadmin create’ing for the projects to come.
Today, I’ll present a code that I’ve just written to parse out the compatible paths from a tree. In my case, the three is the list of subgroups with indexes that one can acquire using Bilbao Crystallographic Server’s marvelous SUBGROUPGRAPH tool. For low indexes, you can already have SUBGROUPGRAPH do this for you by specifying your supergroup G, subgroup H and the index [G:H]. For example, for G=136, H=14 and [G:H] = 8, you can have the following paths drawn:
But as I said, things can get rough when you have a high index. In that case, we can (and we will) try to parse and analyze the paths tree which is outputted by SUBGROUPGRAPH when no index is designated, e.g., for G=136, H=14:
Where you see the maximal subgroups for the involved groups and their corresponding indexes. So, to go from P4_2/mnm (#136) to P2_1/c (#14) with index 8, we can follow the path:
P4_2/mnm (#136) –[2]–> Cmmm (#65) –[2]–> Pmna (#53) –[2]–> P2_1/c (#14)
with the related indexes given in the brackets totaling to 2x2x2 = 8.
The problem I had today was to find possible paths going from G=Fd-3m (#227) to H=Cc (#9) with index [G:H] = 192.
I first parsed the subgroup list into an array with their indexes, then followed the possible paths.
<?PHP
/* Emre S. Tasci <exxx.txxxx@ehu.es> *
* By parsing the possible subgroup paths obtained from
* SUBGROUPGRAPH, finds the paths that are compatible
* with the given terminal super & sub groups and the
* designated index.
*
* 06/06/11 */
# Input Data === 0 ====================================================
$start_sup = 227;
$crit_index = 192;
$end_sub = 9;
$arr = file("data.txt"); # A sample of data.txt for the case of #227->#9
# is included at the end of the code.
# Input Data === 1 ====================================================
$arr_subs = Array();
$arr_labels = Array();
foreach($arr as $line)
{
$auxarr = preg_split("/[ \t]+/",trim($line));
#echo $auxarr[2]."\n";
$sup = $auxarr[1];
$suplabel = $auxarr[2];
$arrsubs = split(";",$auxarr[3]);
$arr_labels[$sup] = $suplabel;
foreach($arrsubs as $subind)
{
$auxarr2 = Array();
preg_match("/([0-9]+)\[([0-9]+)\]/",$subind,$auxarr2);
#print_r($auxarr2);
$sub = $auxarr2[1];
$index = $auxarr2[2];
$arr_subs[$sup][$index][] = $sub;
}
}
#print_r($arr_subs);
# Start exploring
$start_sup = sprintf("%03d",$start_sup);
$end_sub = sprintf("%03d",$end_sub);
$arr_paths = Array();
$arr_paths[$start_sup] = 1;
findpath($start_sup);
function findpath($current_sup)
{
global $arr_paths,$arr_subs,$crit_index,$end_sub;
#echo "*".$current_sup."*\n";
if(strrpos($current_sup,"-") !== FALSE)
{
$current_sup_last_sup = substr($current_sup,strrpos($current_sup,"-")+1);
$current_sup_last_sup = substr($current_sup_last_sup,0,strpos($current_sup_last_sup,"["));
}
else
$current_sup_last_sup = $current_sup;
#echo "*".$current_sup_last_sup."*\n";
$subindexes = array_keys($arr_subs[$current_sup_last_sup]);
foreach($subindexes as $subindex)
{
$subs = $arr_subs[$current_sup_last_sup][$subindex];
echo $subindex.": ".join(",",$subs)."\n";
foreach($subs as $sub)
{
echo $current_sup."-".$sub."[".$subindex."]"."\n";
$arr_paths[$current_sup."-".$sub."[".$subindex."]"] = $arr_paths[$current_sup] * $subindex;
if($arr_paths[$current_sup."-".$sub."[".$subindex."]"]<$crit_index)
findpath($current_sup."-".$sub."[".$subindex."]");
elseif($arr_paths[$current_sup."-".$sub."[".$subindex."]"] == $crit_index && $sub == $end_sub)
echo "\nPath Found: ".$current_sup."-".$sub."[".$subindex."]"."\n\n";
}
}
#print_r($subindexes);
}
print_r($arr_paths);
/*
data.txt:
1 227 Fd-3m 141[3];166[4];203[2];216[2]
2 220 I-43d 122[3];161[4]
3 219 F-43c 120[3];161[4];218[4]
4 218 P-43n 112[3];161[4];220[4]
5 217 I-43m 121[3];160[4];215[2];218[2]
6 216 F-43m 119[3];160[4];215[4]
7 215 P-43m 111[3];160[4];216[2];217[4];219[2]
8 203 Fd-3 070[3]
9 167 R-3c 015[3];161[2];165[3];167[4];167[5];167[7]
10 166 R-3m 012[3];160[2];164[3];166[2];166[4];166[5];166[7];167[2]
11 165 P-3c1 015[3];158[2];163[3];165[3];165[4];165[5];165[7]
12 164 P-3m1 012[3];156[2];162[3];164[2];164[3];164[4];164[5];164[7];165[2]
13 163 P-31c 015[3];159[2];163[3];163[4];163[5];163[7];165[3];167[3]
14 162 P-31m 012[3];157[2];162[2];162[3];162[4];162[5];162[7];163[2];164[3];166[3]
15 161 R3c 009[3];158[3];161[4];161[5];161[7]
16 160 R3m 008[3];156[3];160[2];160[4];160[5];160[7];161[2]
17 159 P31c 009[3];158[3];159[3];159[4];159[5];159[7];161[3]
18 158 P3c1 009[3];158[3];158[4];158[5];158[7];159[3]
19 157 P31m 008[3];156[3];157[2];157[3];157[4];157[5];157[7];159[2];160[3]
20 156 P3m1 008[3];156[2];156[3];156[4];156[5];156[7];157[3];158[2]
21 141 I41/amd 070[2];074[2];088[2];109[2];119[2];122[2];141[3];141[5];141[7];141[9]
22 122 I-42d 043[2];122[3];122[5];122[7];122[9]
23 121 I-42m 042[2];111[2];112[2];113[2];114[2];121[3];121[5];121[7];121[9]
24 120 I-4c2 045[2];116[2];117[2];120[3];120[5];120[7];120[9]
25 119 I-4m2 044[2];115[2];118[2];119[3];119[5];119[7];119[9]
26 118 P-4n2 034[2];118[3];118[5];118[7];118[9];122[2]
27 117 P-4b2 032[2];117[2];117[3];117[5];117[7];117[9];118[2]
28 116 P-4c2 027[2];112[2];114[2];116[3];116[5];116[7];116[9]
29 115 P-4m2 025[2];111[2];113[2];115[2];115[3];115[5];115[7];115[9];116[2];121[2]
30 114 P-421c 037[2];114[3];114[5];114[7];114[9]
31 113 P-421m 035[2];113[2];113[3];113[5];113[7];113[9];114[2]
32 112 P-42c 037[2];112[3];112[5];112[7];112[9];116[2];118[2]
33 111 P-42m 035[2];111[2];111[3];111[5];111[7];111[9];112[2];115[2];117[2];119[2];120[2]
34 109 I41md 043[2];044[2];109[3];109[5];109[7];109[9]
35 088 I41/a 015[2];088[3];088[5];088[7];088[9]
36 074 Imma 012[2];015[2];044[2];046[2];051[2];052[2];053[2];062[2];074[3];074[5];074[7]
37 070 Fddd 015[2];043[2];070[3];070[5];070[7]
38 064 Cmce 012[2];014[2];015[2];036[2];039[2];041[2];053[2];054[2];055[2];056[2];057[2];060[2];061[2];062[2];064[3];064[5];064[7]
39 063 Cmcm 011[2];012[2];015[2];036[2];038[2];040[2];051[2];052[2];057[2];058[2];059[2];060[2];062[2];063[3];063[5];063[7]
40 062 Pnma 011[2];014[2];026[2];031[2];033[2];062[3];062[5];062[7]
41 061 Pbca 014[2];029[2];061[3];061[5];061[7]
42 060 Pbcn 013[2];014[2];029[2];030[2];033[2];060[3];060[5];060[7]
43 059 Pmmn 011[2];013[2];025[2];031[2];056[2];059[2];059[3];059[5];059[7];062[2]
44 058 Pnnm 010[2];014[2];031[2];034[2];058[3];058[5];058[7]
45 057 Pbcm 011[2];013[2];014[2];026[2];028[2];029[2];057[2];057[3];057[5];057[7];060[2];061[2];062[2]
46 056 Pccn 013[2];014[2];027[2];033[2];056[3];056[5];056[7]
47 055 Pbam 010[2];014[2];026[2];032[2];055[2];055[3];055[5];055[7];058[2];062[2]
48 054 Pcca 013[2];014[2];027[2];029[2];032[2];052[2];054[2];054[3];054[5];054[7];056[2];060[2]
49 053 Pmna 010[2];013[2];014[2];028[2];030[2];031[2];052[2];053[2];053[3];053[5];053[7];058[2];060[2]
50 052 Pnna 013[2];014[2];030[2];033[2];034[2];052[3];052[5];052[7]
51 051 Pmma 010[2];011[2];013[2];025[2];026[2];028[2];051[2];051[3];051[5];051[7];053[2];054[2];055[2];057[2];059[2];063[2];064[2]
52 046 Ima2 008[2];009[2];026[2];028[2];030[2];033[2];046[3];046[5];046[7]
53 045 Iba2 009[2];027[2];029[2];032[2];045[3];045[5];045[7]
54 044 Imm2 008[2];025[2];031[2];034[2];044[3];044[5];044[7]
55 043 Fdd2 009[2];043[3];043[5];043[7]
56 042 Fmm2 008[2];035[2];036[2];037[2];038[2];039[2];040[2];041[2];042[3];042[5];042[7]
57 041 Aea2 007[2];009[2];029[2];030[2];032[2];033[2];041[3];041[5];041[7]
58 040 Ama2 006[2];009[2];028[2];031[2];033[2];034[2];040[3];040[5];040[7]
59 039 Aem2 007[2];008[2];026[2];027[2];028[2];029[2];039[2];039[3];039[5];039[7];041[2];045[2];046[2]
60 038 Amm2 006[2];008[2];025[2];026[2];030[2];031[2];038[2];038[3];038[5];038[7];040[2];044[2];046[2]
61 037 Ccc2 009[2];027[2];030[2];034[2];037[3];037[5];037[7]
62 036 Cmc21 008[2];009[2];026[2];029[2];031[2];033[2];036[3];036[5];036[7]
63 035 Cmm2 008[2];025[2];028[2];032[2];035[2];035[3];035[5];035[7];036[2];037[2];044[2];045[2];046[2]
64 034 Pnn2 007[2];034[3];034[5];034[7];043[2]
65 033 Pna21 007[2];033[3];033[5];033[7]
66 032 Pba2 007[2];032[2];032[3];032[5];032[7];033[2];034[2]
67 031 Pmn21 006[2];007[2];031[2];031[3];031[5];031[7];033[2]
68 030 Pnc2 007[2];030[2];030[3];030[5];030[7];034[2]
69 029 Pca21 007[2];029[2];029[3];029[5];029[7];033[2]
70 028 Pma2 006[2];007[2];028[2];028[3];028[5];028[7];029[2];030[2];031[2];032[2];040[2];041[2]
71 027 Pcc2 007[2];027[2];027[3];027[5];027[7];030[2];037[2]
72 026 Pmc21 006[2];007[2];026[2];026[3];026[5];026[7];029[2];031[2];036[2]
73 025 Pmm2 006[2];025[2];025[3];025[5];025[7];026[2];027[2];028[2];035[2];038[2];039[2];042[2]
74 015 C2/c 009[2];013[2];014[2];015[3];015[5];015[7]
75 014 P21/c 007[2];014[2];014[3];014[5];014[7]
76 013 P2/c 007[2];013[2];013[3];013[5];013[7];014[2];015[2]
77 012 C2/m 008[2];010[2];011[2];012[2];012[3];012[5];012[7];013[2];014[2];015[2]
78 011 P21/m 006[2];011[2];011[3];011[5];011[7];014[2]
79 010 P2/m 006[2];010[2];010[3];010[5];010[7];011[2];012[2];013[2]
80 008 Cm 006[2];007[2];008[2];008[3];008[5];008[7];009[2]
81 007 Pc 007[2];007[3];007[5];007[7];009[2]
82 006 Pm 006[2];006[3];006[5];006[7];007[2];008[2]
83 009 Cc 007[2];009[3];009[5];009[7]
*/
When run, the code lists all the paths with the matching ones preceded by “Path Found”, so if you grep wrt it, you’ll have, for this example:
sururi@husniya:/xxx$ php find_trpath_for_index.php |grep Path|sed "s:Path Found\: ::"
227-141[3]-070[2]-015[2]-009[2]-007[2]-007[2]-009[2]
227-141[3]-070[2]-015[2]-013[2]-007[2]-007[2]-009[2]
227-141[3]-070[2]-015[2]-013[2]-013[2]-007[2]-009[2]
227-141[3]-070[2]-015[2]-013[2]-013[2]-015[2]-009[2]
227-141[3]-070[2]-015[2]-013[2]-014[2]-007[2]-009[2]
227-141[3]-070[2]-015[2]-014[2]-007[2]-007[2]-009[2]
227-141[3]-070[2]-015[2]-014[2]-014[2]-007[2]-009[2]
227-141[3]-070[2]-043[2]-009[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-008[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-008[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-012[2]-008[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-008[2]-008[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-008[2]-008[2]-008[2]-009[2]
227-141[3]-074[2]-012[2]-008[2]-009[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-010[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-010[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-012[2]-010[2]-012[2]-008[2]-009[2]
227-141[3]-074[2]-012[2]-010[2]-012[2]-015[2]-009[2]
227-141[3]-074[2]-012[2]-010[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-010[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-012[2]-011[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-011[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-012[2]-011[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-012[2]-008[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-012[2]-008[2]-008[2]-009[2]
227-141[3]-074[2]-012[2]-012[2]-012[2]-008[2]-009[2]
227-141[3]-074[2]-012[2]-012[2]-012[2]-015[2]-009[2]
227-141[3]-074[2]-012[2]-012[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-012[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-012[2]-012[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-013[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-013[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-013[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-012[2]-013[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-014[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-014[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-015[2]-009[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-015[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-012[2]-015[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-012[2]-015[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-015[2]-009[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-015[2]-013[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-015[2]-013[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-015[2]-013[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-015[2]-013[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-015[2]-014[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-015[2]-014[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-008[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-008[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-008[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-008[2]-008[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-008[2]-008[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-008[2]-009[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-027[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-027[2]-037[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-035[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-035[2]-036[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-035[2]-037[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-035[2]-045[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-035[2]-046[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-038[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-038[2]-040[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-038[2]-046[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-039[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-039[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-039[2]-041[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-039[2]-045[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-039[2]-046[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-042[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-042[2]-036[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-042[2]-037[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-042[2]-040[2]-009[2]
227-141[3]-074[2]-044[2]-025[2]-042[2]-041[2]-009[2]
227-141[3]-074[2]-044[2]-031[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-031[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-044[2]-031[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-031[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-031[2]-033[2]-007[2]-009[2]
227-141[3]-074[2]-044[2]-034[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-008[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-008[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-046[2]-008[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-008[2]-008[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-008[2]-008[2]-008[2]-009[2]
227-141[3]-074[2]-046[2]-008[2]-009[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-009[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-026[2]-036[2]-008[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-032[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-028[2]-041[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-030[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-030[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-030[2]-034[2]-007[2]-009[2]
227-141[3]-074[2]-046[2]-030[2]-034[2]-043[2]-009[2]
227-141[3]-074[2]-046[2]-033[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-010[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-010[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-010[2]-012[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-010[2]-012[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-010[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-010[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-011[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-011[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-011[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-013[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-013[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-013[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-013[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-027[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-027[2]-037[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-035[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-035[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-035[2]-037[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-035[2]-045[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-035[2]-046[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-038[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-038[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-038[2]-046[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-039[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-039[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-039[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-039[2]-045[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-039[2]-046[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-042[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-042[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-042[2]-037[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-042[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-025[2]-042[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-026[2]-036[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-032[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-028[2]-041[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-063[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-063[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-063[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-064[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-064[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-051[2]-064[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-053[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-054[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-054[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-054[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-054[2]-027[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-054[2]-027[2]-037[2]-009[2]
227-141[3]-074[2]-051[2]-054[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-054[2]-032[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-055[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-055[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-055[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-055[2]-032[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-057[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-059[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-059[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-059[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-063[2]-012[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-063[2]-012[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-063[2]-036[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-063[2]-038[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-063[2]-038[2]-040[2]-009[2]
227-141[3]-074[2]-051[2]-063[2]-038[2]-046[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-012[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-012[2]-015[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-036[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-039[2]-007[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-039[2]-008[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-039[2]-041[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-039[2]-045[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-039[2]-046[2]-009[2]
227-141[3]-074[2]-051[2]-064[2]-041[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-013[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-013[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-013[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-052[2]-013[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-014[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-014[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-030[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-030[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-030[2]-034[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-030[2]-034[2]-043[2]-009[2]
227-141[3]-074[2]-052[2]-033[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-052[2]-034[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-010[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-010[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-053[2]-010[2]-012[2]-008[2]-009[2]
227-141[3]-074[2]-053[2]-010[2]-012[2]-015[2]-009[2]
227-141[3]-074[2]-053[2]-010[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-010[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-053[2]-013[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-013[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-013[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-053[2]-013[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-014[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-014[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-032[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-028[2]-041[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-030[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-030[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-030[2]-034[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-030[2]-034[2]-043[2]-009[2]
227-141[3]-074[2]-053[2]-031[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-031[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-053[2]-031[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-031[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-031[2]-033[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-052[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-052[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-053[2]-052[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-052[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-052[2]-033[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-052[2]-034[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-052[2]-034[2]-043[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-028[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-028[2]-040[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-028[2]-041[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-053[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-058[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-058[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-058[2]-034[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-058[2]-034[2]-043[2]-009[2]
227-141[3]-074[2]-053[2]-060[2]-013[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-060[2]-013[2]-015[2]-009[2]
227-141[3]-074[2]-053[2]-060[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-060[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-060[2]-030[2]-007[2]-009[2]
227-141[3]-074[2]-053[2]-060[2]-033[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-011[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-011[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-062[2]-011[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-014[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-014[2]-014[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-026[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-026[2]-036[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-029[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-026[2]-036[2]-008[2]-009[2]
227-141[3]-074[2]-062[2]-031[2]-006[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-031[2]-006[2]-008[2]-009[2]
227-141[3]-074[2]-062[2]-031[2]-007[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-031[2]-031[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-031[2]-033[2]-007[2]-009[2]
227-141[3]-074[2]-062[2]-033[2]-007[2]-007[2]-009[2]
227-141[3]-088[2]-015[2]-009[2]-007[2]-007[2]-009[2]
227-141[3]-088[2]-015[2]-013[2]-007[2]-007[2]-009[2]
227-141[3]-088[2]-015[2]-013[2]-013[2]-007[2]-009[2]
227-141[3]-088[2]-015[2]-013[2]-013[2]-015[2]-009[2]
227-141[3]-088[2]-015[2]-013[2]-014[2]-007[2]-009[2]
227-141[3]-088[2]-015[2]-014[2]-007[2]-007[2]-009[2]
227-141[3]-088[2]-015[2]-014[2]-014[2]-007[2]-009[2]
227-141[3]-109[2]-043[2]-009[2]-007[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-008[2]-006[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-008[2]-006[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-008[2]-007[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-008[2]-008[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-008[2]-008[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-008[2]-009[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-006[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-006[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-026[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-026[2]-036[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-027[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-027[2]-037[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-028[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-028[2]-040[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-028[2]-041[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-035[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-035[2]-036[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-035[2]-037[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-035[2]-045[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-035[2]-046[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-038[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-038[2]-040[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-038[2]-046[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-039[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-039[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-039[2]-041[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-039[2]-045[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-039[2]-046[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-042[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-042[2]-036[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-042[2]-037[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-042[2]-040[2]-009[2]
227-141[3]-109[2]-044[2]-025[2]-042[2]-041[2]-009[2]
227-141[3]-109[2]-044[2]-031[2]-006[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-031[2]-006[2]-008[2]-009[2]
227-141[3]-109[2]-044[2]-031[2]-007[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-031[2]-031[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-031[2]-033[2]-007[2]-009[2]
227-141[3]-109[2]-044[2]-034[2]-007[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-008[2]-006[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-008[2]-006[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-008[2]-007[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-008[2]-008[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-008[2]-008[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-008[2]-009[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-006[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-006[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-026[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-026[2]-036[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-027[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-027[2]-037[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-028[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-028[2]-040[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-028[2]-041[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-035[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-035[2]-036[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-035[2]-037[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-035[2]-045[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-035[2]-046[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-038[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-038[2]-040[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-038[2]-046[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-039[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-039[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-039[2]-041[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-039[2]-045[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-039[2]-046[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-042[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-042[2]-036[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-042[2]-037[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-042[2]-040[2]-009[2]
227-141[3]-119[2]-044[2]-025[2]-042[2]-041[2]-009[2]
227-141[3]-119[2]-044[2]-031[2]-006[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-031[2]-006[2]-008[2]-009[2]
227-141[3]-119[2]-044[2]-031[2]-007[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-031[2]-031[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-031[2]-033[2]-007[2]-009[2]
227-141[3]-119[2]-044[2]-034[2]-007[2]-007[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-006[2]-007[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-006[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-026[2]-007[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-026[2]-036[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-027[2]-007[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-027[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-028[2]-007[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-028[2]-040[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-028[2]-041[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-035[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-035[2]-036[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-035[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-035[2]-045[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-035[2]-046[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-038[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-038[2]-040[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-038[2]-046[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-039[2]-007[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-039[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-039[2]-041[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-039[2]-045[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-039[2]-046[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-042[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-042[2]-036[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-042[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-042[2]-040[2]-009[2]
227-141[3]-119[2]-115[2]-025[2]-042[2]-041[2]-009[2]
227-141[3]-119[2]-115[2]-111[2]-035[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-111[2]-035[2]-036[2]-009[2]
227-141[3]-119[2]-115[2]-111[2]-035[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-111[2]-035[2]-045[2]-009[2]
227-141[3]-119[2]-115[2]-111[2]-035[2]-046[2]-009[2]
227-141[3]-119[2]-115[2]-111[2]-112[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-111[2]-120[2]-045[2]-009[2]
227-141[3]-119[2]-115[2]-113[2]-035[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-113[2]-035[2]-036[2]-009[2]
227-141[3]-119[2]-115[2]-113[2]-035[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-113[2]-035[2]-045[2]-009[2]
227-141[3]-119[2]-115[2]-113[2]-035[2]-046[2]-009[2]
227-141[3]-119[2]-115[2]-113[2]-114[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-116[2]-027[2]-007[2]-009[2]
227-141[3]-119[2]-115[2]-116[2]-027[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-116[2]-112[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-116[2]-114[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-121[2]-042[2]-008[2]-009[2]
227-141[3]-119[2]-115[2]-121[2]-042[2]-036[2]-009[2]
227-141[3]-119[2]-115[2]-121[2]-042[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-121[2]-042[2]-040[2]-009[2]
227-141[3]-119[2]-115[2]-121[2]-042[2]-041[2]-009[2]
227-141[3]-119[2]-115[2]-121[2]-112[2]-037[2]-009[2]
227-141[3]-119[2]-115[2]-121[2]-114[2]-037[2]-009[2]
227-141[3]-119[2]-118[2]-034[2]-007[2]-007[2]-009[2]
227-141[3]-122[2]-043[2]-009[2]-007[2]-007[2]-009[2]
227-166[4]-012[3]-008[2]-006[2]-007[2]-009[2]
227-166[4]-012[3]-008[2]-006[2]-008[2]-009[2]
227-166[4]-012[3]-008[2]-007[2]-007[2]-009[2]
227-166[4]-012[3]-008[2]-008[2]-007[2]-009[2]
227-166[4]-012[3]-008[2]-008[2]-008[2]-009[2]
227-166[4]-012[3]-008[2]-009[2]-007[2]-009[2]
227-166[4]-012[3]-010[2]-006[2]-007[2]-009[2]
227-166[4]-012[3]-010[2]-006[2]-008[2]-009[2]
227-166[4]-012[3]-010[2]-012[2]-008[2]-009[2]
227-166[4]-012[3]-010[2]-012[2]-015[2]-009[2]
227-166[4]-012[3]-010[2]-013[2]-007[2]-009[2]
227-166[4]-012[3]-010[2]-013[2]-015[2]-009[2]
227-166[4]-012[3]-011[2]-006[2]-007[2]-009[2]
227-166[4]-012[3]-011[2]-006[2]-008[2]-009[2]
227-166[4]-012[3]-011[2]-014[2]-007[2]-009[2]
227-166[4]-012[3]-012[2]-008[2]-007[2]-009[2]
227-166[4]-012[3]-012[2]-008[2]-008[2]-009[2]
227-166[4]-012[3]-012[2]-012[2]-008[2]-009[2]
227-166[4]-012[3]-012[2]-012[2]-015[2]-009[2]
227-166[4]-012[3]-012[2]-013[2]-007[2]-009[2]
227-166[4]-012[3]-012[2]-013[2]-015[2]-009[2]
227-166[4]-012[3]-012[2]-014[2]-007[2]-009[2]
227-166[4]-012[3]-013[2]-007[2]-007[2]-009[2]
227-166[4]-012[3]-013[2]-013[2]-007[2]-009[2]
227-166[4]-012[3]-013[2]-013[2]-015[2]-009[2]
227-166[4]-012[3]-013[2]-014[2]-007[2]-009[2]
227-166[4]-012[3]-014[2]-007[2]-007[2]-009[2]
227-166[4]-012[3]-014[2]-014[2]-007[2]-009[2]
227-166[4]-012[3]-015[2]-009[2]-007[2]-009[2]
227-166[4]-012[3]-015[2]-013[2]-007[2]-009[2]
227-166[4]-012[3]-015[2]-013[2]-015[2]-009[2]
227-166[4]-012[3]-015[2]-014[2]-007[2]-009[2]
227-166[4]-160[2]-008[3]-006[2]-007[2]-009[2]
227-166[4]-160[2]-008[3]-006[2]-008[2]-009[2]
227-166[4]-160[2]-008[3]-007[2]-007[2]-009[2]
227-166[4]-160[2]-008[3]-008[2]-007[2]-009[2]
227-166[4]-160[2]-008[3]-008[2]-008[2]-009[2]
227-166[4]-160[2]-008[3]-009[2]-007[2]-009[2]
227-166[4]-160[2]-160[2]-008[3]-007[2]-009[2]
227-166[4]-160[2]-160[2]-008[3]-008[2]-009[2]
227-166[4]-160[2]-160[2]-160[2]-008[3]-009[2]
227-166[4]-160[2]-160[2]-160[2]-161[2]-009[3]
227-166[4]-160[2]-161[2]-009[3]-007[2]-009[2]
227-166[4]-160[2]-161[2]-161[4]-009[3]
227-166[4]-160[2]-160[4]-008[3]-009[2]
227-166[4]-160[2]-160[4]-161[2]-009[3]
227-166[4]-166[2]-012[3]-008[2]-007[2]-009[2]
227-166[4]-166[2]-012[3]-008[2]-008[2]-009[2]
227-166[4]-166[2]-012[3]-012[2]-008[2]-009[2]
227-166[4]-166[2]-012[3]-012[2]-015[2]-009[2]
227-166[4]-166[2]-012[3]-013[2]-007[2]-009[2]
227-166[4]-166[2]-012[3]-013[2]-015[2]-009[2]
227-166[4]-166[2]-012[3]-014[2]-007[2]-009[2]
227-166[4]-166[2]-160[2]-008[3]-007[2]-009[2]
227-166[4]-166[2]-160[2]-008[3]-008[2]-009[2]
227-166[4]-166[2]-160[2]-160[2]-008[3]-009[2]
227-166[4]-166[2]-160[2]-160[2]-161[2]-009[3]
227-166[4]-166[2]-166[2]-012[3]-008[2]-009[2]
227-166[4]-166[2]-166[2]-012[3]-015[2]-009[2]
227-166[4]-166[2]-166[2]-160[2]-008[3]-009[2]
227-166[4]-166[2]-166[2]-160[2]-161[2]-009[3]
227-166[4]-166[2]-166[2]-167[2]-015[3]-009[2]
227-166[4]-166[2]-166[2]-167[2]-161[2]-009[3]
227-166[4]-167[2]-015[3]-009[2]-007[2]-009[2]
227-166[4]-167[2]-015[3]-013[2]-007[2]-009[2]
227-166[4]-167[2]-015[3]-013[2]-015[2]-009[2]
227-166[4]-167[2]-015[3]-014[2]-007[2]-009[2]
227-166[4]-167[2]-161[2]-009[3]-007[2]-009[2]
227-166[4]-167[2]-161[2]-161[4]-009[3]
227-166[4]-167[2]-167[4]-015[3]-009[2]
227-166[4]-167[2]-167[4]-161[2]-009[3]
227-166[4]-166[4]-012[3]-008[2]-009[2]
227-166[4]-166[4]-012[3]-015[2]-009[2]
227-166[4]-166[4]-160[2]-008[3]-009[2]
227-166[4]-166[4]-160[2]-161[2]-009[3]
227-166[4]-166[4]-167[2]-015[3]-009[2]
227-166[4]-166[4]-167[2]-161[2]-009[3]
227-203[2]-070[3]-015[2]-009[2]-007[2]-007[2]-009[2]
227-203[2]-070[3]-015[2]-013[2]-007[2]-007[2]-009[2]
227-203[2]-070[3]-015[2]-013[2]-013[2]-007[2]-009[2]
227-203[2]-070[3]-015[2]-013[2]-013[2]-015[2]-009[2]
227-203[2]-070[3]-015[2]-013[2]-014[2]-007[2]-009[2]
227-203[2]-070[3]-015[2]-014[2]-007[2]-007[2]-009[2]
227-203[2]-070[3]-015[2]-014[2]-014[2]-007[2]-009[2]
227-203[2]-070[3]-043[2]-009[2]-007[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-008[2]-006[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-008[2]-006[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-008[2]-007[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-008[2]-008[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-008[2]-008[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-008[2]-009[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-006[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-006[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-026[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-026[2]-036[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-027[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-027[2]-037[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-028[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-028[2]-040[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-028[2]-041[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-035[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-035[2]-036[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-035[2]-037[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-035[2]-045[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-035[2]-046[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-038[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-038[2]-040[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-038[2]-046[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-039[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-039[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-039[2]-041[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-039[2]-045[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-039[2]-046[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-042[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-042[2]-036[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-042[2]-037[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-042[2]-040[2]-009[2]
227-216[2]-119[3]-044[2]-025[2]-042[2]-041[2]-009[2]
227-216[2]-119[3]-044[2]-031[2]-006[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-031[2]-006[2]-008[2]-009[2]
227-216[2]-119[3]-044[2]-031[2]-007[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-031[2]-031[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-031[2]-033[2]-007[2]-009[2]
227-216[2]-119[3]-044[2]-034[2]-007[2]-007[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-006[2]-007[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-006[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-026[2]-007[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-026[2]-036[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-027[2]-007[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-027[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-028[2]-007[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-028[2]-040[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-028[2]-041[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-035[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-035[2]-036[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-035[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-035[2]-045[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-035[2]-046[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-038[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-038[2]-040[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-038[2]-046[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-039[2]-007[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-039[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-039[2]-041[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-039[2]-045[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-039[2]-046[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-042[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-042[2]-036[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-042[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-042[2]-040[2]-009[2]
227-216[2]-119[3]-115[2]-025[2]-042[2]-041[2]-009[2]
227-216[2]-119[3]-115[2]-111[2]-035[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-111[2]-035[2]-036[2]-009[2]
227-216[2]-119[3]-115[2]-111[2]-035[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-111[2]-035[2]-045[2]-009[2]
227-216[2]-119[3]-115[2]-111[2]-035[2]-046[2]-009[2]
227-216[2]-119[3]-115[2]-111[2]-112[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-111[2]-120[2]-045[2]-009[2]
227-216[2]-119[3]-115[2]-113[2]-035[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-113[2]-035[2]-036[2]-009[2]
227-216[2]-119[3]-115[2]-113[2]-035[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-113[2]-035[2]-045[2]-009[2]
227-216[2]-119[3]-115[2]-113[2]-035[2]-046[2]-009[2]
227-216[2]-119[3]-115[2]-113[2]-114[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-116[2]-027[2]-007[2]-009[2]
227-216[2]-119[3]-115[2]-116[2]-027[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-116[2]-112[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-116[2]-114[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-121[2]-042[2]-008[2]-009[2]
227-216[2]-119[3]-115[2]-121[2]-042[2]-036[2]-009[2]
227-216[2]-119[3]-115[2]-121[2]-042[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-121[2]-042[2]-040[2]-009[2]
227-216[2]-119[3]-115[2]-121[2]-042[2]-041[2]-009[2]
227-216[2]-119[3]-115[2]-121[2]-112[2]-037[2]-009[2]
227-216[2]-119[3]-115[2]-121[2]-114[2]-037[2]-009[2]
227-216[2]-119[3]-118[2]-034[2]-007[2]-007[2]-009[2]
227-216[2]-160[4]-008[3]-006[2]-007[2]-009[2]
227-216[2]-160[4]-008[3]-006[2]-008[2]-009[2]
227-216[2]-160[4]-008[3]-007[2]-007[2]-009[2]
227-216[2]-160[4]-008[3]-008[2]-007[2]-009[2]
227-216[2]-160[4]-008[3]-008[2]-008[2]-009[2]
227-216[2]-160[4]-008[3]-009[2]-007[2]-009[2]
227-216[2]-160[4]-160[2]-008[3]-007[2]-009[2]
227-216[2]-160[4]-160[2]-008[3]-008[2]-009[2]
227-216[2]-160[4]-160[2]-160[2]-008[3]-009[2]
227-216[2]-160[4]-160[2]-160[2]-161[2]-009[3]
227-216[2]-160[4]-161[2]-009[3]-007[2]-009[2]
227-216[2]-160[4]-161[2]-161[4]-009[3]
227-216[2]-160[4]-160[4]-008[3]-009[2]
227-216[2]-160[4]-160[4]-161[2]-009[3]
227-216[2]-215[4]-111[3]-035[2]-008[2]-009[2]
227-216[2]-215[4]-111[3]-035[2]-036[2]-009[2]
227-216[2]-215[4]-111[3]-035[2]-037[2]-009[2]
227-216[2]-215[4]-111[3]-035[2]-045[2]-009[2]
227-216[2]-215[4]-111[3]-035[2]-046[2]-009[2]
227-216[2]-215[4]-111[3]-112[2]-037[2]-009[2]
227-216[2]-215[4]-111[3]-120[2]-045[2]-009[2]
227-216[2]-215[4]-160[4]-008[3]-009[2]
227-216[2]-215[4]-160[4]-161[2]-009[3]
227-216[2]-215[4]-219[2]-120[3]-045[2]-009[2]
227-216[2]-215[4]-219[2]-161[4]-009[3]
]]>I started writing stories long before I started writing scientific articles and had this story of mine, “The Most Beautiful Sun of Our Earth” written some 9 years ago (in Turkish). Last week, the idea to re-write/revise it in English and submit it to the Futures occurred to me, so on an insomniatic night (or whatever is the correct English term), I sat down in front of the word processor and pulled it out.
The original story had me and -my then girlfriend- Bengü as the protagonists but I felt that it would be too selfish so, I “borrowed” my dear friend (my battalion as we call it), Andy and his girlfriend Serene for the part. The following day I asked Andy’s permission for the inclusion and also his help to edit this English version of the story, which he kindly accepted.
At the end, I submitted the story to Nature and from the tone of this entry most probably you’ll think that it was accepted for publication and hoorraay and all that but no, yesterday, I got a reply stating that, it wasn’t accepted and that’s it.
The good thing coming out of this is that, Andy actually liked the story very much and this story being one of my favorite SF stories, I’m really happy for it. So, if you’re still interested, you can get a copy of it, completely free from any possible obligations to Nature publishing were it to have been otherwise.. 8)
Emre S. Tasci – “The Most Beautiful Sun of Our Earth”
]]>Recently, while searching recent publications on liquid alloys, I’ve found a very good and resourceful article on liquid Na-K alloys by Dr. J.-F. Wax (“Molecular dynamics study of the static structure of liquid Na-K alloys”, Physica B 403 (2008) 4241-48 | doi:10.1016/j.physb.2008.09.014). In this paper, among other properties, he had also presented the partial pair distribution functions. Partial-pair distribution functions is the averaged probability that you’ll encounter an atom in r distance starting from another atom. In solid phase, due to the symmetry of the crystal structure, you have discrete values instead of a continuous distribution and everything is crystal clear but I guess that’s the price we pay dealing with liquids. 8)
Partial-pair distribution functions are very useful in the sense that they let you derive the bond-length between constituent atom species, corresponding to the first minimums of the plots. Hence, one can see the neighborhood ranges by just looking at the plots.
I’ve sent a mail to Dr. Wax explaining my research topics and my interest and asking for his research data to which he very kindly and generously replied by providing his data.
He had sent me a big file with numerous configurations (coordinates of the atoms) following each other. The first -and the simplest- task was to split the file into seperate configuration files via this bash script:
#!/bin/bash
# Script Name: split_POS4.sh
# Splits the NaK conf files compilation POS4_DAT into sepearate
# conf files.
# Emre S. Tasci <e.tasci@tudelft.nl>
# 01/09/09
linestart=2
p="p"
for (( i = 1; i <= 1000; i++ ))
do
echo $i;
lineend=$(expr $linestart + 2047)
confname=$(printf "%04d" $i)
echo -e "3\n2048" > $confname
sed -n "$linestart,$(echo $lineend$p)" POS4_DAT >> $confname
linestart=$(expr $lineend + 1)
cat $confname | qhull i TO $confname.hull
done
As you can check in the line before the “done” closer, I’ve -proudly- used Qhull software to calculate the convex hull of each configuration. A convex hull is the -hopefully minimum- shape that covers all your system. So, for each configuration I now had 2 files: “xxxx” (where xxxx is the id/number of the configuration set) storing the coordinates (preceded by 3 and 2048, corresponding to dimension and number of atoms information for the qhull to parse & process) and “xxxx.hull” file, generated by the qhull, containing the vertices list of each facet (preceded by the total number of facets).
A facet is (in 3D) the plane formed by the vertice points. Imagine a cube, where an atom lies at each corner, 3 in the inside. So, we can say that our system consists of 11 atoms and the minimal shape that has edges selected from system’s atoms and covering the whole is a cube, with the edges being the 8 atoms at the corners. Now, add an atom floating over the face at the top. Now the convex hull that wraps our new system would no longer be a 6-faced, 8-edged cube but a 9-faced, 9-edged polygon. I need to know the convex hull geometry in order to correctly calculate the coordination number distributions (details will follow shortly).
The aforementioned two files look like:
sururi@husniya hull_calcs $ head 0001
3
2048
0.95278770E+00 0.13334565E+02 0.13376155E+02
0.13025256E+02 0.87618381E+00 0.20168993E+01
0.12745038E+01 0.14068998E-01 0.22887323E+01
0.13066590E+02 0.20788591E+01 0.58119183E+01
0.10468218E+00 0.58523640E-01 0.64288678E+01
0.12839412E+02 0.13117012E+02 0.79093881E+01
0.11918105E+01 0.12163854E+02 0.10270868E+02
0.12673985E+02 0.11642538E+02 0.10514597E+02
sururi@husniya hull_calcs $ head 0001.hull
180
568 1127 1776
1992 104 1551
956 449 1026
632 391 1026
391 632 1543
391 956 1026
966 632 1026
966 15 1039
632 966 1039
To calculate the coordination number distribution, you need to pick each atom and then simply count the other atoms within the cut-off distance depending on your atom’s and the neighboring atom’s species. After you do the counting for every atom, you average the results and voilà ! there is your CN distribution. But here’s the catch: In order to be able to count properly, you must make sure that, your reference atom’s cutoff radii stay within your system — otherwise you’ll be undercounting. Suppose your reference atom is one of the atoms at the corners of the cube: It will (generally/approximately) lack 7/8 of the neighbouring atoms it would normally have if it was to be the center atom of the specimen/sample. Remember that we are studying a continuous liquid and trying to model the system via ensembles of samples. So we should omit these atoms which has neighborhoods outside of our systems. Since neighborhoods are defined through bond-length, we should only include the atoms with a distance to the outside of at least cut-off radii. The “outside” is defined as the region outside the convex hull and it’s designated by the facets (planes).
From the xxxx.hull file, we have the ids of the vertice atoms and their coordinates are stored in the xxxx file. To define the plane equation
we will use the three vertice points, say p1, p2 and p3. From these 3 points, we calculate the normal n as:
then, the plane equation (hence the coefficients a,b,c,d) can be derived by comparing the following equation:
Here, (x0,y0,z0) are just the coordinates of any of the vertice points. And finally, the distance D between a point p1(x1,y1,z1) and a plane (a,b,c,d) is given by:
Since vector operations are involved, I thought it’s better to do this job in Octave and wrote the following script that calculates the coefficients of the facets and then sets down to find out the points whose distances to all the facets are greater than the given cut-off radius:
sururi@husniya trunk $ cat octave_find_atoms_within.m
DummyA = 3;
%conf_index=292;
% Call from the command line like :
% octave -q --eval "Rcut=1.65;conf_index=293;path='$(pwd)';" octave_find_atoms_within.m
% hence specifying the conf_index at the command line.
% Emre S. Tasci <e.tasci@tudelft.nl> *
%
% From the positions and vertices file, calculates the plane
% equations (stored in "facet_coefficients") and then
% filters the atoms with respect to their distances to these
% facets. Writes the output to
% "hull_calcs/xxxx_insiders.txt" file
% 03/09/09 */
conf_index_name = num2str(sprintf("%04d",conf_index));;
%clock
clear facet_coefficients;
facet_coefficients = [];
global p v facet_coefficients
fp = fopen(strcat(path,"/hull_calcs/",conf_index_name));
k = fscanf(fp,"%d",2); % 1st is dim (=3), 2nd is number of atoms (=2048)
p = fscanf(fp,"%f",[k(1),k(2)]);
p = p';
fclose(fp);
%p = load("/tmp/del/delco");
fp = fopen(strcat(path,"/hull_calcs/",conf_index_name,".hull"));
k = fscanf(fp,"%d",1);% number of facets
v = fscanf(fp,"%d",[3,k]);
v = v';
fclose(fp);
%v = load("/tmp/del/delver");
function [a,b,c,d,e] = getPlaneCoefficients(facet_index)
global p v
% Get the 3 vertices
p1 = p(v(facet_index,1)+1,:); % The +1s come from the fact
p2 = p(v(facet_index,2)+1,:); % that qhull enumeration begins
p3 = p(v(facet_index,3)+1,:); % from 0 while Octave's from 1
%printf("%d: %f %f %f %f %f %f %f %f %f\n",facet_index,p1,p2,p3);
% Calculate the normal
n = cross((p2-p1),(p3-p1));
% Calculate the coefficients of the plane :
% ax + by + cz +d = 0
a = n(1);
b = n(2);
c = n(3);
d = -(dot(n,p1));
e = norm(n);
if(e==0)
printf("%f\n",p1);
printf("%f\n",p2);
printf("%f\n",p3);
printf("%d\n",facet_index);
endif
endfunction
function dist = distanceToPlane(point_index,facet_index)
global p facet_coefficients
k = facet_coefficients(facet_index,:);
n = [k(1) k(2) k(3)];
dist = abs(dot(n,p(point_index,:)) + k(4)) / k(5);
endfunction
for i = [1:size(v)(1)]
[a,b,c,d,e] = getPlaneCoefficients(i);
facet_coefficients = [ facet_coefficients ; [a,b,c,d,e]];
endfor
%Rcut = 1.65; % Defined from command line calling
% Find the points that are not closer than Rcut
inside_atoms = [];
for point = [1:size(p)(1)]
%for point = [1:100]
inside=true;
for facet = [1:size(v)(1)]
%for facet = [1:10]
%printf("%d - %d : %f\n",point,facet,dist);
if (distanceToPlane(point,facet)<Rcut)
inside = false;
break;
endif
endfor
if(inside)
inside_atoms = [inside_atoms point-1];
endif
endfor
fp = fopen(strcat(path,"/hull_calcs/",conf_index_name,"_insiders.txt"),"w");
fprintf(fp,"%d\n",inside_atoms);
fclose(fp);
%clock
<?PHP
/* Emre S. Tasci <e.tasci@tudelft.nl> *
*
* parse_configuration_data.php
*
* Written in order to parse the configuation files
* obtained from J.F. Wax in order to calculate the
* Coordination Number distribution. Each configuration
* file contain 2048 atoms' coordinates starting at the
* 3rd line. There is also a convex hull file generated
* using qhull (http://www.qhull.com) that contains the
* indexes (starting from 0) of the atoms that form the
* vertices of the convex hull. Finally there is the
* IP.dat file that identifies each atom (-1 for K; 1 for
* Na -- the second column is the relative masses).
*
* 02/09/09 */
/*
# if present, remove the previous run's results file
if(file_exists("results.txt"))
unlink("results.txt");
*/
error_reporting(E_ALL ^ E_NOTICE);
if(!file_exists("results.txt"))
{
$fp = fopen("results.txt","w");
fwrite($fp,"CONF#"."\tNa".implode("\tNa",range(1,20))."\tK".implode("\tK",range(1,20))."\n");
fclose($fp);
}
# Support for command line variable passing:
for($i=1;$i<sizeof($_SERVER["argv"]);$i++)
{
list($var0,$val0) = explode("=",$_SERVER["argv"][$i]);
$_GET[$var0] = $val0;
}
if($_GET["configuration"])
calculate($_GET["configuration"]);
else
exit("Please specify a configuration number [php parse_configuration_data.php configuration=xxx]\n");
function calculate($conf_index)
{
# Define the atoms array
$arr_atoms = Array();
# Get the information from the files
parse_types($arr_atoms);
$refs = get_vertices($conf_index);
parse_coordinates($arr_atoms,$conf_index);
# Define the Rcut-off values (obtained from partial parid distribution minimums)
$Rcut_arr = Array();
$Rscale = 3.828; # To convert reduced distances to Angstrom (if needed)
$Rscale = 1; # To convert reduced distances to Angstrom
$Rcut_arr["Na"]["Na"] = 1.38 * $Rscale;
$Rcut_arr["Na"]["K"] = 1.65 * $Rscale;
$Rcut_arr["K"]["Na"] = 1.65 * $Rscale;
$Rcut_arr["K"]["K"] = 1.52 * $Rscale;
$Rcut = $Rcut_arr["Na"]["K"]; # Taking the maximum one as the Rcut for safety
/*
# We have everything we need, so proceed to check
# if an atom is at safe distance wrt to the vertices
# Define all the ids
$arr_test_ids = range(0,2047);
# Subtract the ref atoms' ids
$arr_test_ids = array_diff($arr_test_ids, $refs);
sort($arr_test_ids);
//print_r($arr_test_ids);
$arr_passed_atom_ids = Array();
# Check each atom against refs
for($id=0, $size=sizeof($arr_test_ids); $id<$size; $id++)
{
//echo $id."\n";
$arr_atoms[$arr_test_ids[$id]]["in"] = TRUE;
foreach ($refs as $ref)
{
$distance = distance($arr_atoms[$arr_test_ids[$id]],$arr_atoms[$ref]);
//echo "\t".$ref."\t".$distance."\n";
if($distance < $Rcut)
{
$arr_atoms[$arr_test_ids[$id]]["in"] = FALSE;
break;
}
}
if($arr_atoms[$arr_test_ids[$id]]["in"] == TRUE)
$arr_passed_atom_ids[] = $arr_test_ids[$id];
}
*/
# Run the octave script file to calculate the inside atoms
if(!file_exists("hull_calcs/".sprintf("%04d",$conf_index)."_insiders.txt"))
exec("octave -q --eval \"Rcut=".$Rcut.";conf_index=".$conf_index.";path='$(pwd)';\" octave_find_atoms_within.m");
# Read the file generated by Octave
$arr_passed_atom_ids = file("hull_calcs/".sprintf("%04d",$conf_index)."_insiders.txt",FILE_IGNORE_NEW_LINES);
$arr_test_ids = range(0,2047);
$arr_test_ids = array_diff($arr_test_ids, $refs);
sort($arr_test_ids);
for($i=0, $size=sizeof($arr_test_ids);$i<$size;$i++)
$arr_atoms[$arr_test_ids[$i]]["in"]=FALSE;
# Begin checking every passed atom
foreach($arr_passed_atom_ids as $passed_atom_id)
{
$arr_atoms[$passed_atom_id]["in"] = TRUE;
for($i=0, $size=sizeof($arr_atoms); $i<$size; $i++)
{
$distance = distance($arr_atoms[$passed_atom_id],$arr_atoms[$i]);
//echo $passed_atom_id."\t---\t".$i."\n";
if($distance < $Rcut_arr[$arr_atoms[$passed_atom_id]["id"]][$arr_atoms[$i]["id"]] && $distance>0.001)
$arr_atoms[$passed_atom_id]["neighbour_count"] += 1;
}
}
# Do the binning
$CN_Na = Array();
$CN_K = Array();
for($i=0, $size=sizeof($arr_atoms); $i<$size; $i++)
{
if(array_key_exists("neighbour_count",$arr_atoms[$i]))
{
${"CN_".$arr_atoms[$i]['id']}[$arr_atoms[$i]["neighbour_count"]] += 1;
}
}
ksort($CN_Na);
ksort($CN_K);
//print_r($arr_atoms);
//print_r($CN_Na);
//print_r($CN_K);
# Report the results
$filename = "results/".sprintf("%04d",$conf_index)."_results.txt";
$fp = fopen($filename,"w");
fwrite($fp, "### Atoms array ###\n");
fwrite($fp,var_export($arr_atoms,TRUE)."\n");
fwrite($fp, "\n### CN Distribution for Na ###\n");
fwrite($fp,var_export($CN_Na,TRUE)."\n");
fwrite($fp, "\n### CN Distribution for K ###\n");
fwrite($fp,var_export($CN_K,TRUE)."\n");
# Percentage calculation:
$sum_Na = array_sum($CN_Na);
$sum_K = array_sum($CN_K);
foreach($CN_Na as $key=>$value)
$CN_Na[$key] = $value * 100 / $sum_Na;
foreach($CN_K as $key=>$value)
$CN_K[$key] = $value * 100 / $sum_K;
//print_r($CN_Na);
//print_r($CN_K);
fwrite($fp, "\n### CN Distribution (Percentage) for Na ###\n");
fwrite($fp,var_export($CN_Na,TRUE)."\n");
fwrite($fp, "\n### CN Distribution (Percentage) for K ###\n");
fwrite($fp,var_export($CN_K,TRUE)."\n");
fclose($fp);
# Write the summary to the results file
$fp = fopen("results.txt","a");
for($i=1;$i<=20;$i++)
{
if(!array_key_exists($i,$CN_Na))
$CN_Na[$i] = 0;
if(!array_key_exists($i,$CN_K))
$CN_K[$i] = 0;
}
ksort($CN_Na);
ksort($CN_K);
fwrite($fp,sprintf("%04d",$conf_index)."\t".implode("\t",$CN_Na)."\t".implode("\t",$CN_K)."\n");
fclose($fp);
}
function parse_types(&$arr_atoms)
{
# Parse the types
$i = 0;
$fp = fopen("IP.DAT", "r");
while(!feof($fp))
{
$line = fgets($fp,4096);
$auxarr = explode(" ",$line);
if($auxarr[0]==-1)
$arr_atoms[$i]["id"] = "Na";
else
$arr_atoms[$i]["id"] = "K";
$i++;
}
fclose($fp);
array_pop($arr_atoms);
return 0;
}
function get_vertices($conf_index)
{
$arr_refs = Array();
# Get the ids of the vertices of the convex hull
$filename = "hull_calcs/".sprintf("%04d",$conf_index).".hull";
$fp = fopen($filename, "r");
# Bypass the first line
$line = fgets($fp,4096);
while(!feof($fp))
{
$line = fgets($fp,4096);
$auxarr = explode(" ",$line);
array_pop($auxarr);
$arr_refs = array_merge($arr_refs, $auxarr);
}
fclose($fp);
// $arr_refs = array_unique($arr_refs); # This doesn't lose the keys
$arr_refs = array_keys(array_count_values($arr_refs)); # But this does.
return $arr_refs;
}
function parse_coordinates(&$arr_atoms, $conf_index)
{
# Parse the coordinates
$i = 0;
$filename = "hull_calcs/".sprintf("%04d",$conf_index);
$fp = fopen($filename, "r");
# Bypass the first two lines
$line = fgets($fp,4096);
$line = fgets($fp,4096);
while(!feof($fp))
{
$line = fgets($fp,4096);
$arr_atoms[$i]["coords"] = explode(" ",$line);
array_shift($arr_atoms[$i]["coords"]);
$i++;
}
fclose($fp);
array_pop($arr_atoms);
return 0;
}
function distance(&$atom1,&$atom2)
{
# Calculate the distance between the two atoms
$x1 = $atom1["coords"][0];
$x2 = $atom2["coords"][0];
$y1 = $atom1["coords"][1];
$y2 = $atom2["coords"][1];
$z1 = $atom1["coords"][2];
$z2 = $atom2["coords"][2];
return sqrt(pow($x1-$x2,2) + pow($y1-$y2,2) + pow($z1-$z2,2));
}
?>
<?PHP
/* Emre S. Tasci <e.tasci@tudelft.nl> *
* Parses the results/xxxx_results.txt file and generates
* an XYZ file such that the atoms are labeled as the
* vertice atoms, close-vertice atoms and inside atoms..
* 02/09/09 */
# Support for command line variable passing:
for($i=1;$i<sizeof($_SERVER["argv"]);$i++)
{
list($var0,$val0) = explode("=",$_SERVER["argv"][$i]);
$_GET[$var0] = $val0;
}
if($_GET["configuration"])
$conf_index = $_GET["configuration"];
else
exit("Please specify a configuration number [php plotter.php configuration=xxx]\n");
$filename = sprintf("%04d",$conf_index);
# Get the atom array from the results file. Since results file
# contains other arrays as well, we slice the file using sed for the
# relevant part
$last_arr_line = exec('grep "### CN Distribution" results/'.$filename.'_results.txt -n -m1|sed "s:^\([0-9]\+\).*:\1:gi"');
exec('sed -n "2,'.($last_arr_line-1).'p" results/'.$filename.'_results.txt > system_tmp_array_dump_file');
$str=file_get_contents('system_tmp_array_dump_file');
$atom_arr=eval('return '.$str.';');
unlink('system_tmp_array_dump_file');
# Now that we have the coordinates, we itarate over the atoms to check
# the characteristic and designate them in the XYZ file.
$fp = fopen("system_".$filename.".xyz","w");
fwrite($fp,sizeof($atom_arr)."\n");
fwrite($fp,"Generated by plotter.php\n");
for($i=0, $size=sizeof($atom_arr);$i<$size;$i++)
{
if(!array_key_exists("in",$atom_arr[$i]))
$atomlabel = "Mg";#Vertices
elseif($atom_arr[$i]["in"]==TRUE)
$atomlabel = $atom_arr[$i]["id"];#Inside atoms
else
$atomlabel = "O";#Close-vertice atoms
fwrite($fp,$atomlabel."\t".implode("\t",$atom_arr[$i]["coords"]));
}
fclose($fp);
?>
The big orange atoms are the omitted vertice atoms; the small red ones are the atoms critically close to the facets and hence also omitted ones. You can see the purple and yellow atoms which have passed the test sitting happilly inside, with the counting done using them.. 8)
]]>Anyway, recently it occured to me to write another script that would parse the output of the POSCAR2findsymm output and use that output to construct a CIF file that contained the equivalent sites, etc.. So, ladies and gentlemen, here it is (drums roll)…
Example: Code in action
Using the same POSCAR_Pd3S file as in the POSCAR2Cif entry… Feeding it to POSCAR2findsymm and then applying findsymm2cif on it
sururi@husniya tmp $ python /code/POSCAR2findsym/POSCAR2findsym.py POSCAR_Pd3S | php /code/vaspie/findsymm2cif.php speclist=Pd,S
_cell_length_a 7.21915
_cell_length_b 5.56683
_cell_length_c 7.68685
_cell_angle_alpha 90.00000
_cell_angle_beta 90.00000
_cell_angle_gamma 90.00000
_symmetry_Int_Tables_number 40
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
Pd01 Pd 0.75000 0.17834 0.31139
Pd02 Pd 0.75000 0.17845 0.69126
Pd03 Pd 0.50000 0.00000 0.00132
S01 S 0.75000 0.81592 0.50152
Code
#!/usr/bin/php
<?PHP
//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;
}
# ============/code/toolbox/commandline.inc.php========================================
if($_GET["h"]||$_GET["help"])
{
$help = <<<lol
* Script name: findsymmetry2cif.php *
* Converts the specified FINDSYM output to CIF format *
* More Information on FINDSYM : http://stokes.byu.edu/findsym.html *
* *
* (If you have access) More information on script: *
* http://dutsm2017.stm.tudelft.nl/mediawiki/index.php?title=Findsymm2cif
* *
* Emre S. Tasci <e.tasci@tudelft.nl> *
* 23/05/09 *
Usage:
f : input file (Default = stdin)
o : output file (Default = stdout)
identical : if set (identical=1), then the output will be the
transformation matrix and origin shift applied to the coordinates
so that the given coordinates will be regained.
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 findsymmetry2cif.php f=POSCAR_RhBi4 speclist=Bi,Rh identical=1
lol;
echo $help."\n";
exit;
}
$input = "php://stdin";
if($_GET["f"]) $input=$_GET["f"];
$outfile = "php://stdout";
if($_GET["o"]) $outfile=$_GET["o"];
$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"]);
//if(file_exists("findsymm2cif_php.tmp"))
// unlink("findsymm2cif_php.tmp");
if($input == "php://stdin")
{
$lol = file("php://stdin");
$fp1 = fopen("findsymm2cif_php2.tmp","w");
foreach($lol as $line)
fwrite($fp1,$line);
fclose($fp1);
$input = "findsymm2cif_php2.tmp";
}
exec('lstart=`grep -n "Type of each atom:" '.$input.'|sed "s:^\([0-9]\+\).*:\1:gi"`;lstart=`expr $lstart + 1`;lfinish=`grep -n "Position of each atom (dimensionless coordinates)" '.$input.'|sed "s:^\([0-9]\+\).*:\1:gi"`;lfinish=`expr $lfinish - 1`;cat '.$input.' |sed -n "`echo $lstart`,`echo $lfinish`p" | sed \':a;N;$!ba;s/\n/ /g\' > findsymm2cif_php.tmp');
exec('grep "Number of atoms in unit cell:" '.$input.' -A1 | tail -n1 >> findsymm2cif_php.tmp');
exec('grep "Space Group" '.$input.' | awk \'{print $3}\' >> findsymm2cif_php.tmp');
exec('grep "Origin at" '.$input.' | awk \'{print $3 "\t" $4 "\t" $5}\' >> findsymm2cif_php.tmp');
exec('grep "Vectors a,b,c:" '.$input.' -A3 | tail -n3 >> findsymm2cif_php.tmp');
exec('grep "Values of a,b,c,alpha,beta,gamma:" '.$input.' -A1 | tail -n1 >> findsymm2cif_php.tmp');
exec('grep "Wyckoff" '.$input.' -A1 -n | grep "^[0-9]\+-[0-9]"| sed "s:\(^[0-9]\+\)-[0-9]\+\(.*\):\1\t\2:gi" >> findsymm2cif_php.tmp');
if($input == "findsymm2cif_php2.tmp")
{
$input = "php://stdin";
//unlink("findsymm2cif_php2.tmp");
}
//exec("sh findsymm2cif.sh findsymm.out",$file);
$file = file("findsymm2cif_php.tmp");
$type_atoms = split("[ \t]+",trim($file[0]));
$specie_count = Array();
$k = 0;
$specie_count[$k] = 1;
for($i=1; $i < sizeof($type_atoms); $i++)
{
if($type_atoms[$i] != $type_atoms[$i-1])
$k++;
$specie_count[$k]++;
}
for($i=1; $i<sizeof($specie_count); $i++)
$specie_count[$i] += $specie_count[$i-1];
//print_r($specie_count);
$numatoms = trim($file[1]);
$sgno = trim($file[2]);
$origin_shift = split("[ \t]+",trim($file[3]));
$tr_x = split("[ \t]+",trim($file[4]));
$tr_y = split("[ \t]+",trim($file[5]));
$tr_z = split("[ \t]+",trim($file[6]));
$params = split("[ \t]+",trim($file[7]));
$specie = Array();
for($i = 8;$i<sizeof($file);$i++)
$specie[] = split("[ \t]+",trim($file[$i]));
$totatoms = 0;
for($i = 0; $i < sizeof($specie)-1; $i++)
{
$specie[$i][0] = $specie[$i+1][0] - $specie[$i][0] - 1;
$totatoms += $specie[$i][0];
}
$specie[sizeof($specie)-1][0] = $numatoms - $totatoms;
//print_r($specie);
//print_r($specie_count);
$atoms_counted = 0;
for($i = 0; $i < sizeof($specie); $i++)
{
# Going over atoms
$atoms_counted += $specie[$i][0];
for($j = 0; $j<sizeof($specie_count); $j++)
{
if($atoms_counted <= $specie_count[$j])
{
$specie[$i][0] = $species_type_template[$j];
break;
}
}
}
//print_r($specie);
if($_GET["identical"])
{
# Calculate the atom positions corresponding with the given POSCAR
$trmatrix=Array();
$trmatrix[0][0] = $tr_x[0];
$trmatrix[0][1] = $tr_x[1];
$trmatrix[0][2] = $tr_x[2];
$trmatrix[1][0] = $tr_y[0];
$trmatrix[1][1] = $tr_y[1];
$trmatrix[1][2] = $tr_y[2];
$trmatrix[2][0] = $tr_z[0];
$trmatrix[2][1] = $tr_z[1];
$trmatrix[2][2] = $tr_z[2];
for($i=0;$i < sizeof($specie); $i++)
{
$x = $specie[$i][1];
$y = $specie[$i][2];
$z = $specie[$i][3];
$xx = $x*$trmatrix[0][0] + $y*$trmatrix[1][0] + $z*$trmatrix[2][0];
$yy = $x*$trmatrix[0][1] + $y*$trmatrix[1][1] + $z*$trmatrix[2][1];
$zz = $x*$trmatrix[0][2] + $y*$trmatrix[1][2] + $z*$trmatrix[2][2];
$xx += $origin_shift[0];
$yy += $origin_shift[1];
$zz += $origin_shift[2];
$specie[$i][1] = $xx;
$specie[$i][2] = $yy;
$specie[$i][3] = $zz;
for($j=1;$j<4;$j++)
{
while($specie[$i][$j]>=1.0)
$specie[$i][$j]-=1.0;
while($specie[$i][$j]<0.0)
$specie[$i][$j]+=1.0;
}
}
}
$fp = fopen($outfile,"w");
fwrite($fp,"_cell_length_a"."\t".$params[0]."\n");
fwrite($fp,"_cell_length_b"."\t".$params[1]."\n");
fwrite($fp,"_cell_length_c"."\t".$params[2]."\n");
fwrite($fp,"_cell_angle_alpha"."\t".$params[3]."\n");
fwrite($fp,"_cell_angle_beta"."\t".$params[4]."\n");
fwrite($fp,"_cell_angle_gamma"."\t".$params[5]."\n");
fwrite($fp,"_symmetry_Int_Tables_number"."\t".$sgno."\n\n");
fwrite($fp, "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");
//print_r ($specie);
$k = 0;
$specie_type = $specie[0][0];
for($i=0;$i<sizeof($specie);$i++)
{
if($specie[$i][0] == $specie_type)
$k++;
else
{
$k = 1;
$specie_type = $specie[$i][0];
}
fwrite($fp, $specie[$i][0].sprintf("%02d",$k)."\t".$specie[$i][0]."\t".sprintf("%8.5f",$specie[$i][1])."\t".sprintf("%8.5f",$specie[$i][2])."\t".sprintf("%8.5f",$specie[$i][3])."\n");
}
fclose($fp);
?>
]]>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
]]>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);
?>
]]>#!/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
]]>Meeting with Dr. Villars
The biggest event that occured to me was meeting with Dr. Pierre Villars, the person behind the Pauling File database which I have benefitted and still benefitting from. When I had registered to the “MITS 2008 – International Symposium on Materials Database”, an event hosted by the NIMS Materials Database, I had no idea that Dr. Villars would also be attending to the event, so it was a real surprise learning afterwards that I’d have a chance to meet him in person.
Dr. Villars is a very kind and a cheerful person, I had a greatest time knowing him and afterwards I had the privilege of visiting him in his home office in the lovely town of Vitznau, Switzerland where I enjoyed the excellent hospitality of the Villars family. Dr. Villars helped me generously in providing the data in my research project for finding transformation toughening materials through datamining which I’m indebted for.
The MITS 2008 – Int’l Symposium on Materials Database
I was frequently using the databases on the NIMS’ website and had the honour to participate in the last year’s event together with my dear supervisor Dr. Marcel Sluiter. On that occasion, I presented two talks – a materials based one, and another about database query&retrieval methods- titled “Inference from Various Material Properties Using Mutual Information and Clustering Methods” and “Proposal of the XML Specification as a Standard for Materials Database Query and Result Retrieval”, respectively.
I have an ever growing interest on the Japanese Culture since I was a youngster and due to this interest, my expectations from this trip to Japan was high enough for me to be afraid of the fact that I may have disappointed since it was next-to-impossible for anything to fulfill the expectations I had built upon. But, amazingly, I was proved false: the experience I had in Japan were far more superior to anything I had anticipated. The hospitality of the organizing commitee was incredible: everything had been taken care of and everyone was as kind as one could be. I had the chance to meet Dr. Villars as I’ve already wrote about above, and also Dr. Xu, Dr. Vanderah, Dr. Dayal and Dr. Sims and had a great time. It was also due to this meeting, I met a dear friend.
Ongoing Projects
During my time here, I still don’t have anything published and that’s the number one issue that’s bothering me. Apart from my research on the transformation toughening materials, I worked for a time on a question regarding to a related phenomenon, namely about the martensitic transformation mechanism on some transition-metal alloy types. We were able to come up with a nice equation and adequate coefficients that caught the calculation results but couldn’t manage to transfer it to -so to say- a nicely put theory and that doomed the research to being an “ad-hoc” approach. I’ll be returning to it at the first chance I get, so currently, it’s been put on ice (and now remembering how, in the heat of the research I was afraid that someone else would come up with the formula and it would be too late for me.. sad, so sad..). One other (side) project was related to approximate a specific amorphous state in regards with some basis structures obtained using datamining. This project was really exciting and fun to work on because I got to use different methods together: first some datamining, then DFT calculations on the obtained structures and finally thermodynamics to extrapolate to higher temperatures. We’re at the stage of revising the draft at the moment (gratefully!).
Life in the Netherlands & at the TUDelft
It has been almost one and a half year since I started my postdoc here in TUDelft. When we first came here, my daughter was still hadn’t started speaking (yet she had been practicing walking) and now she approximately spent the same amount of time in the Netherlands as she had spent in Turkey (and thankfully, she can speak and also running to everywhere! 8).
When I had come here, I was a nano-physicist, experienced in MD simulations and -independently- a coder when the need arouse. Now, I’m more of a materials guy and coding for the purpose.. 8) I enjoy the wide variety of the jobs & disciplines people around me are working on and definitely benefitting from the works of my friends in our group. I couldn’t ask for more, really. (OK, maybe I’d agree if somebody would come up with an idea for less windy and more sunny weather here in Delft. That’s the one thing I couldn’t still get used to. Delft is a marvelous city with bikes and kind people everywhere, canals full of lillies in the spring, historical buildings everywhere, Amsterdam, Den Haag and Rotterdam just a couple of train stations away and if you wish, Antwerpen, Brussels and Brugges (Belgium) is your next-door neighbour but this dark and closed weather can sometimes be a little mood lowering.) TUDelft is really a wonderful place to be with all the resources and connections that makes you feel you can learn anything and meet anyone that is relevant to the work you’re doing without worrying about the cost. And you’re being kept informed about the events taking place about you through the frequent meetings in your group, in your section, in your department, in your field. I’ve learned a lot and I’m really happy I was accepted here.
So that’s about the things I’ve enjoyed experiencing and living through.. And I promise, I’ll be posting more regularly and frequently from now on..
Yours the chatty.
]]>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:
What the code does:
It has 3 seperate processes. Let’s designate them A, B & C.
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 "<img src=\"/tex/".$file."\" title=\"".$formula."\"><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…
]]>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.
]]>Dr. Villars had already mentioned Miedema’s passing at a young age but still it was really tragic to find such a sad detail on a scientific article.
I checked the web for a Miedema biography but couldn’t find any. The only information I could find was from H. Bakker’s book which I’m quoting the preface below:
PREFACE
Andries Miedema started developing his rather unconventional model in the sixties as a professor at the University of Amsterdam and continued this work from 1971 at Philips Research Laboratories. At first he encountered scepticism in these environments, where scientists were expecting all solutions from quantum- mechanics. Of course, Miedema did not deny that eventually quantum mechanics could produce (almost) exact solutions of problems in solid-state physics and materials science, but the all-day reality he experienced was, – and still is -, that most practical problems that are encountered by researchers are too complex to allow a solution by application of the Schr6dinger equation. Besides, he believed in simplicity of nature in that sense that elements could be characterised by only a few fundamental quantities, whereby it should be possible to describe alloying behaviour. On the basis of a comparison of his estimates of the heat of formation of intermetallic compounds with ample experimental data, he was gradually able to overcome the scepticism of his colleagues and succeeded in convincing the scientific world that his model made sense. This recognition became distinct by the award of the Hewlett Packard prize in 1980 and the Hume-Rothery decoration in 1981. At that time, Miedema himself and many others were successfully using and extending the model to various questions, of course realising that the numbers obtained by the model have the character of estimates, but estimates that could not be obtained in a different fast way. It opened new perspectives in materials science. Maybe the power of the model is best illustrated by a discussion, the present author had with Jim Davenport at Brookhaven National Laboratory in 1987. Dr Davenport showed him a result that was obtained by a one-month run on a Gray 1 supercomputer. This numerical result turned out to be only a few kilo Joules different from the outcome by Miedema’s model, obtained within a few minutes by use of a pocket calculator. Andries Miedema passed away in 1992 at a much too young age of 59 years. His death came as a bad shock not only for the scientific community in The Netherlands, but for scientists all over the world. However, for his family, friends and all the people, who knew him well, it is at least some consolation that he lives on by his model, which is being used widely and which is now part of many student programmes in physics, chemistry and materials science.
It is the aim of the present review to make the reader familiar with the application of the model, rather than to go in-depth into the details of the underlying concepts. After studying the text, the reader should be able, with the aid of a number of tables from Miedema’s papers and his book, given in the appendix, to make estimates of desired quantities and maybe even to extend the model to problems that, so far, were not handled by others. Beside, the reader should be aware of the fact that not all applications could be given in this review. For obvious reasons only part of all existing applications are reported and the reader will find further results of using the model in Miedema’s and co-worker’s papers, book and in literature.
Dr Hans Bakker is an associate professor at the Van der Waals-Zeeman Institute of the University of Amsterdam. In 1970 Andries Miedema was the second supervisor of his thesis on tracer diffusion. He inspired him to subsequent work on diffusion and defects in VIII-IIIA compounds at a time, when, – as Bakker’s talk was recently introduced at a conference -, these materials were not yet relevant’ as they are now. Was it again a foreseeing view of Andries Miedema? Was it not Miedema’s standpoint that intermetallic compounds would become important in the future?
Rather sceptic about the model as Bakker was at first like his colleagues, his scepticism passed into enthusiasm, when he began using the model with success in many research problems. In 1984 Bakker’s group started as one of the first in the world ball milling experiments. Surprisingly this seemingly crude technique turned out to be able to induce well-defined, non-equilibrium states and transformations in solids. Miedema’s model appeared again to be very helpful in understanding and predicting those phenomena.
Hans Bakker
Mauro Magini
Preface from Enthalpies in Alloys : Miedema’s Semi-Emprical Model
H. Bakker
Trans Tech Publications 1998 USA
ISSN 1422-3597
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!
]]>