The ignorance of West Bengal Govt. about fundamental science

A proposal of a billion dollar science experimental facility to build a electron cyclotron in India was put forward in 2011 by the international committee headed by SINP (saha institute of nuclear physics,kolkata) a premier institute in this field.The machine Indian Synchrotron for Materials and Energy Research (ISMER) would be one of the largest in the world after USA,Germany,Japan,France. But only iff west bengal govt offered a 100 Acre land within the span of 50 kms from city and airport for easy transport for physicists across borders . But this is too much to ask for the government who believes in immediate benefit from investment and insists that it will be better to give the land for building AIIMS( medical school) than GOD-KNOWS-WHAT electron synchrotron. This chaps have not even heard of one before, although,  one similar technology is running right in front of their nose at VECC (variable energy cyclotron centre).

Bengal has been the hotspot for fundamental Physics education for long time. It has developed into one of the leading center for making of good scientists starting from Prof. Satyendra Nath Bose, the man who realized the famous “Bose particle”, a now fundamental part of physics to recently Dr. Ashok Sen for his outstanding research in fundamental string theory. What these two men have in common is they hailed from a place called the “intellectual” capital of India:West Bengal and have devoted their lives trying to understand the fundamental constituents of nature.

But the same Bengal, with due respect, is the home of another kind of intellectuals, the “wanna be” intellectuals, who quickly climb up ladder of high politics, who I guess cannot fathom the scale of development and progress a science project can build.

Coming to the points why this should be build:

1. Every nation must strive to achieve something hard/majestic, that is what sets them apart, remember them, just like Maharajahs are remembered by their wonderful monuments which sometimes took decades to complete, the modern day monuments are the Science and Technological marvels. Constructing this will take 10 years almost and employ thousands of people and also advance technology manifold.

2.Kolkata is home of many great scientist, and having this in homeland will make us even more proud.

3.Wiki gives the following uses for Science and Technology apart from fundamental physics:

Bengal have now been producing some software technicians like one in each hundred, only to brain drain to other state or even country doing some poor quality work for foreign demands. This project will employ enthusiast engineers in thousands. I have seen people from Bengal holding high posts in other states, intelligent people, scientists, artists fleeing this land, because of lack of opportunity, lack of cooperation from government, lack of understanding amongst those who take the decisions. Lets bring in a change, a change so that many people can really build something, a thing that will be remembered in history. Lets aim for long term benefit.

http://articles.timesofindia.indiatimes.com/2013-06-17/kolkata/40026493_1_sinp-deutsches-elektronen-synchrotron-saha-institute

http://www.telegraphindia.com/1130301/jsp/nation/story_16619819.jsp#.UgITxOGOWBs

http://www.telegraphindia.com/1130617/jsp/calcutta/story_17016076.jsp#.UgITzeGOWBs

Advertisements

Fast Fourier Transform (FFT) examples

The Discrete Fourier Transform(DFT) is defined by X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi \frac{k}{N} n} latex and the inverse fourier transform is defined as x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot e^{+i 2 \pi \frac{k}{N} n} latex

In [2]:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import time

'''
	Define two functions to compute fourier transforms using two different methods
	1.	dft(x,inverse=False,timer=False): for INVERSE: inverse=true
						  for timed: timer=true
	2.	timedfft(x,inverse=False,timer):for INVERSE: inverse=true
						  for timed: timer=true
	3.	func(x): defines function at x
'''

def func(x):
    f=10.0 
    val=np.sin(2*np.pi*f*x)
    return val

'''
   	Computes the discrete fourier transform /inverse fourier transform
	using basic summation convention.
	OPTIONS:
		inverse: computes inverse fourier transform
		timer:	 prints time required for computation
'''
def dft(x, inverse = False, timer = False) :

    t = time.clock()							#get current time to calculate the computation time
    N = len(x)
    if inverse:
        inv=1
    else:
        inv=-1

    X =[0] * N								#make an array of size N

    for k in xrange(N) :						#compute the required sum
        for n in xrange(N) :
            X[k] += x[n] * np.e**(inv * 2j * np.pi * k * n / N)
        if inverse :
            X[k] /= N

    t = time.clock() - t						#get delta time
    if timer :
        print "Computed DFT of size",N,"in",t,"sec."
    return X

