--- trunk/dp/dp.c	2003/07/09 22:14:06	586
+++ trunk/dp/dp.c	2003/07/10 15:42:41	587
@@ -1,7 +1,7 @@
 #include<stdio.h>
 #include<math.h>
 #include<stdlib.h>
-#define max_sites 10000
+#define max_sites 15000
 #define MYSEED 8973247
 
 int nsites, xlat, ylat;
@@ -15,94 +15,102 @@ int main()
 
 int main()
 {
-	void getmag();
-	double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb);
-	double toterg();
-	void adjust();
-	void mcmove();
-	void store(FILE *outFile);
+  void getmag();
+  double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb);
+  double toterg();
+  void adjust();
+  void mcmove();
+  void store(FILE *outFile);
+  
+  double twopi, nnd, myran, kb;
+  double ent;
+  int icycl, ncycl, imove, nmoves, nsamp;
+  int nx, ny, which, i, j;
+  
+  FILE *outFile;
+  outFile=fopen("dp_file","w");
+  
+  srand48(MYSEED);
+  
+  // These two parameters set the size of the system:
 
-	double twopi, nnd, myran, kb;
-	double ent;
-	int icycl, ncycl, imove, nmoves, nsamp;
-	int nx, ny, which, i, j;
-	
-	FILE *outFile;
-	outFile=fopen("dp_file","w");
-	
-	srand48(MYSEED);
-	
-	domsizex = 231.0;
-	domsizey = 231.0;
-	nnd = 7.0;
-	strength = 10.0;
-	ncycl = 100000000;
-	nmoves = 100;
-	nsamp = 200;
-	twopi = 2.0 * M_PI;
-	dtheta = 0.1;
-	kb = 0.0019872198;
-	beta = 1.0 / (kb * 300.0);
-	z0 = 0.0;
-	deltaz = 0.1;
-	deltaphi = 0.1;
-	kz = kb;
-	ktheta = kb;
-	theta0 = M_PI / 2.0;
+  nnd = 7.0;
+  ny = 100;
 
-	nx = (int) (2.5 * domsizey / sqrt(3.0) / nnd) + 1;
-	ny = (int) (domsizex / nnd) + 1;
+  // we want the domains to be roughly square:
 
-	which = 0; 
-	xlat = 0;
-	for(i=1;i<=nx;i++)
-	{
-		if ((i*sqrt(3.0)*nnd)<domsizex) 
-		{
-			xlat = xlat + 2;
-			ylat = 0;
-			for(j=1;j<=ny;j++)
-			{
-				if((j*nnd)<domsizey)
-				{
-					which = which + 1;
-					x[which] = i * sqrt(3.0) * nnd;
-					y[which] = j * nnd;
-					myran = drand48();
-					phi[which] = myran * twopi;
-					myran = drand48();
-					theta[which] = acos(2.0*myran-1.0);
-					myran = drand48();
-					z[which]=z0 + (2.0*myran-1.0)*deltaz;
-					mu[which] = strength;
-	
-					which = which + 1;
-					x[which] = nnd* sqrt(3.0) * (2.0*i+1.0) / 2.0;
-					y[which] = nnd* ( 2.0*j+1.0 ) / 2.0;
-					myran = drand48();
-					phi[which] = myran * twopi;
-					myran = drand48();
-					theta[which] = acos(2.0*myran-1.0);
-					myran = drand48();
-					z[which] = z0 + (2.0*myran-1.0) * deltaz;
-					mu[which] = strength;
-					ylat = ylat + 2;
-				}
-			}
-		}
-	}
+  nx = (int) ((double)ny / sqrt(3.0));
 
-	nsites= which;
 
-	en=toterg();
-	printf("\nthe initial energy is : \t%f\n\n",en);
+  // periodic box boundaries:
 
-	attemp = 0;
-	nacc = 0;
-	
-	adjust();
+  domsizex = sqrt(3.0) * nx * nnd;
+  domsizey = ny * nnd;
 
