Introduction 2

Topics to cover:

Data manipulation. List and Arrays. Solving Matrix Equations. Curve fitting by least squares method. Numerical Extrapolation.

Refer to Chapter 3, DeVries for certain parts of the topics discussed.

data={{0,-1.1},{1,2.1},{2,3.1},{3,3.98},{4,5.2},{6,7.22}};
lpdata=ListPlot[data,Joined->False,Mesh->All];
fitdata=Fit[data, {1,x},x]
a=fitdata[[1]];
b=fitdata[[2]];
const=0.023142857142859113;slope=1.27257;
y[x_]:=const + slope*x;
py=Plot[y[x],{x,0,6}];
Show[lpdata,py]

data={{0,-1.1},{1,2.1},{2,3.1},{3,3.98},{4,5.2},{6,7.22}};
lpdata=ListPlot[data,Joined->False,Mesh->All];
fitdata=Fit[data, {1,x,x^2},x];
a2=fitdata[[1]];
b2=fitdata[[2]];
b2=Coefficient[b2,x];
c2=fitdata[[3]];
c2=Coefficient[c2,x^2];
y2[x_]:=a2+ b2*x+c2*x^2;
py2:=Plot[y2[x],{x,0,6}];
Show[lpdata,py,py2]

{a,b}
{a2,b2,c2}

b

-0.05853266888150617` (x^2)

NSolve[y[x]==0,x]

Data manipulation: Interpolation

data={{0,-1.1},{1,2.1},{2,3.1},{3,3.98},{4,5.2},{6,7.22}};
y[i_]:=data[[i]][[2]];
ty=Table[y[i],{i,1,6}];
f=ListInterpolation[ty,{{0,6}}];
lpdata=ListPlot[data];
pf=Plot[f[x],{x,0,6}];
Show[lpdata,pf]
xroot=x/.FindRoot[f[x]==0,{x,0}];
f[xroot]

sf=SplineFit[data, Bezier];
sf

f[6]

y[4]

(* table *)
t={a,b,c,d};
et[i_]:=t[[i]];
Table[et[i],{i,1,4}]

data={{0,-1.1},{1,2.1},{2,3.1},{3,3.98},{4,5.2},{6,7.22}};

Length[data]

A={{11,12,13},{21,22,23},{31,32,33}};
B=MatrixForm[{{11,12,13},{21,22,23},{31,32,33}}];
A
B

X={a,b};Y={0,0};
A={{1,1},{1,1}};
Solve[A.X==Y]

m.{x,y}=={a,b}

Data manipulation: Creating data files and how to post-process them.

Exercise: Write a loop to generate say 10 data files, named 1.dat, 2.dat, ...,10.dat. Each *.dat file should contain a randomly generated ‘0’ or ‘1’. Then write another independent program to scan through the *.dat files in the directory. Read in these *.dat sequentially, begining from 1.dat, then 2.dat, ...,10.dat. Identify which of the *.dat is the first one that contains a ‘1’.

The purpose of the above exercise is to give you an idea how to export your numerical results into properly sequenced *.dat files, and how to postprocess them using simple programming instructions. Many data manipulation in actual situations involved the similar procedures as the above exercise.

Assignment: Write a loop to generate say 10 data files, named 1.dat, 2.dat, ...,10.dat. Each *.dat file should contain a randomly generated ‘0’ or ‘1’. Then write another independent program that determines which *.dat files are the 6-th one (starting from 0.dat.) that contains a ‘0’. You program should be smart enough to warn you if there is no such *.dat file exist.

Visualation of xyz data: Download the zipped file at http://www2.fizik.usm.my/tlyoon/teaching/ZCE111/1213SEM1/notes/Introduction/xyzdata/xyzdata.zip
Unzip it. Use mathematica to import the file named “1.dat”. This data file contains the positions of many particles at the timestep t =1, hence the name of the data file “1.dat”. Similarly, data n.dat contains the positions of all the atoms at timestep n. 

Q1: How many atoms are there in the data file?
Q2. Visualise the 3D configuration of these atoms using mathematica. The command to use: ListPointPlot3D[dat1, BoxRatios -> {1, 1, 1}]ListPointPlot3D[dat1, BoxRatios -> {1, 1, 1}].
Q3. Import the files in numerical order and display them one by one in a sequence. 
Q4. Do this using Manipulate.

SetDirectory[NotebookDirectory[]];
FileNames[];
Import["methane.xyz"]

SetDirectory[NotebookDirectory[]];
FileNames[];
Import["sample.xyz"]

SetDirectory[NotebookDirectory[]];
sample=Import["sample.xyz","Table"]

count=0;
For[i=1,i≤Length[sample],i++,
If[Length[sample[[i]]]==4,
Print[sample[[i]]];
count=count+1;
P[count]=Take[sample[[i]],-3];
Print[P[count]];
];
];
lastcount=count;
tP=Table[P[i],{i,1,lastcount}];
Export["sample0.dat",tP];
FileNames[]

list={"H",0.`,0.`,1.089`};
ele[i_]:=list[[i]];
ele[4]