'''
	fourier transform using built in fast fourier method
'''

def timedfft(x,inverse=False,timer=False):

	t=time.clock()
	if not inverse:X=np.fft.fft(x)
	else:	X=np.fft.ifft(x)
	t=time.clock()-t
	if timer: print "Computed DFT using fft of size",len(x),"in",t,"sec."
	return X

'''
	Computing part: 1. Define a function over a interval
			2. Compute fourier transform using the two methods
			3. Compare time
'''

#fix sample size, this gives the the rate of sampling of the function
sample_size=2**8
L=5.0
x=np.linspace(0,L,sample_size)

#there are four plots 1.original signal 2.fourier transformed signal 3.fast fourier transformed signal 4.Reconstructed signal

plt.subplot(411)
plt.xlabel('t')
plt.ylabel('f($t$)')
plt.xlim([0,L/10])
signal=func(x)
plt.plot(x,signal,label='original signal')

'''
				In [20]: run fourier_basic.py
				Computed DFT using fft of size 256 in 0.0 sec.
				Computed DFT of size 256 in 0.48 sec.
				Computed DFT of size 256 in 0.16 sec.

				In [21]: run fourier_basic.py
				Computed DFT using fft of size 4096 in 0.0 sec.
				Computed DFT of size 4096 in 123.94 sec.
				Computed DFT of size 4096 in 43.18 sec.
'''
#compute timed fft
fourierfft=timedfft(signal,timer=True)

#caluclate discrete fourier transform using normal method
fourier =dft(signal,timer=True)

#get the frequencies
n=x.size
d=x[1]-x[0]
freqsfft=np.fft.fftfreq(n,d)

#get positive index
indices=np.where(freqsfft>=0)

#plot
plt.subplot(412)
plt.xlabel('freq')
plt.ylabel('f($\omega$)')
plt.plot(freqsfft,fourierfft,label='FFT')

plt.subplot(413)
plt.xlabel('freq')
plt.ylabel('f($\omega$)')
plt.plot(freqsfft,fourier,color='green',label='Without FFT')

#inverse fourier transform using dft
invfourier=(dft(fourier,inverse=True,timer=True))

plt.subplot(414)
plt.xlabel('t')
plt.ylabel('f($t$)')
plt.xlim([0,L/10])
plt.plot(x,invfourier,'+-',color='red')
display(gcf())
Computed DFT using fft of size 256 in 0.0 sec.
Computed DFT of size 256 in 0.48 sec.
Computed DFT of size 256 in 0.16 sec.

/numpy/core/numeric.py:235: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

N.B.: The fourier transform of the function is symmertric about 0.
Thus it suffices to display only the positive half,
as the negetive half gives us no more information

In [3]:
clf()
plt.subplot(111)
pi=np.pi
sampled=5*np.sin(2*pi*x)+np.sin(2*pi*10*x)
fourier=np.fft.fft(sampled)
#if you want to print only the positive frequencies
index=np.where(freqsfft>=0)
plt.plot(freqsfft[index],fourier[index])
display(gcf())

In [39]:
clf()
sampled=np.exp(-x**2)
fourier=np.fft.fft(sampled)
plt.plot(freqsfft,fourier)
display(gcf())

In [40]:
clf()
from scipy.special import jn
sampled=jn(2,x)
fourier=np.fft.fft(sampled)
plt.plot(freqsfft,fourier)
display(gcf())

2D Fast Poisson Solver

2D Fast Poisson Solver:

Introduction:

Here we will construct poisson solver for a 2D square boundary using a variation
fast fourier transform , dicrete sine transform.There are 4 types of discrete sine transform,you can read more about them in
wikipedia http://en.wikipedia.org/wiki/Discrete_sine_transform

