Coding – Hex, Bugs and More Physics | Emre S. Tasci http://www.emresururi.com/physics a blog about physics, computation, computational physics and materials... Wed, 15 Jan 2014 13:36:05 +0000 en-US hourly 1 https://wordpress.org/?v=4.9.3 Relations… http://www.emresururi.com/physics/?p=205 http://www.emresururi.com/physics/?p=205#respond Wed, 15 Jan 2014 13:22:23 +0000 http://www.emresururi.com/physics/?p=205 Recently, we needed to schematize the relations between the coefficients of a symmetric tensor matrix such that, beyond the tedious numerical values, the basic relations such as equivalency or (-1) times would be visible to the eye.

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

if you really want to, you can verify that its coefficients satisfy the following properties:
Non-zero Elements:
(1,1), (1,8), (1,9), (1,10),
(2,2), (2,5), (2,6), (2,11), (2,14),
(3,3), (3,6), (3,12), (3,15),
(4,4), (4,7), (4,16),
(5,5), (5,11), (5,14),
(6,6), (6,12), (6,15),
(7,7), (7,16),
(8,8), (8,9), (8,10),
(9,9), (9,10),
(10,10),
(11,11), (11,14),
(12,12), (12,14), (12,15),
(13,13), (14,14), (15,15),
(16,16)

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:

tensor_graph1

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

(don’t forget to restart the http server if you’re working via a web interface instead of cli)

Before going further into the graph production, here is a description of producing a random matrix compatible with the given conditions in Octave:

Constructing a sample matrix 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");

Finding out the relations between the coefficients

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).

Painting on a canvas

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:

circ1

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

circ2

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

circ3

Mapping the rows and cols to x and y

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:
map1

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:
map2

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");

). We’re almost ready aside from the connecting the related entries!

Connecting the dots

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

circles_line

Devil is in the details

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.,

circles_line_better

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!

One picture is worth a thousand words

You can have a sneak peak of the working code via here: http://144.122.31.125/cgi-bin/cryst/programs/tensor_graph/

]]> http://www.emresururi.com/physics/?feed=rss2&p=205 0 Quantum Espresso meets bash : Automated running for various values http://www.emresururi.com/physics/?p=156 http://www.emresururi.com/physics/?p=156#comments Tue, 13 Mar 2012 22:43:51 +0000 http://www.emresururi.com/physics/?p=156 Now that we can include comments in our input files and submit our jobs in a tidied up way, it’s time to automatize the “convergence search” process. For this, we’ll be using the same diamond example of the previous entry.

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 &lt;= 9; i++ ))
do
  echo $i;
done

1
2
3
4
5
6
7
8
9

But to work with rational numbers, things get a little bit complex:
for (( i = 1; i &lt;= 5; i++ ))
do
    val=$(awk "BEGIN {print 3+$i/5 }")
    echo $val;
done

3.2
3.4
3.6
3.8
4

Alas, it is still not very practical to everytime calculate the necessary values to produce our intended range. This is where the mighty Octave once again comes to the rescue: [A:I:B] in Octave defines a range that starts from A, and continues by incrementing by I until B is reached (or exceeded). So, here is our example in Octave:
octave:1&gt; [3.2:0.2:4]
ans =

   3.2000   3.4000   3.6000   3.8000   4.0000

You can generate temporary input scripts for Octave but better yet, you can also feed simple scripts directly to Octave:
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

The “-q” parameter causes Octave to run in “quiet” mode, so we just feed the range via pipe to Octave and will harvest the output. In general, we’d like to keep the “ans =” and the blank lines out, and then, when there are more than 1 line’s worth of output, Octave puts the “Column” information as in:
octave:4&gt; [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

all these nuisances can go away via a simple “grep -v” (with the ‘v’ parameter inverting the filter criteria):
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

Replacement:
Now that we have the variable’s changing values, how are we gonna implement it to the input script? By search and replace. Suppose, we have a template in which the variable’s value is specified by “#REPLACE#”, all we need to do is to “sed” the file:
cat diamond.scf.in | sed "s:#REPLACE#:$a:g"&gt;tmp_$a.in

where “$a” is the current value in the iteration.

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)}'

For each iteration we thus obtain the energy for our $a value, so we put them together in two columns into a data file I’d like to call “E_vs_$a.txt”.

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"&gt;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 &gt;&gt; $logfilename
    echo -e $a "    " $energ
done

echo "Results stored in: "$logfilename

Action:
Modify the sample input from the previous post so that it’s lattice parameter A is designated as something like “A = #REPLACE#”, as in:
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

Then call the qe-opti.sh script with the related parameters:
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 &lt; tmp_3.2000.in &gt; 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 &lt; tmp_3.4000.in &gt; 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 &lt; tmp_3.6000.in &gt; 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 &lt; tmp_3.8000.in &gt; 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 &lt; tmp_4.0000.in &gt; 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 &lt; tmp_4.2000.in &gt; 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

Where the “E_vs_a.txt” contains the … E vs a values!:
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

You can easily plot it using GNUPlot or whatever you want. In my case, I prefer GNUPlot:
sururi@husniya:~/shared/qe/diamond_tutorial$ plot E_vs_a.txt

giving me:

plot_E_vs_a

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

plot_E_vs_a_with_fit

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

]]> http://www.emresururi.com/physics/?feed=rss2&p=156 1 Quantum Espresso meets bash : Primer and tidying up http://www.emresururi.com/physics/?p=145 http://www.emresururi.com/physics/?p=145#comments Tue, 13 Mar 2012 15:24:00 +0000 http://www.emresururi.com/physics/?p=145 In my previous postdoc position, DFT calculations were governing my research field but with my current position at the Bilbao Crystallographic Server, we are more on the theoretical side of things then calculations. Nevertheless, recently, my old ‘flame’ have been kindled and I decided to scrape the rust.

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 &lt;111&gt;    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

so, it’s not the tidiest parameter file around but nevertheless, I’ve included some explanations and it’s easily customizable by simple commenting out, right?…

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

I’ve installed the QE to take advantage of the cluster nodes via MPI, so I submit my jobs via:
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.

]]> http://www.emresururi.com/physics/?feed=rss2&p=145 2 Instant Karma: Updating the clones automatically upon push in Mercurial http://www.emresururi.com/physics/?p=134 http://www.emresururi.com/physics/?p=134#respond Mon, 29 Aug 2011 15:00:27 +0000 http://www.emresururi.com/physics/?p=134 Suppose that I have a web server which consists of two identical computers behind a splitter that distributes the incoming requests made to a shared ip. Furthermore, the central repository is being kept on a different computer and as a final thing, assume that I sometimes use the computer located in my office, and sometimes the computer at my home as indicated in the following diagram:
mercurial_scheme

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:

  1. PC@wherever: hg com
  2. PC@wherever: hg push
  3. PC@wherever: ssh CentralRepo
  4. CentralRepo: hg up
  5. CentralRepo: ssh Server#1
  6. Server#1: hg pull
  7. Server#1: hg up
  8. Server#1: ssh Server#2
  9. Server#2: hg pull
  10. Server#2: hg up

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 &gt; /dev/null

This hook just tells the mercurial to execute the ‘incominghook.sh’ script in the .hg directory (relative to the project’s root dir) whenever there’s an incoming data (i.e., “a push to this -central- repository). So what it should do? It should ‘hg up’ the central repository, for starters. Then, it should ‘ping’ the two servers so that they also pull and update from the -now up-to-date- central repository. This is what the incoming.sh script looks like:
#!/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

Again, if it is not very clear how will the servers get into action by a mere SSH connection, I do urge you to check the aforementioned entry.

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)

So that, when “somebody” connects with the relevant key_id file, the server executes the ‘hgpull_project_x.sh’ (located at /home/sururi/bin/) and returns the output, if any. You can guess that this script actually goes to the project’s location, pulls the changes and then ups it. As is evident from the contents of the ‘hgpull_project_x.sh’ script:
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

As you can learn from the strongly referred previous entry, the corresponding line on central repo’s authorized_keys file is something like:
command="hg-ssh /path/to/project/",no-port-forwarding,
no-X11-forwarding,no-agent-forwarding,no-pty ssh-rsa AAAAB...

and we’re done.

]]> http://www.emresururi.com/physics/?feed=rss2&p=134 0 Accessing Mercurial with limited SSH access using key and hg-ssh http://www.emresururi.com/physics/?p=114 http://www.emresururi.com/physics/?p=114#comments Sun, 28 Aug 2011 20:15:00 +0000 http://www.emresururi.com/physics/?p=114 Today, I wondered about (actually needed) the possibility to be able to limit (and hence connect afterwards) the access to the repository center of mercurial, using SSH.

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      |
+-----------------+

I didn’t specified any password because I’m intending to have it used automatically (in the hook procedures) and I specified the location of the pair files as the current directory with the names: id_rsa_example  & id_rsa_example.pub

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

