trunk/dp/dp.c

#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

int nsites, xlat, ylat;
double beta;
double theta[max_sites], x[max_sites], y[max_sites], z[max_sites], phi[max_sites], mu[max_sites];
double domsizex, domsizey;
double magx, magy, magz, magmag, kz, ktheta, strength, z0, theta0;
double en,dtheta,deltaz,deltaphi;
int attemp, nacc;


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(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);
  //

  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;
  domsizey = ny * nnd;

  strength = 10.0;
  ncycl = 1e7;
  nmoves = 100;
  nsamp = 1e4;
  twopi = 2.0 * M_PI;
  dtheta = 0.1;
  kb = 0.0019872198;
  beta = 1.0 / (kb * 300.0);
  z0 = 0.0;
  deltaz = 0.1;
  deltaphi = 0.1; 
  kz = kb; 
  ktheta = kb; 
  theta0 = M_PI / 2.0;
   
  which = 0; 
  xlat = 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);
    }
      
    while(!feof(inFile))
      {
        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);
            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;
          }
      }
    fclose(inFile);
    printf("%f\t%f\t%f\t%f\t%f\n", x[nsites], y[nsites], z[nsites], phi[nsites], theta[nsites]);
  } else {

    for(i=0; i < nx; i++) {
      
      xlat = xlat + 2;
      ylat = 0;
      
      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;

        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();
  printf("\nthe initial energy is : \t%f\n\n",en);
  
  attemp = 0;
  nacc = 0;
  
  adjust();

  for(icycl=1;icycl<=ncycl;icycl++)
        {
                for(imove=1;imove<=nmoves;imove++)
                {
                        mcmove(stream);
                }

                if((icycl%nsamp) == 0)
                {
                        store(outFile);
                        getmag();
                        printf("mag=%f\t%f\t%f\t%f\n", magx, magy, magz, magmag);
                }

                if( (icycl%(ncycl/5))==0 )
                {
                        printf("\n=====> Done\t%d\tout of\t%d\n", icycl, ncycl);
                        adjust();
                }
        }

        if(ncycl!=0)
        {
                if(attemp!=0)
                {
                        printf("Number of att. to disl. a part. :%d\nsuccess: %d(=%f%)\n\n", attemp, nacc, 100*(double)(nacc)/(double)(attemp));
                }
                ent=toterg();
                if(fabs(ent-en)>1e-6) printf("\n###### PROBLEMS ENERGY ###############\n");
                printf("\nTotal energy end of simulation : %f\nrunning energy : %f\ndifference : %f\n\n", ent, en, ent-en);
        }
        fclose(outFile);
        vslDeleteStream( &stream );
        return 0;
}

double toterg()
{
        double eneri(double, double, double, double, double, int, int);  
        
        double xi, yi, zi, phii, thetai, eni;
        int i, jb;
        double ener=0.0;
        for(i=1;i<=nsites;i++)
        {
                xi = x[i];
                yi = y[i];
                zi = z[i];
                phii = phi[i];
                thetai = theta[i];
                jb = i;
                eni=eneri(xi, yi, zi, phii, thetai, i, jb);
                ener = ener + eni;
        }
        return ener;

}

void getmag( )
{
        double  thetai, phii;
        int i;

        magx = 0.0;
        magy = 0.0;
        magz = 0.0;

        for(i=1;i<=nsites;i++)
        {
                phii = phi[i];
                thetai = theta[i];
                magx = magx + mu[i] * cos(phii) * sin(thetai);
                magy = magy + mu[i] * sin(phii) * sin(thetai);
                magz = magz + mu[i] * cos(thetai);
        }

        magx = magx / (double)(nsites);
        magy = magy / (double)(nsites);
        magz = magz / (double)(nsites);
        magmag = sqrt (magx*magx+magy*magy+magz*magz);
}

double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb)
{

        double uxi, uyi, uzi, dx, dy, dz, r, r2, r3, r5, uxj, uyj, uzj, rcut;
        double udotu, rdotui, rdotuj, vij, pre, vint;
        double eni;
        int j;

        pre = 14.38362;

        rcut = 30.0;

        eni = 0.0;
        uxi = mu[i] * cos(phii) * sin(thetai);
        uyi = mu[i] * sin(phii) * sin(thetai);
        uzi = mu[i] * cos(thetai);

        for(j=jb;j<=nsites;j++)
        {
                if(j!=i)
                {
                  dx = x[j]-xi;                 
                  if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizex*copysign(1.0,dx)*(double)(int)(fabs(dx /domsizex)+0.5);
                  dy = y[j]-yi;
                  if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizey*copysign(1.0,dy)*(double)(int)(fabs(dy /domsizey)+0.5);
                  dz = z[j]-zi;
                  
                  r2 = dx*dx+dy*dy+dz*dz;
                  r = sqrt(r2);
                  if(r<rcut)
                    {
                      r3 = r2*r;
                      r5 = r2*r3;
                      
                      uxj = mu[j]*cos(phi[j]) * sin(theta[j]);
                      uyj = mu[j]*sin(phi[j]) * sin(theta[j]);
                      uzj = mu[j] * cos(theta[j]);
                      udotu = uxi*uxj + uyi*uyj + uzi*uzj;
                      rdotui = dx*uxi + dy*uyi + dz*uzi;
                      rdotuj = dx*uxj + dy*uyj + dz*uzj;
                      
                      vij = pre*(udotu/r3 - 3.0*rdotui*rdotuj/r5);
                      eni = eni + vij;
                    }
                }
        }
        vint = 0.5 * kz * (zi-z0) * (zi-z0) + 0.5 * ktheta * (thetai-theta0) * (thetai-theta0);
        eni = eni + vint;
        return eni;
}

