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 581 by xsun, Wed Jul 9 15:14:05 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 < #define max_sites 10000
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  
18   int nsites, xlat, ylat;
# Line 13 | Line 24 | int attemp, nacc;
24   int attemp, nacc;
25  
26  
27 < int main()
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();
34 <        void store(FILE *outFile);
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(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  
46 <        double twopi, nnd, myran, kb;
47 <        double ent;
48 <        int icycl, ncycl, imove, nmoves, nsamp;
49 <        int nx, ny, which, i, j;
29 <        
30 <        FILE *outFile;
31 <        outFile=fopen("dp_file","w");
32 <        
33 <        srand48(MYSEED);
34 <        
35 <        domsizex = 231.0;
36 <        domsizey = 231.0;
37 <        nnd = 7.0;
38 <        strength = 10.0;
39 <        ncycl = 100000000;
40 <        nmoves = 100;
41 <        nsamp = 200;
42 <        twopi = 2.0 * M_PI;
43 <        dtheta = 0.1;
44 <        kb = 0.0019872198;
45 <        beta = 1.0 / (kb * 300.0);
46 <        z0 = 0.0;
47 <        deltaz = 0.1;
48 <        deltaphi = 0.1;
49 <        kz = kb;
50 <        ktheta = kb;
51 <        theta0 = M_PI / 2.0;
46 >  struct        timeval         now_time_val;
47 >  struct        timezone        time_zone;
48 >  struct        tm              *now_tm;
49 >  time_t                        now;
50  
51 <        nx = (int) (2.5 * domsizey / sqrt(3.0) / nnd) + 1;
52 <        ny = (int) (domsizex / nnd) + 1;
51 >  VSLStreamStatePtr stream;
52 >    
53 >  FILE *outFile;
54 >  outFile=fopen("dp_file1","a");
55 >  
56 >  
57 >  //srand48(MYSEED);
58 >  //
59  
60 <        which = 0;
61 <        xlat = 0;
58 <        for(i=1;i<=nx;i++)
59 <        {
60 <                if ((i*sqrt(3.0)*nnd)<domsizex)
61 <                {
62 <                        xlat = xlat + 2;
63 <                        ylat = 0;
64 <                        for(j=1;j<=ny;j++)
65 <                        {
66 <                                if((j*nnd)<domsizey)
67 <                                {
68 <                                        which = which + 1;
69 <                                        x[which] = i * sqrt(3.0) * nnd;
70 <                                        y[which] = j * nnd;
71 <                                        myran = drand48();
72 <                                        phi[which] = myran * twopi;
73 <                                        myran = drand48();
74 <                                        theta[which] = acos(2.0*myran-1.0);
75 <                                        myran = drand48();
76 <                                        z[which]=z0 + (2.0*myran-1.0)*deltaz;
77 <                                        mu[which] = strength;
78 <        
79 <                                        which = which + 1;
80 <                                        x[which] = nnd* sqrt(3.0) * (2.0*i+1.0) / 2.0;
81 <                                        y[which] = nnd* ( 2.0*j+1.0 ) / 2.0;
82 <                                        myran = drand48();
83 <                                        phi[which] = myran * twopi;
84 <                                        myran = drand48();
85 <                                        theta[which] = acos(2.0*myran-1.0);
86 <                                        myran = drand48();
87 <                                        z[which] = z0 + (2.0*myran-1.0) * deltaz;
88 <                                        mu[which] = strength;
89 <                                        ylat = ylat + 2;
90 <                                }
91 <                        }
92 <                }
93 <        }
60 >  gettimeofday(&now_time_val, &time_zone);  /* get the time now  */
61 >  now = now_time_val.tv_sec;                /*  convert to epoch time */
62  
63 <        nsites= which;
63 >  mySeed = (int) now;
64  
65 <        en=toterg();
66 <        printf("\nthe initial energy is : \t%f\n\n",en);
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 >  
72 >  // we want the domains to be roughly square:
73 >  
74 >  nx = (int) ((double)ny / sqrt(3.0));
75 >  
76 >  
77 >  // periodic box boundaries:
78  
79 <        attemp = 0;
80 <        nacc = 0;
102 <        
103 <        adjust();
79 >  domsizex = sqrt(3.0) * nx * nnd;
80 >  domsizey = ny * nnd;
81  
82 <        for(icycl=1;icycl<=ncycl;icycl++)
82 >  strength = 10.0;
83 >  ncycl = 1e7;
84 >  nmoves = 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;
95 >  theta0 = M_PI / 2.0;
96 >  
97 >  which = 0;
98 >  xlat = 0;
99 >
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);
115 >    }
116 >      
117 >    while(!feof(inFile))
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);
128 >            fscanf(inFile,"%s",s);
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 >  } else {
149 >
150 >    for(i=0; i < nx; i++) {
151 >      
152 >      xlat = xlat + 2;
153 >      ylat = 0;
154 >      
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 >        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();
194 >  printf("\nthe initial energy is : \t%f\n\n",en);
195 >  
196 >  attemp = 0;
197 >  nacc = 0;
198 >  
199 >  adjust();
200 >
201 >  for(icycl=1;icycl<=ncycl;icycl++)
202          {
203                  for(imove=1;imove<=nmoves;imove++)
204                  {
205 <                        mcmove();
205 >                        mcmove(stream);
206                  }
207  
208                  if((icycl%nsamp) == 0)
# Line 134 | Line 230 | int main()
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 204 | Line 301 | double eneri(double xi, double yi, double zi, double p
301          {
302                  if(j!=i)
303                  {
304 <                        dx = x[j]-xi;
305 <                        if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
306 <                        dy = y[j]-yi;
307 <                        if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
308 <                        dz = z[j]-zi;
309 <
310 <                        r2 = dx*dx+dy*dy+dz*dz;
311 <                        r = sqrt(r2);
312 <                        if(r<rcut)
313 <                        {
314 <                                r3 = r2*r;
315 <                                r5 = r2*r3;
316 <
317 <                                uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
318 <                                uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
319 <                                uzj = mu[j] * cos(theta[j]);
320 <                                udotu = uxi*uxj + uyi*uyj + uzi*uzj;
321 <                                rdotui = dx*uxi + dy*uyi + dz*uzi;
322 <                                rdotuj = dx*uxj + dy*uyj + dz*uzj;
323 <
324 <                                vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
325 <                                eni = eni + vij;
326 <                        }
304 >                  dx = x[j]-xi;                
305 >                  if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
306 >                  dy = y[j]-yi;
307 >                  if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
308 >                  dz = z[j]-zi;
309 >                  
310 >                  r2 = dx*dx+dy*dy+dz*dz;
311 >                  r = sqrt(r2);
312 >                  if(r<rcut)
313 >                    {
314 >                      r3 = r2*r;
315 >                      r5 = r2*r3;
316 >                      
317 >                      uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
318 >                      uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
319 >                      uzj = mu[j] * cos(theta[j]);
320 >                      udotu = uxi*uxj + uyi*uyj + uzi*uzj;
321 >                      rdotui = dx*uxi + dy*uyi + dz*uzi;
322 >                      rdotuj = dx*uxj + dy*uyj + dz*uzj;
323 >                      
324 >                      vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
325 >                      eni = eni + vij;
326 >                    }
327                  }
328          }
329          vint = 0.5 * kz * (zi-z0) * (zi-z0) + 0.5 * ktheta * (thetai-theta0) * (thetai-theta0);
# Line 267 | 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 275 | 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