Argon_MD.java import javax.swing.*; // "Argon_MD" simulates the ...

snottybugbearΛογισμικό & κατασκευή λογ/κού

3 Νοε 2013 (πριν από 4 χρόνια και 7 μέρες)

116 εμφανίσεις

Argon_MD.java
import javax.swing.*;
// "Argon_MD" simulates the motion of argon atoms. This is a java
// applet
, and its size should be set to 500x400 pixel
public class Argon_MD
 extends JApplet implements Runnable {
int simStatus=0;
////////////////////////////////////
// enter global declaration here
////////////////////////////////////
long nTime       = 0;          // time step counter
int nParticle    = 25;         // the number of particle
double h         = 0.01;       // time step [pico
‐second]
double boxLength = 6.0;        // contrainer's length [nano
‐meter]
double mass      = 40.0;       // argon's mass in dalton
double LJ_Emin   = 0.9964;     // Lennard
‐Jones
 energy minimum [kJ/mol
]
double LJ_Rmin   = 0.3827;     // Lennard
‐Jones
 equilibrium radius [nano
‐meter]
double kB        = 8.31446E‐3; // Boltzman's constant [kJ/mol
/K]
double v0        = 0.05;       // initial speed of each particle [nm
/ps
]
double x[], y[], vx[], vy[], ax[], ay[]; // coordinate, velocity, and acceleration
int nRDFInterval;                        // number of RDF histogram intervals
double aRDF, bRDF;                       // min
 and max range to probe RDF
double rPointRDF[];                      // radius to probe RDF
double solidRDF[], liquidRDF[], gasRDF[];// RDF for 3 phases
int delay;                               // delay in millisecond
Thread animator;                         // animating thread
JSpinner initSpeedSpinner;               // initial speed selector
AnimationMDPanel animationMDPanel;       // particle animation screen
GraphPanel energyGraph;                  // graph displaying energies
GraphPanel RDFPanel;                     // RDF graphs
// "KE" computes total kinetic
 energy of the system
double KE(double vx[], double vy[]){
double E=0; // calculated kinetic
 energy 
int i;      // loop index
// loop thru every particle
for(i=0; i<nParticle; i++)
E = E + 1/2.0 * mass * (vx[i]*vx[i] + vy[i]*vy[i]);
return E;
}
// "LJ_Energy" computes potential energy due to Lennard
‐Jones
 term
double LJ_Energy(double x[], double y[]){
double E=0;   // calculated energy
double ratio; // LJ_Rmin to atomic distance ratio
int i,j;      // particle index
// loop thru all possible pairs 
for(i=0; i<nParticle; i++)
for(j=0; j<nParticle; j++)
if(i<j){
ratio = LJ_Rmin / Math.sqrt( (x[j]‐x[i])*(x[j]‐x[i])
                            +(y[j]‐y[i])*(y[j]‐y[i]));
E = E + LJ_Emin*(Math.pow(ratio,12) ‐ 2.0*Math.pow(ratio,6));
}
Page 1
Argon_MD.java
return E;
}
// "getNetForce" computes net force on atoms, including the Lennard
‐Jones
 term
void getNetForce(double x[],  double y[],     // current position of atoms
                 double Fx[], double Fy[]){   // returned net forces 
int i,j;              // particle index
double LJ_Fx, LJ_Fy;  // Lennard
‐Jones
 force on particle
double F;             // temporary variable
double ratio;         // LJ_Rmin to atomic distance ratio
// net force for the ith
 particle
for(i=0; i<nParticle; i++){
// set summation to zero
Fx[i] = 0.0;
Fy[i] = 0.0;
// loop thru the "other" particle
for(j=0; j<nParticle; j++)
if(i!=j){
// compute LJ_Rmin to atomic distance ratio
ratio = LJ_Rmin / Math.sqrt( (x[j]‐x[i])*(x[j]‐x[i]) 
                           + (y[j]‐y[i])*(y[j]‐y[i]) );
// compute "scaled" manitude
 of force
F = 12.0 * LJ_Emin / LJ_Rmin / LJ_Rmin 
         * (Math.pow(ratio,14) ‐ Math.pow(ratio,8));
// break into x and y component
LJ_Fx = F * (x[i]‐x[j]);
LJ_Fy = F * (y[i]‐y[j]);
// sum to get net net force
Fx[i] = Fx[i] + LJ_Fx;
Fy[i] = Fy[i] + LJ_Fy;
}
}
}
// "setInitPosition" sets initial position to a 
// uniform grid with spacing determined by "a" parameter
void setInitPosition(double x0, double y0,    // lower left corner coordinate
                     double a,                // spacing [nano
‐meter]
                     double x[], double y[]){ // returned coordinates
int maxN;   // number of particle in a single row
int s=0;    // particle counter
int i,j;    // row and column index
// computer number of particles in a single row
maxN = (int)Math.ceil(Math.sqrt(nParticle));
// loop thru all rows and columns
for(i=0; i<maxN; i++)
for(j=0; j<maxN; j++){
x[s] = x0 + i*a;
y[s] = y0 + j*a;
Page 2
Argon_MD.java
s    = s  + 1;
if(s >= nParticle) return;
}
}
// "setInitVelocity" sets initial velocity with every atom 
// having speed "v0" but more or less random direction
void setInitVelocity(double v0,                 // speed [nano
‐meter/pico
‐sec

                     double vx[], double vy[]){ // returned velocities
int i;   // particle index
// loop thru all particles
for(i=0; i<nParticle; i++){
vx[i] = v0*Math.cos(i*Math.PI/2.0);
vy[i] = v0*Math.sin(i*Math.PI/2.0); 
}
}
// "wallCollision" adjusts velocity of the particle to bounce of walls
void wallCollision(double x[],  double y[],    // current x and y coordinates
                   double vx[], double vy[]){  // adjected
 velocities
int i;  // particle index
// loop thru all particles
for(i=0; i<nParticle; i++){
if(x[i] < 0)         vx[i] =  Math.abs(vx[i]); // left wall
if(y[i] < 0)         vy[i] =  Math.abs(vy[i]); // bottom wall
if(x[i] > boxLength) vx[i] = ‐Math.abs(vx[i]); // right wall
if(y[i] > boxLength) vy[i] = ‐Math.abs(vy[i]); // top wall
}
}
// "histogram" generates the histogram of "data" within the range [a,b]
// split into "n" subintervals. 
void histogram(double data[], long nPoint,  // data and number of points
               double a, double b,          // min
 and max range  
               int n,                       // number of subintervals
               double bin[]){               // returned frequencies
double delta=(b‐a)/n;    // subinterval width
int i;                   // data point index
int k;                   // sub interval index
// reset bin to zero
for(k=0; k<n; k++) bin[k] = 0.0;
for(i=0; i < nPoint; i++){
// compute sub interval the data point belongs to
k = (int)Math.floor((data[i]‐a)/delta);
// add to the bin
if(k>=0 && k<n) bin[k] = bin[k]+1.0;
}
}
// "pairDistanceSet" generates all possible pair distance between particles.
// Hence, the set contains nParticle*(nParticle‐1)/2 number of possible pairs.
void pairDistanceSet(double x[], double y[],  // particle coordinates 
Page 3
Argon_MD.java
                     double r[]){             // returned possible distances
int s;   // pair index
int i,j;  // particle index
s=0;
for(i=0; i<nParticle; i++)
for(j=0; j<nParticle; j++)
if(i>j){
r[s] = Math.sqrt( (x[i]‐x[j])*(x[i]‐x[j]) + (y[i]‐y[j])*(y[i]‐y[j]) );
s    = s+1;
}
}
/////////////////////////////////
// enter constructor code here
/////////////////////////////////
void genesis(){
 x = new double[nParticle];
 y = new double[nParticle];
vx = new double[nParticle];
vy = new double[nParticle];
ax = new double[nParticle];
ay = new double[nParticle];
// RDF profiles
nRDFInterval = 100;                   // number of intervals
aRDF = 0.0; bRDF = 1.5;               // probe from [0,1.5] nm
rPointRDF = new double[nRDFInterval]; // mid point of radius
solidRDF  = new double[nRDFInterval]; // solid phase v0=0.05
liquidRDF = new double[nRDFInterval]; // liquid phase v0=0.2
gasRDF    = new double[nRDFInterval]; // gas phase v0=0.5
for(int i=0; i<nRDFInterval; i++){
rPointRDF[i] = (i+1/2.0)*(bRDF‐aRDF)/nRDFInterval+aRDF;
solidRDF[i] = 0.0;
liquidRDF[i] = 0.0;
gasRDF[i] = 0.0;
}
}
/////////////////////////////
// enter code to reset here
/////////////////////////////
void reset(){
// set initial positions
setInitPosition(2,2,LJ_Rmin, x, y);
// set initial velocities: try 0.05(solid), 0.2(liquid), 0.5(gas)
setInitVelocity(v0, vx, vy);
// reset time counter
nTime = 0;
}
////////////////////////////////
// enter code to update here
////////////////////////////////
void update(){
int n,i;
Page 4
Argon_MD.java
for(n=0; n<10; n++){
// compute acceleration
getNetForce(x, y, ax, ay);
for(i=0; i<nParticle; i++){ ax[i] = ax[i]/mass; ay[i] = ay[i]/mass; }
// Euler
‐Cromer
 update
for(i=0; i<nParticle; i++){
vx[i] = vx[i] + h * ax[i];
vy[i] = vy[i] + h * ay[i];
 x[i] =  x[i] + h * vx[i];
 y[i] =  y[i] + h * vy[i];
}
// bounce off the walls
wallCollision(x, y, vx, vy);
// RDF profiles
double rSet[] = new double[(nParticle)*(nParticle‐1)/2];
double bin[] = new double[nRDFInterval];
pairDistanceSet(x, y, rSet);
histogram(rSet, (nParticle)*(nParticle‐1)/2,
  aRDF, bRDF, nRDFInterval, bin);
for(i=0; i<nRDFInterval; i++){
if(Math.abs(v0‐0.05)<0.001)
solidRDF[i] = (nTime*solidRDF[i]   + bin[i])/(nTime+1.0);
if(Math.abs(v0‐0.2)<0.001) 
liquidRDF[i] = (nTime*liquidRDF[i] + bin[i])/(nTime+1.0);
if(Math.abs(v0‐0.5)<0.001) 
gasRDF[i]    = (nTime*gasRDF[i]    + bin[i])/(nTime+1.0);
}
// time step counter
nTime = nTime + 1;
}
}
public void init(){
genesis();
// allocate three sub children windows
// 1. paramPanel 2. cntrlPanel 3. dispPanel
// at the locations shown below
// |‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|
// |         paramPanel            |
// |‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|
// |                               |
// |                               |
// |         dispPanel             |
// |                               |
// |                               |
// |‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|
// |         cntrPanel             |
// |‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|
JPanel paramPanel = new JPanel();
JPanel cntrPanel = new JPanel();
JPanel dispPanel = new JPanel();
///////////////////////
Page 5
Argon_MD.java
// handle 1. paramPanel
///////////////////////
// initial speed selector
String[] initSpeedSet = {"0.05","0.2","0.5"};
SpinnerListModel initSpeedModel = new SpinnerListModel(initSpeedSet);
initSpeedSpinner = new JSpinner(initSpeedModel);
paramPanel.add(new JLabel("Initial Speed [nm/ps]"));
paramPanel.add(initSpeedSpinner);
///////////////////////
// handle 2. cntrPanel
///////////////////////
JButton playButton = new JButton("Play");
JButton pauseButton = new JButton("Pause");
JButton resetButton  = new JButton("Reset");
playButton.addActionListener(new PlayListener());
pauseButton.addActionListener(new PauseListener());
resetButton.addActionListener(new ResetListener());
cntrPanel.setBackground(Color.darkGray);
cntrPanel.add(playButton);
cntrPanel.add(pauseButton);
cntrPanel.add(resetButton);
//////////////////////
// handle 3. dispPanel
//////////////////////
// This middle panel is divided into 2 columns as shown
// 
// |‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|
// |                   |                 |
// | animationMDPanel  |  analysisPanel  |
// |                   |                 |
// |‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐|
dispPanel.setBackground(Color.white);
// 3.1 declare animation screen
animationMDPanel = new AnimationMDPanel(nParticle, 5, Color.blue,
                        0, boxLength, 0, boxLength,
                        200, 200);
dispPanel.add(animationMDPanel);
// 3.2 declare analysisPanel
JPanel analysisPanel = new JPanel();
analysisPanel.setLayout(new BoxLayout(analysisPanel, BoxLayout.Y_AXIS));
analysisPanel.setBackground(Color.white);
// 3.2.1 basis simulation information
JTextArea simInfo = new JTextArea(2,35);
simInfo.setFont(new Font("Courier New", Font.BOLD, 12));
simInfo.setText("Molecular Dynamics of Argon Atoms\n");
simInfo.append("Euler‐Cromer using h = 0.01 ps");
analysisPanel.add(simInfo);
// 3.2.2 kinetic
, total, potential energy graphs
energyGraph = new GraphPanel(3, 100, 200, 100);
energyGraph.setColor(0,Color.blue);
Page 6
Argon_MD.java
energyGraph.setColor(1,Color.red);
energyGraph.setColor(2,Color.green);
analysisPanel.add(energyGraph);
// 3.2.3 energy graph labels
JTextArea energyInfo = new JTextArea(2,35);
energyInfo.setFont(new Font("Courier New", Font.ITALIC, 12));
energyInfo.setText("Kinetic (blue), Total (red), and\n");
energyInfo.append("Potential Energy (green) [kJ/mol]");
analysisPanel.add(energyInfo);
// 3.2.4 RDF graph screen
RDFPanel = new GraphPanel(3, nRDFInterval, 200, 100);
RDFPanel.setColor(0,Color.blue);
RDFPanel.setColor(1,Color.red);
RDFPanel.setColor(2,Color.green);
analysisPanel.add(RDFPanel);
// 3.2.5 RDF labels
JTextArea RDFInfo = new JTextArea(2,35);
RDFInfo.setFont(new Font("Courier New", Font.ITALIC, 12));
RDFInfo.setText("Simple Radial Distribution Function\n");
RDFInfo.append("Solid (blue), Liquid (red), Gas (green)");
analysisPanel.add(RDFInfo);
// done creating analysisPanel
dispPanel.add(analysisPanel);
// add three sub children windows to frame
getContentPane().add(BorderLayout.CENTER,dispPanel);
getContentPane().add(BorderLayout.NORTH,paramPanel);
getContentPane().add(BorderLayout.SOUTH,cntrPanel);
// compute delay between frame
String str = getParameter("fps");
int fps = (str != null) ? Integer.parseInt(str): 100;
delay = (fps > 0)? 1000/fps : 100;
}
// This method is called when the applet
 becomes visible on screen.
// Create a thread and start it.
public void start(){
animator = new Thread(this);
animator.start();
}
// This method is called by the thread that was created in the
// start method. It does the main animation
public void run() {
// main animation loop
while(Thread.currentThread() == animator){
// load initial speed from GUI
v0 = new Float(initSpeedSpinner.getValue().toString());
// check current status and responds
switch(simStatus){
// we are at "reset" state
case 0:
Page 7
Argon_MD.java
reset();
animationMDPanel.setParticlesCoord(nTime*h, x, y);
initSpeedSpinner.setEnabled(true);  
break;
// we are at "update" state
case 1: 
update(); 
initSpeedSpinner.setEnabled(false);
animationMDPanel.setParticlesCoord(nTime*h, x, y);
// energy graph
double K = KE(vx,vy);
double U = LJ_Energy(x,y);
if(nTime%100==0){
energyGraph.pushValue(0, K);
energyGraph.pushValue(1, K+U);
energyGraph.pushValue(2, U);
}
// RDF graph
RDFPanel.setValues(0,nRDFInterval,rPointRDF,solidRDF);
RDFPanel.setValues(1,nRDFInterval,rPointRDF,liquidRDF);
RDFPanel.setValues(2,nRDFInterval,rPointRDF,gasRDF);
break;
// we are at "pause" state
default:
initSpeedSpinner.setEnabled(false);
}
repaint();
try{
Thread.sleep(delay);
}catch(Exception ex){ }
}
}
class PlayListener implements ActionListener{
public void actionPerformed(ActionEvent ev){
simStatus=1;
}
}
class PauseListener implements ActionListener{
public void actionPerformed(ActionEvent ev){
if(simStatus==1) simStatus=2;
}
}
class ResetListener implements ActionListener{
public void actionPerformed(ActionEvent ev){
simStatus=0;
}
}
}
Page 8
Argon_MD.java
///////////////////////////////////
// Utilities Graphics Objects
///////////////////////////////////
/////////////////////////////////////////////////////////
// "AnimationMDPanel" draws a set of particles on 
// screen to animate the motion of the system. Simply
// call "setParticleCoord" method and follow by "repaint()"
class AnimationMDPanel
 extends JPanel {
int nParticle;       // number of particles
double time;         // simulation time in ps
int radius;          // radius in pixel
Color color;         // color to draw particles
double x[], y[];     // coordinates
double xmin, xmax;   // range of x coordinates
double ymin, ymax;   // range of y coordinates
int width, height;   // preferred width and height in pixel 
public AnimationMDPanel(int nParticle,            // number of particles
                int radius,               // radius in pixel
                Color color,              // color of particles
                double xmin, double xmax, // range of x coordinates
                double ymin, double ymax, // range of y coordinates
                int width, int height){   // width and height in pixel 
// set internal variables
this.nParticle = nParticle;
this.radius = radius;
this.color = color;
this.xmin = xmin;
this.xmax = xmax;
this.ymin = ymin; 
this.ymax = ymax;
this.width = width;
this.height = height;
time = 0.0;
// allocate memory
x = new double[nParticle];
y = new double[nParticle];
setBorder(BorderFactory.createLineBorder(Color.black));
}
// set particles coordinates
public void setParticlesCoord(double time, double x[], double y[]){
this.time = time;
for(int i=0; i < nParticle; i++){
this.x[i] = x[i]; this.y[i] = y[i];
}
}
public Dimension getPreferredSize(){
return new Dimension(width,height);
}
public void paintComponent(Graphics g){
// erase screen with background color
Page 9
Argon_MD.java
g.setColor(Color.white);
g.fillRect(0,0,this.getWidth(),this.getHeight());
// print time
g.setColor(Color.black);
String timeString = String.format("t = %.1f ps", time);
int timeHeight = g.getFontMetrics().getHeight();
g.drawString(timeString,5,timeHeight);
double xScale = this.getWidth() /(xmax‐xmin);
double yScale = this.getHeight()/(ymax‐ymin); 
// draw particles
g.setColor(color);
for(int i=0; i < nParticle; i++)
g.fillOval((int)(x[i]*xScale)‐radius,
   (int)(y[i]*yScale)‐radius, 
   2*radius, 2*radius);
}
} // close AnimationMDPanel
/////////////////////////////////////////////////////////
// "GraphPanel" is an object that makes plotting graphs
// much easier. We can plot multiple data sets on the
// same graph screen, and change their color. Ths
// plotting range such as "ymin
" and "ymax
" are determined
// automatically as the data set is loaded into the object.
/////////////////////////////////////////////////////////
class GraphPanel
 extends JPanel {
int nData;                  // number of data set
int maxPoint;               // maximum number of points
int nPoint[];               // number of points in each set
double x[], y[];            // data
Color color[];              // color for each set
double xmin,xmax,ymin,ymax; // range of graph screen
int width, height;          // preferred width and height
public GraphPanel(int nData, int maxPoint, int width, int height){
this.nData = nData;
this.maxPoint = maxPoint;
nPoint = new int[nData];
x = new double[nData*maxPoint];
y = new double[nData*maxPoint];
color = new Color[nData];
this.width = width;
this.height = height;
xmin=0; xmax=0; ymin=0; ymax=0;
for(int i=0; i<nData; i++){
nPoint[i] = 0;
for(int j=0; j<maxPoint; j++){
x[i*maxPoint+j] = xmin;
y[i*maxPoint+j] = ymin;
}
}
setBorder(BorderFactory.createLineBorder(Color.black));
}
Page 10
Argon_MD.java
public void setColor(int i, Color color){
this.color[i] = color;
}
// push ynew
 value next to the last points
public void pushValue(int i, double ynew){
double xnew;
if(nPoint[i]==0) xnew = 0.0; else xnew = x[i*maxPoint+nPoint[i]‐1]+1.0;
if(nPoint[i]==maxPoint){
// scroll first (maxPoint‐1) to the left
for(int j=0; j<(nPoint[i]‐1); j++){ 
x[i*maxPoint+j] = x[i*maxPoint+j+1]; 
y[i*maxPoint+j] = y[i*maxPoint+j+1];
}
// put new value at the last position
x[i*maxPoint+nPoint[i]‐1] = xnew;
y[i*maxPoint+nPoint[i]‐1] = ynew;
}else{
x[i*maxPoint+nPoint[i]] = xnew;
y[i*maxPoint+nPoint[i]] = ynew;
nPoint[i] = nPoint[i] + 1;
}
}
// load a set of data points in
public void setValues(int i, int nPoint, double x[], double y[]){
this.nPoint[i] = 0;
for(int j=0; j<nPoint && j<maxPoint; j++){
this.x[i*maxPoint+j] = x[j]; this.y[i*maxPoint+j] = y[j];
this.nPoint[i] = this.nPoint[i] + 1;
}
}
public Dimension getPreferredSize(){
return new Dimension(width,height);
}
public void paintComponent(Graphics g){
// erase previous screen
g.setColor(Color.white);
g.fillRect(0,0,this.getWidth(),this.getHeight());
// update xmin
 xmax
 ymin
 ymax
xmin=x[0]; xmax=x[0]; ymin=y[0]; ymax=y[0];
for(int i=0; i<nData; i++)
for(int j=0; j<nPoint[i]; j++){
if(x[i*maxPoint+j]<xmin) xmin=x[i*maxPoint+j];
if(x[i*maxPoint+j]>xmax) xmax=x[i*maxPoint+j];
if(y[i*maxPoint+j]<ymin) ymin=y[i*maxPoint+j]; 
if(y[i*maxPoint+j]>ymax) ymax=y[i*maxPoint+j];
}
if((xmax‐xmin)==0) return;
if((ymax‐ymin)==0) return;
Page 11
Argon_MD.java
double xScale = this.getWidth() /(xmax‐xmin);
double yScale = this.getHeight()/(ymax‐ymin); 
// draw ticks
g.setColor(Color.black);
int ytick = this.getHeight()/3;
String tickString = String.format("(%+.2f)", (ymax‐ymin)/3+ymin);
g.drawLine(0, this.getHeight()‐ytick, 5, this.getHeight()‐ytick);
int tickHeight = g.getFontMetrics().getHeight();
g.drawString(tickString,6,this.getHeight()‐ytick+tickHeight/3);
ytick = this.getHeight()*2/3;
tickString = String.format("(%+.2f)", (ymax‐ymin)*2/3+ymin);
g.drawLine(0, this.getHeight()‐ytick, 5, this.getHeight()‐ytick);
tickHeight = g.getFontMetrics().getHeight();
g.drawString(tickString,6,this.getHeight()‐ytick+tickHeight/3);
// draw line connecting scattered data
for(int i=0; i< nData; i++){
g.setColor(color[i]);
for(int j=0; j<(nPoint[i]‐1); j++){
g.drawLine((int)((x[i*maxPoint+j]  ‐xmin)*xScale),
   this.getHeight()‐(int)((y[i*maxPoint+j]  ‐ymin)*yScale),
               (int)((x[i*maxPoint+j+1]‐xmin)*xScale),
               this.getHeight()‐(int)((y[i*maxPoint+j+1]‐ymin)*yScale));
}
}
}
} // close GraphPanel
Page 12