void adjust()
{
        static int attempp, naccp;
        double dzo, dphio, dthetao, frac;

        if((attemp==0)||(attempp>=attemp))
        {
                naccp = nacc;
                attempp = attemp;
        }
        else
        {
                frac = (double)(nacc-naccp)/(double)(attemp-attempp);
                dthetao = dtheta;
                dzo = deltaz;
                dphio = deltaphi;
                dtheta = dtheta * fabs(frac/0.5);
                deltaz = deltaz * fabs(frac/0.5);
                deltaphi = deltaphi * fabs(frac/0.5);

                if((dtheta/dthetao)>1.5) dtheta = dthetao*1.5;
                if((dtheta/dthetao)<0.5) dtheta = dthetao*0.5;
                if((deltaz/dzo)>1.5) deltaz = dzo * 1.5;
                if((deltaz/dzo)<0.5) deltaz = dzo * 0.5;
                if((deltaphi/dphio)>1.5) deltaphi = dphio * 1.5;
                if((deltaphi/dphio)<0.5) deltaphi = dphio * 0.5;
                printf("Max. displ. set to :%f\t%f\t%f\t(old :%f\t%f\t%f)\nFrac. acc.:%f\nattempts:%d\nsucces:%d\n\n", deltaz, deltaphi, dtheta, dzo, dphio, dthetao, frac, attemp-attempp, nacc - naccp);
                naccp = nacc;
                attempp=attemp;
        }
        
}

void mcmove( VSLStreamStatePtr stream )
{
        int o;
        double eno, zn, phin, thetan, enn, myran;
        int jb;

        attemp = attemp + 1;
        jb = 1;
        //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();
        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        zn = z[o] + (2.0*myran-1.0) * deltaz;
        //myran = drand48();
        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        phin = phi[o] + (2.0*myran-1.0) * deltaphi;
        //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();
        vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
        if(myran<=(exp(-beta*(enn-eno))))
        {
                nacc = nacc + 1;
                en = en + (enn-eno);
                z[o] = zn;
                phi[o] = phin;
                theta[o] = thetan;
        }
}