We will solve the 2D electrostatic equation for square boundary,in

\epsilon_0=1 units

The poisson equation solver in 2D basically uses the matrix representation of the discrete poisson equation.
The poisson equation \nabla^2 u= \rho can be solved using fourier transform technique.The algorithm is basically
from
http://www.math.nus.edu.sg/~bao/teach/ma5233/lect10.ppt

This uses dst function(Type -I) from either fft.fftpack or a seperate
dst package I downloaded from internet called transforms

CODE:

In [2]:
from __future__ import division
import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft,fft2,ifft2,ifft,irfft2,rfft2
import random as random
from mpl_toolkits.mplot3d import Axes3D

Define a function for discrete sin transform using fft algorithm,by symmetrizing the function x

In [3]:
def dst(x,axis=-1):
    """Discrete Sine Transform (DST-I)

    Implemented using 2(N+1)-point FFT
    xsym = r_[0,x,0,-x[::-1]]
    DST = (-imag(fft(xsym))/2)[1:(N+1)]

    adjusted to work over an arbitrary axis for entire n-dim array
    """
    n = len(x.shape)
    N = x.shape[axis]
    slices = [None]*3
    for k in range(3):
        slices[k] = []
        for j in range(n):
            slices[k].append(slice(None))
    newshape = list(x.shape)
    newshape[axis] = 2*(N+1)
    xsym = np.zeros(newshape,np.float)
    slices[0][axis] = slice(1,N+1)
    slices[1][axis] = slice(N+2,None)
    slices[2][axis] = slice(None,None,-1)
    for k in range(3):
        slices[k] = tuple(slices[k])
    xsym[slices[0]] = x
    xsym[slices[1]] = -x[slices[2]]
    DST = fft(xsym,axis=axis)
    #print xtilde
    return (-(DST.imag)/2)[slices[0]]

Define 2 dimensional DST, the idst is same as DST because of symmetry

In [4]:
def dst2(x,axes=(-1,-2)):
    return dst(dst(x,axis=axes[0]),axis=axes[1])

def idst2(x,axes=(-1,-2)):
    return dst(dst(x,axis=axes[0]),axis=axes[1])

Now define the poisson solver in 2D square boundary using the algorithm described above,this method can be
generalised to arbitary retangular boundary

In [5]:
def fft_poisson(f,h):

	m,n=f.shape
	f_bar=np.zeros([n,n])
	u_bar = f_bar			# make of the same shape
	u = u_bar

	f_bar=idst2(f)			# f_bar= fourier transform of f

	f_bar = f_bar * (2/n+1)**2  #Normalize
	#u_bar =np.zeros([n,n])
	pi=np.pi
	lam = np.arange(1,n+1)
	lam = -4/h**2 * (np.sin((lam*pi) / (2*(n + 1))))**2 #$compute $\lambda_x$
	#for rectangular domain add $lambda_y$
	for i in xrange(0,n):
		for j in xrange(0,n):
			u_bar[i,j] = (f_bar[i,j]) / (lam[i] + lam[j])

	u=dst2(u_bar)				#sine transform back
	u= u * (2/(n+1))**2			#normalize ,change for rectangular domain
	return u

Fix the initial parameters,domain of the equation etc.

In [6]:
#set bounds a,b,parameters
a = 0; b = 1;	
alpha=10				#alpha is grid points=2^alpha
n=2**alpha
L=b-a					#length of system

xe =np.linspace(a,b,n); 
ye =np.linspace(a,b,n);
x, y = np.meshgrid(xe,ye)
In [7]:
h=L/(n);			#size 
h2 = h**2;			#h squared
hx2=h2; 			#started with a cube,hx2=hy2
hy2=h2;
f=np.zeros([n,n]);

We put charges in the field.

Putting a random negetive charge in the field alog with a posituve charge

In [8]:
#initial conditions

