MDNewton - ver.0.2	
	MD-Newton Ver.0.2 / Wed Mar 2. 1994

This is a Mathematica package for the two-dimensional classical 
Molecular Dynamics calculation of Lenerd-Jones potential system.  

Features of MDNewton ver.0.2 : 
[1] All of the codes of MDNewton were written in Mathematica programming language. Where we referred to the book IPM[3] to improve the codes.
    
[2] The particle's motions were calculated using Newtonian equation of motion.(see MDSS[1] and CSM[2] for more datails)

Please enjoy the world of MDNewton.

- kazume - Iwate university Japan
nisidate@wriron1.s.kanazawa-u.ac.jp


Refs.
[1] Molecular-Dynamics Simulation of Statistical-Mechanical
    Systems (Proceedings of the International School of 
    Physics "Enrico Fermi" Course XCVII)
     G. Ciccoti and W. G. Hoover
     North-Holland Physics Publishing, a devidion of Elserver 
     Science Publishers B. V., 1986
[2] Computer Simulation Methods in Theoretical Physics 
	(2nd edn) 
	 D. W. Heermann
	 Springer-Verlag 1990
[3] Introduction to Programming with Mathematica
	 Richard J. Gaylord, Samuel N. Kamin, Paul R. Wellin
	 TELOS / Springer-Verlag 1993

(* Parameters setting *)
num=3;       (* num^2 = particles *)
rcutt=2;     (* cutt off length of the forces *)
sigma=0.3;   (* L. J. potential parameter *)

lattice=(Flatten[Table[{ i/#-1/(2#), j/#-1/(2#) },{i,1,#},{j,1,#}]])&;

(* Set up the sq lattice for the atoms inside the sq *)
sq=lattice[num];

(* Check the sq lattice *)
Show[Graphics[{PointSize[0.05],RGBColor[0.5,0,1],
 Map[Point,Partition[Flatten[sq],2]]} ],Axes->True,
 AxesLabel->{"X","Y"},
    PlotRange->{{0,1},{0,1}} ];

(* Preperling the dumy matrix *)
dumy1=Table[0,{num^2}];
dumy2=Table[0,{i,1,num^2},{j,1,num^2}];

(* Baundaly condition function *)
bcx=(pos=#;pos-Floor[pos])&;

(* Nearest neighbour view *)
Rdpos=(nnRxy=#1;
	Table[nnRxy[[i,j]]=bcx[#2[[i]]-#2[[j]]],
		{i,1,num^2},{j,i+1,num^2}];nnRxy)&;

(* Delivated Lennard-Jones potential function *)
phi=( ( (sigma/#)^12 - 1/2 (sigma/#)^6 )/#^2)&;

(* MDNewton *)
Newton=({pos=#1};
   ffx=dumy1; ffy=dumy1; (* initialize the forces *)
   posx=Flatten[Table[pos[[i]],{i,1,2num^2,2}]];
   posy=Flatten[Table[pos[[i]],{i,2,2num^2,2}]];

   Rdx=Rdpos[dumy2,posx]; Rdy=Rdpos[dumy2,posy];
      
	Rdus=(Rdx^2+Rdy^2)^0.5; (* distances *)
   Table[mtlfx=0;mtlfy=0;  (* mutual forces *)
      Table[If[Rdus[[i,j]]<rcutt,uij=phi[Rdus[[i,j]]];
	           ffx[[j]]=ffx[[j]]-uij Rdx[[i,j]];
               ffy[[j]]=ffy[[j]]-uij Rdy[[i,j]]];
               mtlfx=mtlfx+uij Rdx[[i,j]];
               mtlfy=mtlfy+uij Rdy[[i,j]],
      {j,i+1,num^2}];
      ffx[[i]]=ffx[[i]]+mtlfx; 
      ffy[[i]]=ffy[[i]]+mtlfy,
   {i,1,num^2}];
   VX=VX+ffx dt/2; VY=VY+ffy dt/2;  (* velocity *)
   posx=bcx[posx+VX dt]; posy=bcx[posy+VY dt]; (* bc *)
   pos=Flatten[Table[{posx[[i]],posy[[i]]},{i,1,num^2}]])&;

(* running the Newton *)

dt=0.005;     (* time step *)
rcutt=0.4;     (* cutt off parameter *) 
sigma=0.1;     (* L. J. potential parameter *)
positions=lattice[num];  (* initial positions *)

(* Random distributed initial verocity *)
(* Needs["Statistics`ContinuousDistributions.m`"]  *)
<</usr/local/math/Packages/Statistics/ContinuousDistributions.m

SeedRandom[1];
VX=Table[Random[NormalDistribution[0.05,0.003]],{num^2}]; 
VY=Table[Random[NormalDistribution[0.05,0.003]],{num^2}];

(* checking the normal distributed velocity *)
ListPlot[Sort[Abs[Join[VX,VY]]]]

(* Checking the Lennard-Jones potential *)
Plot[phi[x],{x,sigma/2,2sigma}, PlotRange->{-10,10}]

cycle = 100;      (* Snap Shot cycle *)
SnapShots = 10;   (* Number of the pictures *)

Timing[
lis=NestList[Nest[Newton,#,cycle]&,positions,SnapShots];]

Table[
Show[Graphics[{PointSize[0.05],RGBColor[0.5,0,1],
	Map[Point,Partition[Flatten[lis[[i]]],2]]} ],
	Axes->True,AxesLabel->{"X","Y"},
		PlotRange->{{0,1},{0,1}} ],{i,1,SnapShots+1}];