trunk/dp/dp.c

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<string.h>
#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(argc, argv)
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();
  void store(FILE *outFile);
  
  double twopi, nnd, myran, kb;
  double ent;
  int icycl, ncycl, imove, nmoves, nsamp;
  int nx, ny, which, i, j;
  
  FILE *inFile;
  char *endptr, s[100000];
  double ux, uy, uz;
  
  FILE *outFile;
  outFile=fopen("dp_file1","a");
  
  srand48(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( strcmp(argv[1], "-resume") == 0) {
    
    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;
      myran = drand48();
      phi[which] = myran * twopi;
      myran = drand48();
      theta[which] = acos(2.0*myran-1.0);
      myran = drand48();
      z[which]=z0 + (2.0*myran-1.0)*deltaz;
      mu[which] = strength;
      
      which = which + 1;

      printf("%d\n", which);

      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();
      phi[which] = myran * twopi;
      myran = drand48();
      theta[which] = acos(2.0*myran-1.0);
      myran = drand48();
      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();
                }

                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);
        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()
{
        int o;
        double eno, zn, phin, thetan, enn, myran;
        int jb;

        attemp = attemp + 1;
        jb = 1;
        myran = drand48();
        o = (int)(nsites*myran) + 1;
        eno=eneri(x[o], y[o], z[o], phi[o], theta[o], o, jb);
        myran = drand48();
        zn = z[o] + (2.0*myran-1.0) * deltaz;
        myran = drand48();
        phin = phi[o] + (2.0*myran-1.0) * deltaphi;
        myran = drand48();
        thetan = theta[o] + (2.0*myran-1.0)*dtheta;
        enn=eneri(x[o], y[o], zn, phin, thetan, o, jb);
        myran = drand48();
        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:	667
Committed:	Tue Aug 5 21:44:35 2003 UTC (20 years, 11 months ago) by xsun
Content type:	text/plain
File size:	8188 byte(s)
Log Message:	* empty log message *
#	User	Rev	Content
1	xsun	581	#include<stdio.h>
2			#include<math.h>
3			#include<stdlib.h>
4	xsun	665	#include<string.h>
5	xsun	587	#define max_sites 15000
6	xsun	581	#define MYSEED 8973247
7
8			int nsites, xlat, ylat;
9			double beta;
10			double theta[max_sites], x[max_sites], y[max_sites], z[max_sites], phi[max_sites], mu[max_sites];
11			double domsizex, domsizey;
12			double magx, magy, magz, magmag, kz, ktheta, strength, z0, theta0;
13			double en,dtheta,deltaz,deltaphi;
14			int attemp, nacc;
15
16
17	xsun	665	int main(argc, argv)
18			int argc;
19			char *argv[];
20	xsun	581	{
21	xsun	587	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	xsun	665	FILE *inFile;
34			char *endptr, s[100000];
35			double ux, uy, uz;
36
37	xsun	587	FILE *outFile;
38	xsun	665	outFile=fopen("dp_file1","a");
39	xsun	587
40			srand48(MYSEED);
41
42			// These two parameters set the size of the system:
43	xsun	581
44	xsun	587	nnd = 7.0;
45			ny = 100;
46	xsun	581
47	xsun	587	// we want the domains to be roughly square:
48	xsun	581
49	xsun	587	nx = (int) ((double)ny / sqrt(3.0));
50	xsun	581
51
52	xsun	587	// periodic box boundaries:
53	xsun	581
54	xsun	587	domsizex = sqrt(3.0) * nx * nnd;
55			domsizey = ny * nnd;
56	xsun	581
57	xsun	587	strength = 10.0;
58			ncycl = 1e7;
59			nmoves = 100;
60	xsun	667	nsamp = 1e4;
61	xsun	587	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	xsun	665
72	xsun	587	which = 0;
73			xlat = 0;
74
75	xsun	665	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	xsun	587	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.0myran-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.0myran-1.0) deltaz;
147			mu[which] = strength;
148			ylat = ylat + 2;
149			}
150			}
151
152			nsites= which;
153	xsun	665	}
154	xsun	587
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	xsun	581	{
165			for(imove=1;imove<=nmoves;imove++)
166			{
167			mcmove();
168			}
169
170			if((icycl%nsamp) == 0)
171			{
172			store(outFile);
173			getmag();
174			printf("mag=%f\t%f\t%f\t%f\n", magx, magy, magz, magmag);
175			}
176
177			if( (icycl%(ncycl/5))==0 )
178			{
179			printf("\n=====> Done\t%d\tout of\t%d\n", icycl, ncycl);
180			adjust();
181			}
182			}
183
184			if(ncycl!=0)
185			{
186			if(attemp!=0)
187			{
188			printf("Number of att. to disl. a part. :%d\nsuccess: %d(=%f%)\n\n", attemp, nacc, 100*(double)(nacc)/(double)(attemp));
189			}
190			ent=toterg();
191			if(fabs(ent-en)>1e-6) printf("\n###### PROBLEMS ENERGY ###############\n");
192			printf("\nTotal energy end of simulation : %f\nrunning energy : %f\ndifference : %f\n\n", ent, en, ent-en);
193			}
194			fclose(outFile);
195			return 0;
196			}
197
198			double toterg()
199			{
200			double eneri(double, double, double, double, double, int, int);
201
202			double xi, yi, zi, phii, thetai, eni;
203			int i, jb;
204			double ener=0.0;
205			for(i=1;i<=nsites;i++)
206			{
207			xi = x[i];
208			yi = y[i];
209			zi = z[i];
210			phii = phi[i];
211			thetai = theta[i];
212			jb = i;
213			eni=eneri(xi, yi, zi, phii, thetai, i, jb);
214			ener = ener + eni;
215			}
216			return ener;
217
218			}
219
220			void getmag( )
221			{
222			double thetai, phii;
223			int i;
224
225			magx = 0.0;
226			magy = 0.0;
227			magz = 0.0;
228
229			for(i=1;i<=nsites;i++)
230			{
231			phii = phi[i];
232			thetai = theta[i];
233			magx = magx + mu[i] * cos(phii) * sin(thetai);
234			magy = magy + mu[i] * sin(phii) * sin(thetai);
235			magz = magz + mu[i] * cos(thetai);
236			}
237
238			magx = magx / (double)(nsites);
239			magy = magy / (double)(nsites);
240			magz = magz / (double)(nsites);
241			magmag = sqrt (magxmagx+magymagy+magz*magz);
242			}
243
244			double eneri(double xi, double yi, double zi, double phii, double thetai, int i, int jb)
245			{
246
247			double uxi, uyi, uzi, dx, dy, dz, r, r2, r3, r5, uxj, uyj, uzj, rcut;
248			double udotu, rdotui, rdotuj, vij, pre, vint;
249			double eni;
250			int j;
251
252			pre = 14.38362;
253
254			rcut = 30.0;
255
256			eni = 0.0;
257			uxi = mu[i] * cos(phii) * sin(thetai);
258			uyi = mu[i] * sin(phii) * sin(thetai);
259			uzi = mu[i] * cos(thetai);
260
261			for(j=jb;j<=nsites;j++)
262			{
263			if(j!=i)
264			{
265	xsun	587	dx = x[j]-xi;
266			if(fabs(dx)>(domsizex/2.0)) dx = dx - domsizexcopysign(1.0,dx)(double)(int)(fabs(dx /domsizex)+0.5);
267			dy = y[j]-yi;
268			if(fabs(dy)>(domsizey/2.0)) dy = dy - domsizeycopysign(1.0,dy)(double)(int)(fabs(dy /domsizey)+0.5);
269			dz = z[j]-zi;
270
271			r2 = dxdx+dydy+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 = uxiuxj + uyiuyj + uzi*uzj;
282			rdotui = dxuxi + dyuyi + dz*uzi;
283			rdotuj = dxuxj + dyuyj + dz*uzj;
284
285			vij = pre(udotu/r3 - 3.0rdotui*rdotuj/r5);
286			eni = eni + vij;
287			}
288	xsun	581	}
289			}
290			vint = 0.5 * kz * (zi-z0) * (zi-z0) + 0.5 * ktheta * (thetai-theta0) * (thetai-theta0);
291			eni = eni + vint;
292			return eni;
293			}
294
295			void adjust()
296			{
297			static int attempp, naccp;
298			double dzo, dphio, dthetao, frac;
299
300			if((attemp==0)\|\|(attempp>=attemp))
301			{
302			naccp = nacc;
303			attempp = attemp;
304			}
305			else
306			{
307			frac = (double)(nacc-naccp)/(double)(attemp-attempp);
308			dthetao = dtheta;
309			dzo = deltaz;
310			dphio = deltaphi;
311			dtheta = dtheta * fabs(frac/0.5);
312			deltaz = deltaz * fabs(frac/0.5);
313			deltaphi = deltaphi * fabs(frac/0.5);
314
315			if((dtheta/dthetao)>1.5) dtheta = dthetao*1.5;
316			if((dtheta/dthetao)<0.5) dtheta = dthetao*0.5;
317			if((deltaz/dzo)>1.5) deltaz = dzo * 1.5;
318			if((deltaz/dzo)<0.5) deltaz = dzo * 0.5;
319			if((deltaphi/dphio)>1.5) deltaphi = dphio * 1.5;
320			if((deltaphi/dphio)<0.5) deltaphi = dphio * 0.5;
321			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);
322			naccp = nacc;
323			attempp=attemp;
324			}
325
326			}
327
328			void mcmove()
329			{
330			int o;
331			double eno, zn, phin, thetan, enn, myran;
332			int jb;
333
334			attemp = attemp + 1;
335			jb = 1;
336			myran = drand48();
337			o = (int)(nsites*myran) + 1;
338			eno=eneri(x[o], y[o], z[o], phi[o], theta[o], o, jb);
339			myran = drand48();
340			zn = z[o] + (2.0myran-1.0) deltaz;
341			myran = drand48();
342			phin = phi[o] + (2.0myran-1.0) deltaphi;
343			myran = drand48();
344			thetan = theta[o] + (2.0myran-1.0)dtheta;
345			enn=eneri(x[o], y[o], zn, phin, thetan, o, jb);
346			myran = drand48();
347			if(myran<=(exp(-beta*(enn-eno))))
348			{
349			nacc = nacc + 1;
350			en = en + (enn-eno);
351			z[o] = zn;
352			phi[o] = phin;
353			theta[o] = thetan;
354			}
355			}
356
357
358			void store(FILE* outFile)
359			{
360			double uxi, uyi, uzi;
361			int i;
362
363			fprintf(outFile,"%d\n",nsites);
364			fprintf(outFile,"%d\t%d\n",xlat, ylat);
365
366			for(i=1;i<=nsites;i++)
367			{
368			uxi = cos(phi[i])*sin(theta[i]);
369			uyi = sin(phi[i])*sin(theta[i]);
370			uzi = cos(theta[i]);
371			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);
372			}
373
374			}