MethPlot Lib

Another fast line of code, so that I know where to look when I need it.
This time is MatplotLib, with LateX in the labels ;-)


import numpy as np
import pylab as pl
data = np.loadtxt('dijmsd5nm')
pl.rc('text', usetex=True)
pl.rc('font', family='serif')
pl.plot(data[:,0], data[:,1], 'ro')
pl.xlabel (r"$\displaystyle\sum_{n=1}^\infty 1$", fontsize= 12, color= 'black')
pl.ylabel(r'\textbf{time}(s)')

pl.show()




A little bit of explanation: in the folder where I use the Python script there are several data files, named something like "dijmsd5nm, dijmsd7_5nm, dijmsd10nm, dijmsd12_5nm, etc", where the underscore stands for the decimal point, since if I want to take these files in windows I am not allowed to have a point otherwise it can be seen as the beginning of the file extension.
If I don't do all the "re.sub" (the substituiting string characters in the file name) matplotlib plots the graphs as they come, i.e not in the good order 5nm->7.5nm->10nm->etc.
So I have to sort them first and then in the for loop that plots they are in order.

And here the refined super confusing version:


import numpy as np
import pylab as pl
import os
import itertools
import re

plotnum= 0
data=[]
dataname=[]
path= '.'
for i in os.listdir(path):
    if os.path.isfile and 'dij' in i:
        dataname.append(i)
        plotnum+=1
#sortthem!
for i in range(0,plotnum):
    for j in range(i,plotnum):
        nameI=re.sub("dijmsd","", dataname[i])
        nameI=re.sub("_", ".", nameI)
        nameI=float(re.sub("nm","",nameI))
        nameJ=re.sub("dijmsd","", dataname[j])
        nameJ=re.sub("_", ".", nameJ)
        nameJ=float(re.sub("nm","",nameJ))
        supp=dataname[i]
        if(nameJ<nameI):
            dataname[i]=dataname[j]
            dataname[j]=supp

for i in range(0,plotnum):
    data.append(np.loadtxt(dataname[i]))
marker = itertools.cycle(('o', '*', 'v', '^', '>', '<', 's', 'D'))

pl.figure(figsize=(8,6), dpi=180)

pl.rc('text', usetex=True)
pl.rc('font', family='serif')
for i in range(0,plotnum):
    name=re.sub("dijmsd","", dataname[i])
    name=re.sub("_", ".", name)
    pl.scatter(data[i][:,0], data[i][:,1],s=30, c=np.random.rand(3,), marker=next(marker), label=name)
pl.title(r"$\delta_{ij}^2$ per shell")
pl.xlabel (r"Shell")
pl.ylabel(r"$\delta_{ij}^2$ $\left[\r{A} ^2\right]$")
pl.xlim(0.,20.)
pl.ylim(0.02, 0.041)
pl.legend(loc="lower right")

pl.savefig("pictures/exercice_2.png", dpi=180)

pl.show()

BasHack

So i discovered that I basically use this blog myself because I don't remember my own scripts (I am a really confused worker) so here's is another that took me 15 min to remember and re write:


#!/bin/bash
for d in */; do
awk '(NR>1) {print $1 "\t" $3 }' ${d}statistics > dij${d%/}
awk '(NR>1) {print $1 "\t" $7 }' ${d}statistics > kij${d%/}
head -n -2 dij${d%/} > temp; mv temp dij${d%/}
head -n -2 kij${d%/} > temp; mv temp kij${d%/}
done 

SnHATEs and Pythons

Again, convert set of .xyz files into .data files to feed LAMMPS with.
This time I did it in Python and man, that was easier!

#!/usr/bin/env python
import sys, math

try:
infilename = sys.argv[1]; outfilename =sys.argv[2]
except:
print ("Error\nUsage:", sys.argv[0], "infile outfile"); sys.exit(1)

ifile = open(infilename, 'r')
ofile = open(outfilename, 'w')

xlo=0.
xhi=0.
atomtypes= 1
lines=ifile.readlines()
for line in lines:
pair=line.split()
if (line==lines[0]):
N=int(pair[0])

if (len(pair)==4 and pair[0][0]!='#'):
if (float(pair[1])>xhi):
xhi=float(pair[1])
if (float(pair[1])<xlo):
xlo=float(pair[1])



ofile.write("\n")
ofile.write ("%d atoms\n" %N)
ofile.write ("1 atom types\n")
ofile.write ("%g %g xlo xhi\n" %(xlo-20.,xhi+20.))
ofile.write ("%g %g ylo yhi\n" %(xlo-20.,xhi+20.))
ofile.write ("%g %g zlo zhi\n" %(xlo-20.,xhi+20.))
ofile.write("\nAtoms\n\n")

for i in range(2, len(lines)):
ofile.write (str(i-1) +" "+ lines[i])

print (xlo, xhi)

