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 [ ]:

//

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

Electrostatics through simulation

Here are the results of solving the poisson equation in 2D boundary: I wrote this code for our class which will go opensource book soon-involving the introduction of python in physics teaching My solver uses fft(fast fourier trnasform) hence its pretty much fast. Enjoy the Pics.Write to me if you find problems,bugs,or dont agree on a thing : 1.Arbitrary charge distribution:

2.Dipole Charge distribution

3.Quardapole Charge distribution.

4. Unstable twi dipole side by side.