-	for(icycl=1;icycl<=ncycl;icycl++)
+  strength = 10.0;
+  ncycl = 1e7;
+  nmoves = 100;
+  nsamp = 1e4;
+  twopi = 2.0 * M_PI;
+  dtheta = 0.1;
+  kb = 0.0019872198;
+  beta = 1.0 / (kb * 300.0);
+  z0 = 0.0;
+  deltaz = 0.1;
+  deltaphi = 0.1;
+  kz = kb;
+  ktheta = kb;
+  theta0 = M_PI / 2.0;
+    
+  which = 0; 
+  xlat = 0;
+
+  for(i=0; i < nx; i++) {
+
+    xlat = xlat + 2;
+    ylat = 0;
+    
+    for(j=0; j<ny; j++) {     
+      
+      which = which + 1;
+      x[which] = i * sqrt(3.0) * nnd;
+      y[which] = j * nnd;
+      myran = drand48();
+      phi[which] = myran * twopi;
+      myran = drand48();
+      theta[which] = acos(2.0*myran-1.0);
+      myran = drand48();
+      z[which]=z0 + (2.0*myran-1.0)*deltaz;
+      mu[which] = strength;
+      
+      which = which + 1;
+
+      printf("%d\n", which);
+
+      x[which] = nnd* sqrt(3.0) * (2.0*i + 1.0) / 2.0;
+      y[which] = nnd* ( 2.0*j+1.0 ) / 2.0;
+      myran = drand48();
+      phi[which] = myran * twopi;
+      myran = drand48();
+      theta[which] = acos(2.0*myran-1.0);
+      myran = drand48();
+      z[which] = z0 + (2.0*myran-1.0) * deltaz;
+      mu[which] = strength;
+      ylat = ylat + 2;
+    }
+  }
+
+  nsites= which;
+  
+  en=toterg();
+  printf("\nthe initial energy is : \t%f\n\n",en);
+  
+  attemp = 0;
+  nacc = 0;
+  
+  adjust();
+
+  for(icycl=1;icycl<=ncycl;icycl++)
 	{
 		for(imove=1;imove<=nmoves;imove++)
 		{
@@ -204,29 +212,29 @@ double eneri(double xi, double yi, double zi, double p
 	{
 		if(j!=i)
 		{
-			dx = x[j]-xi;
-			if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
-			dy = y[j]-yi;
-			if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
-			dz = z[j]-zi;
-
-			r2 = dx*dx+dy*dy+dz*dz;
-			r = sqrt(r2);
-			if(r<rcut)
-			{
-				r3 = r2*r;
-				r5 = r2*r3;
-
-				uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
-				uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
-				uzj = mu[j] * cos(theta[j]);
-				udotu = uxi*uxj + uyi*uyj + uzi*uzj;
-				rdotui = dx*uxi + dy*uyi + dz*uzi;
-				rdotuj = dx*uxj + dy*uyj + dz*uzj;
-
-				vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
-				eni = eni + vij;
-			}
+		  dx = x[j]-xi;			
+		  if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
+		  dy = y[j]-yi;
+		  if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
+		  dz = z[j]-zi;
+		  
+		  r2 = dx*dx+dy*dy+dz*dz;
+		  r = sqrt(r2);
+		  if(r<rcut)
+		    {
+		      r3 = r2*r;
+		      r5 = r2*r3;
+		      
+		      uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
+		      uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
+		      uzj = mu[j] * cos(theta[j]);
+		      udotu = uxi*uxj + uyi*uyj + uzi*uzj;
+		      rdotui = dx*uxi + dy*uyi + dz*uzi;
+		      rdotuj = dx*uxj + dy*uyj + dz*uzj;
+		      
+		      vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
+		      eni = eni + vij;
+		    }
 		}
 	}
 	vint = 0.5 * kz * (zi-z0) * (zi-z0) + 0.5 * ktheta * (thetai-theta0) * (thetai-theta0);