--- trunk/dp/dp.c	2003/08/05 21:44:35	667
+++ trunk/dp/dp.c	2003/08/13 16:43:53	692
@@ -1,7 +1,17 @@
 #include<stdio.h>
+#include<sys/time.h>
+#include<unistd.h>
+#include<time.h>
 #include<math.h>
 #include<stdlib.h>
 #include<string.h>
+#include "mkl_vsl.h"
+
+#define SEED    1
+#define BRNG    VSL_BRNG_MCG31
+#define METHOD  0
+#define N       1
+
 #define max_sites 15000
 #define MYSEED 8973247
 
@@ -14,41 +24,56 @@ int main(argc, argv)
 int attemp, nacc;
 
 
-int main(argc, argv)
-int argc;
-char *argv[];
+int main(int argc, char *argv[])
 {
   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 mcmove(VSLStreamStatePtr stream);
   void store(FILE *outFile);
   
   double twopi, nnd, myran, kb;
   double ent;
   int icycl, ncycl, imove, nmoves, nsamp;
   int nx, ny, which, i, j;
+  int readfile, mySeed;
   
   FILE *inFile;
   char *endptr, s[100000];
   double ux, uy, uz;
-  
+
+  struct        timeval         now_time_val;
+  struct        timezone        time_zone;
+  struct        tm              *now_tm;
+  time_t                        now;
+
+  VSLStreamStatePtr stream;
+    
   FILE *outFile;
   outFile=fopen("dp_file1","a");
   
-  srand48(MYSEED);
   
-  // These two parameters set the size of the system:
+  //srand48(MYSEED);
+  //
 
+  gettimeofday(&now_time_val, &time_zone);  /* get the time now  */
+  now = now_time_val.tv_sec;                /*  convert to epoch time */
+
+  mySeed = (int) now;
+
+  vslNewStream(&stream, BRNG, mySeed);
+  
+  // These two parameters set the size of the system:
+  
   nnd = 7.0;
   ny = 100;
-
+  
   // we want the domains to be roughly square:
-
+  
   nx = (int) ((double)ny / sqrt(3.0));
-
-
+  
+  
   // periodic box boundaries:
 
   domsizex = sqrt(3.0) * nx * nnd;
@@ -64,92 +89,105 @@ char *argv[];
   beta = 1.0 / (kb * 300.0);
   z0 = 0.0;
   deltaz = 0.1;
-  deltaphi = 0.1;
-  kz = kb;
-  ktheta = kb;
+  deltaphi = 0.1; 
+  kz = kb; 
+  ktheta = kb; 
   theta0 = M_PI / 2.0;
    
   which = 0; 
   xlat = 0;
 
-  if( strcmp(argv[1], "-resume") == 0) {
+  if (argc > 1) {
+    if( strcmp(argv[1], "-resume") == 0) 
+      readfile = 1;
+    else
+      readfile = 0;
+  } else 
+    readfile = 0;
+
     
+  if (readfile == 1) {
+    
     inFile = fopen("./l_con","r");
     if (inFile == NULL){
-	    printf("Error opening file\n");
-	    exit(-1);
+      printf("Error opening file\n");
+      exit(-1);
     }
-    
+      
     while(!feof(inFile))
-    {
-	    fgets(s, 500, inFile);
-	    nsites = atoi(s);
+      {
+	fgets(s, 500, inFile);
+	nsites = atoi(s);
+	fscanf(inFile,"%s",s);
+	xlat= atoi(s);
+	fscanf(inFile,"%s",s);
+	ylat= atoi(s);
+	for(i=1; i<=nsites; i++)
+	  {
 	    fscanf(inFile,"%s",s);
-	    xlat= atoi(s);
 	    fscanf(inFile,"%s",s);
-	    ylat= atoi(s);
-	    for(i=1; i<=nsites; i++)
-	    {
-		    fscanf(inFile,"%s",s);
-		    fscanf(inFile,"%s",s);
-		    x[i] = strtod(s, &endptr);
-		    fscanf(inFile,"%s",s);
-		    y[i] = strtod(s, &endptr);
-		    fscanf(inFile,"%s",s);
-		    z[i] = strtod(s, &endptr);
-		    fscanf(inFile,"%s",s);
-		    phi[i] = strtod(s, &endptr);
-		    fscanf(inFile,"%s",s);
-		    ux = strtod(s, &endptr);
-		    fscanf(inFile,"%s",s);
-		    uy = strtod(s, &endptr);
-		    fscanf(inFile,"%s",s);
-		    uz = strtod(s, &endptr);
-		    theta[i] = acos(uz);
-		    mu[i] = strength;
-	    }
-    }
+	    x[i] = strtod(s, &endptr);
+	    fscanf(inFile,"%s",s);
+	    y[i] = strtod(s, &endptr);
+	    fscanf(inFile,"%s",s);
+	    z[i] = strtod(s, &endptr);
+	    fscanf(inFile,"%s",s);
+	    phi[i] = strtod(s, &endptr);
+	    fscanf(inFile,"%s",s);
+	    ux = strtod(s, &endptr);
+	    fscanf(inFile,"%s",s);
+	    uy = strtod(s, &endptr);
+	    fscanf(inFile,"%s",s);
+	    uz = strtod(s, &endptr);
+	    theta[i] = acos(uz);
+	    mu[i] = strength;
+	  }
+      }
     fclose(inFile);
     printf("%f\t%f\t%f\t%f\t%f\n", x[nsites], y[nsites], z[nsites], phi[nsites], theta[nsites]);
-  }
+  } else {
 
-  else{
-  for(i=0; i < nx; i++) {
-
-    xlat = xlat + 2;
-    ylat = 0;
-    
-    for(j=0; j<ny; j++) {     
+    for(i=0; i < nx; i++) {
       
-      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;
+      xlat = xlat + 2;
+      ylat = 0;
       
-      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;
-    }
-  }
+      for(j=0; j<ny; j++) {     
+	
+	which = which + 1;
+	x[which] = i * sqrt(3.0) * nnd;
+	y[which] = j * nnd;
+	
+	vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
+	//myran = drand48();
+	phi[which] = myran * twopi;
+	//myran = drand48();
+	vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
+	theta[which] = acos(2.0*myran-1.0);
+	//myran = drand48();
+	vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
+	z[which]=z0 + (2.0*myran-1.0)*deltaz;
+	mu[which] = strength;
+	
+	which = which + 1;
 
-  nsites= 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();
+	vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
+	phi[which] = myran * twopi;
+	//myran = drand48();
+	vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
+	theta[which] = acos(2.0*myran-1.0);
+	//myran = drand48();
+	vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
+	z[which] = z0 + (2.0*myran-1.0) * deltaz;
+	mu[which] = strength;
+	ylat = ylat + 2;
+      }
+    }
+    
+    nsites= which;
   }
   
   en=toterg();
@@ -164,7 +202,7 @@ char *argv[];
 	{
 		for(imove=1;imove<=nmoves;imove++)
 		{
-			mcmove();
+			mcmove(stream);
 		}
 
 		if((icycl%nsamp) == 0)
@@ -192,6 +230,7 @@ char *argv[];
 		printf("\nTotal energy end of simulation : %f\nrunning energy : %f\ndifference : %f\n\n", ent, en, ent-en);
 	}
 	fclose(outFile);
+	vslDeleteStream( &stream );
 	return 0;
 }
 
@@ -325,7 +364,7 @@ void mcmove()
 	
 }
 
-void mcmove()
+void mcmove( VSLStreamStatePtr stream )
 {
 	int o;
 	double eno, zn, phin, thetan, enn, myran;
@@ -333,17 +372,22 @@ void mcmove()
 
 	attemp = attemp + 1;
 	jb = 1;
-	myran = drand48();
+	//myran = drand48();
+        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
 	o = (int)(nsites*myran) + 1;
 	eno=eneri(x[o], y[o], z[o], phi[o], theta[o], o, jb);
-	myran = drand48();
+	//myran = drand48();
+        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
 	zn = z[o] + (2.0*myran-1.0) * deltaz;
-	myran = drand48();
+	//myran = drand48();
+        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
 	phin = phi[o] + (2.0*myran-1.0) * deltaphi;
-	myran = drand48();
+	//myran = drand48();
+        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
 	thetan = theta[o] + (2.0*myran-1.0)*dtheta;
 	enn=eneri(x[o], y[o], zn, phin, thetan, o, jb);
-	myran = drand48();
+	//myran = drand48();
+        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
 	if(myran<=(exp(-beta*(enn-eno))))
 	{
 		nacc = nacc + 1;