list={"H",0.`,0.`,1.089`};
Take[list,-3]

dummy=Import["sample0.dat"];
P[1]=dummy[[1]];P[2]=dummy[[2]];
P[3]=dummy[[3]];

Deltaphi=0.1Pi;
T={
{1,0,0},
{0,Cos[Deltaphi],-Sin[Deltaphi]},
{0,Sin[Deltaphi],Cos[Deltaphi]}
};
Pprime[1]=T.P[1];
Pprime[2]=T.P[2];
Pprime[3]=T.P[3];
tPprim=Table[Pprime[i],{i,1,3}]
tP

Graphics3D[Line[tPprim]]

Graphics3D[
Point[tPprim],PlotRange->{{-2,2},{-2,2},{-2,2}}];
Graphics3D[Line[tPprim]];
Graphics3D[
{Line[tPprim],
Point[tPprim],PlotRange->{{-2,2},{-2,2},{-2,2},AspectRatio->1}
}
]

Graphics3D[
Point[tPprim],PlotRange->{{-2,2},{-2,2},{-2,2}}];
Graphics3D[Line[tPprim]];
Graphics3D[
{Line[tP],
Point[tP],PlotRange->{{-2,2},{-2,2},{-2,2},AspectRatio->1}
}
]

P12=P[2]-P[1];P32=P[3]-P[2];P31=P[3]-P[1];
theta13=P12.P32
theta23=P31.P12
theta12=P32.P31

P12prime=Pprime[2]-Pprime[1];P32prime=Pprime[3]-Pprime[2];P31prime=Pprime[3]-Pprime[1];
theta13=P12prime.P32prime
theta23=P31prime.P12prime
theta12=P32prime.P31prime

Export["sample1.dat",tPprim]

Download the file bdsimplified.dat from http://www2.fizik.usm.my/tlyoon/teaching/ZCE111/1213SEM1/notes/Introduction/xyzdata/.This is a structured data file with a format. The first colume is the index assigned to an atom, column 2 is the type of atom, and column 3 is the biding energy of that atoms. For example, in the second row of bdsimplified.dat,

2 3 -5.76691 

means the atoms is identified as atom with index 2, of type 3, and its binding energy is -5.76691. Write a program to calculate the average biding energy for atom type 3.

Download the file bd1.dat from http://www2.fizik.usm.my/tlyoon/teaching/ZCE111/1213SEM1/notes/Introduction/xyzdata/.This is a structured data file with a slightly more complicated format. It keeps the record of the binding energy of the atoms at all time step. Write a program to that plots the average binding energy of atom type 3 as a function of time step.

Homework assignment: Download http://www2.fizik.usm.my/tlyoon/teaching/ZCE111/1213SEM1/notes/Introduction/xyzdata/bd2.dat. .This is a structured data file with a slightly more complicated format. It keeps the record of the binding energy of the atoms at all time step. Write a program to that plots the average binding energy of atom type 1 as a function of time step. Modify your program to plot the average binding energy of atom type 2 as a function of time step. Note that the biding energy for atom type 1 is in column 3, while  the biding energy for atom type 2 is in column 4.

List and Arrays: Some exercises.

Solving Matrix Equations using Mathematica.