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 665 by xsun, Mon Aug 4 23:40:54 2003 UTC

# Line 1 | Line 1
1   #include<stdio.h>
2   #include<math.h>
3   #include<stdlib.h>
4 < #define max_sites 10000
4 > #include<string.h>
5 > #define max_sites 15000
6   #define MYSEED 8973247
7  
8   int nsites, xlat, ylat;
# Line 13 | Line 14 | int attemp, nacc;
14   int attemp, nacc;
15  
16  
17 < int main()
17 > int main(argc, argv)
18 > int argc;
19 > char *argv[];
20   {
21 <        void getmag();
22 <        double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb);
23 <        double toterg();
24 <        void adjust();
25 <        void mcmove();
26 <        void store(FILE *outFile);
21 >  void getmag();
22 >  double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb);
23 >  double toterg();
24 >  void adjust();
25 >  void mcmove();
26 >  void store(FILE *outFile);
27 >  
28 >  double twopi, nnd, myran, kb;
29 >  double ent;
30 >  int icycl, ncycl, imove, nmoves, nsamp;
31 >  int nx, ny, which, i, j;
32 >  
33 >  FILE *inFile;
34 >  char *endptr, s[100000];
35 >  double ux, uy, uz;
36 >  
37 >  FILE *outFile;
38 >  outFile=fopen("dp_file1","a");
39 >  
40 >  srand48(MYSEED);
41 >  
42 >  // These two parameters set the size of the system:
43  
44 <        double twopi, nnd, myran, kb;
45 <        double ent;
27 <        int icycl, ncycl, imove, nmoves, nsamp;
28 <        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;
44 >  nnd = 7.0;
45 >  ny = 100;
46  
47 <        nx = (int) (2.5 * domsizey / sqrt(3.0) / nnd) + 1;
54 <        ny = (int) (domsizex / nnd) + 1;
47 >  // we want the domains to be roughly square:
48  
49 <        which = 0;
57 <        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 <        }
49 >  nx = (int) ((double)ny / sqrt(3.0));
50  
95 –        nsites= which;
51  
52 <        en=toterg();
98 <        printf("\nthe initial energy is : \t%f\n\n",en);
52 >  // periodic box boundaries:
53  
54 <        attemp = 0;
55 <        nacc = 0;
102 <        
103 <        adjust();
54 >  domsizex = sqrt(3.0) * nx * nnd;
55 >  domsizey = ny * nnd;
56  
57 <        for(icycl=1;icycl<=ncycl;icycl++)
57 >  strength = 10.0;
58 >  ncycl = 1e7;
59 >  nmoves = 100;
60 >  nsamp = 100;
61 >  twopi = 2.0 * M_PI;
62 >  dtheta = 0.1;
63 >  kb = 0.0019872198;
64 >  beta = 1.0 / (kb * 300.0);
65 >  z0 = 0.0;
66 >  deltaz = 0.1;
67 >  deltaphi = 0.1;
68 >  kz = kb;
69 >  ktheta = kb;
70 >  theta0 = M_PI / 2.0;
71 >  
72 >  which = 0;
73 >  xlat = 0;
74 >
75 >  if( strcmp(argv[1], "-resume") == 0) {
76 >    
77 >    inFile = fopen("./l_con","r");
78 >    if (inFile == NULL){
79 >            printf("Error opening file\n");
80 >            exit(-1);
81 >    }
82 >    
83 >    while(!feof(inFile))
84 >    {
85 >            fgets(s, 500, inFile);
86 >            nsites = atoi(s);
87 >            fscanf(inFile,"%s",s);
88 >            xlat= atoi(s);
89 >            fscanf(inFile,"%s",s);
90 >            ylat= atoi(s);
91 >            for(i=1; i<=nsites; i++)
92 >            {
93 >                    fscanf(inFile,"%s",s);
94 >                    fscanf(inFile,"%s",s);
95 >                    x[i] = strtod(s, &endptr);
96 >                    fscanf(inFile,"%s",s);
97 >                    y[i] = strtod(s, &endptr);
98 >                    fscanf(inFile,"%s",s);
99 >                    z[i] = strtod(s, &endptr);
100 >                    fscanf(inFile,"%s",s);
101 >                    phi[i] = strtod(s, &endptr);
102 >                    fscanf(inFile,"%s",s);
103 >                    ux = strtod(s, &endptr);
104 >                    fscanf(inFile,"%s",s);
105 >                    uy = strtod(s, &endptr);
106 >                    fscanf(inFile,"%s",s);
107 >                    uz = strtod(s, &endptr);
108 >                    theta[i] = acos(uz);
109 >                    mu[i] = strength;
110 >            }
111 >    }
112 >    fclose(inFile);
113 >    printf("%f\t%f\t%f\t%f\t%f\n", x[nsites], y[nsites], z[nsites], phi[nsites], theta[nsites]);
114 >  }
115 >
116 >  else{
117 >  for(i=0; i < nx; i++) {
118 >
119 >    xlat = xlat + 2;
120 >    ylat = 0;
121 >    
122 >    for(j=0; j<ny; j++) {    
123 >      
124 >      which = which + 1;
125 >      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;
134 >      
135 >      which = which + 1;
136 >
137 >      printf("%d\n", which);
138 >
139 >      x[which] = nnd* sqrt(3.0) * (2.0*i + 1.0) / 2.0;
140 >      y[which] = nnd* ( 2.0*j+1.0 ) / 2.0;
141 >      myran = drand48();
142 >      phi[which] = myran * twopi;
143 >      myran = drand48();
144 >      theta[which] = acos(2.0*myran-1.0);
145 >      myran = drand48();
146 >      z[which] = z0 + (2.0*myran-1.0) * deltaz;
147 >      mu[which] = strength;
148 >      ylat = ylat + 2;
149 >    }
150 >  }
151 >
152 >  nsites= which;
153 >  }
154 >  
155 >  en=toterg();
156 >  printf("\nthe initial energy is : \t%f\n\n",en);
157 >  
158 >  attemp = 0;
159 >  nacc = 0;
160 >  
161 >  adjust();
162 >
163 >  for(icycl=1;icycl<=ncycl;icycl++)
164          {
165                  for(imove=1;imove<=nmoves;imove++)
166                  {
# Line 204 | Line 262 | double eneri(double xi, double yi, double zi, double p
262          {
263                  if(j!=i)
264                  {
265 <                        dx = x[j]-xi;
266 <                        if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
267 <                        dy = y[j]-yi;
268 <                        if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
269 <                        dz = z[j]-zi;
270 <
271 <                        r2 = dx*dx+dy*dy+dz*dz;
272 <                        r = sqrt(r2);
273 <                        if(r<rcut)
274 <                        {
275 <                                r3 = r2*r;
276 <                                r5 = r2*r3;
277 <
278 <                                uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
279 <                                uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
280 <                                uzj = mu[j] * cos(theta[j]);
281 <                                udotu = uxi*uxj + uyi*uyj + uzi*uzj;
282 <                                rdotui = dx*uxi + dy*uyi + dz*uzi;
283 <                                rdotuj = dx*uxj + dy*uyj + dz*uzj;
284 <
285 <                                vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
286 <                                eni = eni + vij;
287 <                        }
265 >                  dx = x[j]-xi;                
266 >                  if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
267 >                  dy = y[j]-yi;
268 >                  if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
269 >                  dz = z[j]-zi;
270 >                  
271 >                  r2 = dx*dx+dy*dy+dz*dz;
272 >                  r = sqrt(r2);
273 >                  if(r<rcut)
274 >                    {
275 >                      r3 = r2*r;
276 >                      r5 = r2*r3;
277 >                      
278 >                      uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
279 >                      uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
280 >                      uzj = mu[j] * cos(theta[j]);
281 >                      udotu = uxi*uxj + uyi*uyj + uzi*uzj;
282 >                      rdotui = dx*uxi + dy*uyi + dz*uzi;
283 >                      rdotuj = dx*uxj + dy*uyj + dz*uzj;
284 >                      
285 >                      vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
286 >                      eni = eni + vij;
287 >                    }
288                  }
289          }
290          vint = 0.5 * kz * (zi-z0) * (zi-z0) + 0.5 * ktheta * (thetai-theta0) * (thetai-theta0);

Diff Legend

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