#f[round(n/2):n,round(n/2):n]=1; #Initial condition
f[round((n+1)/2),round((n+1)/2-10)]=20
#f[round((n+1)/2),round((n+1)/2+10)]=-20
f[random.randint(0,n),random.randint(0,n)]=-10	#put a random charge
In [9]:
nx,ny=np.array(f.shape)-1		 #last index
U=np.zeros([n,n])

set boundary conditions

In [10]:
# BOUNDARY CONDITIONS
#set U(x,y)=g(x,y)
U[0,:]=0
U[nx,:]=0#5e-5
U[:,0]=0
U[:,ny]=0
In [11]:
##homogenize boundary condition

f[1,:]=f[1,:]+U[0,:]/hx2;
f[nx-1,:]=f[nx-1,:]+U[n-1,:]/hx2;

f[:,1]=f[:,1]+U[:,0]/hy2;
f[:,ny-1]=f[:,ny-1]+U[:,n-1]/hy2;

Solve the equation:

In [12]:
U=fft_poisson(f,h)
In [13]:
plt.figure()
#plt.imshow((U),cmap='hot')
plt.contour(U,50)
#plt.contour(f)

plt.show()
In [14]:
display(gcf())

In [15]:
imshow(U)
Out[15]:
<matplotlib.image.AxesImage at 0xb2b8a8c>
In [16]:
display(gcf())

 

In [ ]:

//

Hall Effect- a different outlook for seeing Hall coefficient

Hall effect is an extraordinary phenomena which involves generation of a subsidiary Electric field V_H when a electric field and magnetic field is applied to a conductor. You can get a basic introduction from wiki. If you are not familiar with it

But familiar Hall effect introduces it through a bound rectangular boundary, which limits the current destity in y direction to zero in steady state. In brief the Hall coefficient is defined as R_H=\frac{E_y}{J_x B}=\frac{1}{ne}

where ”J_x” is the current density of the carrier electrons, and E_y is the induced electric field,” n” is the carrier charge density and e is the carrier charge.

But consider the Hall effect in a 2D ring.The electric field applied in the ends if inner and outer ring and the magnetic field in the transverse direction.Then solving the Lorentz force equation along with Drude theory model (eqn 1)we get a 2 coupled equation,which can be solved directly

\frac{dp}{dt}=-e(E+v\times B)-\frac{p}{\tau} —–(1)

which gives from drude theory (see Ascroft Mermin)

\sigma_0 E_x=\omega_c \tau J_y+J_x ----(2) \\    \sigma_0 E_y=-\omega_c \tau J_x+J_y----(3)

Solving the two equations by multiplying the equations successively by p_x,p_y and adding/sutracting we get

\sigma_0 (E\times J)_z=\omega_c \tau J.J—-(4)

and

\sigma_0 (E. J)_{xy}= J.J—-(5)

where \sigma_0 (E\times J)_z=E_xJ_y-E_yJ_x—-(6)

\sigma_0 (E. J)_{xy}=J_xE_x+J_yE_y—-(7)

So we can define tan (\theta_H)=\omega_c\tau

which in common terms is the hall angle

where \omega_c=\frac{eH}{m}

Therefore the relation can be written in a more generally

\frac{tan (\theta_H)}{H}=\frac{\sigma}{ne} \\using \\ \sigma=\frac{ne^2\tau}{m}

The hall coefficient would then be defined by

R_H=\frac{tan (\theta_H)}{\sigma H}=\frac{1}{ne}

The angle between the net electric field and the net current is given by the angle tan (\theta)

This to me is a more general statement for Hall coefficient which reduces to the normal hall coefficient if we consider appropriate boundaries and noticing tan(theta) must be dimensionless and considering E_x<<E_y we get

R_H=\frac{E_y}{J_x H}=\frac{1}{ne}

which is the usual hall-coeficient.

Therefore the definition of Hall coefficient should be the tangent of Hall angle between the net Electric field and current in the conductor per unit applied magnetic field per unit conductivity.

Also this would help to resolve some issues regarding energy conservation in hall effect.

Suggestions are welcome.

Aritra Kundu

PS:here is an interesting article regarding hall effect in circular geometry http://www.jstor.org/stable/983911?seq=4

