ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/dp/dp.c
(Generate patch)

Comparing trunk/dp/dp.c (file contents):
Revision 665 by xsun, Mon Aug 4 23:40:54 2003 UTC vs.
Revision 692 by xsun, Wed Aug 13 16:43:53 2003 UTC

# Line 1 | Line 1
1   #include<stdio.h>
2 + #include<sys/time.h>
3 + #include<unistd.h>
4 + #include<time.h>
5   #include<math.h>
6   #include<stdlib.h>
7   #include<string.h>
8 + #include "mkl_vsl.h"
9 +
10 + #define SEED    1
11 + #define BRNG    VSL_BRNG_MCG31
12 + #define METHOD  0
13 + #define N       1
14 +
15   #define max_sites 15000
16   #define MYSEED 8973247
17  
# Line 14 | Line 24 | int attemp, nacc;
24   int attemp, nacc;
25  
26  
27 < int main(argc, argv)
18 < int argc;
19 < char *argv[];
27 > int main(int argc, char *argv[])
28   {
29    void getmag();
30    double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb);
31    double toterg();
32    void adjust();
33 <  void mcmove();
33 >  void mcmove(VSLStreamStatePtr stream);
34    void store(FILE *outFile);
35    
36    double twopi, nnd, myran, kb;
37    double ent;
38    int icycl, ncycl, imove, nmoves, nsamp;
39    int nx, ny, which, i, j;
40 +  int readfile, mySeed;
41    
42    FILE *inFile;
43    char *endptr, s[100000];
44    double ux, uy, uz;
45 <  
45 >
46 >  struct        timeval         now_time_val;
47 >  struct        timezone        time_zone;
48 >  struct        tm              *now_tm;
49 >  time_t                        now;
50 >
51 >  VSLStreamStatePtr stream;
52 >    
53    FILE *outFile;
54    outFile=fopen("dp_file1","a");
55    
40 –  srand48(MYSEED);
56    
57 <  // These two parameters set the size of the system:
57 >  //srand48(MYSEED);
58 >  //
59  
60 +  gettimeofday(&now_time_val, &time_zone);  /* get the time now  */
61 +  now = now_time_val.tv_sec;                /*  convert to epoch time */
62 +
63 +  mySeed = (int) now;
64 +
65 +  vslNewStream(&stream, BRNG, mySeed);
66 +  
67 +  // These two parameters set the size of the system:
68 +  
69    nnd = 7.0;
70    ny = 100;
71 <
71 >  
72    // we want the domains to be roughly square:
73 <
73 >  
74    nx = (int) ((double)ny / sqrt(3.0));
75 <
76 <
75 >  
76 >  
77    // periodic box boundaries:
78  
79    domsizex = sqrt(3.0) * nx * nnd;
# Line 57 | Line 82 | char *argv[];
82    strength = 10.0;
83    ncycl = 1e7;
84    nmoves = 100;
85 <  nsamp = 100;
85 >  nsamp = 1e4;
86    twopi = 2.0 * M_PI;
87    dtheta = 0.1;
88    kb = 0.0019872198;
89    beta = 1.0 / (kb * 300.0);
90    z0 = 0.0;
91    deltaz = 0.1;
92 <  deltaphi = 0.1;
93 <  kz = kb;
94 <  ktheta = kb;
92 >  deltaphi = 0.1;
93 >  kz = kb;
94 >  ktheta = kb;
95    theta0 = M_PI / 2.0;
96    
97    which = 0;
98    xlat = 0;
99  
100 <  if( strcmp(argv[1], "-resume") == 0) {
100 >  if (argc > 1) {
101 >    if( strcmp(argv[1], "-resume") == 0)
102 >      readfile = 1;
103 >    else
104 >      readfile = 0;
105 >  } else
106 >    readfile = 0;
107 >
108      
109 +  if (readfile == 1) {
110 +    
111      inFile = fopen("./l_con","r");
112      if (inFile == NULL){
113 <            printf("Error opening file\n");
114 <            exit(-1);
113 >      printf("Error opening file\n");
114 >      exit(-1);
115      }
116 <    
116 >      
117      while(!feof(inFile))
118 <    {
119 <            fgets(s, 500, inFile);
120 <            nsites = atoi(s);
118 >      {
119 >        fgets(s, 500, inFile);
120 >        nsites = atoi(s);
121 >        fscanf(inFile,"%s",s);
122 >        xlat= atoi(s);
123 >        fscanf(inFile,"%s",s);
124 >        ylat= atoi(s);
125 >        for(i=1; i<=nsites; i++)
126 >          {
127              fscanf(inFile,"%s",s);
88 –            xlat= atoi(s);
128              fscanf(inFile,"%s",s);
129 <            ylat= atoi(s);
130 <            for(i=1; i<=nsites; i++)
131 <            {
132 <                    fscanf(inFile,"%s",s);
133 <                    fscanf(inFile,"%s",s);
134 <                    x[i] = strtod(s, &endptr);
135 <                    fscanf(inFile,"%s",s);
136 <                    y[i] = strtod(s, &endptr);
137 <                    fscanf(inFile,"%s",s);
138 <                    z[i] = strtod(s, &endptr);
139 <                    fscanf(inFile,"%s",s);
140 <                    phi[i] = strtod(s, &endptr);
141 <                    fscanf(inFile,"%s",s);
142 <                    ux = strtod(s, &endptr);
143 <                    fscanf(inFile,"%s",s);
144 <                    uy = strtod(s, &endptr);
145 <                    fscanf(inFile,"%s",s);
107 <                    uz = strtod(s, &endptr);
108 <                    theta[i] = acos(uz);
109 <                    mu[i] = strength;
110 <            }
111 <    }
129 >            x[i] = strtod(s, &endptr);
130 >            fscanf(inFile,"%s",s);
131 >            y[i] = strtod(s, &endptr);
132 >            fscanf(inFile,"%s",s);
133 >            z[i] = strtod(s, &endptr);
134 >            fscanf(inFile,"%s",s);
135 >            phi[i] = strtod(s, &endptr);
136 >            fscanf(inFile,"%s",s);
137 >            ux = strtod(s, &endptr);
138 >            fscanf(inFile,"%s",s);
139 >            uy = strtod(s, &endptr);
140 >            fscanf(inFile,"%s",s);
141 >            uz = strtod(s, &endptr);
142 >            theta[i] = acos(uz);
143 >            mu[i] = strength;
144 >          }
145 >      }
146      fclose(inFile);
147      printf("%f\t%f\t%f\t%f\t%f\n", x[nsites], y[nsites], z[nsites], phi[nsites], theta[nsites]);
148 <  }
148 >  } else {
149  
150 <  else{
117 <  for(i=0; i < nx; i++) {
118 <
119 <    xlat = xlat + 2;
120 <    ylat = 0;
121 <    
122 <    for(j=0; j<ny; j++) {    
150 >    for(i=0; i < nx; i++) {
151        
152 <      which = which + 1;
153 <      x[which] = i * sqrt(3.0) * nnd;
126 <      y[which] = j * nnd;
127 <      myran = drand48();
128 <      phi[which] = myran * twopi;
129 <      myran = drand48();
130 <      theta[which] = acos(2.0*myran-1.0);
131 <      myran = drand48();
132 <      z[which]=z0 + (2.0*myran-1.0)*deltaz;
133 <      mu[which] = strength;
152 >      xlat = xlat + 2;
153 >      ylat = 0;
154        
155 <      which = which + 1;
156 <
157 <      printf("%d\n", which);
158 <
159 <      x[which] = nnd* sqrt(3.0) * (2.0*i + 1.0) / 2.0;
160 <      y[which] = nnd* ( 2.0*j+1.0 ) / 2.0;
161 <      myran = drand48();
162 <      phi[which] = myran * twopi;
163 <      myran = drand48();
164 <      theta[which] = acos(2.0*myran-1.0);
165 <      myran = drand48();
166 <      z[which] = z0 + (2.0*myran-1.0) * deltaz;
167 <      mu[which] = strength;
168 <      ylat = ylat + 2;
169 <    }
170 <  }
155 >      for(j=0; j<ny; j++) {    
156 >        
157 >        which = which + 1;
158 >        x[which] = i * sqrt(3.0) * nnd;
159 >        y[which] = j * nnd;
160 >        
161 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
162 >        //myran = drand48();
163 >        phi[which] = myran * twopi;
164 >        //myran = drand48();
165 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
166 >        theta[which] = acos(2.0*myran-1.0);
167 >        //myran = drand48();
168 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
169 >        z[which]=z0 + (2.0*myran-1.0)*deltaz;
170 >        mu[which] = strength;
171 >        
172 >        which = which + 1;
173  
174 <  nsites= which;
174 >        x[which] = nnd* sqrt(3.0) * (2.0*i + 1.0) / 2.0;
175 >        y[which] = nnd* ( 2.0*j+1.0 ) / 2.0;
176 >        //myran = drand48();
177 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
178 >        phi[which] = myran * twopi;
179 >        //myran = drand48();
180 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
181 >        theta[which] = acos(2.0*myran-1.0);
182 >        //myran = drand48();
183 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
184 >        z[which] = z0 + (2.0*myran-1.0) * deltaz;
185 >        mu[which] = strength;
186 >        ylat = ylat + 2;
187 >      }
188 >    }
189 >    
190 >    nsites= which;
191    }
192    
193    en=toterg();
# Line 164 | Line 202 | char *argv[];
202          {
203                  for(imove=1;imove<=nmoves;imove++)
204                  {
205 <                        mcmove();
205 >                        mcmove(stream);
206                  }
207  
208                  if((icycl%nsamp) == 0)
# Line 192 | Line 230 | char *argv[];
230                  printf("\nTotal energy end of simulation : %f\nrunning energy : %f\ndifference : %f\n\n", ent, en, ent-en);
231          }
232          fclose(outFile);
233 +        vslDeleteStream( &stream );
234          return 0;
235   }
236  
# Line 325 | Line 364 | void adjust()
364          
365   }
366  
367 < void mcmove()
367 > void mcmove( VSLStreamStatePtr stream )
368   {
369          int o;
370          double eno, zn, phin, thetan, enn, myran;
# Line 333 | Line 372 | void mcmove()
372  
373          attemp = attemp + 1;
374          jb = 1;
375 <        myran = drand48();
375 >        //myran = drand48();
376 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
377          o = (int)(nsites*myran) + 1;
378          eno=eneri(x[o], y[o], z[o], phi[o], theta[o], o, jb);
379 <        myran = drand48();
379 >        //myran = drand48();
380 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
381          zn = z[o] + (2.0*myran-1.0) * deltaz;
382 <        myran = drand48();
382 >        //myran = drand48();
383 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
384          phin = phi[o] + (2.0*myran-1.0) * deltaphi;
385 <        myran = drand48();
385 >        //myran = drand48();
386 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
387          thetan = theta[o] + (2.0*myran-1.0)*dtheta;
388          enn=eneri(x[o], y[o], zn, phin, thetan, o, jb);
389 <        myran = drand48();
389 >        //myran = drand48();
390 >        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
391          if(myran<=(exp(-beta*(enn-eno))))
392          {
393                  nacc = nacc + 1;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines