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 587 by xsun, Thu Jul 10 15:42:41 2003 UTC

# Line 1 | Line 1
1   #include<stdio.h>
2   #include<math.h>
3   #include<stdlib.h>
4 < #define max_sites 10000
4 > #define max_sites 15000
5   #define MYSEED 8973247
6  
7   int nsites, xlat, ylat;
# Line 15 | Line 15 | int main()
15  
16   int main()
17   {
18 <        void getmag();
19 <        double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb);
20 <        double toterg();
21 <        void adjust();
22 <        void mcmove();
23 <        void store(FILE *outFile);
18 >  void getmag();
19 >  double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb);
20 >  double toterg();
21 >  void adjust();
22 >  void mcmove();
23 >  void store(FILE *outFile);
24 >  
25 >  double twopi, nnd, myran, kb;
26 >  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 >  // These two parameters set the size of the system:
36  
37 <        double twopi, nnd, myran, kb;
38 <        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;
37 >  nnd = 7.0;
38 >  ny = 100;
39  
40 <        nx = (int) (2.5 * domsizey / sqrt(3.0) / nnd) + 1;
54 <        ny = (int) (domsizex / nnd) + 1;
40 >  // we want the domains to be roughly square:
41  
42 <        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 <        }
42 >  nx = (int) ((double)ny / sqrt(3.0));
43  
95 –        nsites= which;
44  
45 <        en=toterg();
98 <        printf("\nthe initial energy is : \t%f\n\n",en);
45 >  // periodic box boundaries:
46  
47 <        attemp = 0;
48 <        nacc = 0;
102 <        
103 <        adjust();
47 >  domsizex = sqrt(3.0) * nx * nnd;
48 >  domsizey = ny * nnd;
49  
50 <        for(icycl=1;icycl<=ncycl;icycl++)
50 >  strength = 10.0;
51 >  ncycl = 1e7;
52 >  nmoves = 100;
53 >  nsamp = 1e4;
54 >  twopi = 2.0 * M_PI;
55 >  dtheta = 0.1;
56 >  kb = 0.0019872198;
57 >  beta = 1.0 / (kb * 300.0);
58 >  z0 = 0.0;
59 >  deltaz = 0.1;
60 >  deltaphi = 0.1;
61 >  kz = kb;
62 >  ktheta = kb;
63 >  theta0 = M_PI / 2.0;
64 >    
65 >  which = 0;
66 >  xlat = 0;
67 >
68 >  for(i=0; i < nx; i++) {
69 >
70 >    xlat = xlat + 2;
71 >    ylat = 0;
72 >    
73 >    for(j=0; j<ny; j++) {    
74 >      
75 >      which = which + 1;
76 >      x[which] = i * sqrt(3.0) * nnd;
77 >      y[which] = j * nnd;
78 >      myran = drand48();
79 >      phi[which] = myran * twopi;
80 >      myran = drand48();
81 >      theta[which] = acos(2.0*myran-1.0);
82 >      myran = drand48();
83 >      z[which]=z0 + (2.0*myran-1.0)*deltaz;
84 >      mu[which] = strength;
85 >      
86 >      which = which + 1;
87 >
88 >      printf("%d\n", which);
89 >
90 >      x[which] = nnd* sqrt(3.0) * (2.0*i + 1.0) / 2.0;
91 >      y[which] = nnd* ( 2.0*j+1.0 ) / 2.0;
92 >      myran = drand48();
93 >      phi[which] = myran * twopi;
94 >      myran = drand48();
95 >      theta[which] = acos(2.0*myran-1.0);
96 >      myran = drand48();
97 >      z[which] = z0 + (2.0*myran-1.0) * deltaz;
98 >      mu[which] = strength;
99 >      ylat = ylat + 2;
100 >    }
101 >  }
102 >
103 >  nsites= which;
104 >  
105 >  en=toterg();
106 >  printf("\nthe initial energy is : \t%f\n\n",en);
107 >  
108 >  attemp = 0;
109 >  nacc = 0;
110 >  
111 >  adjust();
112 >
113 >  for(icycl=1;icycl<=ncycl;icycl++)
114          {
115                  for(imove=1;imove<=nmoves;imove++)
116                  {
# Line 204 | Line 212 | double eneri(double xi, double yi, double zi, double p
212          {
213                  if(j!=i)
214                  {
215 <                        dx = x[j]-xi;
216 <                        if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
217 <                        dy = y[j]-yi;
218 <                        if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
219 <                        dz = z[j]-zi;
220 <
221 <                        r2 = dx*dx+dy*dy+dz*dz;
222 <                        r = sqrt(r2);
223 <                        if(r<rcut)
224 <                        {
225 <                                r3 = r2*r;
226 <                                r5 = r2*r3;
227 <
228 <                                uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
229 <                                uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
230 <                                uzj = mu[j] * cos(theta[j]);
231 <                                udotu = uxi*uxj + uyi*uyj + uzi*uzj;
232 <                                rdotui = dx*uxi + dy*uyi + dz*uzi;
233 <                                rdotuj = dx*uxj + dy*uyj + dz*uzj;
234 <
235 <                                vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
236 <                                eni = eni + vij;
237 <                        }
215 >                  dx = x[j]-xi;                
216 >                  if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
217 >                  dy = y[j]-yi;
218 >                  if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
219 >                  dz = z[j]-zi;
220 >                  
221 >                  r2 = dx*dx+dy*dy+dz*dz;
222 >                  r = sqrt(r2);
223 >                  if(r<rcut)
224 >                    {
225 >                      r3 = r2*r;
226 >                      r5 = r2*r3;
227 >                      
228 >                      uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
229 >                      uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
230 >                      uzj = mu[j] * cos(theta[j]);
231 >                      udotu = uxi*uxj + uyi*uyj + uzi*uzj;
232 >                      rdotui = dx*uxi + dy*uyi + dz*uzi;
233 >                      rdotuj = dx*uxj + dy*uyj + dz*uzj;
234 >                      
235 >                      vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
236 >                      eni = eni + vij;
237 >                    }
238                  }
239          }
240          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