Home
|
Course Outline
|
Lectures
|
Homework
|
Files
The following applet, which is generated by the program PivotAlgorithm.java, demonstrates the pivot enrichment algorithm.
The code for this program is based on a fortran code saw.f by Andrea Pelissetto, which I found on the Web. It incorporates some changes and improvements on the simple description of the pivot algorithm on pages 406-407:
Here is the code which generates this applet:
// PivotAlgorithm.java
import comphys.Easy;
import comphys.graphics.*;
public class PivotAlgorithm extends Animation {
int L = 16; // size of lattice in each dimension
int N = 8; // number of bonds
int[] xMonomer, yMonomer; // polymer configuration
int[] xTrial, yTrial; // trial configuration
int polymers; // number of polymers generated
double RSquaredSum; // to compute mean end-to-end distance
int nPivot; // pivot point
int pivots; // successful pivot attempts
void initial () {
// allocate arrays if necessary
if (xMonomer == null || xMonomer.length != N + 1) {
xMonomer = new int[N + 1];
yMonomer = new int[N + 1];
xTrial = new int[N + 1];
yTrial = new int[N + 1];
}
// rod polymer with fixed start position at (N + 1, N + 1)
for (int i = 0; i <= N; i++) {
xMonomer[i] = N + 1 - i;
yMonomer[i] = N + 1;
}
// initialize hash table
initHash();
polymers = 1;
RSquaredSum = 0;
nPivot = -1;
pivots = 0;
data();
}
void data () {
double dx = xMonomer[0] - xMonomer[N];
double dy = yMonomer[0] - yMonomer[N];
RSquaredSum += dx * dx + dy * dy;
}
// Hash table with linear probing to check whether a site is occupied
// Start monomer in polymer is fixed at (N + 1, N + 1) which guarantees
// that the point (0, 0) is empty, so 0 signifies an empty table entry
int hashSize; // size of hash table
int[] xHash, yHash; // hash tables to store monomer (x, y) coords
int hash1, hash2; // hash function = x * hash1 + y * hash2 + 1
void initHash () {
// selected prime numbers
int[] primes = {47, 97, 499, 997, 4999, 9973, 49999};
// hash table size must be > (N - 1), choose it to be at
// least five times bigger to minimize key collisions
for (int i = 0; i < primes.length; i++) {
hashSize = primes[i];
if (hashSize > 5 * N)
break;
}
// allocate tables if necessary
if (xHash == null || xHash.length != hashSize) {
xHash = new int[hashSize];
yHash = new int[hashSize];
}
// zero tables
for (int i = 0; i < hashSize; i++)
xHash[i] = yHash[i] = 0;
// set hash function parameters
hash1 = 37;
hash2 = hash1 * hash1 + 1;
}
boolean occupied (int x, int y) {
// compute hash of (x, y)
int hash = (hash1 * x + hash2 * y) % hashSize;
// linear probe of hash table
boolean found = false;
while (true) {
int xStored = xHash[hash];
if (xStored == 0) {
// insert into (x, y) into table
xHash[hash] = x;
yHash[hash] = y;
break;
} else {
int yStored = yHash[hash];
if (x == xStored && y == yStored) {
found = true;
break;
} else
hash = ++hash % hashSize;
}
}
return found;
}
void pivot () {
// choose a random pivot monomer 0 < nPivot < N
nPivot = 1 + (int) (Math.random() * (N - 1));
int nMove = N - nPivot; // number of monomers to be moved
int nMin = nMove; // smaller of nPivot and nMove
if (nPivot < nMove)
nMin = nPivot;
int xPivot = xMonomer[nPivot];
int yPivot = yMonomer[nPivot];
xTrial[nPivot] = xPivot;
yTrial[nPivot] = yPivot;
// select random square lattice symmetry operation
int symmetryOperation = 1 + (int) (Math.random() * 7);
// set up transformation matrix
int xx = 0, xy = 0, yx = 0, yy = 0;
switch (symmetryOperation) {
case 1: // 90 deg clockwise rotation
xy = 1; yx = -1; break;
case 2: // 180 deg clockwise rotation
xx = -1; yy = -1; break;
case 3: // 270 deg clockwise rotation
xy = -1; yx = 1; break;
case 4: // reflect in x = -y
xy = -1; yx = -1; break;
case 5: // reflect in y = 0
xx = 1; yy = -1; break;
case 6: // reflect in x = y
xy = 1; yx = 1; break;
case 7: // reflect in x = 0
xx = -1; yy = 1; break;
default: // should not occur
xx = 1; yy = 1; break;
}
boolean moveAllowed = true;
// monomers symmetrically located about the pivot
for (int i = 1; i <= nMin; i++) {
// monomers beyond the pivot to be moved
int x = xMonomer[nPivot + i] - xPivot;
int y = yMonomer[nPivot + i] - yPivot;
int temp = xx * x + xy * y + xPivot;
y = yx * x + yy * y + yPivot;
x = temp;
if (occupied(x, y)) {
moveAllowed = false;
break;
}
xTrial[nPivot + i] = x;
yTrial[nPivot + i] = y;
// fixed monomers between start and pivot
x = xMonomer[nPivot - i];
y = yMonomer[nPivot - i];
if (occupied(x, y)) {
moveAllowed = false;
break;
}
xTrial[nPivot - i] = x;
yTrial[nPivot - i] = y;
}
// remaining monomers
if (nMove > nPivot && moveAllowed) {
// monomers beyond pivot to be moved
for (int i = 1; i <= nMove - nPivot; i++) {
int x = xMonomer[nPivot + nMin + i] - xPivot;
int y = yMonomer[nPivot + nMin + i] - yPivot;
int temp = xx * x + xy * y + xPivot;
y = yx * x + yy * y + yPivot;
x = temp;
if (occupied(x, y)) {
moveAllowed = false;
break;
}
xTrial[nPivot + nMin + i] = x;
yTrial[nPivot + nMin + i] = y;
}
} else if (moveAllowed) {
// fixed monomers between start and pivot
for (int i = 1; i <= nPivot - nMove; i++) {
int x = xMonomer[nPivot - nMin - i];
int y = yMonomer[nPivot - nMin - i];
if (occupied(x, y)) {
moveAllowed = false;
break;
}
xTrial[nPivot - nMin - i] = x;
yTrial[nPivot - nMin - i] = y;
}
}
if (moveAllowed) {
// replace polymer by trial polymer
for (int i = 0; i <= N; i++) {
xMonomer[i] = xTrial[i];
yMonomer[i] = yTrial[i];
}
++pivots;
}
++polymers;
data();
// clear hash table
for (int i = 0; i < hashSize; i++)
xHash[i] = 0;
}
class Movie extends Plot {
Movie () {
setSize(400, 400);
setBackground("white");
}
public void paint () {
clear();
double r = 0.15;
double x = xMonomer[0];
double y = yMonomer[0];
setWindow(N + 0.5 - L/2, N + 1.5 + L/2,
N + 0.5 - L/2, N + 1.5 + L/2);
for (int i = 0; i < N; i++) {
double x1 = xMonomer[i + 1];
double y1 = yMonomer[i + 1];
setColor("green");
plotLine(x, y, x1, y1);
if (i == nPivot)
setColor("red");
else
setColor("black");
floodCircle(x - r, x + r, y - r, y + r);
x = x1;
y = y1;
}
setColor("black");
floodCircle(x - r, x + r, y - r, y + r);
}
}
class Output extends Plot {
Output () {
setSize(400, 30);
setBackground("blue");
}
public void paint () {
clear();
setWindow(0, 100, 0, 3);
setColor("white");
plotString("Polymer " + polymers, 5, 1);
plotString(" = " + Easy.format(RSquaredSum / polymers, 4),
40, 1);
plotStringLeft("Pivots "
+ Easy.format(pivots * 100.0 / polymers, 3)
+ " %", 95, 1);
}
}
Movie movie;
Output output;
Reader NReader, LReader;
public void init () {
initial();
add(movie = new Movie());
add(output = new Output());
add(NReader = new Reader("N =", N));
add(LReader = new Reader("L =", L));
addControlPanel();
}
public void step () {
pivot();
movie.repaint();
output.repaint();
}
public void reset () {
N = NReader.readInt();
L = LReader.readInt();
initial();
movie.repaint();
output.repaint();
}
public static void main (String[] args) {
PivotAlgorithm pivotAlgorithm = new PivotAlgorithm();
pivotAlgorithm.frame("Pivot Algorithm Simulation", 430, 600);
}
}
Questions or comments:
phygons@acsu.buffalo.edu