Consider the following 16×16 matrix:

   0.06483   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.65809   0.65809   0.60365   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.10993   0.00000   0.00000   0.37862   0.37862   0.00000   0.00000   0.00000   0.00000   0.86568   0.00000   0.00000   0.42424   0.00000   0.00000
   0.00000   0.00000   0.10993   0.00000   0.00000   0.12949   0.00000   0.00000   0.00000   0.00000   0.00000  -0.86568   0.00000   0.00000   0.42424   0.00000
   0.00000   0.00000   0.00000   0.10259   0.00000   0.00000   0.84469   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.73695
   0.00000   0.37862   0.00000   0.00000   0.01584   0.00000   0.00000   0.00000   0.00000   0.00000   0.84812   0.00000   0.00000   0.87983   0.00000   0.00000
   0.00000   0.37862   0.12949   0.00000   0.00000   0.01584   0.00000   0.00000   0.00000   0.00000   0.00000  -0.84812   0.00000   0.00000   0.87983   0.00000
   0.00000   0.00000   0.00000   0.84469   0.00000   0.00000   0.01168   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.87050
   0.65809   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.88440   0.28986   0.45572   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.65809   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.28986   0.88440   0.45572   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.60365   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.45572   0.45572   0.89801   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.86568   0.00000   0.00000   0.84812   0.00000   0.00000   0.00000   0.00000   0.00000   0.35993   0.00000   0.00000   0.19336   0.00000   0.00000
   0.00000   0.00000  -0.86568   0.00000   0.00000  -0.84812   0.00000   0.00000   0.00000   0.00000   0.00000   0.35993   0.00000   0.62239  -0.62239   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   1.18908   0.00000   0.00000   0.00000
   0.00000   0.42424   0.00000   0.00000   0.87983   0.00000   0.00000   0.00000   0.00000   0.00000   0.19336   0.62239   0.00000   0.54930   0.00000   0.00000
   0.00000   0.00000   0.42424   0.00000   0.00000   0.87983   0.00000   0.00000   0.00000   0.00000   0.00000  -0.62239   0.00000   0.00000   0.54930   0.00000
   0.00000   0.00000   0.00000   0.73695   0.00000   0.00000   0.87050   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.55094

Equality relations:
(1,9) = (1,8)
(3,3) = (2,2)
(2,6) = (2,5)
(3,12) = -(2,11)
(3,15) = (2,14)
(6,6) = (5,5)
(6,12) = -(5,11)
(6,15) = (5,14)
(9,9) = (8,8)
(9,10) = (8,10)
(12,12) = (11,11)
(12,15) = -(12,14)
(15,15) = (14,14)

These kind of relations are sought upon tensor matrices (at least, in group theory 8), kind of postulating the imposed limits and… (have I said “relations”?). So, the sort-of-important-like thingy is to be able to spot them after evaluating some calculations in the computers.

Have you verified yet? Good! Now let’s take a look at this:

So, this is what we’ll be discussing in this post – preparation of such graphs. We’ll be using the PHP module Imagick for ImageMagick functions and objects. On Ubuntu, you can install the package from the command line as:

sudo apt-get install php5-imagick

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

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:

not bad for beginning, right?

The circle function has 4 parameters: the (x,y) coordinates of the center and the coordinates of a point that is on the edge of the circle. Let’s draw another circle, this time an unfilled one below this filled one:

# Declare the object (/holder) 
$image = new Imagick();

 # Define the canvas (size 640x640 pixels with white background)
 $image --->newImage( 640, 640, new ImagickPixel( 'white' ) );

# Define the 1st, filled circle:
$circle = new ImagickDraw();
$circle->setFillColor("#000000");
$circle->circle(320,320,350,320);

# Add (/register) it to the canvas:
$image->drawImage($circle);

# Define the 2nd, unfilled circle:
$circle = new ImagickDraw();
$circle->setStrokeColor("#000000");
$circle->setStrokeWidth(2);
$circle->setFillColor("none");
$circle->circle(320,350,350,350);

# Add (/register) it to the canvas:
$image->drawImage($circle);

# Designate the image type as PNG:
$image->setImageFormat( "png" );

# Export (/write) it to a file "out.png":
$fp = fopen("out.png","w");
fwrite($fp,$image);
fclose($fp);

The “stroke” methods are responsible for the perimeter and the fillcolor setting to “none” ensures that we have an otherwise invisible circle.

So far, so good.. Let’s add a line that passes between the circles from (120,50) to (520,50):

# Define the line:
$line = new ImagickDraw();
$line->setStrokeColor("#000000");
$line->setStrokeWidth(2);
$line->line(120,350,520,350);
$image->drawImage($line);

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:

doesn’t look right, right? And here’s the code that produced it:

# Declare the object (/holder)
$image = new Imagick(); 

# Define the canvas (size 640x640 pixels with white background) 
$image--->newImage( 640, 640, new ImagickPixel( 'white' ) );

# Define the 1st, filled circle:
$dxy = 640 / 16;
for($i=1;$i<=16;$i++)
{
for($j=1;$j<=16;$j++)
     {
         $circle = new ImagickDraw();
         $circle->setFillColor("#000000");
         $circle->circle($dxy*($j-1),$dxy*($i-1),$dxy*($j-1)+$dxy/6,$dxy*($i-1));

         # Add (/register) it to the canvas:
         $image->drawImage($circle);
     }
}

# Designate the image type as PNG:
$image->setImageFormat( "png" );

# Export (/write) it to a file "out.png":
$fp = fopen("out.png","w");
fwrite($fp,$image);
fclose($fp);

What we need to do is to shift it half of the periodicity in the x & y directions:

That’s more like it — this shift was introduced to the code by changing the coordinate parameters of the circle(s):
$circle->circle($dxy*($j-1/2),$dxy*($i-1/2),$dxy*($j-1/2)+$dxy/6,$dxy*($i-1/2));