DL_POLY a quick user guide-Installation to first simulation in Molecular Dynamics

This summer I spend my time in Ecole Centrale Paris, trying to run a computer simulation on NPLIN.I will post about that later, but this was my first introduction to molecular dynamics simulation,  and there I faced many problems.So I thought to write up a bit, in case it might be of help for someone in future.I am posting this , but you sure can find help from the DISCO forums.And yes, read through the manual!!

In this, I explain how we can run DL_POLY in parallel.It really did a 30 minuite simulation in 30 seconds in my 4X2 core processor in laptop

Installation:
1.Download DL_POLY_4 from ftp://ftp.dl.ac.uk/ccp5/DL_POLY/
2.Obtain Liscense by mailing the authors for home/personal/educational use
3.Unzip,Untar,Unencrypt
4.In Linux install open mpi by typing
$ sudo apt-get install libopenmpi-dev openmpi-bin openmpi-doc
$ sudo apt-get install ssh
Generate keys and add to ssh
$ ssh-keygen -t dsa
$ cat id_dsa.pub >> authorized_keys

Test program:
#include <stdio.h>
#include <mpi.h>

int main(int argc, char** argv)

{
int myrank, nprocs;

MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
for (;;){
printf(“Hello from processor %d of %d\n”, myrank, nprocs);
}
MPI_Finalize();
return 0;
}

Testing:
$ mpicc MPI_Hello.c -o MPI_Hello
$ mpiexec -n 5 MPI_Hello         #specify to run on 5 processors,press Ctrl+Z to exit
You can install htop to see if the process is really parallel [sudo apt-get install htop]
Compiling DL_POLY for Laptop Dell XPS 15z
Requires fortran to be installed.[sudo apt-get install gfortran]
1.In the DL_POLY directory>Build>copy “Makefile_MPI” to source directory and rename to “make”
2.In terminal type “sudo make hpc”
3.The process will create a DLPOLY.Z in the execute directory. Which also has a JAVA gui which can be used.(Personally I find using GUI a bit cumbersome and problematic,so I run using terminal)
4. Terminal >browse to DL_POLY execcute directory> type $./gui to invoke gui

Running sample scripts:
Some example and basic examples can be found online in
http://www.ccp5.ac.uk/DL_POLY/TUTORIAL/EXERCISES/

Although they used GUI, I am going to use a mixture of both

DLPOLY requires some necessary files to run a simulation,The CONTROL,CONFIG,FIELD files.I am going to explain them as we go,but first lets see if its working!
Download the files for EXERCISE 1 which finds out the phase transition for KCl.
The CONFIG file lists the system in brief,it gives the coordinates and velocities of all atoms in the system in initial condition.
The CONTROL file is important,it determins the thermodynamics of the system.The temperature, pressure, cutoff, and configuration like NVT,NPT etc.
FIELD defines any other field which cannot be characterised by concrete mathematical formulas

The CONFIG file uses its own self-consistent set of units.

physical quantity symbol unit value
time t0 1×10−12 seconds (picoseconds)
length l0 1×10−10 metres (Angstroms)
mass m0 1.66054×10−27 kilograms (amu)
charge q0 1.60218×10−19 Coulombs (electron charge)
energy E0=m0(l0/t0)2 1.66054×10−23 Joules = 10 J mol−1
pressure p0=E0/l30 1.66054×107 Pascal = 166.054 bar
Planck constant 6.35078E0t0
Boltzmann constant kB 0.83145E0/K

The others can be found from the user manual(Page 7 of DL_POLY_4_USER_MANUAL)

will learn the process through a process through the first example of DL_POLY guide.i.e. phase transition of KCl

http://www.ccp5.ac.uk/DL_POLY/TUTORIAL/EXERCISES/expmt1.html

This exercise studies a well-known phase transition in potassium chloride using constant pressure molecular dynamics. The objective is to develop the best practice in using such algorithms and to learn how phase transitions can be induced, detected and monitored in a simulation.