void store(FILE* outFile)
{
        double uxi, uyi, uzi;
        int i;
        
        fprintf(outFile,"%d\n",nsites);
        fprintf(outFile,"%d\t%d\n",xlat, ylat);
        
        for(i=1;i<=nsites;i++)
        {
                uxi = cos(phi[i])*sin(theta[i]);
                uyi = sin(phi[i])*sin(theta[i]);
                uzi = cos(theta[i]);
                fprintf(outFile,"%s\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n","Ar",x[i],y[i],z[i],phi[i],uxi,uyi,uzi);
        }

}
Revision:	692
Committed:	Wed Aug 13 16:43:53 2003 UTC (20 years, 10 months ago) by xsun
Content type:	text/plain
File size:	9483 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	xsun	581	#include<stdio.h>
2	xsun	692	#include<sys/time.h>
3			#include<unistd.h>
4			#include<time.h>
5	xsun	581	#include<math.h>
6			#include<stdlib.h>
7	xsun	665	#include<string.h>
8	xsun	692	#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	xsun	587	#define max_sites 15000
16	xsun	581	#define MYSEED 8973247
17
18			int nsites, xlat, ylat;
19			double beta;
20			double theta[max_sites], x[max_sites], y[max_sites], z[max_sites], phi[max_sites], mu[max_sites];
21			double domsizex, domsizey;
22			double magx, magy, magz, magmag, kz, ktheta, strength, z0, theta0;
23			double en,dtheta,deltaz,deltaphi;
24			int attemp, nacc;
25
26
27	xsun	692	int main(int argc, char *argv[])
28	xsun	581	{
29	xsun	587	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	xsun	692	void mcmove(VSLStreamStatePtr stream);
34	xsun	587	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	xsun	692	int readfile, mySeed;
41	xsun	587
42	xsun	665	FILE *inFile;
43			char *endptr, s[100000];
44			double ux, uy, uz;
45	xsun	692
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	xsun	587	FILE *outFile;
54	xsun	665	outFile=fopen("dp_file1","a");
55	xsun	587
56
57	xsun	692	//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	xsun	587	// These two parameters set the size of the system:
68	xsun	692
69	xsun	587	nnd = 7.0;
70			ny = 100;
71	xsun	692
72	xsun	587	// we want the domains to be roughly square:
73	xsun	692
74	xsun	587	nx = (int) ((double)ny / sqrt(3.0));
75	xsun	692
76
77	xsun	587	// periodic box boundaries:
78	xsun	581
79	xsun	587	domsizex = sqrt(3.0) * nx * nnd;
80			domsizey = ny * nnd;
81	xsun	581
82	xsun	587	strength = 10.0;
83			ncycl = 1e7;
84			nmoves = 100;
85	xsun	667	nsamp = 1e4;
86	xsun	587	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	xsun	692	deltaphi = 0.1;
93			kz = kb;
94			ktheta = kb;
95	xsun	587	theta0 = M_PI / 2.0;
96	xsun	665
97	xsun	587	which = 0;
98			xlat = 0;
99
100	xsun	692	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	xsun	665
109	xsun	692	if (readfile == 1) {
110
111	xsun	665	inFile = fopen("./l_con","r");
112			if (inFile == NULL){
113	xsun	692	printf("Error opening file\n");
114			exit(-1);
115	xsun	665	}
116	xsun	692
117	xsun	665	while(!feof(inFile))
118	xsun	692	{
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	xsun	665	fscanf(inFile,"%s",s);
128			fscanf(inFile,"%s",s);
129	xsun	692	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	xsun	665	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	xsun	692	} else {
149	xsun	665
150	xsun	692	for(i=0; i < nx; i++) {
151	xsun	587
152	xsun	692	xlat = xlat + 2;
153			ylat = 0;
154	xsun	587
155	xsun	692	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.0myran-1.0)deltaz;
170			mu[which] = strength;
171
172			which = which + 1;
173	xsun	587
174	xsun	692	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.0myran-1.0) deltaz;
185			mu[which] = strength;
186			ylat = ylat + 2;
187			}
188	xsun	587	}
189	xsun	692
190			nsites= which;
191	xsun	587	}
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	xsun	581	{
203			for(imove=1;imove<=nmoves;imove++)
204			{
205	xsun	692	mcmove(stream);
206	xsun	581	}
207
208			if((icycl%nsamp) == 0)
209			{
210			store(outFile);
211			getmag();
212			printf("mag=%f\t%f\t%f\t%f\n", magx, magy, magz, magmag);
213			}
214
215			if( (icycl%(ncycl/5))==0 )
216			{
217			printf("\n=====> Done\t%d\tout of\t%d\n", icycl, ncycl);
218			adjust();
219			}
220			}
221
222			if(ncycl!=0)
223			{
224			if(attemp!=0)
225			{
226			printf("Number of att. to disl. a part. :%d\nsuccess: %d(=%f%)\n\n", attemp, nacc, 100*(double)(nacc)/(double)(attemp));
227			}
228			ent=toterg();
229			if(fabs(ent-en)>1e-6) printf("\n###### PROBLEMS ENERGY ###############\n");
230			printf("\nTotal energy end of simulation : %f\nrunning energy : %f\ndifference : %f\n\n", ent, en, ent-en);
231			}
232			fclose(outFile);
233	xsun	692	vslDeleteStream( &stream );
234	xsun	581	return 0;
235			}
236
237			double toterg()
238			{
239			double eneri(double, double, double, double, double, int, int);
240
241			double xi, yi, zi, phii, thetai, eni;
242			int i, jb;
243			double ener=0.0;
244			for(i=1;i<=nsites;i++)
245			{
246			xi = x[i];
247			yi = y[i];
248			zi = z[i];
249			phii = phi[i];
250			thetai = theta[i];
251			jb = i;
252			eni=eneri(xi, yi, zi, phii, thetai, i, jb);
253			ener = ener + eni;
254			}
255			return ener;
256
257			}
258
259			void getmag( )
260			{
261			double thetai, phii;
262			int i;
263
264			magx = 0.0;
265			magy = 0.0;
266			magz = 0.0;
267
268			for(i=1;i<=nsites;i++)
269			{
270			phii = phi[i];
271			thetai = theta[i];
272			magx = magx + mu[i] * cos(phii) * sin(thetai);
273			magy = magy + mu[i] * sin(phii) * sin(thetai);
274			magz = magz + mu[i] * cos(thetai);
275			}
276
277			magx = magx / (double)(nsites);
278			magy = magy / (double)(nsites);
279			magz = magz / (double)(nsites);
280			magmag = sqrt (magxmagx+magymagy+magz*magz);
281			}
282
283			double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb)
284			{
285
286			double uxi, uyi, uzi, dx, dy, dz, r, r2, r3, r5, uxj, uyj, uzj, rcut;
287			double udotu, rdotui, rdotuj, vij, pre, vint;
288			double eni;
289			int j;
290
291			pre = 14.38362;
292
293			rcut = 30.0;
294
295			eni = 0.0;
296			uxi = mu[i] * cos(phii) * sin(thetai);
297			uyi = mu[i] * sin(phii) * sin(thetai);
298			uzi = mu[i] * cos(thetai);
299
300			for(j=jb;j<=nsites;j++)
301			{
302			if(j!=i)
303			{
304	xsun	587	dx = x[j]-xi;
305			if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizexcopysign(1.0,dx)(double)(int)(fabs(dx /domsizex)+0.5);
306			dy = y[j]-yi;
307			if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizeycopysign(1.0,dy)(double)(int)(fabs(dy /domsizey)+0.5);
308			dz = z[j]-zi;
309
310			r2 = dxdx+dydy+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 = uxiuxj + uyiuyj + uzi*uzj;
321			rdotui = dxuxi + dyuyi + dz*uzi;
322			rdotuj = dxuxj + dyuyj + dz*uzj;
323
324			vij = pre(udotu/r3 - 3.0rdotui*rdotuj/r5);
325			eni = eni + vij;
326			}
327	xsun	581	}
328			}
329			vint = 0.5 * kz * (zi-z0) * (zi-z0) + 0.5 * ktheta * (thetai-theta0) * (thetai-theta0);
330			eni = eni + vint;
331			return eni;
332			}
333
334			void adjust()
335			{
336			static int attempp, naccp;
337			double dzo, dphio, dthetao, frac;
338
339			if((attemp==0)\|\|(attempp>=attemp))
340			{
341			naccp = nacc;
342			attempp = attemp;
343			}
344			else
345			{
346			frac = (double)(nacc-naccp)/(double)(attemp-attempp);
347			dthetao = dtheta;
348			dzo = deltaz;
349			dphio = deltaphi;
350			dtheta = dtheta * fabs(frac/0.5);
351			deltaz = deltaz * fabs(frac/0.5);
352			deltaphi = deltaphi * fabs(frac/0.5);
353
354			if((dtheta/dthetao)>1.5) dtheta = dthetao*1.5;
355			if((dtheta/dthetao)<0.5) dtheta = dthetao*0.5;
356			if((deltaz/dzo)>1.5) deltaz = dzo * 1.5;
357			if((deltaz/dzo)<0.5) deltaz = dzo * 0.5;
358			if((deltaphi/dphio)>1.5) deltaphi = dphio * 1.5;
359			if((deltaphi/dphio)<0.5) deltaphi = dphio * 0.5;
360			printf("Max. displ. set to :%f\t%f\t%f\t(old :%f\t%f\t%f)\nFrac. acc.:%f\nattempts:%d\nsucces:%d\n\n", deltaz, deltaphi, dtheta, dzo, dphio, dthetao, frac, attemp-attempp, nacc - naccp);
361			naccp = nacc;
362			attempp=attemp;
363			}
364
365			}
366
367	xsun	692	void mcmove( VSLStreamStatePtr stream )
368	xsun	581	{
369			int o;
370			double eno, zn, phin, thetan, enn, myran;
371			int jb;
372
373			attemp = attemp + 1;
374			jb = 1;
375	xsun	692	//myran = drand48();
376			vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
377	xsun	581	o = (int)(nsites*myran) + 1;
378			eno=eneri(x[o], y[o], z[o], phi[o], theta[o], o, jb);
379	xsun	692	//myran = drand48();
380			vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
381	xsun	581	zn = z[o] + (2.0myran-1.0) deltaz;
382	xsun	692	//myran = drand48();
383			vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
384	xsun	581	phin = phi[o] + (2.0myran-1.0) deltaphi;
385	xsun	692	//myran = drand48();
386			vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
387	xsun	581	thetan = theta[o] + (2.0myran-1.0)dtheta;
388			enn=eneri(x[o], y[o], zn, phin, thetan, o, jb);
389	xsun	692	//myran = drand48();
390			vdRngUniform( METHOD, stream, N, &myran, 0.0, 1.0);
391	xsun	581	if(myran<=(exp(-beta*(enn-eno))))
392			{
393			nacc = nacc + 1;
394			en = en + (enn-eno);
395			z[o] = zn;
396			phi[o] = phin;
397			theta[o] = thetan;
398			}
399			}
400
401
402			void store(FILE* outFile)
403			{
404			double uxi, uyi, uzi;
405			int i;
406
407			fprintf(outFile,"%d\n",nsites);
408			fprintf(outFile,"%d\t%d\n",xlat, ylat);
409
410			for(i=1;i<=nsites;i++)
411			{
412			uxi = cos(phi[i])*sin(theta[i]);
413			uyi = sin(phi[i])*sin(theta[i]);
414			uzi = cos(theta[i]);
415			fprintf(outFile,"%s\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n","Ar",x[i],y[i],z[i],phi[i],uxi,uyi,uzi);
416			}
417
418			}