And from this moment on -hopefully- you should be able to connect without the remote asking for your password (since we’re not using the default place & filename for the key, i.e. ~/.ssh/id_rsa and ~/.ssh/id_rsa.pub, this newly generated key should be included with the “-i” (identity file) option as in: ssh -i /tmp/tmp/id_rsa_example yourusername:remote)… but, that wasn’t exactly what we wanted. We are thinking of a scenario where you’d like to (semi-)freely distribute the SSH access to the central repository to your developers while making sure that they wouldn’t (couldn’t) do something nasty while they are connected to the remote computer.

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

By adding the “command” option to the beginning of this line would make sure that whenever somebody with the corresponding key connects, that command is run automatically, the results passed back and the connection is terminated afterwards. So, far example. modifying the entry upstairs as:
command="date" ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAnq5rMLnoab+2F28g/nb58RBENWtX395TuyDFsYkalGaZxrziwDoau/wglkU19DbcAVKgw0p6lMEIuh2iALOppRzxrTgFFhJkL1dxzkugbbPEoSWyfrj9FivzpnxHWgRHQApQeWUBOZhroDTURwfqcyC9SW020CR57jLWfgw+idqwtCu+ZBYmEyHSJcZIH2mWXLrUQ8OalxCFVaLKL50Lpc7V8XJPs+Pg6MPVgfDUqMdjrGkAF7j4viOHTjDWP1h4Ngim70dOeyxWtuqbCbxM4APTShaqET42sj1jHxL2m1dJzXX8s/gEdN0O09hZPhI6rlC+ANWIdJ1vJfODMXWaQ== sururi@local== sururi@vala

would cause the following behaviour from the remote:
sururi@local:~$ ssh -i /tmp/tmp/id_rsa_example sururi@remote
Sun Aug 28 22:05:54 CEST 2011
Connection to remote closed.

I’m a bit tired now writing in full details, so will skip some obvious things from now on. Mercurial has an SSH-wrapper called ‘hg-ssh’ exactly for this purpose, and you can use it by specifying the paths of the repositories as command arguments as in:
command="hg-ssh /path/to/repository" ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAQEAnq5rMLnoab+2F28g/nb58RBENWtX395TuyDFsYkalGaZxrziwDoau/wglkU19DbcAVKgw0p6lMEIuh2iALOppRzxrTgFFhJkL1dxzkugbbPEoSWyfrj9FivzpnxHWgRHQApQeWUBOZhroDTURwfqcyC9SW020CR57jLWfgw+idqwtCu+ZBYmEyHSJcZIH2mWXLrUQ8OalxCFVaLKL50Lpc7V8XJPs+Pg6MPVgfDUqMdjrGkAF7j4viOHTjDWP1h4Ngim70dOeyxWtuqbCbxM4APTShaqET42sj1jHxL2m1dJzXX8s/gEdN0O09hZPhI6rlC+ANWIdJ1vJfODMXWaQ== sururi@local

(and for more security, I would suggest you include the other options as well, such as:)
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

and you can pull now on the local computer from the central repository at the remote computer with a command like:
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

]]> http://www.emresururi.com/physics/?feed=rss2&p=114 2 from SVN to Mercurial http://www.emresururi.com/physics/?p=112 http://www.emresururi.com/physics/?p=112#respond Tue, 23 Aug 2011 14:28:07 +0000 http://www.emresururi.com/physics/?p=112 I’ve been using SVN for revision control system for some time. I certainly benefited from it to organize or fall-back to a stable version (“the good old times”) when I was in trouble while I was computing solo but there were two things that always upset me while using it:

  • Commiting meant publishing: This issue might be worked around when you have your trunk (we actually have the trunk for the ‘most recent’ status and a ‘released’ branch for the tested and stable programs of that trunk) and each of the developers play in their branches until it looks OK and hence merge it to trunk (and hopefully eventually make it to the ‘released’ branch), but the problem is, most of the time I’m in the middle of the coding process (in my branch) and it would be nice if I didn’t affect the main repository everytime I made a commit (it’s really not very nice) – in other words, even in my branch I would like to save while not making them available. The answer to this lies of course in distributed revision control systems (dvcs).
  • Merging: Merging is hell in SVN when you have many developers. Been there, suffered that.

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.

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=112 0     Exploring through the paths of subgroups, collecting indexes on the way… http://www.emresururi.com/physics/?p=102 http://www.emresururi.com/physics/?p=102#respond Mon, 06 Jun 2011 14:46:23 +0000 http://www.emresururi.com/physics/?p=102 Having a personal wiki has its pros and cons with the pros being obvious but, on the other hand, one doesn’t feel like reposting what he has already committed into the corresponding wikipedia entry (and since reposting means filtering out the sensitive information and generalizing it, it is double the effort). So this fact, more or less, explains my absence of activity in this blog for the last 1.5 years. But then, there is one negative aspect of personal wikipedia “blogging” – it’s like putting everything under the rug : you classify the data/information and then you forget that it is there. That’s when “actual” blogging (like this one) comes handy. So, sorry & welcome back..

    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:

    subgroupgraph_tempfile963

    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:

    maxsubs_subggraph

    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]

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=102 0     Coordination Numbers http://www.emresururi.com/physics/?p=91 http://www.emresururi.com/physics/?p=91#respond Thu, 03 Sep 2009 15:31:10 +0000 http://www.emresururi.com/physics/?p=91 Recently, we are working on liquid alloys and for our work, we need the coordination number distribution. It’s kind of a gray area when you talk about coordination number of atoms in a liquid but it’s better than nothing to understand the properties and even more. I’d like to write a log on the tools I’ve recently coded/used in order to refer to in the future if a similar situation occurs. Also, maybe there are some others conducting similar research so, the tools introduced here may be helpful to them.

    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

    (The compilation was of 1000 configurations, each with 2048 atoms)

    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

    The sets of 3 numbers in the xxxx.hull file is the ids/numbers of the atoms forming the edges of the facets, i.e., they are the vertices.

    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

    This script is runned within the following php code which then takes the filtered atoms and does the counting using them:
    <?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));
}

?>

    The code above generates a “results/xxxx_results.txt” file containing all the information obtained in arrays format and also appends to “results.txt” file the relevant configuration file’s CN distribution summary. The systems can be visualized using the output generated by the following plotter.php script:
    <?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);
?>

    which generates a “system_xxxx.xyz” file, ready to be viewed in a standard molecule viewing software. It designates the vertice atoms (of the convex hull, remember? 8)as “Mg” and the atoms close to them as “O” for easy-viewing purpose. A sample parsed file looks like this:

    The filtered case of a sample Na-K system

    The filtered case of a sample Na-K system

    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)

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=91 0     POSCAR2Cif with symmetries discovered http://www.emresururi.com/physics/?p=90 http://www.emresururi.com/physics/?p=90#respond Thu, 04 Jun 2009 08:16:57 +0000 http://www.emresururi.com/physics/?p=90 In previous posts, I had submitted two codes:

    • POSCAR2Cif : that converts a given POSCAR file to CIF, so to speak, ruthlessly, i.e. just as is, without checking for symmetries or anything.
    • POSCAR2findsymm : another converter that takes a POSCAR file, and prepares an input file such that Harold Stokes’ findsym code from the ISOTROPY package would proceed it and deduce the symmetry information. If findsym executable is not accessible in your system (or if you run the script with the “t” flag set to on), it would just print the input to the screen and you could use it to fill the input form of the code’s web implementation, instead

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

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=90 0     Some unit cell functions for Octave http://www.emresururi.com/physics/?p=89 http://www.emresururi.com/physics/?p=89#respond Sat, 23 May 2009 16:43:56 +0000 http://www.emresururi.com/physics/?p=89 I’ve written some functions to aid in checking the transformations as well as switching between the parameter ↔ vector representations.

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

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

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

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

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

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

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

    Example

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

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=89 0     POSCAR2Cif http://www.emresururi.com/physics/?p=88 http://www.emresururi.com/physics/?p=88#comments Sat, 23 May 2009 16:26:12 +0000 http://www.emresururi.com/physics/?p=88 Converts POSCAR files to CIF format. Uses Octave (latvec2par) to convert the unit cell vectors to lattice cell parameters. It is assumed (compulsory) that the atom coordinates in the POSCAR file are in fractional (direct) coordinates.

    Parameters

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

    Example

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

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

    Code

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

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=88 1     POSCAR2findsymm http://www.emresururi.com/physics/?p=87 http://www.emresururi.com/physics/?p=87#comments Tue, 24 Mar 2009 09:17:37 +0000 http://www.emresururi.com/physics/?p=87 A very simple code I wrote while studying the Python language that fetches the information from a VASP POSCAR file and prepares (and feeds) the input for Harold Stokes’ findsym (you should actually download the executable included in the ISOTROPY package instead of the referred web version) program to find the unit cell symmetry.

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

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

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

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

   Accepts alternative POSCAR filename for argument
   (default = POSCAR)

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

                                             10/03/09 */