The CONFIG file:

DL_POLY: Potassium Chloride Test Case                                           #heading
2         3      2000    0.5000000000E-02

!levcfg(0-2,2 :position,velocity,forces included in file)

! imcon(periodic boundary condition:0,1,2,3,6)

!megatm(total number of particles)

!crystallographic information for periodic boundary(x,y,z component of a,b,c cell vector)

18.6960000000        0.0000000000        0.0000000000
0.0000000000       18.6960000000        0.0000000000
0.0000000000        0.0000000000       18.6960000000
!since levcfg = 2 each record will contain 3 data position,veocity,force
K+               1                            #atom name, index
-7.608595309        -7.897790000        -7.892053559    #position
1.056610291        -1.218664448         3.345828610    #velocity
-1979.558687         739.7961625         1027.996603    #force
K+               2
-4.678562643        -4.509858530        -7.813252872
-0.7078383910         1.220785303       -0.9440783885
212.7300659        -4290.721728         1363.866117
K+               3
-7.821414265        -4.635443539        -4.732164540
2.637614561        0.5778767520E-01    -1.704765568
150.3309776        -812.6932914         1429.413120
K+               4
-4.577735234        -7.808222424        -4.540282303
4.950411107         3.576253827        0.2276058039
-782.7309730         1322.553499        -205.1240756
Cl-              5
-7.810782981        -4.845250366        -7.815286768
1.643324971       -0.6479194620         3.393130560
407.2394357         1675.508460         748.2506938

The details about preparing this file with other packages(Athen) will be discussed
Next we prepare the CONTROL file:

**–begins here–<We will at first discuss only relevant parameters>–**
DL_POLY: Potassium Chloride Test Case    !heading
temperature      600.0                !system temperature in Kelvin
pressure         144.0                !system pressure
steps             4000                !maximum no of steps,2000 is enough
equilibration      500                  !time for the system to equlibrate
timestep         0.005
multiple             1
print               10
stats               10
scale               10
restart scale
rdf                 10             !write radial distribution function at 10 timescale interval
print rdf
cutoff             7.0                !define cutoff for interactions to 7 Angstorms
rvdw               7.0
delr               1.0
ewald precision  1.e-6
ensemble nst hoover 0.1 1.0            !Ensemble
trajectory 20 30 0                !write history file
job time          3600                !set job time
close time         100
finish                        !control file ends

FIELD FILE:

DL_POLY: Potassium Chloride Test Case

UNITS eV    !energy units

MOLECULES 1    ! number of molecule type

Potassium Chloride    !Molecule name

nummols 27        !number of molecules

atoms 8        ! no of atoms in this molecular type

K+           39.0983       0.994    4

Cl-          35.4530      -0.994    4

finish

VDW 3            !number of van der wall interactions

K+      K+      buck      3796.9      0.2603      124.90

Cl-     Cl-     buck      1227.2      0.3214      124.90

K+      Cl-     buck      4117.9      0.3048      0.0000

CLOSE

RUNNING DLPOLY:

To run the program make sure the CONFIG,CONTROL,FIELD files are there in the execute folder, from terminal invoke dlpoly

./dlpoly.z

 

This will run the simulation and write the output files

ANALYZING OUTPUTS:

 

  1. We will visualize the resulting system.This can be done through GUI by plotting the REVCON file(which is the system at the end of simulation)
    You can also visualize the CONFIG file.
    Another better alternative of visualizing the files is to convert the CONFIG file to XYZ files and plotting them on VMD

  2. Checking the Diffusion coefficient in OUTPUT file, if it is less than 10-5then it is liquid state.A value of zero suggests it to be solid state.

  3. Now we will analyze the output from the simulation using the RDFDAT file that was created.From GUI plot the RDF -> select the atom types K+ and Cl-

    This will extract the RDF for K+ and Cl- from the RDFDAT file and create a separate XY file.For various simulation parameters you can create different RDF files and merge them together and plot so that you can compare where the system is heading to and draw good conclusions.I have written a small script in python to do that.