addiction to GNU pot

Since I do not have Origin (the scientific software, not the EA Steam-copy-paste) I tried to use alternatives like SciDavis or QtPlot. After struggling for hours to get what I needed, I went back to the roots using GNUPlot.

I hereby post a little script that I used to print the same function with different parameters. The function is a LogNormal distribution and I think that, for a fast visualization, my script is more then enough for everyone.
It's not that easy to find all the commands together in a single example and I don't want  to waste any more hours looking around.

Comments start with '#' character.



f(x,m,s,d)=1/((x-d)*s*sqrt(2*pi))*exp(-(log(x-d)-m)**2/(2*s*s))
#the function I want to plot with parameters m, s, d and variable x

set autoscale

set border 3 front 
set tics nomirror out scale 0.75
#these two lines remove the top and right borders, leaving only the 2 xy axes

set yrange [0:0.22]
set xrange [0:35]
#set the min-max range of the two axes x and yù

set ytics 0.025
set xtics 2
#set the spacing on axes, the tics

set ylabel "p(D)"
set xlabel "E.C.D. [nm]"
#sets the label to the axes

plot f(x, 2.38, 0.38, 0) notitle lt -1 lw 2
#plots the 1st function without graph legend, with lines of color -1 (black) and line width of 2 (thicker line than default)

replot f(x, 1.60, 0.46, 1.8129) notitle lt 1 lw 2
#plots the 2nd graph without erasing the first and changing only the color to 2 (red)

pause -1
#this is needed if you call gnuplot via terminal in order to have the image not disappears after few moments

Skullbasher and Scriptbasher

Then comes the time when you need to use a program (PowDog) that takes as input only file with format .qpow and you only have the .xyz files.
Moreover this program wants as inline command the atom element. So you script something with bash that transforms .xyz into .qpow and then starts the program N times with different inline command each time.

Mission accomplished.

#!/bin/bash 

for f in *.xyz; do 

   awk '{t=$1; $1=$2; $2=$3; $3=$4; $4=t; print }' $f > $f.qpowd 

done 

for f in *.qpowd; do 

#the "${f%%-*}" only takes the first 2 characters of the file's name 

   ./powdog -i $f -F form.factor -e "${f%%-*}" 

done

I sing the parsing song

One simple format to store your atomic positions is .xyz format.
The file starts with a number (N) that tells how many atoms (and how many rows) will follow.
The remaining file is made of 4 columns with type of atom and 3 cordinates. An example is the following:

4

Pd 1 1 2.
Pd 3 1.4 1.
Pd 1 4.3 11
Pd 4 1.2 3.

Now, for every simulation or elaboration, i need to read these kind of files and, preferably, create 4 arrays one of strings and the other 3 of float.
The best thing would be to have a separate function outside the main, in order to call it when needed.
So the values must be passed by reference and not by value, using C pointers. 
The problems arise because of the clumsy way C handles strings, which basically are arrays of char.
Moreover, I do not expect the atom type to be a string of more than 5 characters, so what I need to do when passing by reference my array of strings (which is a Nx6 matrix of char) is to send a triple pointer to the function that is, actually, a double pointer to a string and everything blows up without a careful use of & and *.
Here's the result.
That's it.  It took time to get it so simply and now I am happy and I wanted to share this in order to avoid a lot of complications to those interested.

1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <string.h>
4. 
   
5. 
   
6. void  parserxyz(char *fileneeded, double **atomx, double **atomy, double **atomz, char (**atomtype)[5], unsigned int *N)
7. {
8. FILE * fp;
9. unsigned int i;
11. fp = fopen(fileneeded, "r");
12. if (fp == NULL) {
13.  exit(EXIT_FAILURE);
14. }
15. fscanf(fp,"%u", N);
17. *atomx = (double*) malloc(sizeof(double)*(*N));
18. *atomy = (double*) malloc(sizeof(double)*(*N));
19. *atomz = (double*) malloc(sizeof(double)*(*N));
20. *atomtype= malloc(5*sizeof(char) * (*N));
22. for (i=0; i<(*N); ++i)
23. {
25. fscanf(fp,"%s %lf %lf %lf", (*atomtype)[i],  &(*atomx)[i], &(*atomy)[i], &(*atomz)[i] ); 
26. }
28. 
    
29. fclose(fp);
30. 
    
31. }

Introduction

I started my PhD.
I want to share what is going on and I want to remember what I have done.
I am basically a really bad student with a master degree in physics and now I work at the Engineering Department of my University.
English is not my first language and I will try my best.

My research proposal is all about working with Molecular Dynamics to simulate atomic and molecular systems.
Ok, nothing fancy like black holes or M-Theory: I am just not smart enough to make a living  out of theoretical physics (I was a theoretical physicist once, now I remain only a very mathematical oriented engineer).

Right, I have nothing more to say. Let's get started.