"""

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

force_output = False

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

f = open(filename,'r')

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

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

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

initial_guess = "P" # For the centering

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

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

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

# print pos_vec

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

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

    Example Usage:

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

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

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=87 1     LaTeX to PNG and all that equations.. http://www.emresururi.com/physics/?p=85 http://www.emresururi.com/physics/?p=85#respond Mon, 23 Mar 2009 14:32:29 +0000 http://www.emresururi.com/physics/?p=85 After migrating this blog to a new server, I had some difficulties in writing up the equations. It was due to the fact that on this server, I wasn’t able to execute the necessary commands to parse the latex code and convert it to PNG format.

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

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

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

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

    The code is as follows:

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

    To apply this:

    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…

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=85 0     Two bash scripts for VASP http://www.emresururi.com/physics/?p=84 http://www.emresururi.com/physics/?p=84#comments Mon, 20 Oct 2008 12:56:20 +0000 http://www.emresururi.com/physics/?p=84 It’s of great importance to have sufficient number of k-points when doing DFT calculations. The most practical (but CPU time costly) way is to check energies for different number of k-points until some kind of "convergence" is reached.

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

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

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

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

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

    #!/bin/bash

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

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

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

    # run the vasp
    vasp > vasp.out

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

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

    After it ends, you should see something like

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

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

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

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

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

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

    ./energy_extract_from_OUTCAR.sh > energies.txt

    And then, obtaining the last column via awk

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

    Now you can also easily plot the difference file.

    Hope this scripts will be of any help.

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=84 1     Miedema et al.’s Enthalpy code — 25 years after.. http://www.emresururi.com/physics/?p=81 http://www.emresururi.com/physics/?p=81#comments Mon, 28 Jul 2008 13:34:28 +0000 http://www.emresururi.com/physics/?p=81 Yes, it’s been 25 years since A.K. Niessen, A.R. Miedema et al.’s "Model Predictions for the Enthalpy of Formation of Transition Metal Alloys II" titled article (CALPHAD 7(1) 51-70 1983) was published. It can be thought as a sequel to Miedema et al.’s 1979 (Calphad 1 341-359 1979) and 1980 dated (Physica 100 1-28 1980) papers with some update on the data as well as a computer code to calculate the enthalpies of formation written in ALGOL.

    I have been currently working on these papers and by the way ported the code presented in CALPHAD 1983 to FORTRAN and here it is:

     

    C      THE IMPLEMENTATION OF THE ALGOL CODE FROM
C      A.K. Niessen, F.R. de Boer, R. Boom, P.F. de Chatel
C      W.C.M. Mattens and A.R. Miedema
C      CALPHAD Vol.7 No.1 pp. 51-70, 1983
C      Model Predictions for the Enthalpy of Formation of Transition 
C        Metal Alloys II"
C
C      by Emre S. Tasci, TUDelft (2008)

      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
C      DOUBLE PRECISION va,aa,p,r,x,dn,dg,ng,pq,xa,xm,csa,csm
C      DOUBLE PRECISION cas,cms,fam,dph,vm,am,fma,var,vmr
C      DOUBLE PRECISION vas, vms, xas, xms, xva,xvm, dgl,pqrs, pqrl
      INTEGER arrout, z1,z2,a,m,n,w,dh,ga,gm,gal,gml,nr,dHtrans
      LOGICAL detp, detr, bool
      CHARACTER(LEN=5) Symbol
      DIMENSION arrout(15)
      DIMENSION Phi(75), acf(75), dHtrans(75), Nws113(75), Rcf(75),
     + V213(75), xxx(9), Symbol(75)
      COMMON /CDH/ x, dh
      COMMON /CDPH/ dph, Phi, Nws113, ng, r, Rcf, p, pq, dn, pqrs, pqrl

      OPEN(UNIT=3,FILE='inpdata.txt',STATUS='OLD')
      DO I=1,75
        READ(3,*)nr,Symbol(nr),Phi(nr),Nws113(nr),V213(nr),
     +  acf(nr),Rcf(nr),dHtrans(nr)
C       WRITE(*,*)nr,Symbol(nr),Phi(nr),Nws113(nr),V213(nr),
C     + acf(nr),Rcf(nr),dHtrans(nr)
      END DO
      CLOSE(UNIT=3)

      xxx(1) = 0.001
      xxx(2) = 1.0/6.0
      xxx(3) = 1.0/4.0
      xxx(4) = 1.0/3.0
      xxx(5) = 1.0/2.0
      xxx(6) = 2.0/3.0
      xxx(7) = 3.0/4.0
      xxx(8) = 5.0/6.0
      xxx(9) = 0.999

      DO z1 = 1,46 
          WRITE(*,200)z1,Symbol(z1),Phi(z1),Nws113(z1)**3,
     +      V213(z1)**(3.0/2),dHtrans(z1)
200       FORMAT(I2," ",A5," Phi: ",F5.2,"V  Nws: ",F5.2,
     +     "d.u.  Vmole: "
     +     ,F5.2,"cm3  DeltaHtrans:",I3,"kJ/mole",/)
          WRITE(*,250)"M","AM5","AM3","AM2","AM",
     +     "MA2","MA3","MA5","AinM","AM","MinA"
250       FORMAT(A,"      ",9(A4,"  "),A4,/)
          DO z2 = 1, 75
             m = z2
             a = z1
             CALL PQRSL(a,m)

             va = V213(a)
             aa = acf(a)
             ga = dHtrans(a)
             gal = ga

             vm = V213(m)
             am = acf(m)
             gm = dHtrans(m)

             IF (m.EQ.67.OR.m.EQ.68) THEN
                 gml = 0
             ELSE
                 gml = gm
             END IF

             n = 1

             DO I = 1,9
              xa = xxx(I)
              IF(xa.EQ.0.001.OR.xa.EQ.0.999) THEN
                  bool = .TRUE.
              ELSE
                  bool = .FALSE.
              END IF

              xm = 1 - xa
              xva = xa * va
              xvm = xm * vm

              dgl = xa * gal + xm * gml
              dg = xa * ga + xm * gm
              csa = xva / (xva + xvm)

              csm = 1 - csa
              fam = csm * (1 + 8 * csa * csa * csm * csm)
              fma = fam * csa / csm

              var = va * (1 + aa * csm * dph)
              vmr = vm * (1 - am * csa * dph)

              vas = va * (1 + aa * fam * dph)
              vms = vm * (1 - am * fma * dph)

              xar = xa * var
              xas = xa * vas
              xmr = xm * vmr
              xms = xm * vms

              cas = xas / (xas+xms)
              cms = 1 - cas
              csa = xar / (xar + xmr)
              csm = 1 - csa

              fam = cms * (1 + 8 * cas * cas * cms * cms)
              fma = fam * cas / cms

              IF(bool) THEN
                  IF(n.EQ.1) THEN 
                      GOTO 292
                  ELSE
                      GOTO 392
                  END IF
              END IF
                  IF(xa.EQ.0.5) THEN
                      x = csm * xar * p * pqrl + dgl
                      CALL MAXI
                      arrout(n+5) = dh
                  END IF
                  x = fam * xas * p * pqrs + dg
                  CALL MAXI
                  arrout(n) = dh
                  arrout(11) = gml
                  GOTO 492
292               x = (csm * xar * p * pqrl + dgl) / xa
                  CALL MAXI
                  arrout(n) = dh
                  arrout(11) = gml
                  GOTO 492
392               x = (csm * xar * p * pqrl + dgl) / xm
                  CALL MAXI
                  arrout(n) = dh
                  arrout(12) = gal
492               n = n + 1
             END DO
C            DO J=1,15
C              WRITE(*,500) J, arrout(J)
C            END DO
             WRITE(*,600)Symbol(m),arrout(2),arrout(3),arrout(4),
     +        arrout(5),
     +        arrout(6),arrout(7),arrout(8),arrout(1),arrout(10),
     +        arrout(9)
             IF(Symbol(m).EQ."Ni".OR.Symbol(m).EQ."Pd".OR.Symbol(m).EQ.
     +       "Lu".OR.Symbol(m).EQ."Pt".OR.Symbol(m).EQ."Pu".OR.
     +       Symbol(m).EQ."Au".OR.Symbol(m).EQ."H".OR.Symbol(m).EQ.
     +       "Cs".OR.Symbol(m).EQ."Mg".OR.Symbol(m).EQ."Ba".OR.
     +       Symbol(m).EQ."Hg".OR.Symbol(m).EQ."Tl".OR.Symbol(m).EQ.
     +       "Pb".OR.Symbol(m).EQ."Bi")WRITE(*,*)""
          END DO
      WRITE(*,"(2A,/)")"-----------------------------------",
     + "----------------------------"
      END DO

C500   FORMAT("O[",I3,"] = ",I50)
600   FORMAT(A5,"  ",9(I4,"  "),I4)

      STOP
      END

      SUBROUTINE MAXI
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      INTEGER dh
      COMMON /CDH/ x, dh
      IF(x.LT.-999.OR.x.GT.999) THEN
        dh = 999
      ELSE
        dh = NINT(x)
      END IF
C      WRITE(*,*)"x:",dh
      RETURN
      END

      SUBROUTINE PQRSL(a,b)
      IMPLICIT DOUBLE PRECISION(A-H,N-Z)
      INTEGER a,b
      LOGICAL detp,detr
      DIMENSION Phi(75), acf(75), dHtrans(75), Nws113(75), Rcf(75),
     + V213(75)
      COMMON /CDPH/ dph, Phi, Nws113, ng, r, Rcf, p, pq, dn, pqrs, pqrl

          dph = Phi(a) - Phi(b)
          dn = Nws113(a)-Nws113(b)
          ng = (1/Nws113(a) + 1/Nws113(b))/2

          IF (detr(a).EQV.detr(b)) THEN
              r = 0
          ELSE
              r = Rcf(a)*Rcf(b)
          END IF

          IF (detp(a).AND.detp(b)) THEN
              p = 1.15 * 12.35
C             BOTH TRUE
          ELSE IF (detp(a).EQV.detp(b)) THEN
              p = 12.35 / 1.15
C             BOTH FALSE
          ELSE
              p = 12.35
          END IF
          p = p / ng
          pq = -dph * dph + 9.4 * dn * dn
          pqrs = pq - r
          pqrl = pq - 0.73 * r

      RETURN
      END

      LOGICAL FUNCTION  detp(z)
      INTEGER z
      IF (z.LT.47.AND.(z.NE.23.AND.z.NE.31)) THEN
          detp = .true.
      ELSE
          detp = .false.
      END IF
      RETURN
      END

      LOGICAL FUNCTION  detr(z)
      INTEGER z
      IF (z.LT.47.OR.(z.GT.54.AND.z.LT.58)) THEN
          detr = .true.
      ELSE
          detr = .false.
      END IF
      RETURN
      END

    Even though ALGOL doesn’t exist anymore, being a highly semantic language, it is still used in forms of pseudo-code to present algorithms and because of this reason, it didn’t take much to implement the code in FORTRAN. You can download the datafile from here. And also if you don’t want to run the code, here is the result file. And, as you can see, I haven’t trimmed or commented the code since I’ll feed it the updated values and move on forward, sorry for that!

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=81 1     Sub-groups of space groups http://www.emresururi.com/physics/?p=78 http://www.emresururi.com/physics/?p=78#comments Thu, 03 Jul 2008 13:06:41 +0000 http://www.emresururi.com/physics/?p=78 ISOTROPY is a highly effective program written by the guru of applications of Group Theory on Phase Transitions, Harold T. Stokes (For those who are interested in the topic, I recommend his Introduction to Isotropy Subgroups and Displacive Phase Transitions titled paper(co-authored by Dorian M. Hatch) and Structures and phase transitions in perovskites – a group-theoretical approach paper (co-authored by Christopher J. Howard).

    I’m pretty new and noobie in Group Theory – it was one of the subjects I found myself always running away. Anyway, I wanted to limit my search for some specific phase transition search in the database. Say, if an A-B binary has been reported in S1 and S2 structures, I didn’t want to go all the way looking for possible transition mechanisms if S1->S2 isn’t allowed by the group theory at all!

    In the meantime, I’ve begun experimenting with the ISOTROPY program, yet I’m still at the very beginning. You can find detailed information and an online version in the related website, but to give an example, suppose we’d like to know about the subgroups of spacegroup 123 (P4/mmm). We enter the following commands:

    VALUE PARENT 123
    SHOW PARENT
    SHOW SUBGROUP
    SHOW LATTICE
    DISPLAY ISOTROPY

    and it displays a list of which a portion is presented below:

    Parent     Lattice Subgroup   Lattice
    123 P4/mmm q       123 P4/mmm q
    123 P4/mmm q       47 Pmmm    o
    123 P4/mmm q       83 P4/m    q
    123 P4/mmm q       65 Cmmm    o-b
    123 P4/mmm q       10 P2/m    m
    123 P4/mmm q       12 C2/m    m-b
    123 P4/mmm q       2 P-1      t
    123 P4/mmm q       89 P422    q
    123 P4/mmm q       111 P-42m  q
    123 P4/mmm q       99 P4mm    q
    123 P4/mmm q       115 P-4m2  q
    123 P4/mmm q       25 Pmm2    o
    123 P4/mmm q       38 Amm2    o-b
    123 P4/mmm q       6 Pm       m
    123 P4/mmm q       123 P4/mmm q
    123 P4/mmm q       127 P4/mbm q
    123 P4/mmm q       127 P4/mbm q

    You may have noticed that, there are some recurring lines – this is due to some other properties that we have not specifically selected for display, i.e., it would display lines only containing "123 P4/mmm" if we hadn’t specifically asked for subgroup and lattice information (in analogy with the degeneracies in QM).

    I needed a table containing all the subgroups of each space-group but there were two problems:

    1.  As far as I know, ISOTROPY does not accept an input file where you can pass the commands you’ve prepared via a script/macro. So, I was afraid that I would type "VALUE PARENT 1 – DISPLAY ISOTROPY – VALUE PARENT 2 – DISPLAY ISOTROPY – …" (by the way, ISOTROPY has a clever interpreter so that you can just type the first letters of commands and parameters as long as they are not dubious).

    2. Again, as far as I know, you can’t direct the output to a file, meaning you have to copy the output manually and paste into a file (x230 in my case).

    Luckily, it turned out that:

    1. You could not supply an input file, but you could instead, paste the contents of such an input file and ISOTROPY would process it line by line.

    2. Stokes had included a very useful and helpful property that all the sessions was recorded automatically in the "iso.log" file.

    So, all I had to do was

    1. Write a script that would produce the query to retrieve the subgroups of all 230 space-groups.

    2. Paste the outcome of this script to ISOTROPY.

    3. Write another script that would parse the "iso.log" file and split these reported subgroups per spacegroup (and while I’m at it, remove the duplicate lines and order the subgroups).

    Maybe you won’t believe me but I am also not happy to paste lines and lines of code to this blog’s entries. I will deal it soon but in the mean time, please bear with them a little longer.. Ath the moment, I will describe the process and will include them at the end.

    * build the query via query_builder.php
    * feed this query to ISOTROPY
    * split the space groups from the "iso.log" via split_spacegroups.php
    * run sortall.php which :
        * runs the sortkillduplines.php on each of these files
        * sorts again via the sort.pl wrt to the 4th col (space number of the subgroup)
        * removes the blank lines
    * And that’s it..

    Here comes the codes (php+perl)

    query_builder.php:

    <?
 // Generates the query that will be "pasted" into the ISOTROPY prompt
 // The query will be stored in isotropy_query.txt
$fp = fopen("isotropy_query.txt","w");
 fwrite($fp,"PAGE 1000\nSHOW SUBGROUP\nSHOW LATTICE\nSHOW PARENT\n");
 for($i=1;$i<231;$i++)
     fwrite($fp,"VALUE PARENT $i\nDISPLAY ISOTROPY\n");
 fwrite($fp,"\n\n");
 fclose($fp);
 ?>


    sortkillduplines:

    <?
// PERL sorts the mentioned file (via $_GET["file"]) and deletes duplicate lines
require("commandline.inc.php");
if($_GET["file"]) $file = $_GET["file"];
 else
 {
     $in = fopen(’php://stdin’, ‘r’);
     $fw = fopen("sortkilldup_tmp.tmp","w");
     while(!feof($in))
     {
         fwrite($fw, rtrim(fgets($in, 4096))."\n");
     }
     fclose($fw);
     $file = "sortkilldup_tmp.tmp";
 }
 if($_GET["seperator"])$f_seperator = TRUE;//Means that an extra seperator has been placed for this order purpose and is seperated from the text by the parameter defined with seperator=xxx.
if($_GET["n"]||$_GET["numeric"]||$_GET["num"])$f_numeric = TRUE; // If it is TRUE than it is assumed that the first word of each line is a number and will be sort accordingly.
 // Reference : http://www.stonehenge.com/merlyn/UnixReview/col06.html
if($f_numeric) $exec = exec("perl -e ‘print sort numerically <>;\nsub numerically { \$a <=> \$b; }’ ".$file." > sorttemp.txt");
 else $exec = exec("perl -e ‘print sort <>’ ".$file." > sorttemp.txt");
 $fi = file("sorttemp.txt");
 unlink("sorttemp.txt");
 $curlinecont = "wrwerwer";
 if(!$f_seperator)
 {
     for($i=0; $i<sizeof($fi); $i++)
     {
         $line = rtrim($fi[$i]);
         if($curlinecont!=$line)
         {
             echo $line."\n";
             $curlinecont = $line;
         }
     }
 }
 else
 {
     $frs = fopen("sorttemp2.txt","w"); // will capture the dup killed output in this file and resort it later.
     for($i=0; $i<sizeof($fi); $i++)
     {
         $lineorg = rtrim($fi[$i]);
         $linearr = explode($_GET["seperator"],$lineorg);
         if($curlinecont!=$linearr[0])
         {
             fwrite($frs, $linearr[1]."\n");
             $curlinecont = $linearr[0];
         }
     }
     fclose($frs);
     // resort the file 
     passthru("perl -e ‘print sort <>’ sorttemp2.txt");
     unlink("sorttemp2.txt");
 }
   // unlink("sortkilldup_tmp.tmp");
 ?>


    split_spacegroups.php:

    <?
 // splits the isotropy.log file generated by the query : 
 //      PAGE 1000
 //      SHOW SUBGROUP
 //      SHOW LATTICE
 //      SHOW PARENT
 //      VALUE PARENT 1
 //      DISPLAY ISOTROPY
 //      â€¦
 //      VALUE PARENT 230
 //      DISPLAY ISOTROPY
 // with respect to the parent groups
$file = fopen("iso.log", "r") or exit("Unable to open file!");
 //Output a line of the file until the end is reached
 $cur = 0;
 $fp = fopen("dummy.txt","w"); // we need to open a dummy file in order to include the fclose() in the iteration
 while(!feof($file))
 {
     // Read the line
     $line = fgets($file);
    // Check if it is a spacegroup information entry by looking at the first word
     // - if this word is a number, than process that line
     if(ereg("^([0-9]+)",$line,$num))
     {
         $num = $num[0];
         if($cur!=$num)
         {
             // We have a new space group, so create its file:
             fclose($fp); // First close the previous space group’s file
             $filename = sprintf("%03d",$num)."_sg.txt";
             $fp = fopen($filename,"w");
             $cur = $num; // Set the current check data to the sg number
         }
         fwrite($fp,$line);
     }
}
 fclose($file);
 unlink("dummy.txt");
?>


    sortall.php:

    <?
 // sorts & kills dupes in all the spacegroup files
for($num=1;$num<231;$num++)
 {
     $filename = sprintf("%03d",$num)."_sg.txt";
     $exec = exec("php ../../toolbox/sortkillduplines.php file=".$filename." > temp");
     $exec = exec("./sort.pl temp > temp2");
     $exec = exec("perl -pi -e \"s/^\\n//\" < temp2 > ".$filename);
 }
 unlink("temp");
 unlink("temp2");
?>

    or you can directly skip all these steps and download the generated spacegroup files from here. 8)

     

    Lattice values are given in Schoenflies Notation where it corresponds to (with Pearson notation):

    T : triclinic (AP)
    M : prmitive monoclinic (MP)
    M-B : base-centered monoclinic (MC)
    O : primitive orthorhombic (OP)
    O-B : base-centered orthorhombic (OC)
    O-V : body-centered orthorhombic (OI)
    O-F : face-centered orthorhombic (OF)
    Q : primitive tetragonal (TP)
    Q-V : body-centered tetragonal (TI)
    RH : trigonal (HR)
    H : hexagonal (HP)
    C : primitive cubic (CP)
    C-F : face-centered cubic (CF)
    C-V : body-centered cubic (CI)

    Quick reference on Bravais Lattices and info on how/why they are called after Bravais but not Frankenheim..

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=78 1     The marriage of heaven (common elements) and hell (code writing code). http://www.emresururi.com/physics/?p=76 http://www.emresururi.com/physics/?p=76#respond Thu, 26 Jun 2008 10:57:09 +0000 http://www.emresururi.com/physics/?p=76 So, as Blake had said, "Energy is Eternal Delight." or in other words, let’s merge the preceding two entries to write a code for a function that produces MySQL queries that will fetch the common elements involved with the given elements:

    function fetch_elements($el,$order,$total)
{
    // produces the query to select the elements that have binaries with the $el
    // $order and $total are used to adjust the tab positions and designating the sets
 
    $abc = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
    $tab = str_repeat("\t",($total-$order));
    $lom = <<<LOM
$tab SELECT symbol_A AS symb, val1 AS val$abc[$order]
$tab FROM `dbl014`
$tab WHERE symbol_B = "$el"
$tab GROUP BY Symbol_A
$tab UNION
$tab SELECT symbol_B AS symb, val1 AS val$abc[$order]
$tab FROM `dbl014`
$tab WHERE (
$tab     symbol_A = "$el"
$tab     AND symbol_B != ""
$tab )
$tab GROUP BY Symbol_B
LOM;
    return $lom;
}
 
function merge_queries($str1,$str2,$or1,$or2,$tab)
{
    // merges two node (or node group) queries (formed via fetch_elements() function)
    // using the INNER JOIN function
    return "$tab SELECT *\n$tab FROM (\n".$str1."\n$tab ) AS $or1\n$tab INNER JOIN (\n".$str2."\n$tab ) AS $or2\n$tab USING (symb)\n$tab ";
}
 
function fetch_common_elements($ar_el)
{
    // produces the query that selects the common elements that have binaries with all of the elements given in the $ar_el array (using their symbols)
    $total = sizeof($ar_el);
    $abc = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
    //$tab = str_repeat("\t",($total-$order-1));
    $str[0] = fetch_elements($ar_el[0],0,$total-1);
    $str[1] = fetch_elements($ar_el[1],1,$total);
    $stri = merge_queries($str[0],$str[1],"A","B",str_repeat("\t",($total-2)));
    for($i=2;$i<$total;$i++)
    {
        $str[$i] = fetch_elements($ar_el[$i],$i,$total);
        $stri = merge_queries($str[$i],$stri,$abc[$i],$abc[($i+1)],str_repeat("\t",($total-$i-1)));
    }
    return $stri;
}

    Now, we can call it with, for example:
    echo fetch_common_elements(Array("Ga","Fe","Gd","Y","Ti","Ge","Ce"));
    Which will output:

     SELECT *
 FROM (
         SELECT symbol_A AS symb, val1 AS valG
         FROM `dbl014`
         WHERE symbol_B = "Ce"
         GROUP BY Symbol_A
         UNION
         SELECT symbol_B AS symb, val1 AS valG
         FROM `dbl014`
         WHERE (
             symbol_A = "Ce"
             AND symbol_B != ""
         )
         GROUP BY Symbol_B
 ) AS G
 INNER JOIN (
         SELECT *
         FROM (
                 SELECT symbol_A AS symb, val1 AS valF
                 FROM `dbl014`
                 WHERE symbol_B = "Ge"
                 GROUP BY Symbol_A
                 UNION
                 SELECT symbol_B AS symb, val1 AS valF
                 FROM `dbl014`
                 WHERE (
                     symbol_A = "Ge"
                     AND symbol_B != ""
                 )
                 GROUP BY Symbol_B
         ) AS F
         INNER JOIN (
                 SELECT *
                 FROM (
                         SELECT symbol_A AS symb, val1 AS valE
                         FROM `dbl014`
                         WHERE symbol_B = "Ti"
                         GROUP BY Symbol_A
                         UNION
                         SELECT symbol_B AS symb, val1 AS valE
                         FROM `dbl014`
                         WHERE (
                             symbol_A = "Ti"
                             AND symbol_B != ""
                         )
                         GROUP BY Symbol_B
                 ) AS E
                 INNER JOIN (
                         SELECT *
                         FROM (
                                 SELECT symbol_A AS symb, val1 AS valD
                                 FROM `dbl014`
                                 WHERE symbol_B = "Y"
                                 GROUP BY Symbol_A
                                 UNION
                                 SELECT symbol_B AS symb, val1 AS valD
                                 FROM `dbl014`
                                 WHERE (
                                     symbol_A = "Y"
                                     AND symbol_B != ""
                                 )
                                 GROUP BY Symbol_B
                         ) AS D
                         INNER JOIN (
                                 SELECT *
                                 FROM (
                                         SELECT symbol_A AS symb, val1 AS valC
                                         FROM `dbl014`
                                         WHERE symbol_B = "Gd"
                                         GROUP BY Symbol_A
                                         UNION
                                         SELECT symbol_B AS symb, val1 AS valC
                                         FROM `dbl014`
                                         WHERE (
                                             symbol_A = "Gd"
                                             AND symbol_B != ""
                                         )
                                         GROUP BY Symbol_B
                                 ) AS C
                                 INNER JOIN (
                                         SELECT *
                                         FROM (
                                                 SELECT symbol_A AS symb, val1 AS valA
                                                 FROM `dbl014`
                                                 WHERE symbol_B = "Ga"
                                                 GROUP BY Symbol_A
                                                 UNION
                                                 SELECT symbol_B AS symb, val1 AS valA
                                                 FROM `dbl014`
                                                 WHERE (
                                                     symbol_A = "Ga"
                                                     AND symbol_B != ""
                                                 )
                                                 GROUP BY Symbol_B
                                         ) AS A
                                         INNER JOIN (
                                                 SELECT symbol_A AS symb, val1 AS valB
                                                 FROM `dbl014`
                                                 WHERE symbol_B = "Fe"
                                                 GROUP BY Symbol_A
                                                 UNION
                                                 SELECT symbol_B AS symb, val1 AS valB
                                                 FROM `dbl014`
                                                 WHERE (
                                                     symbol_A = "Fe"
                                                     AND symbol_B != ""
                                                 )
                                                 GROUP BY Symbol_B
                                         ) AS B
                                         USING (symb)
 
                                 ) AS D
                                 USING (symb)
 
                         ) AS E
                         USING (symb)
 
                 ) AS F
                 USING (symb)
 
         ) AS G
         USING (symb)
 
 ) AS H
 USING (symb)

    and when fed to the database will yield:

     

     thank y’all.. 8)

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=76 0     A code to write code (as in a war to end -all- wars – 8P) http://www.emresururi.com/physics/?p=74 http://www.emresururi.com/physics/?p=74#respond Thu, 19 Jun 2008 15:03:53 +0000 http://www.emresururi.com/physics/?p=74 Lately, I’ve been studying on the nice and efficient algorithm of Dominik Endres and he’s been very helpful in replying to my relentless mails about the algorithm and the software. Now, if you have checked the algorithm, there is a part where iteration is needed to be done over the possible bin boundary distributions. Via a smart workaround he avoids this and uses a recursive equation.

    For my needs, I want to follow the brute-force approach. To requote the situation, suppose we have K=4 number of classes, let’s label these as 0,1,2,3. There are 3 different ways to group these classes into 2 (M=2-1=1) bins obeying the following rules:

    * There is always at least one classification per bin.

    * Bins are adjascent to each other.

    * Bins can not overlap.

    The number of configurations are given by the formula:

    N=(K-1)! / ( (K-M-1)! M!)   (still haven’t brought back my latex renderer, sorry for that!)

     

    and these 3 possible configurations for K=4, M=1 are:

    |—-|—-|—-|

    0            3  : Meaning, the first bin only covers the 0 class while the second bin covers 1, 2, 3.

         1       3 : Meaning, the first bin covers 0, 1 while the second bin covers 2, 3.

             2   3 : Meaning, the first bin covers 0, 1, 2 while the second bin covers 3.

     

    Now, think that you are writing a code to get this values, it is easy when you pseudo-code it:

    for k_0 = 0 : K-1-M

       for k_1 = k_0+1 : K-M

          …

             for k_(K-3) = k_(K-4)+1 : K-4

                for k_(K-2) = k_(K-3)+1 : K-3

                   print("k_0 – k_1 – … – k_(K-2) – K-1")

                endfor

             endfor

          …

       endfor

    endfor

    easier said than done because the number of the k_ variables and the for loops are changing with M.

    So after some trying in Octave for recursive loops with the help of the feval() function, I gave up and wrote a PHP code that writes an Octave code and it was really fun! 8)

    <?
    //Support for command line variable passing:
    for($i=1;$i<sizeof($_SERVER["argv"]);$i++)
    {
        list($var0,$val0) = explode("=",$_SERVER["argv"][$i]);
        $_GET[$var0] = $val0;
    }

    $K = $_GET["K"];
    $M = $_GET["M"];

    if(!$_GET["K"])$_GET["K"] = 10;
    if(!$_GET["M"])$_GET["M"] = 2;
    if(!$_GET["f"])$_GET["f"] = "lol.m";

    $fp = fopen($_GET["f"],"w");
    fwrite($fp, "count=0;\n");

    fwrite($fp, str_repeat("\t",$m)."for k0 = 0:".($K-$M-1)."\n");

    for($m=1;$m<$M;$m++)
    {
        fwrite($fp, str_repeat("\t",$m)."for k".$m." = k".($m-1)."+1:".($K-$M-1+$m)."\n");
    }

    fwrite($fp, str_repeat("\t",$m)."printf(\"".str_repeat("%d – ",$M)."%d\\n\",");

    for($i=0;$i<$m;$i++)
        fwrite($fp, "k".$i.",");
    fwrite($fp,($K-1).")\n");

    fwrite($fp, str_repeat("\t",$m)."count++;\n");

    for($m=$M-1;$m>=0;$m–)
    {
        fwrite($fp, str_repeat("\t",$m)."endfor\n");
    }

    fwrite($fp,"if(count!=factorial(".$K."-1)/(factorial(".$K."-".$M."-1)*factorial(".$M.")))\n\tprintf(\"Houston, we have a problem..\");\nendif\n");

    //fwrite($fp,"printf(\"count : %d\\nfactorial : %d\\n\",count,factorial(".$K."-1)/(factorial(".$K."-".$M."-1)*factorial(".$M.")));\n");

    fclose($fp);

    system("octave ".$_GET["f"]);
    ?>

    Now, if you call this code with for example:

    > php genfor.php K=4 M=1

    it produces and runs this Octave file lol.m:

    count=0;
    for k0 = 0:2
        printf("%d – %d\n",k0,3)
        count++;
    endfor
    if(count!=factorial(4-1)/(factorial(4-1-1)*factorial(1)))
        printf("Houston, we have a problem..");
    endif

    which outputs:

    0 – 3
    1 – 3
    2 – 3

    And you get to keep your hands clean 8)

     

    Later addition (27/6/2008) You can easily modify the code to also give you possible permutations for a set with K items and groups of M items (where ordering is not important). For example, for K=5, say {0,1,2,3,4} and M=4:

    0-1-2-3
    0-1-2-4
    0-1-3-4
    0-2-3-4
    1-2-3-4

    5 (=K!/(M!(K-M)!)) different sets.

    The modifications we make to the code are:

    * Add an optional parameter that will tell the code that, we want permutation instead of binning:

    if($_GET["perm"]||$_GET["p"]||$_GET["P"])
    {
        $f_perm = TRUE;
        $K+=1;
    }

    via the f_perm flag, we can check the type of the operation anywhere within the code. K is incremented because, previously, we were setting the end point fixed so that the last bin would always include the last point for the bins to cover all the space. But for permutationing we are not bounded by this fixication so, we have an additional freedom.

    * Previously, we were just including the last point in the print() statement, we should fix this if we are permutating:

    $str = "";

    for($i=0;$i<$m;$i++)
        $str .= "k".$i.",";
    if(!$f_perm) fwrite($fp,$str.($K-1).")\n");
    else fwrite($fp, substr($str,0,-1).")\n");

    The reason for instead of now writing directly to the file in the iteration we just store it in a variable (str) is to deal with the trailing comma if we are permuting, as simple as that!

    and the difference between before and after the modifications is:

    sururi@dutsm0175 Octave $ diff genfor.php genfor2.php
    10a11,16
    > if($_GET["perm"]||$_GET["p"]||$_GET["P"])
    > {
    >     $f_perm = TRUE;
    >     $K+=1;
    > }
    >
    21c27,30
    < fwrite($fp, str_repeat("\t",$m)."printf(\"".str_repeat("%d – ",$M)."%d\\n\",");
    —
    > if(!$f_perm)fwrite($fp, str_repeat("\t",$m)."printf(\"".str_repeat("%d – ",$M)."%d\\n\",");
    > else fwrite($fp, str_repeat("\t",$m)."printf(\"".str_repeat("%d – ",($M-1))."%d\\n\",");
    >
    > $str = "";
    24,25c33,35
    <     fwrite($fp, "k".$i.",");
    < fwrite($fp,($K-1).")\n");
    —
    >     $str .= "k".$i.",";
    > if(!$f_perm) fwrite($fp,$str.($K-1).")\n");
    > else fwrite($fp, substr($str,0,-1).")\n");

    for the really lazy people (like myself), here is the whole deal! 8)

    <?
     //Support for command line variable passing:
    for($i=1;$i<sizeof($_SERVER["argv"]);$i++)
    {
        list($var0,$val0) = explode("=",$_SERVER["argv"][$i]);
        $_GET[$var0] = $val0;
    }


    $K = $_GET["K"];
    $M = $_GET["M"];

    if(!$_GET["K"])$_GET["K"] = 10;
    if(!$_GET["M"])$_GET["M"] = 2;
    if(!$_GET["f"])$_GET["f"] = "lol.m";

    if($_GET["perm"]||$_GET["p"]||$_GET["P"])
    {
        $f_perm = TRUE;
        $K+=1;
    }

    $fp = fopen($_GET["f"],"w");
    fwrite($fp, "count=0;\n");

    fwrite($fp, str_repeat("\t",$m)."for k0 = 0:".($K-$M-1)."\n");

    for($m=1;$m<$M;$m++)
    {
        fwrite($fp, str_repeat("\t",$m)."for k".$m." = k".($m-1)."+1:".($K-$M-1+$m)."\n");
    }

    if(!$f_perm)fwrite($fp, str_repeat("\t",$m)."printf(\"".str_repeat("%d – ",$M)."%d\\n\",");
    else fwrite($fp, str_repeat("\t",$m)."printf(\"".str_repeat("%d – ",($M-1))."%d\\n\",");

    $str = "";

    for($i=0;$i<$m;$i++)
        $str .= "k".$i.",";
    if(!$f_perm) fwrite($fp,$str.($K-1).")\n");
    else fwrite($fp, substr($str,0,-1).")\n");

    fwrite($fp, str_repeat("\t",$m)."count++;\n");

    for($m=$M-1;$m>=0;$m–)
    {
        fwrite($fp, str_repeat("\t",$m)."endfor\n");
    }

    fwrite($fp,"if(count!=factorial(".$K."-1)/(factorial(".$K."-".$M."-1)*factorial(".$M.")))\n\tprintf(\"Houston, we have a problem..\");\nendif\n");

    fwrite($fp,"printf(\"count : %d\\nfactorial : %d\\n\",count,factorial(".$K."-1)/(factorial(".$K."-".$M."-1)*factorial(".$M.")));\n");

    fclose($fp);

    system("octave ".$_GET["f"]);
    ?>

    After the modifications,

    php genfor.php K=5 M=3 p=1

    will produce

    0 – 1 – 2
    0 – 1 – 3
    0 – 1 – 4
    0 – 2 – 3
    0 – 2 – 4
    0 – 3 – 4
    1 – 2 – 3
    1 – 2 – 4
    1 – 3 – 4
    2 – 3 – 4

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=74 0     Specification of an extensible and portable file format for electronic structure and crystallographic data http://www.emresururi.com/physics/?p=67 http://www.emresururi.com/physics/?p=67#respond Tue, 06 May 2008 14:09:17 +0000 http://www.emresururi.com/physics/?p=67 Specification of an extensible and portable file format for electronic structure and crystallographic data

    X. Gonzea, C.-O. Almbladh, A. Cucca, D. Caliste, C. Freysoldt, M.A.L. Marques, V. Olevano, Y. Pouillon and M.J. Verstraet

    Comput. Mater. Sci. (2008) doi:10.1016/j.commatsci.2008.02.023

    Abstract

    In order to allow different software applications, in constant evolution, to interact and exchange data, flexible file formats are needed. A file format specification for different types of content has been elaborated to allow communication of data for the software developed within the European Network of Excellence “NANOQUANTA”, focusing on first-principles calculations of materials and nanosystems. It might be used by other software as well, and is described here in detail. The format relies on the NetCDF binary input/output library, already used in many different scientific communities, that provides flexibility as well as portability across languages and platforms. Thanks to NetCDF, the content can be accessed by keywords, ensuring the file format is extensible and backward compatible.

    Keywords

    Electronic structure; File format standardization; Crystallographic datafiles; Density datafiles; Wavefunctions datafiles; NetCDF

    PACS classification codes

    61.68.+n; 63.20.dk; 71.15.Ap

    Problems with the binary representations:

    1. lack of portability between big-endian and little-endian platforms (and vice-versa), or between 32-bit and 64-bit platforms;
    2. difficulties to read the files written by F77/90 codes from C/C++ software (and vice-versa)
    3. lack of extensibility, as one file produced for one version of the software might not be readable by a past/forthcoming version

    (Network Common Data Form) library solves these issues. (Alas, it is binary and in the examples presented in  http://www.unidata.ucar.edu/software/netcdf/examples/files.html, it seems like the CDF representations are larger in size than those of the CDL metadata representations)

    The idea of standardization of file formats is not new in the electronic structure community [X. Gonze, G. Zerah, K.W. Jakobsen, and K. Hinsen. Psi-k Newsletter 55, February 2003, pp. 129-134. URL: <http://psi-k.dl.ac.uk>] (This article is a must read article, by the way. It’s a very insightful and concerning overview of the good natured developments in the atomistic/molecular software domain. It is purifying in the attempt to produce something good..). However, it proved difficult to achieve without a formal organization, gathering code developers involved in different software project, with a sufficient incentive to realize effective file exchange between such software.

    HDF (Hierarchical Data Format) is also an alternative. NetCDF is simpler to use if the data formats are flat, while HDF has definite advantages if one is dealing with hierarchical formats. Typically, we will need to describe many different multi-dimensional arrays of real or complex numbers, for which NetCDF is an adequate tool.

    Although a data specification might be presented irrespective of the library actually used for the implementation, such a freedom might lead to implementations using incompatible file formats (like NetCDF and XML for instance). This possibility would in effect block the expected standardization gain. Thus, as a part of the standardization, we require the future implementations of our specification to rely on the NetCDF library only. (Alternative to independent schemas is presented as a mandatory format which is against the whole idea of development. The reason for NetCDF being incompatible with XML lies solely on the inflexibility/inextensibility of the prior format. Although such a readily defined and built format is adventageous considering the huge numerical data and parameters the ab inito software uses, it is very immobile and unnecessary for compound and element data.)

    The compact representation brought by NetCDF can be by-passed in favour of text encoding (very agreable given the usage purposes: XML type extensible schema usage is much more adequate for structure/composition properties). We are aware of several standardization efforts [J. Junquera, M. Verstraete, X. Gonze, unpublished] and [J. Mortensen, F. Jollet, unpublished. Also check Minutes of the discussion on XML format for PAW setups], relying on XML, which emphasize addressing by content to represent such atomic data.

    4 types of data types are distinguished:

    • (A) The actual numerical data (which defines whether a file contains wavefunctions, a density, etc), for which a name must have been agreed in the specification.
    • (B) The auxiliary data that is mandatory to make proper usage of the actual numerical data of A-type. The name and description of this auxiliary information is also agreed.
    • (C) The auxiliary data that is not mandatory to make proper usage of the A-type numerical data, but for which a name and description has been agreed in the specification.
    • (D) Other data, typically code-dependent, whose availability might help the use of the file for a specific code. The name of these variables should be different from the names chosen for agreed variables of A–C types. Such type D data might even be redundant with type A–C data.

    The NetCDF interface adapts the dimension ordering to the programming language used. The notation here is C-like, i.e. row-major storage, the last index varying the fastest. In FORTRAN, a similar memory mapping is obtained by reversing the order of the indices. (So, the ordering/reverse-ordering is handled by the interface/library)

    Concluding Remarks

    We presented the specifications for a file format relying on the NetCDF I/O library with a content related to electronic structure and crystallographic data. This specification takes advantage of all the interesting properties of NetCDF-based files, in particular portability and extensibility. It is designed for both serial and distributed usage, although the latter characteristics was not presented here.

    Several software in the Nanoquanta and ETSF [15] context can produce or read this file format: ABINIT [16] and [17], DP [18], GWST, SELF [19], V_Sim [20]. In order to further encourage its use, a library of Fortran routines [5] has been set up, and is available under the GNU LGPL licence.

    Additional Information:

    Announcement for a past event
    ( http://www.tddft.org/pipermail/fsatom/2003-February/000004.html )

    CECAM – psi-k – SIMU joint Tutorial

    1) Title : Software solutions for data exchange and code gluing.

    Location : Lyon Dates : 8-10 october, 2003

    Purpose : In this tutorial, we will teach software tools and standards that have recently emerged in view of the exchange of data (text and binary) and gluing of codes : (1) Python, as scripting language, its interfaces with C and FORTRAN ; (2) XML, a standard for representing structured data in text files ; (3) netCDF, a library and file format for the exchange and storage of binary data, and its interfaces with C, Fortran, and Python

    Organizers : X. Gonze gonze@pcpm.ucl.ac.be

    K. Hinsen hinsen@cnrs-orleans.fr

    2) Scientific content

    Recent discussions, related to the CECAM workshop on "Open Source Software for Microscopic Simulations", June 19-21, 2002, to the GRID concept (http://www.gridcomputing.com), as well as to the future Integrated Infrastructure Initiative proposal linked to the European psi-k network (http://psi-k.dl.ac.uk), have made clear that one challenge for the coming years is the ability to establish standards for accessing codes, transferring data between codes, testing codes against each other, and become able to "glue" them (this being facilitated by the Free Software concept).

    In the present tutorial, we would like to teach three "software solutions" to face this challenge : Python, XML and netCDF.

    Python is now the de facto "scripting langage" standard in the computational physics and chemistry community. XML (eXtended Markup Language) is a framework for building mark-up languages, allowing to set-up self-describing documents, readable by humans and machines. netCDF allows binary files to be portable accross platforms. It is not our aim to cover all possible solutions to the above-mentioned challenges (e.g. PERL, Tcl, or HDF), but these three have proven suitable for atomic-scale simulations, in the framework of leading projects like CAMPOS (http://www.fysik.dtu.dk/campos), MMTK (http://dirac.cnrs-orleans.fr/MMTK), and GROMACS (http://www.gromacs.org). Other software projects like ABINIT (http://www.abinit.org) and PWSCF (http://www.pwscf.org – in the DEMOCRITOS context), among others, have made clear their interest for these. All of these software solutions can be used without having to buy a licence.

    Tentative program of the tutorial. Lectures in the morning, hands-on training in the afternoon.

    1st day ——- 2h Python basics 1h Interface : Python/C or FORTRAN 1h XML basics Afternoon Training with Python, and interfaces with C and FORTRAN

    2nd day ——- 2h Python : object oriented (+ an application to GUI and Tk) 1h Interface : Python/XML 1h Interface : XML + C or FORTRAN Afternoon Training with XML + interfaces

    3rd day ——- 1h Python : numerical 1h netCDF basics 1h Interface : netCDF/Python 1h Interface : netCDF/C or FORTRAN Afternoon Training with netCDF + interfaces

    3) List of lecturers

    K. Hinsen (Orleans, France), organizer X. Gonze (Louvain-la-Neuve, Belgium), organizer K. Jakobsen (Lyngby, Denmark), instructor J. Schiotz (Lyngby, Denmark), instructor J. Van Der Spoel (Groningen, The Netherlands), instructor M. van Loewis (Berlin, Germany), instructor

    4) Number of participants : around 20, Most of the participants should be PhD students, postdoc or young permanent scientists, involved in code development. It is assumed that the attendants have a good knowledge of UNIX, and C or FORTRAN.

    Our budget will allow contributing to travel and local expenses of up to 20 participants.

    XML and NetCDF:

    from Specification of file formats for NANOQUANTA/ETSF Specification – www.etsf.eu/research/software/nq_specff_v2.1_final.pdf

    Section 1. General considerations concerning the present file format specifications.

    One has to consider separately the set of data to be included in each of different types of files, from their representation. Concerning the latter, one encounters simple text files, binary files, XML-structured files, NetCDF files, etc … It was already decided previously (Nanoquanta meeting Maratea Sept. 2004) to evolve towards formats that deal appropriately with the self- description issue, i.e. XML and NetCDF. The inherent flexibility of these representations will also allow to evolve specific versions of each type of files progressively, and refine earlier working proposals. The same direction has been adopted by several groups of code developers that we know of.

    Information on NetCDF and XML can be obtained from the official Web sites,

    http://www.unidata.ucar.edu/software/netcdf/ and

    http://www.w3.org/XML/

    There are numerous other presentations of these formats on the Web, or in books.

    The elaboration of file formats based on NetCDF has advanced a lot during the Louvain-la- Neuve mini-workshop. There has been also some remarks about XML.

    Concerning XML :

    (A) The XML format is most adapted for the structured representation of relatively small quantity of data, as it is not compressed.

    (B) It is a very flexible format, but hard to read in Fortran (no problem in C, C++ or Python). Recently, Alberto Garcia has set up a XMLF90 library of routines to read XML from Fortran. http://lcdx00.wm.lc.ehu.es/~wdpgaara/xml/index.html Other efforts exists in this direction http://nn-online.org/code/xml/

    Concerning NetCDF

    • (A) Several groups of developers inside NQ have already a good experience of using it, for the representation of binary data (large files).
    • (B) Although there is no clear advantage of NetCDF compared to HDF (another possibility for large binary files), this experience inside the NQ network is the main reason for preferring it. By the way, NetCDF and HDF are willing to merge (this might take a few years, though).
    • (C) File size limitations of NetCDF exist, see appendix D, but should be overcome in the future.

    Thanks to the flexibility of NetCDF, the content of a NetCDF file format suitable for use for NQ softwares might be of four different types :

    (1) The actual numerical data (that defines a file for wavefunctions, or a density file, etc …), whose NetCDF description would have been agreed.

    (2) The auxiliary data that are mandatory to make proper usage of the actual numerical data. The NetCDF description of these auxiliary data should also be agreed.

    (3) The auxiliary data that are not mandatory, but whose NetCDF description has been agreed, in a larger context.

    (4) Other data, typically code-dependent, whose existence might help the use of the file for a specific code.

    References:
    [5] URL: <http://www.etsf.eu/index.php?page=tools>.

    [15] URL: <http://www.etsf.eu/>.

    [16] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J.-Y. Raty and D.C. Allan, Comput. Mater. Sci. 25 (2002), pp. 478–492.

    [17] X. Gonze, G.-M. Rignanese, M. Verstraete, J.-M. Beuken, Y. Pouillon, R. Caracas, F. Jollet, M. Torrent, G. Zérah, M. Mikami, Ph. Ghosez, M. Veithen, J.-Y. Raty, V. Olevano, F. Bruneval, L. Reining, R. Godby, G. Onida, D.R. Hamann and D.C. Allan, Zeit. Kristall. 220 (2005), pp. 558–562.

    [18] URL: <http://dp-code.org>.

    [19] URL: <http://www.bethe-salpeter.org>.

    [20] URL: http://www-drfmc.cea.fr/sp2m/L_Sim/V_Sim/index.en.html.

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=67 0     Element Arrays http://www.emresururi.com/physics/?p=65 http://www.emresururi.com/physics/?p=65#respond Thu, 10 Apr 2008 09:05:18 +0000 http://www.emresururi.com/physics/?p=65 For future reference when coding, I thought it would be useful to include these arrays. All of the arrays are ordered with respect to their Atomic Numbers.

     

    Atomic Numbers Array

    "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106", "107", "108", "109", "110", "111", "112", "113", "114", "115", "116", "117", "118"

     

    Symbols Array

    "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Uub", "Uut", "Uuq", "Uup", "Uuh", "Uus", "Uuo"

     

    Names Array

    "Hydrogen", "Helium", "Lithium", "Beryllium", "Boron", "Carbon", "Nitrogen", "Oxygen", "Fluorine", "Neon", "Sodium", "Magnesium", "Aluminum", "Silicon", "Phosphorus", "Sulfur", "Chlorine", "Argon", "Potassium", "Calcium", "Scandium", "Titanium", "Vanadium", "Chromium", "Manganese", "Iron", "Cobalt", "Nickel", "Copper", "Zinc", "Gallium", "Germanium", "Arsenic", "Selenium", "Bromine", "Krypton", "Rubidium", "Strontium", "Yttrium", "Zirconium", "Niobium", "Molybdenum", "Technetium", "Ruthenium", "Rhodium", "Palladium", "Silver", "Cadmium", "Indium", "Tin", "Antimony", "Tellurium", "Iodine", "Xenon", "Cesium", "Barium", "Lanthanum", "Cerium", "Praseodymium", "Neodymium", "Promethium", "Samarium", "Europium", "Gadolinium", "Terbium", "Dysprosium", "Holmium", "Erbium", "Thulium", "Ytterbium", "Lutetium", "Hafnium", "Tantalum", "Tungsten", "Rhenium", "Osmium", "Iridium", "Platinum", "Gold", "Mercury", "Thallium", "Lead", "Bismuth", "Polonium", "Astatine", "Radon", "Francium", "Radium", "Actinium", "Thorium", "Protactinium", "Uranium", "Neptunium", "Plutonium", "Americium", "Curium", "Berkelium", "Californium", "Einsteinium", "Fermium", "Mendelevium", "Nobelium", "Lawrencium", "Rutherfordium", "Dubnium", "Seaborgium", "Bohrium", "Hassium", "Meitnerium", "Darmstadtium", "Roentgenium", "Ununbium", "Ununtrium", "Ununquadium", "Ununpentium", "Ununhexium", "Ununseptium", "Ununoctium"

     

    Pettifor’s Mendeleev Numbers Array
    D.G. Pettifor, Mater. Sci. Tech. 4 (1988) 2480 | D.G. Pettifor, Bonding and Structure of Molecules and Solids, Oxford University Press (1995) p.13

    "103", "1", "12", "77", "86", "95", "100", "101", "102", "2", "11", "73", "80", "85", "90", "94", "99", "3", "10", "16", "19", "51", "54", "57", "60", "61", "64", "67", "72", "76", "81", "84", "89", "93", "98", "4", "9", "15", "25", "49", "53", "56", "59", "62", "65", "69", "71", "75", "79", "83", "88", "92", "97", "5", "8", "14", "33", "32", "31", "30", "29", "28", "18", "27", "26", "24", "23", "22", "21", "17", "20", "50", "52", "55", "58", "63", "66", "68", "70", "74", "78", "82", "87", "91", "96", "6", "7", "13", "48", "47", "46", "45", "44", "43", "42", "41", "40", "39", "38", "37", "36", "35", "34", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"

     

    Villars’ Mendeleev Numbers Array
    P. Villars et al., J. Alloys Compd. 367 (2004) 167] – doi:10.1016/j.jallcom.2003.08.060

    "92", "98", "1", "67", "72", "77", "82", "87", "93", "99", "2", "68", "73", "78", "83", "88", "94", "100", "3", "7", "11", "43", "46", "49", "52", "55", "58", "61", "64", "69", "74", "79", "84", "89", "95", "101", "4", "8", "12", "44", "47", "50", "53", "56", "59", "62", "65", "70", "75", "80", "85", "90", "96", "102", "5", "9", "13", "15", "17", "19", "21", "23", "25", "27", "29", "31", "33", "35", "37", "39", "41", "45", "48", "51", "54", "57", "60", "63", "66", "71", "76", "81", "86", "91", "97", "103", "6", "10", "14", "16", "18", "20", "22", "24", "26", "28", "30", "32", "34", "36", "38", "40", "42", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0"

    You can also download the SQL version for these values from here.

    Hope it helps..

    ]]>     http://www.emresururi.com/physics/?feed=rss2&p=65 0