# Import modules

import sys

import matplotlib.pyplot as plt

import numpy as np

”’

to be invoked as

plot_rdf <option={-m,-p}> <no of files> <list of files>

option -m merges the given number of files using single X cordinates and

y coordinate.Option -p plots the given file

eg

plot_rdf -m 2 RDF_100.XY RDF_14.XY

merges the two files

”’

#import pylab as pl

option=sys.argv[1]

no_of_options=int(sys.argv[2])

index=3

l_index=index+no_of_options

#process data

if(option==’-m’):

   data=np.loadtxt(sys.argv[index])

   shape_data=data.shape

   while(index<l_index-1):

       index=index+1

       data2=np.loadtxt(sys.argv[index])[:,1].reshape(shape_data[0],1)

       data=np.append(data,data2,1)

   f=open(‘RDF_Merged.XY’,’w’)

#TODO: write names of field and parameters in first column

   for row in data:

       X=”

       for element in row:

           X=X+’ ‘+str(element)

       f.write(X+’\n’)

   f.close()

elif (option==’-p’):

   # Open the data

   datafile = sys.argv[index]

   #f = open(datafile, ‘r’)

   data = np.loadtxt(datafile)

   shape_data=data.shape

   i=1

   #plot all Y’s

   while(i<shape_data[1]):

       plt.plot(data[:,0],data[:,i])

       i=i+1

   # Get axis labels

   xaxis_label = raw_input(‘X-axis label > ‘)

   yaxis_label = raw_input(‘Y-axis label > ‘)

   # Show the plot

   plt.ylabel(yaxis_label)

   plt.xlabel(xaxis_label)

   plt.show()

#Stop

The RDF is the measure of how the probability deviates from an uniform distribution.From the radial distribution function a number of parameters can be determined.We can also identify the phase transition from the change in pattern in the RDF

Picture stolen from Jacek Dziedzic,Gdansk University of Technology,Computational nanotechnology

Diracs Puzzle

I will post here with some musings I came across another blog:
http://www.minustwofish.com/Here is a question that was made to Dirac, the great genius !Here is it

Three fishermen come back from the sea. Each collapses in their respective tent. Fisherman #1 wakes up and decides to get his share of the catch. He counts the fish, realizes it is not a number divisible by three, throws away one fish to the sea correcting the situation, and takes a third of the remaining fish into his tent. Fisherman #2 wakes up later and also decides he is going to get his share of the catch. Fisherman #2 wants a third of the fish he sees. It is not a multiple of three, but he throws away one fish and takes a third of the fish and goes to sleep. Fisherman #3 wakes up after, and does the same: he throws away one fish, takes a third of the fish, and goes into his tent.
What is the smallest number of fish for which this would happen?

Ah its simple isnt it?

Guess what was Dirac’s answer -2! no wonder he came up with the idea of filled negative electrons in Dirac sea

Quantum Harmonic Oscillator -Simulations

Heya…For the past couple of weeks I have been doing simulations for quantum systems, like the harmonic oscillator(mainly) and square well..as a part of the course..We produced some really nice graphics (using python) for multiple harmonic wells.
In this post I will be sharing some of the simulations..It can give a good insight to the physics of quantum mechanics.
Plots of wavefunctions for  multiple harmonic wells

Some useful points to be noted
  1. The degeneracy in each levels increase as the number of wells is increased.(remember the lattice structure?)
  2. For large number of wells, it begins to show band structure.
  3. The separation between the wells determine the magnitude of their interaction

formation of band structure

Random Walker – Python Simulation

SyntaxHighlighter.config.bloggerMode=true; SyntaxHighlighter.config.clipboardSwf=’http://alexgorbatchev.com/pub/sh/current/scripts/clipboard.swf&#8217;; SyntaxHighlighter.all(); I havent posted for long ..have been quite busy..here is some very simple code to animate random walk I just love seeing how the random walk simulates…experiment with it..and freely redistribute..;) Recently Im learing python simulations…so this was what i did while reading statistical mechanics by richel