So, depending on the value (zero/non-zero) of an element of our tensor, we can put it as a big circle (r=$dxy/6) or a small one (r=$dxy/12). In addition, if it’s negative, we can have it drawn as an unfilled one (

$circle->setStrokeColor("#000000");
$circle->setStrokeWidth(2);
$circle->setFillColor("none");

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

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

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/

How to Prepare an Input File for Surface Calculations

Emre S. Tasci

Description

Tools

Procedure

Initial structural data

Figure 1 VESTA’s “Unit cell” interface

Figure 2 VESTA’s “Structure parameters” interface

Figure 3 Fm3m platinum carbide

Constructing the Supercell

Figure 4 Transformation matrix for the supercell

Figure 5 A 3x3x3 supercell

,

Figure 6 Designating the boundaries for the construction of the supercell.

Figure 7 (a)Single unit cell (conventional); (b) 3x3x3 supercell

Lattice Planes

etasci@metu.edu.tr

Figure 8 Selecting the (111) lattice plane

The plane’s distance to the origin can be specified in terms of Å or interplane distance, as well as selecting 3 or more atoms lying on the plane (multiple selections can be realised by holding down SHIFT while clicking on the atoms) and then clicking on the “Calculate the best plane for the selected atoms” button to have the plane passing from those positions to appear. This is shown in Figure 9↓, done within a 3x3x3 supercell.

Figure 9 Selecting the (111) lattice plane in the 3x3x3 supercell

Figure 10 The new unit cell marked out by the lattice plane “boundaries”

Trimming the new unit cell

Figure 11 Filtered new unit cell with the original cell’s lattice vectors

Transforming into the new unit cell

Figure 12 The reference points for the new axe1s

A Word of caution

Figure 13 The transformed structure with missing atoms

Figure 14 TRANSTRU input form

Final Touch

Figure 15 Prepared surface with vacuum

Figure 16 New unit cell with periodicity applied

Figure 17 Alternative surface with the Pt atoms on the top (Left: prior to vacuum introduction; Center: with vacuum; Right: tiled in periodicity)

Alternative Ways

Conclusion

Acknowledgements

Appendix: Walkthru for PtC (111) Surface

1. Open VESTA, (“File→New Structure”)
   1. “Unit Cell”:
      1. Space Group: 225 (F m -3 m)
      2. Lattice parameter a: 4.5 Å
   2. “Structure Parameters”:
      1. New→Symbol:C; Label:C; x,y,z=0.0
      2. New→Symbol:Pt; Label:Pt; x,y,z=0.5
   3. “Boundary..”
      1. x(min)=y(min)=z(min)=-2
      2. x(max)=y(max)=z(max)=5
   4. “Edit→Lattice Planes…”
      1. (hkl)=(111); d(Å)=9.09327
      2. (hkl)=(111); d(Å)=1.29904
      3. (hkl)=(-110); d(Å)=3.18198
      4. (hkl)=(-110); d(Å)=-3.18198
      5. (hkl)=(11-2); d(Å)=3.67423
      6. (hkl)=(11-2); d(Å)=-4.59279
      (compare with Figure 10↑)
   5. Select (“Select tool” – shortcut “s” key) and delete (shortcut “Del” key) all the atoms lying out of the area designated by the lattice planes.
   6. Select an “origin-atom” of a corner and the 3 “axe-atoms” in the edge of each direction of the unit cell, write down their fractional coordinates
      1. o(1,0,-1/2)
      2. a(0,1,-1/2)
      3. b(1/2,-1/2,1/2)
      4. c(2,1,1/2)
   7. Subtract the origin-atom coordinates from the rest and construct the transformation matrix accordingly
      P = ⎡⎢⎢⎢⎣ − 1 − 121 1 − 121 011⎤⎥⎥⎥⎦
   8. Transform the initial unit cell with this matrix via TRANSTRU (http://www.cryst.ehu.es/cryst/transtru.html) (Figure 14↑)
   9. Save the resulting (low symmetry) structure as CIF, open it in VESTA, export it to VASP, select “Cartesian coordinates”
   10. Open the VASP file in an editor, increase the c-lattice parameter to a big value (e.g. 25.000) to introduce vacuum. Save it and open it in VESTA.

References

[1] Mois Ilia Aroyo, Juan Manuel Perez-Mato, Cesar Capillas, Eli Kroumova, Svetoslav Ivantchev, Gotzon Madariaga, Asen Kirov, and Hans Wondratschek. Bilbao crystallographic server: I. databases and crystallographic computing programs. Zeitschrift für Kristallographie, 221(1):1527, 01 2006. doi: 10.1524/zkri.2006.221.1.15.

[2] Mois I. Aroyo, Asen Kirov, Cesar Capillas, J. M. Perez-Mato, and Hans Wondratschek. Bilbao crystallographic server. ii. representations of crystallographic point groups and space groups. Acta Crystallographica Section A, 62(2):115128, Mar 2006. http://dx.doi.org/10.1107/S0108767305040286

[3] M I Aroyo, J M Perez-Mato, D Orobengoa, E Tasci, G De La Flor, and A Kirov. Crystallography online: Bilbao crystallographic server. Bulg Chem Commun, 43(2):183197, 2011.

[4] S. R. Hall, F. H. Allen, and I. D. Brown. The crystallographic information le (cif): a new standard archive file for crystallography. Acta Crystallographica Section A, 47(6):655685, Nov 1991.

[5] G. Kresse and J. Furthmüller. Ecient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B, 54:1116911186, Oct 1996.

[6] Koichi Momma and Fujio Izumi. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. Journal of Applied Crystallography, 44(6):12721276, Dec 2011.

[7] A Zaoui and M Ferhat. Dynamical stability and high pressure phases of platinum carbide. Solid State Communications, 151(12):867869, 6 2011.

[8] Saulius Grazulis, Adriana Daskevic, Andrius Merkys, Daniel Chateigner, Luca Lutterotti, Miguel Quirós, Nadezhda R. Serebryanaya, Peter Moeck, Robert T. Downs, and Armel Le Bail. Crystallography open database (cod): an open-access collection of crystal structures and platform for world-wide collaboration. Nucleic Acids Research, 40(D1):D420D427, 2012.

[9] Alec Belsky, Mariette Hellenbrandt, Vicky Lynn Karen, and Peter Luksch. New developments in the inorganic crystal structure database (icsd): accessibility in support of materials research and design. Acta Crystallographica Section B, 58(3 Part 1):364369, Jun 2002.

[10] Springer, editor. The Landolt-Boernstein Database (http://www.springermaterials.com/navigation/). Springer Materials.

[11] Yongsheng Zhang. First-principles statistical mechanics approach to step decoration at solid surfaces . PhD thesis, Frei Universitat Berlin, 2008.

❬✶❪ ▼✐s ■❧✐❛ ❆r2✱ ❏❛ ▼❛❡❧ P❡r❡3✲▼❛t✱ ❈❡s❛r ❈❛✐❧❧❛s✱ ❊❧✐ ❑r❛✱ ❙❡ts❧❛ ■❛t❝❡✱ ●t3 ▼❛❞❛r✐❛❛✱

❆s❡ ❑✐r✱ ❛❞ ❍❛s ❲❞r❛ts❝❡❦ ✐❧❜❛ ❝r2st❛❧❧r❛✐❝ s❡r❡r✿ ■ ❞❛t❛❜❛s❡s ❛❞ ❝r2st❛❧❧r❛✐❝ ❝t✐

rr❛s ❩❡✐ts❝r✐❢t ❢ür ❑r✐st❛❧❧r❛✐❡✱ ✷✷✶✭✶✮✿✶✺✕✷✼✱ ✵✶ ✷✵✵✻ ❞✐✿ ✶✵✶✺✷✹3❦r✐✷✵✵✻✷✷✶✶✶✺

❬✷❪ ▼✐s ■ ❆r2✱ ❆s❡ ❑✐r✱ ❈❡s❛r ❈❛✐❧❧❛s✱ ❏ ▼ P❡r❡3✲▼❛t✱ ❛❞ ❍❛s ❲❞r❛ts❝❡❦ ✐❧❜❛ ❝r2st❛❧❧r❛✐❝

s❡r❡r ✐✐ r❡r❡s❡t❛t✐s ❢ ❝r2st❛❧❧r❛✐❝ ✐t rs ❛❞ s❛❝❡ rs ❆❝t❛ ❈r2st❛❧❧r❛✐❝❛ ❙❡❝t✐ ❆✱

✻✷✭✷✮✿✶✶✺✕✶✷✽✱ ▼❛r ✷✵✵✻

❬✸❪ ▼ ■ ❆r2✱ ❏ ▼ P❡r❡3✲▼❛t✱ ❉ ❖r❜❡❛✱ ❊ ❚❛s❝✐✱ ● ❉❡ ▲❛ ❋❧r✱ ❛❞ ❆ ❑✐r ❈r2st❛❧❧r❛2 ❧✐❡✿ ✐❧❜❛

❝r2st❛❧❧r❛✐❝ s❡r❡r ❧ ❈❡ ❈✱ ✹✸✭✷✮✿✶✽✸✕✶✾✼✱ ✷✵✶✶

❬✹❪ ❙ ❘ ❍❛❧❧✱ ❋ ❍ ❆❧❧❡✱ ❛❞ ■ ❉ r✇ ❚❡ ❝r2st❛❧❧r❛✐❝ ✐❢r❛t✐ ✜❧❡ ✭❝✐❢✮✿ ❛ ❡✇ st❛❞❛r❞ ❛r❝✐❡ ✜❧❡

❢r ❝r2st❛❧❧r❛2 ❆❝t❛ ❈r2st❛❧❧r❛✐❝❛ ❙❡❝t✐ ❆✱ ✹✼✭✻✮✿✻✺✺✕✻✽✺✱ ◆ ✶✾✾✶❍❖❲ ❚❖ P❘❊P❆❘❊ ❆◆ ■◆P❯❚ ❋■▲❊ ❋❖❘ ❙❯❘❋❆❈❊ ❈❆▲❈❯▲❆❚■❖◆❙

✶✻

❬✺❪ ● ❑r❡ss❡ ❛❞ ❏ ❋rtü❧❧❡r ❊✣❝✐❡t ✐t❡r❛t✐❡ s❝❡❡s ❢r ❛❜ ✐✐t✐ tt❛❧✲❡❡r2 ❝❛❧❝❧❛t✐s s✐ ❛ ❧❛❡✲✇❛❡

❜❛s✐s s❡t P2s ❘❡ ✱ ✺✹✿✶✶✶✻✾✕✶✶✶✽✻✱ ❖❝t ✶✾✾✻

❬✻❪ ❑✐❝✐ ▼❛ ❛❞ ❋❥✐ ■3✐ ❱❊❙❚❆✸ ❢r tr❡❡✲❞✐❡s✐❛❧ ✐s❛❧✐3❛t✐ ❢ ❝r2st❛❧✱ ❧❡tr✐❝ ❛❞ r❧✲

2 ❞❛t❛ ❏r❛❧ ❢ ❆❧✐❡❞ ❈r2st❛❧❧r❛2✱ ✹✹✭✻✮✿✶✷✼✷✕✶✷✼✻✱ ❉❡❝ ✷✵✶✶

❬✼❪ ❆ ❩❛✐ ❛❞ ▼ ❋❡r❛t ❉2❛✐❝❛❧ st❛❜✐❧✐t2 ❛❞ ✐ r❡ssr❡ ❛s❡s ❢ ❧❛t✐ ❝❛r❜✐❞❡ ❙❧✐❞ ❙t❛t❡ ❈✐✲

❝❛t✐s✱ ✶✺✶✭✶✷✮✿✽✻✼✕✽✻✾✱ ✻ ✷✵✶✶

❬✽❪ ❙❛❧✐s ●r❛➸❧✐s✱ ❆❞r✐❛❛ ❉❛➨❦❡✐✱ ❆❞r✐s ▼❡r❦2s✱ ❉❛✐❡❧ ❈❛t❡✐❡r✱ ▲❝❛ ▲tt❡rtt✐✱ ▼✐❡❧ ◗✐rós✱

◆❛❞❡3❞❛ ❘ ❙❡r❡❜r2❛❛2❛✱ P❡t❡r ▼❡❝❦✱ ❘❜❡rt ❚ ❉✇s✱ ❛❞ ❆r❡❧ ▲❡ ❛✐❧ ❈r2st❛❧❧r❛2 ❡ ❞❛t❛✲

❜❛s❡ ✭❝❞✮✿ ❛ ❡✲❛❝❝❡ss ❝❧❧❡❝t✐ ❢ ❝r2st❛❧ str❝tr❡s ❛❞ ❧❛t❢r ❢r ✇r❧❞✲✇✐❞❡ ❝❧❧❛❜r❛t✐ ◆❝❧❡✐❝

❆❝✐❞s ❘❡s❡❛r❝✱ ✹✵✭❉✶✮✿❉✹✷✵✕❉✹✷✼✱ ✷✵✶✷

❬✾❪ ❆❧❡❝ ❡❧s❦2✱ ▼❛r✐❡tt❡ ❍❡❧❧❡❜r❛❞t✱ ❱✐❝❦2 ▲2 ❑❛r❡✱ ❛❞ P❡t❡r ▲❦s❝ ◆❡✇ ❞❡❡❧❡ts ✐ t❡ ✐r❛✐❝

❝r2st❛❧ str❝tr❡ ❞❛t❛❜❛s❡ ✭✐❝s❞✮✿ ❛❝❝❡ss✐❜✐❧✐t2 ✐ srt ❢ ❛t❡r✐❛❧s r❡s❡❛r❝ ❛❞ ❞❡s✐ ❆❝t❛ ❈r2st❛❧❧r❛✐❝❛

❙❡❝t✐ ✱ ✺✽✭✸ P❛rt ✶✮✿✸✻✹✕✸✻✾✱ ❏ ✷✵✵✷

❬✶✵❪ ❙r✐❡r✱ ❡❞✐tr ❚❡ ▲❛❞❧t✲❡rst❡✐ ❉❛t❛❜❛s❡ ✭tt✿✇✇✇sr✐❡r❛t❡r✐❛❧s❝❛✐❛t✐✮ ❙r✐❡r

▼❛t❡r✐❛❧s

❬✶✶❪ ❨s❡ ❩❛ ❋✐rst✲r✐❝✐❧❡s st❛t✐st✐❝❛❧ ❡❝❛✐❝s ❛r❛❝ t st❡ ❞❡❝r❛t✐ ❛t s❧✐❞ sr❢❛❝❡s P❉ t❡s✐s✱

❋r❡✐ ❯✐❡rs✐t❛t ❡r❧✐✱ ✷✵✵✽

Suppose that we want to calculate the energies corresponding to various values of the lattice parameter a, ranging from 3.2 to 4, incrementing by 2: 3.2, 3.4, … , 4.0

Incrementation:
In bash, you can define a for loop to handle this:

for (( i = 1; i <= 9; i++ ))
do
  echo $i;
done

1
2
3
4
5
6
7
8
9

for (( i = 1; i <= 5; i++ ))
do
    val=$(awk "BEGIN {print 3+$i/5 }")
    echo $val;
done

3.2
3.4
3.6
3.8
4

octave:1> [3.2:0.2:4]
ans =

   3.2000   3.4000   3.6000   3.8000   4.0000

sururi@husniya:~/shared/qe/diamond_tutorial$ echo "[3.2:0.2:4]"|octave -q
ans =

   3.2000   3.4000   3.6000   3.8000   4.0000

octave:4> [3:0.1:4]
ans =

 Columns 1 through 9:

   3.0000   3.1000   3.2000   3.3000   3.4000   3.5000   3.6000   3.7000   3.8000

 Columns 10 and 11:

   3.9000   4.0000

sururi@vala:/tmp$ echo "[3.2:0.1:4]"|octave -q|sed -n "2,1000p"| grep -v "Column"

   3.2000   3.3000   3.4000   3.5000   3.6000   3.7000   3.8000   3.9000

   4.0000

sururi@vala:/tmp$ echo "[3.2:0.1:4]"|octave -q|sed -n "2,1000p"| grep -v "Column"|grep -v -e "^$"
   3.2000   3.3000   3.4000   3.5000   3.6000   3.7000   3.8000   3.9000
   4.0000

cat diamond.scf.in | sed "s:#REPLACE#:$a:g">tmp_$a.in

Parsing the Energy:
We submit this job via mpirun (in our case via mpi-pw.x of the previous post) and parse the output for the final energy:

grep -e ! tmp_$a.out|tail -n1|awk '{print $(NF-1)}'

The Script:
Putting all these together, here is the script in its full glory:

#!/bin/bash

# Takes in the template Quantum Espresso input file, runs it iteratively
# changing the #REPLACE# parameters iteratively with the range elements,
# grepping the energy and putting it along with the current value of the
# #REPLACE# parameter.
#
# Emre S. Tasci, 03/2012

# Example: Suppose that the file "diamond.scf.in" contains the following line:
#   A = #REPLACE#
# Calling "qe-opti.sh diamond.scf.in a 3.2 0.2 4.2" would cause that line to be modified as
#   A = 3.2
# and submitted, followed by
#   A = 3.4, 3.6, .., 4.2 in the sequential runs.
#
# at the end, having constructed the E_vs_a.txt file.

if [ $# -ne 5 ]
then
    echo "Usage: qe-opti.sh infile_template.in variable_name initial_value increment     final_value"
    echo "Example: qe-opti.sh diamond.scf.in a 3.2 0.2 4.2"
    echo -e "\nNeeds octave to be accessible via \"octave\" command for the range to     work"
    echo -e "\nEmre S. Tasci, 03/2012"
    exit
fi

logext=".txt"
logfilename=E_vs_$2$logext
rm -f $logfilename

range=$(echo "[$3:$4:$5]"|octave -q|sed -n "2,1000p"|grep -vr "^$"|grep -v "Column")
for a in $range
do
    # echo $a
    cat $1|sed "s:#REPLACE#:$a:g">tmp_$a.in
    /home/sururi/bin/mpi-pw.x 4 tmp_$a.in
energ=$(grep -e ! tmp_$a.out|tail -n1|awk '{print $(NF-1)}')
    # echo $energ
    echo -e $a "    " $energ >> $logfilename
    echo -e $a "    " $energ
done

echo "Results stored in: "$logfilename

sururi@bebop:~/shared/qe/diamond_tutorial$ grep "REPLACE" -B2 -A2 diamond.scf.in
    # a,b,c in Angstrom =======
    #A = 3.56712
    A = #REPLACE#
    B = 2.49
    C = 2.49

sururi@bebop:~/shared/qe/diamond_tutorial$ qe-opti.sh diamond.scf.in a 3.2 0.2 4.2
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.2000.in > tmp_3.2000.out

Start : Tue Mar 13 23:21:07 CET 2012
Finish: Tue Mar 13 23:21:08 CET 2012
Elapsed time:0:00:01

3.2000      -22.55859979
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.4000.in > tmp_3.4000.out

Start : Tue Mar 13 23:21:08 CET 2012
Finish: Tue Mar 13 23:21:09 CET 2012
Elapsed time:0:00:01

3.4000      -22.67717405
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.6000.in > tmp_3.6000.out

Start : Tue Mar 13 23:21:09 CET 2012
Finish: Tue Mar 13 23:21:10 CET 2012
Elapsed time:0:00:01

3.6000      -22.69824702
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_3.8000.in > tmp_3.8000.out

Start : Tue Mar 13 23:21:10 CET 2012
Finish: Tue Mar 13 23:21:11 CET 2012
Elapsed time:0:00:01

3.8000      -22.65805473
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_4.0000.in > tmp_4.0000.out

Start : Tue Mar 13 23:21:11 CET 2012
Finish: Tue Mar 13 23:21:13 CET 2012
Elapsed time:0:00:02

4.0000      -22.57776354
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < tmp_4.2000.in > tmp_4.2000.out

Start : Tue Mar 13 23:21:13 CET 2012
Finish: Tue Mar 13 23:21:14 CET 2012
Elapsed time:0:00:01

4.2000      -22.47538891
Results stored in: E_vs_a.txt

sururi@bebop:~/shared/qe/diamond_tutorial$ cat E_vs_a.txt
3.2000      -22.55859979
3.4000      -22.67717405
3.6000      -22.69824702
3.8000      -22.65805473
4.0000      -22.57776354
4.2000      -22.47538891

sururi@husniya:~/shared/qe/diamond_tutorial$ plot E_vs_a.txt

Other than the ‘plot’, I’ve also written another one, fit for data sets just like these where you have a parabol-like data and you’d like to try fitting it and here is what it does (automatically):

sururi@husniya:~/shared/qe/diamond_tutorial$ plotfitparabol.sh E_vs_a.txt
func    f(x)=0.674197*x*x-4.881275*x-13.855234
min, f(min)    3.620066 -22.690502

Here are the source codes of the ‘plot’ and ‘plotfitparabol.sh’ scripts. Doei!

sururi@mois:/vala/bin$ cat plot
#!/bin/bash 
gnuplot -persist <<END 
set term postscript enhanced color
set output "$1.ps"
plot "$1" $2 $3

set term wxt
replot
END
sururi@mois:/vala/bin$ cat plotfitparabol.sh 
#!/bin/bash 

# Usage : plotfitparabol.sh <datafile>

# Reads the data points from the given file
# then uses Octave's "polyfit" function to calculate the coefficients 
# for the quadratic fit and also the minimum (maximum) value.
# Afterwards, prepares a GNUPlot script to display:
#    the data set (scatter)
#    fit function
#    the extremum point's value
#
# Written with Quantum Espresso optimizing procedures in mind.

 
if [ $# -ne 1 ]
then
    echo "Usage: plotfitparabol.sh <datafile>"
    echo "Example: plotfitparabol.sh E_vs_A.txt"
    echo -e "\nNeeds octave to be accessible via \"octave\" command to work"
    echo -e "\nEmre S. Tasci <emre.tasci@ehu.es>, 03/2012"
    exit
fi

octave_script="data=load(\"$1\");m=polyfit(data(:,1),data(:,2),2);min=-m(2)*0.5/m(1);fmin=min*min*m(1)+min*m(2)+m(3);printf(\"f(x)=%f*x*x%+f*x%+f|%f\t%f\",m(1),m(2),m(3),min,fmin)"
oct=$(echo $octave_script|octave -q)
#echo $oct

func=${oct/|*/ }
minfmin=${oct/*|/ }
echo -e "func \t" $func
echo -e "min, f(min) \t" $minfmin

echo '$func
set term postscript enhanced color
set label "$1" at graph 0.03, graph 0.94
set output "$1.ps"
plot "$1" title "data", f(x) title "$func", "<echo '$minfmin'" title "$minfmin"

set term wxt
replot'

gnuplot -persist <<END 
$func
set term postscript enhanced color
set label "$1" at graph 0.03, graph 0.94
set output "$1.ps"
plot "$1" title "data", f(x) title "$func", "<echo '$minfmin'" title "$minfmin"

set term wxt
replot
END

As a change, I selected and installed the Quantum Espresso (QE) package instead of VASP that I was using in the past, main reason QE being open source.

As a newbie to QE and long-gone returner to DFT, I started to practice by mostly ye olde trial/error approach.

The first thing that surprised me the lack of a commenting option in regard to QE input files! (Actually, I’m almost sure that it surely is supported but I have yet to find!). Here’s my sample input file:

&control
    title='Diamond Relaxing',

    # Calculation Type ==========
    calculation='scf'
    #calculation='relax',
    #calculation='vc-relax'
    #calculation='cp'

    # Restart ===================
    restart_mode='from_scratch',
    #restart_mode='restart' 

    # Read/Write units ==========
    #ndr = 51,
    #ndw = 51,

    #nstep = 100,
    #nstep = 1000,
    nstep = 2000,

    # max_seconds = 600 # jobs stops after max_seconds CPU time.

    iprint = 10, # band energies are written every iprint iterations
    #tstress = .TRUE., # calculate stress. It is set to .TRUE. automatically if
                       # calculation='vc-md' or 'vc-relax'
    #tprnfor = .TRUE., # print forces. Set to .TRUE. if calculation='relax','md','vc-md'
    #dt = 5.0d0, # time step for molecular dynamics, in Rydberg atomic units
                 # (1 a.u.=4.8378 * 10^-17 s : beware, the CP code use
                 #  Hartree atomic units, half that much!!!)

    prefix='diamond_rel', # prepended to input/output filenames
    pseudo_dir = '/usr/share/espresso-4.3.2/pseudo/',
    outdir='/home/sururi/shared/tmp/'
    etot_conv_thr = 1.d-6 # convergence threshold on total energy (a.u) for ionic minimization
    forc_conv_thr = 1.d-3 # convergence threshold on forces (a.u) for ionic minimization
 /
 &system
##  ibrav
##    0          "free", see above                 not used
##    1          cubic P (sc)                      not used
##    2          cubic F (fcc)                     not used
##    3          cubic I (bcc)                     not used
##    4          Hexagonal and Trigonal P        celldm(3)=c/a
##    5          Trigonal R, 3fold axis c        celldm(4)=cos(alpha)
##   -5          Trigonal R, 3fold axis <111>    celldm(4)=cos(alpha)
##    6          Tetragonal P (st)               celldm(3)=c/a
##    7          Tetragonal I (bct)              celldm(3)=c/a
##    8          Orthorhombic P                  celldm(2)=b/a,celldm(3)=c/a
##    9          Orthorhombic base-centered(bco) celldm(2)=b/a,celldm(3)=c/a
##   10          Orthorhombic face-centered      celldm(2)=b/a,celldm(3)=c/a
##   11          Orthorhombic body-centered      celldm(2)=b/a,celldm(3)=c/a
##   12          Monoclinic P, unique axis c     celldm(2)=b/a,celldm(3)=c/a,
##                                               celldm(4)=cos(ab)
##  -12          Monoclinic P, unique axis b     celldm(2)=b/a,celldm(3)=c/a,
##                                               celldm(5)=cos(ac)
##   13          Monoclinic base-centered        celldm(2)=b/a,celldm(3)=c/a,
##                                               celldm(4)=cos(ab)
##   14          Triclinic                       celldm(2)= b/a,
##                                               celldm(3)= c/a,
##                                               celldm(4)= cos(bc),
##                                               celldm(5)= cos(ac),
##                                               celldm(6)= cos(ab
    ibrav= 2,

# Cell dimensions ========================
    # a,b,c in Angstrom =======
    A = 3.56712
    B = 2.49
    C = 2.49
    cosAB=0.5
    cosAC=0.5
    cosBC=0.5

    ## OR ##
    ## celldm(1:6) in Bohr =====
    ## 1 Bohr =  0.529177249 A
    ## 1 A = 1.889725989 Bohr
    # celldm(1) = 6.7409 # alat
    # celldm(2) = 4.7054 # b/a
    # celldm(3) = 4.7054 # c/a
    # celldm(4) = 0.5 # cos(ab)
    # celldm(5) = 0.5 # cos(ac)
    # celldm(6) = 0.5 # cos(ab)

   nat=2,
   ntyp= 1,

   ecutwfc = 40, # kinetic energy cutoff (Ry) for wavefunctions
   # ecutrho = 160 # kinetic energy cutoff (Ry) for charge density and potential.
                   # For norm-conserving pseudopotential you should stick to the
                   # default value
   # nbnd = 6 # number of electronic states (bands) to be calculated.
              # Default:
              # for an insulator, nbnd = number of valence bands (nbnd = # of electrons /2);
              # for a metal, 20% more (minimum 4 more)
 /
 &electrons
   emass = 400.d0
   emass_cutoff = 2.5d0

   #electron_dynamics = 'sd' # CP specific
   #electron_dynamics = 'bfgs' # CP specific

  conv_thr = 1D-6 # Convergence threshold for selfconsistency

 /
 &ions
    #ion_dynamics = 'none'
    #ion_dynamics = 'bfgs'
    #ion_dynamics = 'sd'

    # === damping === 0 ==
    # ion_dynamics = 'damp'
    # ion_damping = 0.2
    # ion_velocities = 'zero' # CP specific -- initial ionic velocities
    # === damping === 1 ==

    #ion_nstepe = 10 # CP specific -- number of electronic steps per ionic step. (Def: 1)
 /
 &cell
    cell_dynamics = 'none',
    #cell_dynamics = 'bfgs', # ion_dynamics must be 'bfgs' too
    #cell_dofree = 'xyz'
    #cell_factor = 3
    press = 0.0d,
    # press_conv_thr = 0.05 # Convergence threshold on the pressure for variable cell relaxation
 /
ATOMIC_SPECIES
 C  12.0107 C.pz-vbc.UPF
ATOMIC_POSITIONS
 C 0.0 0.0 0.0
 C 0.25 0.25 0.25
K_POINTS automatic
   4 4 4 1 1 1

And then the problem of (yet undiscovered) comments in the QE input files hits us. But with bash, it’s pretty easy. So, if you save the above contents to a file called, say, “diamond.scf.in” and filter it via grep and sed, you’re done:

sururi@husniya:~/shared/qe/diamond_tutorial$ cat diamond.scf.in | grep -v -E "^[\t ]*#" | sed "s:\([^#]*\)#.*:\1:"
&control
    title='Diamond Relaxing',

    calculation='scf'

    restart_mode='from_scratch',

    nstep = 2000,

    iprint = 10, 

    prefix='diamond_rel',
    pseudo_dir = '/usr/share/espresso-4.3.2/pseudo/',
    outdir='/home/sururi/shared/tmp/'
    etot_conv_thr = 1.d-6
    forc_conv_thr = 1.d-3
 /
 &system
    ibrav= 2,

    A = 3.56712
    B = 2.49
    C = 2.49
    cosAB=0.5
    cosAC=0.5
    cosBC=0.5

   nat=2,
   ntyp= 1,

   ecutwfc = 40,
 /
 &electrons
   emass = 400.d0
   emass_cutoff = 2.5d0

  conv_thr = 1D-6 

 /
 &ions

 /
 &cell
    cell_dynamics = 'none',
    press = 0.0d,
 /
ATOMIC_SPECIES
 C  12.0107 C.pz-vbc.UPF
ATOMIC_POSITIONS
 C 0.0 0.0 0.0
 C 0.25 0.25 0.25
K_POINTS automatic
   4 4 4 1 1 1

mpirun -n #_of_processes /usr/share/espresso-4.3.2/bin/pw.x < input_file > output_file

which I guess is more or less similar to what you are doing to run your jobs.

So, combining this “comment filtering” and job submission, I’ve written the following bash script (“mpi-pw.x“) to automatize (and keep log of) the job:

sururi@cowboy:~$ cat /home/sururi/bin/mpi-pw.x 
#!/bin/bash

# Script to submit Quantum Espresso jobs via mpi
# EST, Fri Mar  9 17:02:10 CET 2012

function timer()
{
    if [[ $# -eq 0 ]]; then
        echo $(date '+%s')
    else
        local  stime=$1
        etime=$(date '+%s')

        if [[ -z "$stime" ]]; then stime=$etime; fi

        dt=$((etime - stime))
        ds=$((dt % 60))
        dm=$(((dt / 60) % 60))
        dh=$((dt / 3600))
        printf '%d:%02d:%02d' $dh $dm $ds
    fi
}

if [ $# -ne 2 ]
then
    echo "Usage: mpi-pw.x num_of_processes infile.in"
    exit
fi

job=$2
job_wo_ext=${job%.*}

outext=".out"
logext=".log"
echo -e "Running: mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < $2 > $job_wo_ext$outext\n"
echo -e "Running: mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < $2 > $job_wo_ext$outext\n">$job_wo_ext$logext
echo "Start : "$(date)
echo "Start : "$(date)>>$job_wo_ext$logext

# Take out the comments from the input file
#grep -v -E "^[\t ]*#" $2 > in_wo_comments.in.tmp
grep -v -E "^[\t ]*#" $2 | sed "s:\([^#]*\)#.*:\1:" > in_wo_comments.in.tmp

t=$(timer)
mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < in_wo_comments.in.tmp > $job_wo_ext$outext
echo "Finish: "$(date)
echo "Finish: "$(date)>>$job_wo_ext$logext
printf 'Elapsed time:%s\n\n' $(timer $t)
printf 'Elapsed time:%s\n\n' $(timer $t) >>$job_wo_ext$logext

rm -f in_wo_comments.in.tmp

(I’ve picked the “timer” function from Mitch Frazier’s entry in LinuxJournal)

What it does is, really simple: it strips out the extension of the input file ($2 – the second argument in calling); filters out the comments, constructing the temporary input file “in_wo_comments.in.tmp”; starts the timer; submits the job using the number of processes passed in the calling ($1 – first argument) and directs the output to a file whose name is the extension stripped input filename + “.out”; when the job is (this way or that) done, stops the timer and calculates the elapsed time via the ‘timer’ function. During all this time, it prints these information both on screen and to a file called extension stripped input filename + “.log”.

I have a very similar another script called “mpi-cp.x” that I’ve forked for Carr-Parinelli jobs (just search & replace all the occurrences of “mpi-pw.x” with “mpi-cp.x” and you’re done 8).

And here is the script in action:

sururi@cowboy:~/shared/qe/diamond_tutorial$ ll
total 16
drwxr-xr-x 2 sururi sururi 4096 2012-03-13 15:59 ./
drwxr-xr-x 7 sururi sururi 4096 2012-03-13 15:10 ../
-rw-r--r-- 1 sururi sururi 4954 2012-03-13 15:59 diamond.scf.in

sururi@cowboy:~/shared/qe/diamond_tutorial$ mpi-pw.x 4 diamond.scf.in 
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < diamond.scf.in > diamond.scf.out

Start : Tue Mar 13 16:18:20 CET 2012
Finish: Tue Mar 13 16:18:21 CET 2012
Elapsed time:0:00:01

sururi@cowboy:~/shared/qe/diamond_tutorial$ grep -e ! diamond.scf.out 
!    total energy              =     -22.70063050 Ry

sururi@cowboy:~/shared/qe/diamond_tutorial$ cat diamond.scf.log 
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < diamond.scf.in > diamond.scf.out

Start : Tue Mar 13 16:18:20 CET 2012
Finish: Tue Mar 13 16:18:21 CET 2012
Elapsed time:0:00:01

Now that we have the mpi-pw.x and mpi-cp.x, our next entry in the series will be the automated search for optimized values.

Recently, while searching recent publications on liquid alloys, I’ve found a very good and resourceful article on liquid Na-K alloys by Dr. J.-F. Wax (“Molecular dynamics study of the static structure of liquid Na-K alloys”, Physica B 403 (2008) 4241-48 | doi:10.1016/j.physb.2008.09.014). In this paper, among other properties, he had also presented the partial pair distribution functions. Partial-pair distribution functions is the averaged probability that you’ll encounter an atom in r distance starting from another atom. In solid phase, due to the symmetry of the crystal structure, you have discrete values instead of a continuous distribution and everything is crystal clear but I guess that’s the price we pay dealing with liquids. 8)

Partial-pair distribution functions are very useful in the sense that they let you derive the bond-length between constituent atom species, corresponding to the first minimums of the plots. Hence, one can see the neighborhood ranges by just looking at the plots.

I’ve sent a mail to Dr. Wax explaining my research topics and my interest and asking for his research data to which he very kindly and generously replied by providing his data.

He had sent me a big file with numerous configurations (coordinates of the atoms) following each other. The first -and the simplest- task was to split the file into seperate configuration files via this bash script:

#!/bin/bash
# Script Name: split_POS4.sh
# Splits the NaK conf files compilation POS4_DAT into sepearate
# conf files.
# Emre S. Tasci <e.tasci@tudelft.nl>
#                                              01/09/09

linestart=2
p="p"
for (( i = 1; i <= 1000; i++ ))
do
    echo $i;
    lineend=$(expr $linestart + 2047)
    confname=$(printf "%04d" $i)
    echo -e "3\n2048" > $confname
    sed -n "$linestart,$(echo $lineend$p)" POS4_DAT >> $confname
    linestart=$(expr $lineend + 1)
    cat $confname | qhull i TO $confname.hull
done

As you can check in the line before the “done” closer, I’ve -proudly- used Qhull software to calculate the convex hull of each configuration. A convex hull is the -hopefully minimum- shape that covers all your system. So, for each configuration I now had 2 files: “xxxx” (where xxxx is the id/number of the configuration set) storing the coordinates (preceded by 3 and 2048, corresponding to dimension and number of atoms information for the qhull to parse & process) and “xxxx.hull” file, generated by the qhull, containing the vertices list of each facet (preceded by the total number of facets).

A facet is (in 3D) the plane formed by the vertice points. Imagine a cube, where an atom lies at each corner, 3 in the inside. So, we can say that our system consists of 11 atoms and the minimal shape that has edges selected from system’s atoms and covering the whole is a cube, with the edges being the 8 atoms at the corners. Now, add an atom floating over the face at the top. Now the convex hull that wraps our new system would no longer be a 6-faced, 8-edged cube but a 9-faced, 9-edged polygon. I need to know the convex hull geometry in order to correctly calculate the coordination number distributions (details will follow shortly).

The aforementioned two files look like:

sururi@husniya hull_calcs $ head 0001
3
2048
 0.95278770E+00 0.13334565E+02 0.13376155E+02
 0.13025256E+02 0.87618381E+00 0.20168993E+01
 0.12745038E+01 0.14068998E-01 0.22887323E+01
 0.13066590E+02 0.20788591E+01 0.58119183E+01
 0.10468218E+00 0.58523640E-01 0.64288678E+01
 0.12839412E+02 0.13117012E+02 0.79093881E+01
 0.11918105E+01 0.12163854E+02 0.10270868E+02
 0.12673985E+02 0.11642538E+02 0.10514597E+02
sururi@husniya hull_calcs $ head 0001.hull
180
568 1127 1776
1992 104 1551
956 449 1026
632 391 1026
391 632 1543
391 956 1026
966 632 1026
966 15 1039
632 966 1039

To calculate the coordination number distribution, you need to pick each atom and then simply count the other atoms within the cut-off distance depending on your atom’s and the neighboring atom’s species. After you do the counting for every atom, you average the results and voilà! there is your CN distribution. But here’s the catch: In order to be able to count properly, you must make sure that, your reference atom’s cutoff radii stay within your system — otherwise you’ll be undercounting. Suppose your reference atom is one of the atoms at the corners of the cube: It will (generally/approximately) lack 7/8 of the neighbouring atoms it would normally have if it was to be the center atom of the specimen/sample. Remember that we are studying a continuous liquid and trying to model the system via ensembles of samples. So we should omit these atoms which has neighborhoods outside of our systems. Since neighborhoods are defined through bond-length, we should only include the atoms with a distance to the outside of at least cut-off radii. The “outside” is defined as the region outside the convex hull and it’s designated by the facets (planes).

From the xxxx.hull file, we have the ids of the vertice atoms and their coordinates are stored in the xxxx file. To define the plane equation

we will use the three vertice points, say p1, p2 and p3. From these 3 points, we calculate the normal n as:

then, the plane equation (hence the coefficients a,b,c,d) can be derived by comparing the following equation:

Here, (x0,y0,z0) are just the coordinates of any of the vertice points. And finally, the distance D between a point p1(x1,y1,z1) and a plane (a,b,c,d) is given by:

Since vector operations are involved, I thought it’s better to do this job in Octave and wrote the following script that calculates the coefficients of the facets and then sets down to find out the points whose distances to all the facets are greater than the given cut-off radius:

sururi@husniya trunk $ cat octave_find_atoms_within.m
DummyA = 3;
%conf_index=292;

% Call from the command line like :
%   octave -q --eval "Rcut=1.65;conf_index=293;path='$(pwd)';" octave_find_atoms_within.m
% hence specifying the conf_index at the command line.

% Emre S. Tasci <e.tasci@tudelft.nl>                    *
%
% From the positions and vertices file, calculates the plane
% equations (stored in "facet_coefficients") and then
% filters the atoms with respect to their distances to these
% facets. Writes the output to
% "hull_calcs/xxxx_insiders.txt" file
%                                              03/09/09 */

conf_index_name = num2str(sprintf("%04d",conf_index));;
%clock
clear facet_coefficients;
facet_coefficients = [];
global p v facet_coefficients
fp = fopen(strcat(path,"/hull_calcs/",conf_index_name));
k = fscanf(fp,"%d",2); % 1st is dim (=3), 2nd is number of atoms (=2048)
p = fscanf(fp,"%f",[k(1),k(2)]);
p = p';
fclose(fp);
%p = load("/tmp/del/delco");
fp = fopen(strcat(path,"/hull_calcs/",conf_index_name,".hull"));
k = fscanf(fp,"%d",1);% number of facets
v = fscanf(fp,"%d",[3,k]);
v = v';
fclose(fp);
%v = load("/tmp/del/delver");

function [a,b,c,d,e] = getPlaneCoefficients(facet_index)
    global p v
    % Get the 3 vertices
    p1 = p(v(facet_index,1)+1,:); % The +1s come from the fact 
    p2 = p(v(facet_index,2)+1,:); % that qhull enumeration begins
    p3 = p(v(facet_index,3)+1,:); % from 0 while Octave's from 1

    %printf("%d: %f %f %f %f %f %f %f %f %f\n",facet_index,p1,p2,p3);

    % Calculate the normal
    n = cross((p2-p1),(p3-p1));

    % Calculate the coefficients of the plane :
    % ax + by + cz +d = 0
    a = n(1);
    b = n(2);
    c = n(3);
    d = -(dot(n,p1));
    e = norm(n);
    if(e==0)
        printf("%f\n",p1);
        printf("%f\n",p2);
        printf("%f\n",p3);
        printf("%d\n",facet_index);
    endif
endfunction

function dist = distanceToPlane(point_index,facet_index)
    global p facet_coefficients
    k = facet_coefficients(facet_index,:);
    n = [k(1) k(2) k(3)];
    dist = abs(dot(n,p(point_index,:)) + k(4)) / k(5);
endfunction

for i = [1:size(v)(1)]
    [a,b,c,d,e] = getPlaneCoefficients(i);
    facet_coefficients = [ facet_coefficients ; [a,b,c,d,e]];
endfor

%Rcut = 1.65; % Defined from command line calling
% Find the points that are not closer than Rcut
inside_atoms = [];
for point = [1:size(p)(1)]
%for point = [1:100]
    inside=true;
    for facet = [1:size(v)(1)]
    %for facet = [1:10]
        %printf("%d - %d : %f\n",point,facet,dist);
        if (distanceToPlane(point,facet)<Rcut)
            inside = false;
            break;
        endif
    endfor
    if(inside)
        inside_atoms = [inside_atoms point-1];
    endif
endfor

fp = fopen(strcat(path,"/hull_calcs/",conf_index_name,"_insiders.txt"),"w");
fprintf(fp,"%d\n",inside_atoms);
fclose(fp);
%clock

<?PHP
/* Emre S. Tasci <e.tasci@tudelft.nl>                    *
 *
 * parse_configuration_data.php
 *
 * Written in order to parse the configuation files
 * obtained from J.F. Wax in order to calculate the 
 * Coordination Number distribution. Each configuration 
 * file contain 2048 atoms' coordinates starting at the 
 * 3rd line. There is also a convex hull file generated 
 * using qhull (http://www.qhull.com) that contains the 
 * indexes (starting from 0) of the atoms that form the 
 * vertices of the convex hull. Finally there is the 
 * IP.dat file that identifies each atom (-1 for K; 1 for 
 * Na -- the second column is the relative masses).
 *
 *                                              02/09/09 */

/*
# if present, remove the previous run's results file
if(file_exists("results.txt"))
    unlink("results.txt");
 */

error_reporting(E_ALL ^ E_NOTICE);

if(!file_exists("results.txt"))
{
    $fp = fopen("results.txt","w");
    fwrite($fp,"CONF#"."\tNa".implode("\tNa",range(1,20))."\tK".implode("\tK",range(1,20))."\n");
    fclose($fp);
}

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

if($_GET["configuration"])
    calculate($_GET["configuration"]);
else
    exit("Please specify a configuration number [php parse_configuration_data.php configuration=xxx]\n");

function calculate($conf_index)
{
    # Define the atoms array
    $arr_atoms = Array();

    # Get the information from the files
    parse_types($arr_atoms);
    $refs = get_vertices($conf_index);
    parse_coordinates($arr_atoms,$conf_index);

    # Define the Rcut-off values (obtained from partial parid distribution minimums)
    $Rcut_arr = Array();
    $Rscale = 3.828; # To convert reduced distances to Angstrom (if needed)
    $Rscale = 1; # To convert reduced distances to Angstrom
    $Rcut_arr["Na"]["Na"] = 1.38 * $Rscale;
    $Rcut_arr["Na"]["K"]  = 1.65 * $Rscale;
    $Rcut_arr["K"]["Na"]  = 1.65 * $Rscale;
    $Rcut_arr["K"]["K"]   = 1.52 * $Rscale;
    $Rcut = $Rcut_arr["Na"]["K"]; # Taking the maximum one as the Rcut for safety

    /*
    # We have everything we need, so proceed to check 
    # if an atom is at safe distance wrt to the vertices

    # Define all the ids
    $arr_test_ids = range(0,2047);
    # Subtract the ref atoms' ids
    $arr_test_ids = array_diff($arr_test_ids, $refs);
    sort($arr_test_ids);
    //print_r($arr_test_ids);

    $arr_passed_atom_ids = Array();
    # Check each atom against refs
    for($id=0, $size=sizeof($arr_test_ids); $id<$size; $id++)
    {
        //echo $id."\n";
        $arr_atoms[$arr_test_ids[$id]]["in"] = TRUE;
        foreach ($refs as $ref)
        {
            $distance = distance($arr_atoms[$arr_test_ids[$id]],$arr_atoms[$ref]);
            //echo "\t".$ref."\t".$distance."\n";
            if($distance < $Rcut)
            {
                $arr_atoms[$arr_test_ids[$id]]["in"] = FALSE;
                break;
            }
        }
        if($arr_atoms[$arr_test_ids[$id]]["in"] == TRUE)
            $arr_passed_atom_ids[] = $arr_test_ids[$id];
    }
     */

    # Run the octave script file to calculate the inside atoms
    if(!file_exists("hull_calcs/".sprintf("%04d",$conf_index)."_insiders.txt"))
    exec("octave -q --eval \"Rcut=".$Rcut.";conf_index=".$conf_index.";path='$(pwd)';\" octave_find_atoms_within.m");

    # Read the file generated by Octave
    $arr_passed_atom_ids = file("hull_calcs/".sprintf("%04d",$conf_index)."_insiders.txt",FILE_IGNORE_NEW_LINES);

    $arr_test_ids = range(0,2047);
    $arr_test_ids = array_diff($arr_test_ids, $refs);
    sort($arr_test_ids);
    for($i=0, $size=sizeof($arr_test_ids);$i<$size;$i++)
        $arr_atoms[$arr_test_ids[$i]]["in"]=FALSE;

    # Begin checking every passed atom
    foreach($arr_passed_atom_ids as $passed_atom_id)
    {
        $arr_atoms[$passed_atom_id]["in"] = TRUE;
        for($i=0, $size=sizeof($arr_atoms); $i<$size; $i++)
        {
            $distance = distance($arr_atoms[$passed_atom_id],$arr_atoms[$i]);
            //echo $passed_atom_id."\t---\t".$i."\n";
            if($distance < $Rcut_arr[$arr_atoms[$passed_atom_id]["id"]][$arr_atoms[$i]["id"]] && $distance>0.001)
                $arr_atoms[$passed_atom_id]["neighbour_count"] += 1;
        }
    }

    # Do the binning
    $CN_Na = Array();
    $CN_K  = Array();
    for($i=0, $size=sizeof($arr_atoms); $i<$size; $i++)
    {
        if(array_key_exists("neighbour_count",$arr_atoms[$i]))
        {
            ${"CN_".$arr_atoms[$i]['id']}[$arr_atoms[$i]["neighbour_count"]] += 1;
        }
    }

    ksort($CN_Na);
    ksort($CN_K);

    //print_r($arr_atoms);
    //print_r($CN_Na);
    //print_r($CN_K);

    # Report the results
    $filename = "results/".sprintf("%04d",$conf_index)."_results.txt";
    $fp = fopen($filename,"w");
    fwrite($fp, "### Atoms array ###\n");
    fwrite($fp,var_export($arr_atoms,TRUE)."\n");
    fwrite($fp, "\n### CN Distribution for Na ###\n");
    fwrite($fp,var_export($CN_Na,TRUE)."\n");
    fwrite($fp, "\n### CN Distribution for K ###\n");
    fwrite($fp,var_export($CN_K,TRUE)."\n");

    # Percentage calculation:
    $sum_Na = array_sum($CN_Na);
    $sum_K  = array_sum($CN_K);
    foreach($CN_Na as $key=>$value)
        $CN_Na[$key] = $value * 100 / $sum_Na;
    foreach($CN_K  as $key=>$value)
        $CN_K[$key]  = $value * 100 / $sum_K;
    //print_r($CN_Na);
    //print_r($CN_K);
    fwrite($fp, "\n### CN Distribution (Percentage) for Na ###\n");
    fwrite($fp,var_export($CN_Na,TRUE)."\n");
    fwrite($fp, "\n### CN Distribution (Percentage) for K ###\n");
    fwrite($fp,var_export($CN_K,TRUE)."\n");
    fclose($fp);

    # Write the summary to the results file
    $fp = fopen("results.txt","a");
    for($i=1;$i<=20;$i++)
    {
        if(!array_key_exists($i,$CN_Na))
            $CN_Na[$i] = 0;
        if(!array_key_exists($i,$CN_K))
            $CN_K[$i] = 0;
    }
    ksort($CN_Na);
    ksort($CN_K);
    fwrite($fp,sprintf("%04d",$conf_index)."\t".implode("\t",$CN_Na)."\t".implode("\t",$CN_K)."\n");
    fclose($fp);

}

function parse_types(&$arr_atoms)
{
    # Parse the types
    $i = 0;
    $fp = fopen("IP.DAT", "r");
    while(!feof($fp))
    {
        $line = fgets($fp,4096);
        $auxarr = explode(" ",$line);
        if($auxarr[0]==-1)
            $arr_atoms[$i]["id"] = "Na";
        else
            $arr_atoms[$i]["id"] = "K";
        $i++;
    }
    fclose($fp);
    array_pop($arr_atoms);
    return 0;
}

function get_vertices($conf_index)
{
    $arr_refs = Array();

    # Get the ids of the vertices of the convex hull
    $filename = "hull_calcs/".sprintf("%04d",$conf_index).".hull";
    $fp = fopen($filename, "r");

    # Bypass the first line
    $line = fgets($fp,4096);

    while(!feof($fp))
    {
        $line = fgets($fp,4096);
        $auxarr = explode(" ",$line);
        array_pop($auxarr);
        $arr_refs = array_merge($arr_refs, $auxarr);
    }
    fclose($fp);
    // $arr_refs = array_unique($arr_refs); # This doesn't lose the keys
    $arr_refs = array_keys(array_count_values($arr_refs)); # But this does.
    return $arr_refs;
}

function parse_coordinates(&$arr_atoms, $conf_index)
{
    # Parse the coordinates
    $i = 0;
    $filename = "hull_calcs/".sprintf("%04d",$conf_index);
    $fp = fopen($filename, "r");

    # Bypass the first two lines
    $line = fgets($fp,4096);
    $line = fgets($fp,4096);

    while(!feof($fp))
    {
        $line = fgets($fp,4096);
        $arr_atoms[$i]["coords"] = explode(" ",$line);
        array_shift($arr_atoms[$i]["coords"]);
        $i++;
    }
    fclose($fp);
    array_pop($arr_atoms);
    return 0;
}

function distance(&$atom1,&$atom2)
{
    # Calculate the distance between the two atoms
    $x1 = $atom1["coords"][0];
    $x2 = $atom2["coords"][0];
    $y1 = $atom1["coords"][1];
    $y2 = $atom2["coords"][1];
    $z1 = $atom1["coords"][2];
    $z2 = $atom2["coords"][2];
    return sqrt(pow($x1-$x2,2) + pow($y1-$y2,2) + pow($z1-$z2,2));
}

?>

<?PHP
/* Emre S. Tasci <e.tasci@tudelft.nl>                    *
 * Parses the results/xxxx_results.txt file and generates 
 * an XYZ file such that the atoms are labeled as the 
 * vertice atoms, close-vertice atoms and inside atoms..
 *                                              02/09/09 */

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

if($_GET["configuration"])
    $conf_index = $_GET["configuration"];
else 
    exit("Please specify a configuration number [php plotter.php configuration=xxx]\n");

$filename = sprintf("%04d",$conf_index);

# Get the atom array from the results file. Since results file 
# contains other arrays as well, we slice the file using sed for the 
# relevant part
$last_arr_line = exec('grep "### CN Distribution" results/'.$filename.'_results.txt -n -m1|sed "s:^\([0-9]\+\).*:\1:gi"');
exec('sed -n "2,'.($last_arr_line-1).'p" results/'.$filename.'_results.txt > system_tmp_array_dump_file');
$str=file_get_contents('system_tmp_array_dump_file');
$atom_arr=eval('return '.$str.';');
unlink('system_tmp_array_dump_file');

# Now that we have the coordinates, we itarate over the atoms to check 
# the characteristic and designate them in the XYZ file.
$fp = fopen("system_".$filename.".xyz","w");
fwrite($fp,sizeof($atom_arr)."\n");
fwrite($fp,"Generated by plotter.php\n");
for($i=0, $size=sizeof($atom_arr);$i<$size;$i++)
{
    if(!array_key_exists("in",$atom_arr[$i]))
        $atomlabel = "Mg";#Vertices
    elseif($atom_arr[$i]["in"]==TRUE)
        $atomlabel = $atom_arr[$i]["id"];#Inside atoms
    else
        $atomlabel = "O";#Close-vertice atoms
    fwrite($fp,$atomlabel."\t".implode("\t",$atom_arr[$i]["coords"]));
}
fclose($fp);
?>

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

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

Parameters

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

Example

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

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

Code

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

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

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

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

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

The code is as follows:

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

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

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

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

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

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.

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

sururi@dutsm0175 tmp $ cat booklet.sh
#!/bin/bash
#
# call with file.pdf
file=$1
filebase=$(basename $file .pdf)
pdftops $file output.ps
psbook output.ps tmp.ps
pstops "4:0L@.7(21cm,0)+1L@.7(21cm,14.85cm),2R@.7(0,29.7cm)+3R@.7(0,14.85cm)" tmp.ps > ${filebase}-booklet.ps
rm -f output.ps tmp.ps
echo "Converting back to pdf …"
ps2pdf ${filebase}-booklet.ps
rm -f ${filebase}-booklet.ps
sururi@dutsm0175 tmp $

There may be some problems when your path includes spaces (even if you escape them) and the first page may be in the wrong side (and inverted) but these are really minor annoyances compared to what this script achieves.

Usage is pretty simple:

booklet.sh <input_file_name.pdf>

and the output file will be input_file_name-booklet.pdf

Btw, this will be the first post under the new category Tools.