# Random_walker.py
# Uses turtle class to simulate random walk through pseudo random generator
# external interrupt is enabled ,such that one can give external impacts
# to the random walker and see if that small external interrupts change random walk pattern
#1.give brisk shocks with up /right/left keys
#2.try combinations
# Aritra Kundu
# arsenous.blogspot.com


import turtle
import numpy as np
turtle.setup(400,500) # determine the window size
wn = turtle.Screen() # get a reference to the window
wn.title("Handling keypresses!") # Change the window title
tess = turtle.Turtle() # Create our favorite turtle

# The next four functions are our "event handlers".
def h1():
tess.forward(30)

def h2():
tess.left(45)

def h3():
tess.right(45)

def h4():
wn.bye() # Close down the turtle window

#These lines "wire up" keypresses to the handlers we've defined.
wn.onkey(h1, "Up")
wn.onkey(h2, "Left")
wn.onkey(h3, "Right")
wn.onkey(h4, "q")

wn.listen()

#listen responses from keyboard
N=1000 # no of steps
for i in xrange(1,N):
x=np.random.randint(-15,15) #step size for random walker,generate integer
theta=np.random.randint(0,360) # rotation angle
tess.forward(x)
tess.right(theta)

turtle.clear()
#wn.mainloop()

The Physics of Love –part 1

The Beginning:

Finding a Girlfriend

The drake equation,the LOVE equation

The Drake equation describes the possibility of finding extraterrestrial intelligent life in outer space,who can possibly contact us using radio signals.It was founded by astrobiologist Dr.Frank Drake in 1961,Emeritus Professor of Astronomy and Astrophysics at the  University of California,Santa Cruz.You can read more about Drake equation at wikipedia
But what has a equation describing possibility of exo life got to do with finding a girlfriend.Well it turns out fitting in right parameters will give a rough estimate according to your needs.This is a variation of the Secretary problem. Although it has low success rate (37 %)..well its good to have something than nothing…The drake equation :
G=R.fp.ne.fl.fi.fe.L
G=no of civilization capable of interstellar communication
R=rate of formation of stars
f =fraction of stars that has planets
ne=planets similar to earth
fl=fraction of earth like planets supporting life of any kind
fi=fraction when intelligent life develops
L=length of time such civilizations survive

Some other people some other constraints,well that merely reduces the possibility to find the number of civilizations.

The drake equation can be modified to represent the finding of a compatible match for two species who settled on earth long back,a guy(from mars) with a girl(from venus).I call it the “love theory”.First done in 1999.[see link 1]

In love theory:
G=R.fw.fb .fa.fu..fs 
G=no of potential girls for your match
R=total population of India(if you dont want to have girlfiends from abroad!!)
fw=fraction of the population who are of opposite sex(in my case i will consider female)
fb=fraction of them in your state(bengal,in my case)
fa=fraction of women with appropriate age
fu=fraction of them having university education
fs=and still harder,fraction of them whom i find attractive.
R=1.21billion
fw=.486
fb=0.08
fa=.24 [estimate 63.4/3,wiki]
fu=.30
fs=.2

summing up these are realy low valus of finding a girlfriend.
G=1.21*.0005 billion
=.0006774 billion=677400 people all over india.
this excludes the fact that the woman will find you attractive too.To incoparate this we can (approximately) just square the all the terms except R in expression as the same situation applies to the girl(well the coefficients for girls should be multiplied in fact ).Thus the net total chance decreases.
G=R.(fw.fb .fa.fu..fs .L)2
which gives an estimate of 302 matches.Thus in lifetime we have 302 girls who can potentially be our life partners….pretty low ? but this excludes the fact that the girls may be already committed…

But still we find one of them right there waiting for us…

hope you enjoyed this..more information on this can be found by searching “why i dont have a girl friend”
N.B.:Ever since a boy starts brooding over his incapability and failure in love life,and his futile search over years(for a grad student),he concentrated more in his alcoholic work,and even brings down the probability of finding a new one.
References: