MD++: My Documents/Codes/MD++.svn/trunk/src/ising.cpp Source File

00001 // Last Modified : Fri Sep  5 18:24:11 2008
00002 
00003 #include "ising.h"
00004 
00005 void IsingFrame::initparser()
00006 {
00007     char s[100];
00008 
00009     /* Parser */
00010     bindvar("myname",myname,STRING);
00011     bindvar("command",command,STRING);
00012     bindvar("input",input,DOUBLE);
00013     bindvar("NX",&_NX,INT);
00014     bindvar("NY",&_NY,INT);
00015     bindvar("NZ",&_NZ,INT);
00016     bindvar("nsample",&_NSAMPLE,INT);
00017     bindvar("nsuccess",&_NSUCCESS,INT);
00018 
00019     /* Simulation setting */
00020     bindvar("kBT",&kBT,DOUBLE);
00021     bindvar("J",&J,DOUBLE);
00022     bindvar("H",&H,DOUBLE);
00023     
00024     bindvar("totalsteps",&totalsteps,INT);
00025     bindvar("savepropfreq",&savepropfreq,INT);
00026     bindvar("randseed",&_RANDSEED,INT);
00027 
00028     /* Forward Flux Sampling (FFS) parameters */
00029     bindvar("curstep",&curstep,INT);
00030     bindvar("continue_curstep",&continue_curstep,INT);
00031     bindvar("N_lgst_cluster",&N_lgst_cluster,INT);
00032     bindvar("N_lgst_avg",&N_lgst_avg,INT);
00033     bindvar("saveFFScn",&saveFFScn,INT);
00034     bindvar("FFSoption",&FFSoption,INT);
00035     bindvar("FFSpruning",&FFSpruning,INT);
00036     bindvar("FFScn_weight",&FFScn_weight,DOUBLE);
00037     bindvar("FFS_Pp",&FFS_Pp,DOUBLE);
00038     bindvar("FFS0_check",&FFS0_check,INT);
00039     bindvar("lambda_A",&lambda_A,INT);
00040     bindvar("lambda_B",&lambda_B,INT);
00041     bindvar("FFScurstep",&FFScurstep,INT);
00042     bindvar("FFSautoend",&FFSautoend,INT);
00043 
00044     /* File input and output */
00045     bindvar("savecn",&savecn,INT);
00046     bindvar("saveprop",&saveprop,INT);
00047     bindvar("savecnfreq",&savecnfreq,INT);
00048     bindvar("savepropfreq",&savepropfreq,INT);
00049     bindvar("printfreq",&printfreq,INT);
00050     bindvar("calcryfreq",&calcryfreq,INT); 
00051     bindvar("filecounter",&filecounter,INT);
00052     bindvar("FFSfilecounter",&FFSfilecounter,INT);
00053     bindvar("incnfile",incnfile,STRING);
00054     bindvar("finalcnfile",finalcnfile,STRING);
00055     bindvar("outpropfile",outpropfile,STRING);
00056     bindvar("intercnfile",intercnfile,STRING);
00057     bindvar("FFScnfile",FFScnfile,STRING);
00058     bindvar("zipfiles",&zipfiles,INT);
00059 
00060     /* Visualization */
00061     bindvar("win_width",&win_width,INT);
00062     bindvar("win_height",&win_height,INT);
00063     bindvar("plotfreq",&plotfreq,INT);
00064     bindvar("atomradius",atomradius,DOUBLE);
00065 
00066     bindvar("backgroundcolor",backgroundcolor,STRING);
00067     bindvar("rotateangles",rotateangles,DOUBLE);
00068     for(int i=0;i<MAXSPECIES;i++)
00069     {
00070         sprintf(s,"atomcolor%d",i);
00071         bindvar(s,atomcolor[i],STRING);
00072     }
00073     bindvar("atomcolor",atomcolor[0],STRING);
00074     for(int i=0;i<MAXCOLORS;i++)
00075     {
00076         sprintf(s,"color%02d",i);
00077         bindvar(s,colornames[i],STRING);
00078     }
00079 }    
00080 
00081 int IsingFrame::exec(char *name)
00082 {
00083     if(Organizer::exec(name)==0) return 0;
00084 
00085     /* Parser */
00086     bindcommand(name,"runcommand",runcommand());
00087 
00088     bindcommand(name,"initspin",initspin());    
00089     bindcommand(name,"initpath",initpath((int)input[0]));
00090 
00091     /* Monte Carlo simulation */
00092     bindcommand(name,"MCrun",MCrun());    
00093     bindcommand(name,"srand",srand((unsigned int)_RANDSEED));
00094     bindcommand(name,"srand48",srand48((unsigned int)_RANDSEED));
00095     bindcommand(name,"srandbytime",srand((unsigned int)time(NULL)));
00096     bindcommand(name,"srand48bytime",srand48((unsigned int)time(NULL)));
00097 
00098     bindcommand(name,"copyPathtoS",copyPathtoS((int)input[0],(int)input[1]));
00099     bindcommand(name,"SampleMCpath",SampleMCpath(input[0],input[1],(int)input[2],(int)input[3]));
00100     bindcommand(name,"WalkonChain",WalkonChain(input[0],input[1],(int)input[2]));
00101     bindcommand(name,"ComputeSuccessRate",_NSUCCESS=ComputeSuccessRate(_NSAMPLE,input[0],input[1],
00102                                                                        (int)input[2]));
00103     bindcommand(name,"AnalyzeConfigs",AnalyzeConfigs((int)input[0],(int)input[1],_NSAMPLE,
00104                                                      input[2],input[3],(int)input[4]));
00105     /* FFS implementation */
00106     bindcommand(name,"FFSprocess",FFSprocess());
00107     bindcommand(name,"calcrystalorder",calcrystalorder());
00108     bindcommand(name,"allocDFS",allocDFS());    
00109     bindcommand(name,"assign_Lam",assign_Lam());
00110             
00111     /* File input and output */
00112     bindcommand(name,"writecn",writefinalcnfile());
00113     bindcommand(name,"readcn",readcn());
00114     bindcommand(name,"openintercnfile",openintercnfile());
00115     bindcommand(name,"openFFScnfile",openFFScnfile());
00116     bindcommand(name,"writeintercn",writeintercnfile());
00117     bindcommand(name,"writeFFScn",writeFFScnfile());
00118     bindcommand(name,"setfilecounter",setfilecounter());
00119     bindcommand(name,"setFFSfilecounter",setFFSfilecounter());
00120     bindcommand(name,"openpropfile",openpropfile());
00121     bindcommand(name,"writepath",writepath((int)input[0]));
00122     bindcommand(name,"readpath",readpath((int)input[0]));
00123 
00124     /* Visualization */
00125     bindcommand(name,"openwin",openwin());
00126     bindcommand(name,"plot",plot());
00127     bindcommand(name,"alloccolors",alloccolors());
00128     bindcommand(name,"alloccolorsX",alloccolorsX());
00129     bindcommand(name,"rotate",rotate());
00130     bindcommand(name,"saverot",saverot());
00131     bindcommand(name,"reversergb",win->reversergb());
00132     bindcommand(name,"wintogglepause",wintogglepause());
00133     
00134     return -1;
00135 }
00136 
00137 void IsingFrame::runcommand()
00138 {
00139     char extcomm[400];
00140     
00141     LFile::SubHomeDir(command,command);
00142     INFO("run:  "<<command);
00143     if(nolog)
00144         system(command);
00145     else
00146     {
00147         if (ncom==0) system("rm -f B*.log");
00148         strcpy(extcomm,command);
00149         strcpy(extcomm+strlen(extcomm)," > B");
00150         sprintf(extcomm+strlen(extcomm),"%d.log",ncom);
00151         system(extcomm); 
00152         sprintf(extcomm,"echo -- '\n'B%d.log: \"%s\"'\n' >> B.log",ncom,command);
00153         system(extcomm);
00154         sprintf(extcomm,"cat B%d.log >> B.log",ncom);
00155         system(extcomm);
00156         ncom++;
00157     }
00158 }
00159 
00160 void IsingFrame::initvars()
00161 {
00162     initparser();
00163     strcpy(command,"echo Hello World");
00164     
00165     strcpy(incnfile,"spin.cn");
00166     strcpy(intercnfile,"inter.cn");
00167     strcpy(FFScnfile,"FFS.cn");
00168     strcpy(finalcnfile,"final.cn");
00169     strcpy(outpropfile,"prop.out");
00170 
00171     strcpy(bondcolor,"red");
00172     strcpy(highlightcolor,"purple");
00173 
00174     for(int i=0;i<MAXSPECIES;i++)
00175     {
00176         atomradius[i]=0.1;
00177         strcpy(atomcolor[i],"orange");
00178     }
00179     for(int i=0;i<MAXCOLORS;i++)
00180     {
00181         strcpy(colornames[i],"gray50");
00182     }
00183 
00184 }
00185 
00186 #ifdef _USETCL
00187 /********************************************************************/
00188 /* Initialize Tcl parser */
00189 /********************************************************************/
00190 int IsingFrame::Tcl_AppInit(Tcl_Interp *interp)
00191 {
00192     /*
00193      * Tcl_Init reads init.tcl from the Tcl script library.
00194      */
00195     if (Tcl_Init(interp) == TCL_ERROR) {
00196         return TCL_ERROR;
00197     }
00198     /*
00199      * Register application-specific commands.
00200      */
00201     Tcl_CreateCommand(interp, "MD++", MD_Cmd,
00202                       (ClientData)NULL, (Tcl_CmdDeleteProc *)NULL);
00203 
00204     Tcl_CreateCommand(interp, "MD++_Set", MD_SetVar,
00205                       (ClientData)NULL, (Tcl_CmdDeleteProc *)NULL);
00206 
00207     Tcl_CreateCommand(interp, "MD++_Get", MD_GetVar,
00208                       (ClientData)NULL, (Tcl_CmdDeleteProc *)NULL);
00209 
00210     /*
00211       Random_Init(interp);
00212       Blob_Init(interp);
00213     */
00214     return TCL_OK;
00215 }
00216 
00217 int Tcl_Main_Wrapper(int argc, char *argv[])
00218 {
00219     int i, pid;
00220     char oldargv[1000], newargv[1000], cmd[1000], c, *p, *filename;
00221     FILE *oldfile, *newfile;
00222 
00223     pid=getpid();
00224     for(i=1;i<argc;i++)
00225     {
00226         strcpy(oldargv,argv[i]);
00227         /* identify suffix .script */
00228         p = strrchr(oldargv,'.');
00229         if(p!=NULL) p = strstr(p,".script");
00230         if(p!=NULL)
00231         {
00232             /* identify file name */
00233             filename = strrchr(oldargv,'/');
00234             if(filename==NULL)
00235                 filename = oldargv;
00236             else
00237                 filename ++;
00238 
00239             p[0] = 0; sprintf(newargv,"/tmp/%s.tmp%d.%s",filename,pid,"tcl");
00240             p[0] = '.';
00241 
00242             INFO_Printf("replace: %s -> %s\n",oldargv,newargv);
00243 
00244             strcpy(argv[i],newargv);
00245 
00246             /* create new file */
00247             oldfile = fopen(oldargv,"r");
00248             newfile = fopen(newargv,"w");
00249 
00250             if(oldfile==NULL) ERROR("old input file open failure");
00251             if(newfile==NULL) ERROR("new input file open failure");
00252             fprintf(newfile,"# -*-shell-script-*-\n");
00253             fprintf(newfile,"#TCL script file created from %s\n\n",oldargv);
00254             fprintf(newfile,"MD++ {\n");
00255             for(;;)
00256             {
00257                 c=getc(oldfile);
00258                 if(c==EOF) break;
00259                 putc(c,newfile);
00260             }
00261             fprintf(newfile,"\n}\n");
00262             fclose(oldfile);
00263             fclose(newfile);
00264         }
00265     }
00266     INFO_Printf("\n");
00267 
00268     setenv("TCL_LIBRARY","/usr/share/tcl8.4",0);
00269     Tcl_Main(argc, argv, Tcl_AppInit);
00270 
00271     /* remove tmp files */
00272     sprintf(cmd,"rm -f /tmp/*.tmp%d.tcl\n",pid);
00273     //INFO_Printf(cmd);
00274     system(cmd);
00275 
00276     return TCL_OK;
00277 }
00278 #endif
00279 
00280 
00281 
00282 /* Memory Allocation */
00283 void IsingFrame::Alloc()
00284 {
00285     Ntot = _NX*_NY*_NZ;
00286 
00287     Realloc(_S,int,Ntot);
00288     memset(_S,-1,sizeof(int)*Ntot);
00289 
00290     Realloc(_S0,int,Ntot);
00291     memset(_S0,-1,sizeof(int)*Ntot);    
00292 }
00293 
00294 void IsingFrame::copyStoS0()
00295 {
00296     memcpy(_S0,_S,sizeof(int)*Ntot);
00297 }
00298 
00299 void IsingFrame::copyS0toS()
00300 {
00301     memcpy(_S,_S0,sizeof(int)*Ntot);
00302 }
00303 
00304 
00305 /* initialization */
00306 void IsingFrame::initspin()
00307 {
00308     int i, j, k, size;
00309 
00310     size = _NX*_NY*_NZ;
00311     
00312     Alloc();
00313 
00314     if (input[0] == -1)
00315     {
00316         memset(_S,-1,sizeof(int)*size);
00317     }
00318     else if (input[0] == 1)
00319     {
00320         memset(_S,1,sizeof(int)*size);
00321     }
00322     else if (input[0] == 0)
00323     {
00324         /* randomize */
00325         for(i=0;i<_NX;i++)
00326             for(j=0;j<_NY;j++)
00327                 for(k=0;k<_NZ;k++)
00328                 {
00329                     _S[i*_NY*_NZ+j*_NZ+k]=2*(drand48()>0.5)-1;
00330                 }
00331     }
00332 }
00333 
00334 /********************************************************************/
00335 /* Monte Carlo Simulation */
00336 /********************************************************************/
00337 void IsingFrame::MCrun()
00338 {
00339 #define MAXNN 6
00340     int i, j, k, ip, im, jp, jm, kp, km;
00341     int snn, acc, iter, key, Ntemp, Ntemp2;
00342     double Rp[2*MAXNN+1], Rm[2*MAXNN+1], R;
00343     int Sum; double avg_temp;
00344     
00345     Sum=0;
00346     /* precompute flip rate */
00347     for(i=-MAXNN;i<=MAXNN;i++)
00348     {
00349         Rp[i+MAXNN] = exp(-2*(J*i+H)/kBT);
00350         Rm[i+MAXNN] = exp( 2*(J*i+H)/kBT);
00351     }
00352 
00353     INFO_Printf("Ntot = %d\n",Ntot);
00354     /* compute total spin */
00355     Stot = 0;
00356     for(i=0;i<Ntot;i++) Stot+=_S[i];
00357 
00358     /* continue step implmented Sep/01/2008, Seunghwa */
00359     if(continue_curstep) step0 = curstep;
00360     else step0 = 0;
00361     
00362     /* FFS algorithm, works only if automatic end is initiated */
00363     if (FFSautoend==1) {
00364         Ntemp = (int) floor(0.5*Ntot);
00365         Ntemp2= (int) _Lam_array[FFSoption*2];
00366         if ( Ntemp < Ntemp2 ) {
00367             _Lam_array[FFSoption]=Ntemp;
00368         } else {
00369             _Lam_array[FFSoption]=Ntemp2;
00370         }
00371     } 
00372 
00373     for(curstep=step0;curstep<(step0+totalsteps);curstep++)
00374     { 
00375         key = 0;
00376         for(iter=0;iter<Ntot;iter++)
00377         {
00378             while(win!=NULL)
00379             {
00380                 if(win->IsPaused()) sleep(1);
00381                 else break;
00382             }
00383 
00384         /* simulation step */
00385         i = (int) floor(drand48()*_NX);
00386         j = (int) floor(drand48()*_NY);
00387 
00388         ip = (i+1)%_NX;
00389         im = (i-1+_NX)%_NX;
00390         jp = (j+1)%_NY;
00391         jm = (j-1+_NY)%_NY;
00392 
00393         if (_NZ==1)
00394         {
00395             snn = _S[ip*_NY+j]+_S[im*_NY+j]+_S[i*_NY+jp]+_S[i*_NY+jm];
00396             if (_S[i*_NY+j]>0)
00397                 R = Rp[snn+MAXNN];
00398             else
00399                 R = Rm[snn+MAXNN];
00400             acc = (drand48()<=R);
00401             if (acc)
00402             {
00403                 _S[i*_NY+j] = -_S[i*_NY+j];
00404                 Stot += 2*_S[i*_NY+j];
00405             }
00406         }
00407         else if (_NZ>1)
00408         {
00409             k = (int) floor(drand48()*_NZ);
00410             kp = (k+1)%_NZ;
00411             km = (k-1+_NZ)%_NZ;
00412             snn = _S[ip*_NY*_NZ+j*_NZ+k] + _S[im*_NY*_NZ+j*_NZ+k]
00413                 + _S[i*_NY*_NZ+jp*_NZ+k] + _S[i*_NY*_NZ+jm*_NZ+k]
00414                 + _S[i*_NY*_NZ+j*_NZ+kp] + _S[i*_NY*_NZ+j*_NZ+km];
00415 
00416             if (_S[i*_NY*_NZ+j*_NZ+k]>0)
00417                 R = Rp[snn+MAXNN];
00418             else
00419                 R = Rm[snn+MAXNN];
00420             acc = (drand48()<=R);
00421             if (acc)
00422             {
00423                 _S[i*_NY*_NZ+j*_NZ+k] = -_S[i*_NY*_NZ+j*_NZ+k];
00424                 Stot += 2*_S[i*_NY*_NZ+j*_NZ+k];
00425             }
00426         }
00427         
00428         
00429         if((saveFFScn==1)&&(curstep%calcryfreq==0)&&(iter==Ntot-1))
00430             calcrystalorder();
00431 
00432         if((curstep%printfreq==0)&&(iter==Ntot-1))
00433         {
00434             INFO_Printf("curstep = %6d Stot = %6d (%6d/%6d up) Nc = %6d Check =%d\n",curstep,Stot,(Stot+Ntot)/2,Ntot,N_lgst_cluster,FFS0_check);
00435         }
00436 
00437         if( saveFFScn == 1 && curstep!=0 && curstep!=step0 && (iter==Ntot-1))
00438         {
00439             FFScurstep=curstep;
00440             key=FFSprocess(); 
00441             if (FFSoption==0 && N_lgst_cluster > lambda_B)
00442                 initspin();
00443         }
00444         
00445         if ((saveprop)&&(iter==Ntot-1))
00446         {
00447             if((curstep%savepropfreq)==0) pf.write(this);
00448         }
00449             if (key==1) break;
00450         if((savecn)&&(iter==Ntot-1))
00451             if((curstep%savecnfreq)==0&&(curstep!=0))
00452                 intercn.write(this);
00453         
00454         if(iter==Ntot-1) {
00455             Sum=Sum+N_lgst_cluster; 
00456             winplot();
00457         }
00458         }
00459         if (key==1) break;
00460     }
00461     avg_temp=Sum/totalsteps;
00462     N_lgst_avg=(int) avg_temp;
00463 }
00464 
00465 /* compute energy of current config */
00466 double IsingFrame::Energy()
00467 {
00468     int i, j, k, ip, im, jp, jm, kp, km;
00469     int snn, snntot;
00470     
00471     Stot = 0;
00472     for(i=0;i<Ntot;i++)
00473         Stot+=_S[i];
00474     
00475     if(_NZ==1)
00476     {
00477         snntot = 0;
00478         for(i=0;i<_NX;i++)
00479           for(j=0;j<_NY;j++)
00480           {
00481               ip = (i+1)%_NX;
00482               im = (i-1+_NX)%_NX;
00483               jp = (j+1)%_NY;
00484               jm = (j-1+_NY)%_NY;
00485               
00486               snn = _S[ip*_NY+j]+_S[im*_NY+j]+_S[i*_NY+jp]+_S[i*_NY+jm];
00487               snntot += snn*_S[i*_NY+j];
00488           }
00489     }
00490     else
00491     {
00492         snntot = 0;
00493         for(i=0;i<_NX;i++)
00494           for(j=0;j<_NY;j++)
00495               for(k=0;k<_NZ;k++)
00496               {
00497                   ip = (i+1)%_NX;
00498                   im = (i-1+_NX)%_NX;
00499                   jp = (j+1)%_NY;
00500                   jm = (j-1+_NY)%_NY;
00501                   kp = (k+1)%_NZ;
00502                   km = (k-1+_NZ)%_NZ;
00503               
00504                   snn = _S[ip*_NY*_NZ+j*_NZ+k] + _S[im*_NY*_NZ+j*_NZ+k]
00505                       + _S[i*_NY*_NZ+jp*_NZ+k] + _S[i*_NY*_NZ+jm*_NZ+k]
00506                       + _S[i*_NY*_NZ+j*_NZ+kp] + _S[i*_NY*_NZ+j*_NZ+km];
00507                   
00508                   snntot += snn*_S[i*_NY*_NZ+j*_NZ+k];
00509           }
00510     }
00511     E = -snntot*J*0.5-Stot*H;
00512     return E;
00513 }
00514 
00515 
00516 /********************************************************************/
00517 /* Path sampling */
00518 /********************************************************************/
00519 void Path::Alloc(int nc,int nx,int ny,int nz)
00520 {
00521     maxlen = nc;
00522     _NX = nx; _NY = ny; _NZ = nz; 
00523     Ntot = _NX*_NY*_NZ;
00524     
00525     Realloc(_S0,int,Ntot);
00526     memset(_S0,-1,sizeof(int)*Ntot);
00527 
00528     Realloc(ix,int,maxlen);
00529     Realloc(iy,int,maxlen);
00530     Realloc(iz,int,maxlen);
00531     Realloc(newS,int,maxlen);
00532     Realloc(newStot,int,maxlen);    
00533 }
00534 
00535 void Path::SetS0(int *s)
00536 {
00537     memcpy(_S0,s,sizeof(int)*Ntot);
00538 }
00539 
00540 void Path::Free()
00541 {
00542     free(_S0);
00543     free(ix);
00544     free(iy);
00545     free(iz);
00546     free(newS);
00547     free(newStot);
00548 }
00549 
00550 void Path::Append(int i,int j,int k,int s)
00551 {
00552     if(len==maxlen)
00553     {
00554         /* expand memory allocation */
00555         maxlen = maxlen + blocksize;
00556         Realloc(ix,int,maxlen);
00557         Realloc(iy,int,maxlen);
00558         Realloc(iz,int,maxlen);
00559         Realloc(newS,int,maxlen);
00560         Realloc(newStot,int,maxlen);    
00561     }
00562     else if(len>maxlen)
00563     {
00564         ERROR("Path: memory out of bound");
00565     }
00566 
00567     /* Append entry at the end of list */
00568     ix[len]=i;
00569     iy[len]=j;
00570     iz[len]=k;
00571     newS[len]=s;
00572     /* assuming that the new spin always flips sign */
00573     if(len>0) newStot[len]=newStot[len-1]+s*2;
00574     else newStot[len]=S0tot+s*2;
00575 
00576     len ++;
00577 }
00578 
00579 /* reading and writing paths */
00580 int IsingFrame::readpath(int pn)
00581 {
00582     class Path *p;
00583     int i, j, nc, nx, ny, nz, ci, cj, ck, cs, newStot;
00584     char *buffer; char *pp, *q; 
00585     char fname[300];
00586 
00587     p = NULL;
00588     if(pn==1)      p=&pathA;
00589     else if(pn==2) p=&pathB;
00590     else if(pn==3) p=&pathC;
00591     else ERROR("initpath: pn="<<pn<<" not allowed");
00592 
00593     /* read file: incnfile */
00594     LFile::SubHomeDir(incnfile,fname);
00595     LFile::LoadToString(fname,buffer,0);
00596 
00597     pp=buffer;
00598     
00599     q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00600     sscanf(q, "%d", &nc);
00601 
00602     q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00603     sscanf(q, "%d", &nx);
00604 
00605     q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00606     sscanf(q, "%d", &ny);
00607 
00608     q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00609     sscanf(q, "%d", &nz);
00610 
00611     if((nc>p->maxlen)|(nx!=p->_NX)|(ny!=p->_NY)|(nz!=p->_NZ))
00612         p->Alloc(nc,nx,ny,nz);
00613 
00614     /* read matrix _S0 */
00615     newStot = 0;
00616     for(i=0;i<p->Ntot;i++)
00617     {
00618         q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00619         sscanf(q, "%d", p->_S0+i);
00620         newStot += p->_S0[i];
00621     }
00622 
00623     /* read incremental changes */
00624     p->len = 0;
00625     for(j=0;j<p->len;j++)
00626     {
00627         q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00628         sscanf(q, "%d",&ci);
00629         q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00630         sscanf(q, "%d",&cj);
00631         q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00632         sscanf(q, "%d",&ck);
00633         q=pp; pp=strchr(pp,'\n'); if(pp) *(char *)pp++=0;
00634         sscanf(q, "%d",&cs);
00635         
00636         p->Append(ci,cj,ck,cs);
00637 
00638         newStot += p->newS[i]*2; /* assuming spin always flips in each step of path*/
00639         p->newStot[i] = newStot;        
00640     }
00641     Free(buffer);
00642     DUMP("readpath finished");
00643     return 0;
00644 }
00645 
00646 
00647 int IsingFrame::writepath(int pn)
00648 {
00649     class Path *p;
00650     int i;
00651     FILE *fp;
00652 
00653     p = NULL;
00654     if(pn==1)      p=&pathA;
00655     else if(pn==2) p=&pathB;
00656     else if(pn==3) p=&pathC;
00657     else ERROR("writepath: pn="<<pn<<" not allowed");
00658 
00659     /* read file: finalcnfile */
00660     fp=fopen(finalcnfile,"w");
00661 
00662     fprintf(fp,"%d\n",p->len);
00663     fprintf(fp,"%d\n",p->_NX);
00664     fprintf(fp,"%d\n",p->_NY);
00665     fprintf(fp,"%d\n",p->_NZ);
00666 
00667     /* write spin matrix */
00668     for(i=0;i<p->Ntot;i++)
00669     {
00670         fprintf(fp,"%d\n", p->_S0[i]);
00671     }
00672 
00673     for(i=0;i<p->len;i++)
00674     {
00675         fprintf(fp,"%d\n%d\n%d\n%d\n",p->ix[i],p->iy[i],p->iz[i],p->newS[i]);
00676     }
00677     DUMP("writepath finished");
00678     return 0;
00679 }
00680 
00681 void IsingFrame::initpath(int pn)
00682 {
00683     class Path *p;
00684     int i, j, k, ix, iy, iz, *a;
00685     
00686     p = NULL;
00687     if(pn==1)      p=&pathA;
00688     else if(pn==2) p=&pathB;
00689     else if(pn==3) p=&pathC;
00690     else ERROR("initpath: pn="<<pn<<" not allowed");
00691 
00692     /* create a rand path from all spin down (-1) to all spin up (+1) */
00693     p->Alloc(Ntot,_NX,_NY,_NZ);
00694 
00695     /* set initial spin to point down (-1) */
00696     memset(p->_S0,-1,sizeof(int)*Ntot);
00697     p->S0tot = -Ntot;
00698     
00699     /* create random permutation array */
00700     a = (int*) malloc(sizeof(int)*Ntot);
00701     for(i=0;i<Ntot;i++) a[i]=i;
00702     for(i=0;i<Ntot;i++)
00703     {
00704         j=(int)floor((Ntot-i)*drand48());
00705         k = a[i];
00706         a[i] = a[i+j];
00707         a[i+j] = k;
00708     }
00709 
00710     /* start builing the path from the beginning */
00711     p->len = 0;
00712     for(i=0;i<Ntot;i++)
00713     {
00714         iz = a[i]%_NZ;
00715         iy = ((a[i]-iz)/_NZ)%_NY;
00716         ix = (a[i]-iz-iy*_NZ)/_NX;
00717         p->Append(ix,iy,iz,1);
00718     }
00719 
00720     free(a);
00721 }
00722 
00723 void IsingFrame::copyPathtoS(int pn,int nstep)
00724 {
00725     class Path *p;
00726     int i, id;
00727 
00728     p = NULL;
00729     if(pn==1)      p=&pathA;
00730     else if(pn==2) p=&pathB;
00731     else if(pn==3) p=&pathC;
00732     else ERROR("copyPathtoS: pn="<<pn<<" not allowed");
00733 
00734     if((_NX!=p->_NX)|(_NY!=p->_NY)|(_NZ!=p->_NZ))
00735         ERROR("copyPathtoS: dimension mismatch!");
00736     
00737     memcpy(_S,p->_S0,sizeof(int)*Ntot);
00738 
00739     for(i=0;i<nstep;i++)
00740     {
00741         id = (p->ix[i])*_NY*_NZ+(p->iy[i])*_NZ+(p->iz[i]);
00742         _S[id] = p->newS[i];
00743     }
00744     //winplot();
00745 }
00746 
00747 
00748 /* sampling of Monte Carlo path */
00749 void IsingFrame::SampleMCpath(double smin, double smax,int nmax,int pn)
00750 {  /* starting from _S (current spin) configuartion
00751     *  do MC simulation, record spin flips into path
00752     *  terminate if fraction of spin up < smin, or > smax
00753     *  or if maximum iteration exceeds nmax
00754     */
00755     class Path *p;
00756     int i, j, ip, im, jp, jm, snn, acc, iter;
00757     double Rp[9], Rm[9], R, frac_spin_up;
00758 
00759     p = NULL;
00760     if(pn==1)      p=&pathA;
00761     else if(pn==2) p=&pathB;
00762     else if(pn==3) p=&pathC;
00763     else ERROR("SampleMCpath: pn="<<pn<<" not allowed");
00764 
00765     //INFO("SampleMCpath");
00766     
00767     p->SetS0(_S);
00768     p->len = 0;
00769 
00770     /* precompute flip rate */
00771     for(i=-4;i<=4;i++)
00772     {
00773         Rp[i+4] = exp(-2*J*(i+H)/kBT);
00774         Rm[i+4] = exp( 2*J*(i+H)/kBT);
00775     }
00776 
00777     /* compute total spin */
00778     Stot = 0;
00779     for(i=0;i<Ntot;i++)
00780         Stot+=_S[i];
00781     p->S0tot = Stot;
00782     
00783      for(iter=0;iter<nmax;iter++)
00784      {
00785          while(win!=NULL)
00786          {
00787              if(win->IsPaused()) sleep(1);
00788              else break;
00789          }
00790 
00791          /* simulation step */
00792          i = (int) floor(drand48()*_NX);
00793          j = (int) floor(drand48()*_NY);
00794 
00795          ip = (i+1)%_NX;
00796          im = (i-1+_NX)%_NX;
00797          jp = (j+1)%_NY;
00798          jm = (j-1+_NY)%_NY;
00799          
00800          snn = _S[ip*_NY+j]+_S[im*_NY+j]+_S[i*_NY+jp]+_S[i*_NY+jm];
00801          if (_S[i*_NY+j]>0)
00802              R = Rp[snn+4];
00803          else
00804              R = Rm[snn+4];
00805          
00806          acc = (drand48()<=R);
00807          if (acc)
00808          {
00809              _S[i*_NY+j] = -_S[i*_NY+j];
00810              Stot += 2*_S[i*_NY+j];
00811              p->Append(i,j,0,_S[i*_NY+j]);
00812          }
00813 
00814          frac_spin_up = (Stot+Ntot)/2.0/Ntot;
00815 
00816          if((frac_spin_up<smin)|(frac_spin_up>smax))
00817              break;
00818          
00819 //         if(iter%printfreq==0)
00820 //            INFO_Printf("iter = %6d Stot = %6d (%d/%d up)\n",
00821 //                        iter,Stot,(Stot+Ntot)/2,Ntot);
00822 
00823          //if(iter%plotfreq==0) winplot();
00824      }
00825 //     INFO_Printf("SampleMCpath: Stot %d - %d\n",Stot,p->newStot[p->len-1]);     
00826 //     INFO_Printf("SampleMCpath: num spin up %d\n",(Stot+Ntot)/2);     
00827 }
00828 
00829 void IsingFrame::copypath(int pdes, int psrc)
00830 {
00831     /* copy path psrc to pdes */
00832     class Path *pd, *ps;
00833 
00834     pd = NULL;
00835     if(pdes==1)      pd=&pathA;
00836     else if(pdes==2) pd=&pathB;
00837     else if(pdes==3) pd=&pathC;
00838     else ERROR("copypath: pdes="<<pdes<<" not allowed");
00839     
00840     ps = NULL;
00841     if(psrc==1)      ps=&pathA;
00842     else if(psrc==2) ps=&pathB;
00843     else if(psrc==3) ps=&pathC;
00844     else ERROR("copypath: psrc="<<psrc<<" not allowed");
00845 
00846     pd->Alloc(ps->maxlen,ps->_NX,ps->_NY,ps->_NZ);
00847 
00848     pd->len = ps->len;
00849     pd->Ntot = ps->Ntot;
00850     pd->S0tot = ps->S0tot;
00851     
00852     memcpy(pd->_S0,ps->_S0,sizeof(int)*(ps->Ntot));
00853     memcpy(pd->ix,ps->ix,sizeof(int)*(ps->maxlen));
00854     memcpy(pd->iy,ps->iy,sizeof(int)*(ps->maxlen));
00855     memcpy(pd->iz,ps->iz,sizeof(int)*(ps->maxlen));
00856     memcpy(pd->newS,ps->newS,sizeof(int)*(ps->maxlen));
00857     memcpy(pd->newStot,ps->newStot,sizeof(int)*(ps->maxlen));
00858 
00859 }
00860 
00861     
00862 void IsingFrame::cutpath(int pn, int istart, int newlen)
00863 {
00864     /* only retain segments from istart to iend-1 */
00865     class Path *p;
00866     int i, id;
00867     
00868     p = NULL;
00869     if(pn==1)      p=&pathA;
00870     else if(pn==2) p=&pathB;
00871     else if(pn==3) p=&pathC;
00872     else ERROR("cutpath: pn="<<pn<<" not allowed");
00873 
00874     if((istart<0)||(newlen<0)||(newlen>p->len))
00875         ERROR("cutpath: istart="<<istart<<"  newlen="<<newlen);
00876 
00877     /* walk along the path to istart-1 */
00878     for(i=0;i<istart;i++)
00879     {
00880         id = (p->ix[i])*_NY*_NZ+(p->iy[i])*_NZ+(p->iz[i]);
00881         p->_S0[id] = p->newS[i];
00882         p->S0tot = p->newStot[i];
00883     }
00884 
00885     /* move path segment up front */
00886     memmove(p->ix,p->ix+istart,sizeof(int)*(newlen));
00887     memmove(p->iy,p->iy+istart,sizeof(int)*(newlen));
00888     memmove(p->iz,p->iz+istart,sizeof(int)*(newlen));
00889     memmove(p->newS,p->newS+istart,sizeof(int)*(newlen));
00890     memmove(p->newStot,p->newStot+istart,sizeof(int)*(newlen));
00891             
00892 
00893     memset(p->ix+newlen,0,sizeof(int)*(p->maxlen-newlen));
00894     memset(p->iy+newlen,0,sizeof(int)*(p->maxlen-newlen));
00895     memset(p->iz+newlen,0,sizeof(int)*(p->maxlen-newlen));
00896     memset(p->newS+newlen,0,sizeof(int)*(p->maxlen-newlen));
00897     memset(p->newStot+newlen,0,sizeof(int)*(p->maxlen-newlen));
00898     
00899     p->len = newlen;
00900 }
00901 
00902 void IsingFrame::mergepath(int phead, int ptail)
00903 {
00904     /* merge ptail to the end of phead */
00905     class Path *ph, *pt;
00906     int i;
00907     
00908     ph = NULL;
00909     if(phead==1)      ph=&pathA;
00910     else if(phead==2) ph=&pathB;
00911     else if(phead==3) ph=&pathC;
00912     else ERROR("mergepath: phead="<<phead<<" not allowed");
00913     
00914     pt = NULL;
00915     if(ptail==1)      pt=&pathA;
00916     else if(ptail==2) pt=&pathB;
00917     else if(ptail==3) pt=&pathC;
00918     else ERROR("mergepath: ptail="<<ptail<<" not allowed");
00919 
00920     if((ph->_NX!=pt->_NX)||(ph->_NY!=pt->_NY)||(ph->_NZ!=pt->_NZ))
00921         ERROR("mergepath: dimension mismatch");
00922 
00923     if(ph->newStot[ph->len-1]!=pt->S0tot)
00924     {
00925         INFO_Printf("ph Stot_end = %d pt S0tot = %d",ph->newStot[ph->len-1],pt->S0tot);
00926         ERROR("mergepath: total spin mismatch");
00927     }
00928     for(i=0;i<pt->len;i++)
00929         ph->Append(pt->ix[i],pt->iy[i],pt->iz[i],pt->newS[i]);
00930 
00931 }
00932 
00933 void IsingFrame::reversepath(int pn)
00934 {
00935     /* reverse direction of path */
00936     class Path *p;
00937     int i, k, id, len;
00938 
00939     p = NULL;
00940     if(pn==1)      p=&pathA;
00941     else if(pn==2) p=&pathB;
00942     else if(pn==3) p=&pathC;
00943     else ERROR("cutpath: pn="<<pn<<" not allowed");
00944     
00945     /* walk along the path to the end */
00946     len = p->len;
00947     for(i=0;i<len;i++)
00948     {
00949         id = (p->ix[i])*_NY*_NZ+(p->iy[i])*_NZ+(p->iz[i]);
00950         p->_S0[id] = p->newS[i];
00951         p->S0tot = p->newStot[i];        
00952     }
00953 
00954     //k=0;
00955     //for(i=0;i<p->Ntot;i++) k+=p->_S0[i];
00956     //INFO_Printf("reversepath: S0 = %d %d\n",k,p->S0tot);
00957     
00958     /* swap the arrays */
00959     for(i=0;i<len/2;i++)
00960     {
00961         k=p->ix[i]; p->ix[i]=p->ix[len-i-1]; p->ix[len-i-1]=k;
00962         k=p->iy[i]; p->iy[i]=p->iy[len-i-1]; p->iy[len-i-1]=k;
00963         k=p->iz[i]; p->iz[i]=p->iz[len-i-1]; p->iz[len-i-1]=k;
00964         k=p->newS[i]; p->newS[i]=p->newS[len-i-1]; p->newS[len-i-1]=k;
00965     }
00966 
00967     /* change the sign of newS */
00968     for(i=0;i<len;i++)
00969         p->newS[i] = -p->newS[i];
00970 
00971     /* update newStot */
00972     if(len>0) p->newStot[0]=p->S0tot+p->newS[0]*2;
00973     for(i=1;i<len;i++)
00974     {
00975         p->newStot[i]=p->newStot[i-1]+p->newS[i]*2;
00976     }
00977     //INFO_Printf("reversepath: Stot_end = %d\n",p->newStot[len-1]);
00978 }
00979 
00980 int IsingFrame::ComputeSuccessRate(int nsample,double smin,double smax,int nmax)
00981 {
00982     int nsucc, isample;
00983     double frac_spin_up;
00984 
00985     copyStoS0();
00986     nsucc = 0;
00987 
00988     for(isample=0;isample<nsample;isample++)
00989     {
00990             
00991         SampleMCpath(smin, smax, nmax, 3); /* resample path C starting from S */
00992         frac_spin_up = (pathC.newStot[pathC.len-1]+Ntot)/2.0/Ntot;
00993             
00994         if(frac_spin_up<=smin)
00995         {
00996             /* failure */
00997         }
00998         else if(frac_spin_up>=smax)
00999         {
01000             /* success */
01001             nsucc ++;
01002         }
01003         else
01004         {
01005             /* cannot tell */
01006             WARNING("ComputeSuccessRate: nmax too short to tell");
01007         }
01008         copyS0toS();
01009      }
01010      INFO_Printf("ComputeSuccessRate: nsample = %d nsuccess = %d sratio = %.5f\n",
01011                  nsample,nsucc,(1.0*nsucc)/nsample);
01012      return nsucc;
01013 }
01014     
01015 void IsingFrame::AnalyzeConfigs(int nstart,int nend,int nsample,double smin, double smax,int nmax)
01016 {
01017     int iconfig, isample, nsucc;
01018     double frac_spin_up;
01019     
01020     initpath(3); /* initialize path C */
01021     for(iconfig=nstart;iconfig<=nend;iconfig++)
01022     {
01023         /* read in cn (spin) file */
01024         sprintf(incnfile,"inter%04d.cn",iconfig);
01025         readcn();  copyStoS0();
01026         Energy();
01027         INFO_Printf("Analyzing %s  E = %e  Stot = %d Nup = %d\n",incnfile,E,Stot,(Stot+Ntot)/2);
01028 
01029         curstep = 0; winplot();
01030         
01031         nsucc = 0;
01032         for(isample=0;isample<nsample;isample++)
01033         {
01034             /* read in cn (spin) file again */
01035             if(isample>0) copyS0toS();
01036             
01037             SampleMCpath(smin, smax, nmax, 3); /* resample path C starting from S */
01038             frac_spin_up = (pathC.newStot[pathC.len-1]+Ntot)/2.0/Ntot;
01039             
01040             if(frac_spin_up<=smin)
01041             {
01042                 /* failure */
01043             }
01044             else if(frac_spin_up>=smax)
01045             {
01046                 /* success */
01047                 nsucc ++;
01048             }
01049             else
01050             {
01051                 /* cannot tell */
01052                 WARNING("AnalyzeConfig: nmax too short to tell");
01053             }
01054             if(isample%printfreq==0) INFO_Printf("nsucc = %d / %d\n",nsucc,isample+1);
01055         }
01056         INFO_Printf("Analyzing %s  E = %e  Stot = %d nsucc = %d nsample = %d sratio = %.5f\n",
01057                     incnfile,E,Stot,nsucc,nsample,(1.0*nsucc)/nsample);
01058     }
01059 }
01060 
01061 
01062 void IsingFrame::WalkonChain(double smin,double smax,int nmax)
01063 {
01064     int len, slide, dslide;
01065     double frac_spin_up;
01066     
01067     /* randomize path A */
01068     initpath(1);
01069     initpath(2);
01070     initpath(3);
01071 
01072     len = pathA.len;
01073     slide = (int) floor(drand48()*len);
01074     for(curstep=0;curstep<totalsteps;curstep++)
01075     {
01076         //slide = 3500;
01077         //slide = 2000;
01078     
01079 
01080         copyPathtoS(1, slide);             /* move S along path A */
01081         winplot();                 /* plot S before path sampling */
01082         if(savecn && (curstep%savecnfreq==0) && (curstep!=0) )
01083         {
01084             _NSUCCESS = ComputeSuccessRate(_NSAMPLE,smin,smax,nmax);
01085             intercn.write(this);   /* save saddle conf as intercn */
01086         }
01087 
01088         SampleMCpath(smin, smax, nmax, 3); /* resample path C */
01089         //winplot();                 /* plot S after path sampling */
01090 
01091         copypath(2, 1);               /* save path A to path B */
01092     
01093         frac_spin_up = (pathC.newStot[pathC.len-1]+Ntot)/2.0/Ntot;
01094 
01095         //INFO_Printf("Num spin up = %d\n",(pathC.newStot[pathC.len-1]+Ntot)/2);
01096     
01097         INFO_Printf("WalkonChain: curstep = %6d slide = %6d  len = %7d",
01098                     curstep,slide,len);
01099         if(frac_spin_up<=smin)
01100         { /* reach state A (failure) */
01101             //INFO("spin down");
01102             INFO_Printf(" (-) down %.3f\n",frac_spin_up);
01103             reversepath(1);           /* reverse path A */
01104             cutpath(1, 0, len-slide); /* truncate path A */
01105             mergepath(1, 3);          /* merge path A and C */
01106             reversepath(1);           /* reverse path A */
01107             slide = pathC.len;
01108             //dslide = (int) floor(drand48()*(len-slide)*0.3);
01109             dslide = (int) floor(drand48()*Ntot);
01110             slide += dslide;
01111             len = pathA.len;
01112             if(slide>len) slide-=dslide;
01113         }
01114         else if(frac_spin_up>=smax)
01115         { /* reach state B (success) */
01116             //INFO("spin up");
01117             INFO_Printf(" (+) up   %.3f\n",frac_spin_up);            
01118             cutpath(1, 0, slide);  /* truncate path A from 0 to slide */
01119             mergepath(1, 3);       /* merge path A and C */
01120             //slide -= (int) floor(drand48()*slide*0.3);
01121             dslide = (int) floor(drand48()*Ntot);
01122             slide -= dslide;
01123             len = pathA.len;
01124             if(slide<0) slide+=dslide;
01125         }
01126         else
01127         {
01128             //ERROR("WalkonChain: frac_spin_up="<<frac_spin_up);
01129             /* maximum iteration exceeded, try again */
01130             INFO_Printf(" try again %.3f\n",frac_spin_up);            
01131             slide = (int) floor(drand48()*len);
01132         }
01133 
01134         if(savecn && (curstep%savecnfreq==0) && (curstep!=0) )
01135         {
01136             INFO("save inter-path-A.pcn");
01137             strcpy(finalcnfile,"inter-path-A.pcn");
01138             writepath(1);
01139         }
01140     }       
01141 }
01142 
01143 
01144 /********************************************************************/
01145 /* Forward Flux Sampling (FFS) algorithm
01146 /********************************************************************/
01147 int IsingFrame::FFSprocess()
01148 {
01149     int i;
01150     if (FFSoption==0) 
01151     {
01152         if (N_lgst_cluster >= _Lam_array[0] && FFS0_check == 0 )
01153         {
01154             FFScn.write(this,zipfiles,true);
01155             FFS0_check=1;
01156         }
01157         if (FFS0_check==1 && N_lgst_cluster <= lambda_A )
01158             FFS0_check=0;
01159         return 0;
01160     }
01161     else
01162     {
01163         if (N_lgst_cluster >= _Lam_array[FFSoption])
01164         {
01165             FFS0_check = 1;
01166             FFScn.write(this,zipfiles,true);
01167             return 1;
01168         }
01169         
01170         if ( FFSpruning > 0 )
01171         {
01172             for (i=FFSoption;i>1;i--)
01173             {
01174                 if ( N_lgst_cluster<=_Lam_array[i-2] && _Lam_check[i-2]==0 )
01175                 {
01176                     if (drand48()>FFS_Pp)
01177                     {
01178                         FFScn_weight=FFScn_weight/(1-FFS_Pp);
01179                         _Lam_check[i-2]=1;
01180                     }
01181                     else
01182                     {
01183                         FFS0_check = -1;
01184                     }
01185                 }
01186                 if ( N_lgst_cluster <= lambda_A ) FFS0_check=-1;
01187             }
01188             if ( N_lgst_cluster <= lambda_A ) FFS0_check=-1;
01189         }
01190         else
01191         {
01192             if ( FFScurstep > step0 + totalsteps || N_lgst_cluster <= lambda_A )
01193                 FFS0_check = -1;
01194         }
01195         
01196         if (FFS0_check<0) return 1;
01197         else return 0;
01198     }
01199 }
01200 
01201 void IsingFrame::calcrystalorder()
01202 {
01203     int i, Nc_i;
01204     int N_cluster, Size_cluster[Ntot];
01205     N_cluster = 0;
01206     N_lgst_cluster = 0;
01207 
01208     for (i=0;i<Ntot;i++)
01209     {
01210         DFS_mark[i]=0;
01211         Size_cluster[i]=0;
01212     }
01213 
01214     for (i=0;i<Ntot;i++)
01215     {
01216         if (_S[i]>0)
01217         {
01218             N_cluster++;
01219             N_cluster_temp=0;
01220             if (N_cluster == 1)
01221                 Nc_i=i;
01222             DFS(i);
01223             if (N_cluster_temp > N_lgst_cluster)
01224             {
01225                 N_lgst_cluster = N_cluster_temp;
01226                 Nc_i=i;
01227             }
01228         }
01229     }
01230 }
01231 
01232 void IsingFrame::allocDFS()
01233 {
01234     Realloc(DFS_mark,int,Ntot);
01235 }
01236 
01237 void IsingFrame::assign_Lam()
01238 {
01239     int i, n;
01240     n = input[0];
01241     Realloc(_Lam_array,int,n);
01242     Realloc(_Lam_check,int,n);
01243     for (i=1;i<=n;i++)
01244     {
01245         _Lam_array[i-1]=input[i];
01246         _Lam_check[i-1]=0;
01247     }
01248     FFS0_check=1;
01249 }
01250 
01251 void IsingFrame::DFS(int i)
01252 {
01253     int j, ip, im, jp, jm, ind[4];
01254     int icur, jcur;
01255     div_t divresult;
01256 
01257     if (DFS_mark[i]==1 || _S[i]<0)
01258         return;
01259     else
01260     {
01261         DFS_mark[i]=1;
01262         N_cluster_temp++;
01263         divresult = div(i,_NY);
01264         icur = divresult.quot;
01265         jcur = divresult.rem;
01266         ip = (icur+1)%_NX;
01267         im = (icur-1+_NX)%_NX;
01268         jp = (jcur+1)%_NY;
01269         jm = (jcur-1+_NY)%_NY;
01270         ind[0]=ip*_NY+jcur;
01271         ind[1]=im*_NY+jcur;
01272         ind[2]=icur*_NY+jp;
01273         ind[3]=icur*_NY+jm;
01274         for (j=0;j<4;j++)
01275         {
01276             if (_S[ind[j]] < 0 || DFS_mark[ind[j]] ==1 )
01277                 continue;
01278             DFS(ind[j]);
01279         } 
01280     }
01281 
01282 }
01283 
01284 /********************************************************************/
01285 /* File input and output */
01286 /********************************************************************/
01287 int IsingFrame::readcn()
01288 {
01289     initcn.open(incnfile,LFile::O_Read);
01290     return initcn.read(this);
01291 }
01292 
01293 int IsingFrame::setfilecounter()
01294 {
01295     intercn.setcount(filecounter);
01296     return 0;
01297 }
01298 
01299 int IsingFrame::setFFSfilecounter()
01300 {
01301     FFScn.setcount(FFSfilecounter);
01302     return 0;
01303 }
01304 
01305 int IsingFrame::writefinalcnfile(int zip, bool bg)
01306 {
01307     finalcn.open(finalcnfile);
01308     return finalcn.write(this,zip,bg);
01309 }
01310 
01311 void IsingFrame::saveintercn(int step)
01312 {
01313     if(savecn)
01314         if((step%savecnfreq)==0&&(step!=0))
01315             writeintercnfile();
01316 }
01317 
01318 int IsingFrame::writeintercnfile(int zip,bool bg)
01319 {
01320     return intercn.write(this,zip,bg);
01321 }
01322 
01323 int IsingFrame::writeFFScnfile(int zip, bool bg)
01324 {
01325     return FFScn.write(this,zip,bg);
01326 }
01327 
01328 int IsingFrame::openintercnfile()
01329 {
01330     return intercn.open(intercnfile);
01331 }
01332 
01333 int IsingFrame::openFFScnfile()
01334 {
01335     return FFScn.open(FFScnfile);
01336 }
01337 
01338 int IsingFrame::openpropfile()
01339 {
01340     return pf.open(outpropfile);
01341 }
01342 
01343 /* Configuration Files */
01344 char * CNFile::describe()
01345 {
01346     static char tmp[500];
01347     sprintf(tmp,"%s","spin configurations, Format:\n"
01348             "_NX _NY _NZ; s(1,1,1); s(1,1,2); ...; s(_NX,_NY,_NZ)");
01349     return tmp;
01350 }
01351 
01352 int CNFile::writeblock(void *p)
01353 {
01354     int i;
01355     IsingFrame &d=*((IsingFrame *)p);
01356     f->Printf("%10d\n%10d\n%10d\n",d._NX,d._NY,d._NZ);
01357     for(i=0;i<d.Ntot;i++)
01358     {
01359         f->Printf("%3d\n",d._S[i]);
01360     }
01361     f->Printf("%10d\n",d._NSAMPLE);
01362     f->Printf("%10d\n",d._NSUCCESS);
01363     return 0;
01364 }
01365 
01366 int CNFile::readblock(void *p)
01367 {
01368     int i;
01369     char *buffer, *pp, *q;
01370     IsingFrame &d=*((IsingFrame *)p);
01371     
01372     LFile::LoadToString(fname,buffer,0);
01373 
01374     pp=buffer;
01375     sscanf(pp, "%d %d %d", &d._NX, &d._NY, &d._NZ);
01376     INFO("readblock: NX="<<d._NX<<" NY="<<d._NY<<" NZ="<<d._NZ);    
01377     
01378     d.Alloc();
01379     
01380     pp=strchr(pp, '\n');
01381     pp++;
01382     pp=strchr(pp, '\n');
01383     pp++;
01384     pp=strchr(pp, '\n');
01385     pp++;
01386 
01387     d.Stot = 0;
01388     for(i=0;i<d.Ntot;i++)
01389     {
01390         q=pp;
01391         pp=strchr(pp,'\n');
01392         if(pp) *(char *)pp++=0;
01393         sscanf(q, "%d",&(d._S[i]));
01394         d.Stot += d._S[i];
01395     }
01396     q=pp; pp=strchr(pp,'\n'); 
01397     if(pp) 
01398     {
01399        *(char *)pp++=0;
01400        sscanf(q, "%d",&(d._NSAMPLE));
01401        q=pp; pp=strchr(pp,'\n'); 
01402        if(pp)
01403        {
01404           *(char *)pp++=0;
01405           sscanf(q, "%d",&(d._NSUCCESS));
01406        }
01407     }
01408 
01409     Free(buffer);
01410     DUMP("readblock finished");
01411     return 0;
01412 }
01413 
01414 /* Output Property Files */
01415 char * PropFile::describe()
01416 {
01417     static char tmp[500];
01418     sprintf(tmp,"%s","Format:\n"
01419             "curstep total_spin num_spin_up ");
01420     return tmp;
01421 }
01422 
01423 int PropFile::writeentry(void *p)
01424 {
01425     char tmp[500];
01426     IsingFrame &d=*((IsingFrame *)p);
01427     sprintf(tmp,"%10d %10d %10.0f %10d %10d %2d\n",
01428             d.curstep,d.Stot,(d.Stot+d.Ntot)/2.0,d.Ntot,d.N_lgst_cluster,d.FFS0_check);
01429     *f<<tmp;
01430     return 0;
01431 }
01432 
01433 /********************************************************************/
01434 /* Visualization */
01435 /********************************************************************/
01436 
01437 void IsingFrame::openwin()
01438 {
01439     openwindow(win_width,win_height,dirname);
01440 }
01441 
01442 int IsingFrame::openwindow(int w,int h,const char *n)
01443 {
01444     if(win!=NULL)
01445         if(win->alive) return -1;
01446     win=new YWindow(w,h,n,true,true,false);
01447     if(win==NULL) return -1;
01448     else
01449     {
01450         win->setinterval(100);
01451         win->Run();
01452         return 0;
01453     }
01454 }
01455 
01456 void IsingFrame::closewindow() {delete(win);}
01457 
01458 void IsingFrame::winplot()
01459 {
01460     if(win!=NULL)
01461         if(win->alive)
01462             if((curstep%plotfreq)==0)
01463             {
01464                 plot();
01465             }
01466 }
01467 
01468 void IsingFrame::winplot(int step)
01469 {
01470     if(win!=NULL)
01471         if(win->alive)
01472             if((step%plotfreq)==0)
01473             {
01474                 plot();
01475             }
01476 }        
01477 
01478 void IsingFrame::wintogglepause()
01479 {
01480     if(win!=NULL)
01481         if(win->IsAlive())
01482             win->TogglePause();
01483 }
01484 
01485 #include "namecolor.c"
01486 void IsingFrame::alloccolors()
01487 {
01488     int r,g,b;
01489     if(win!=NULL)
01490     {
01491         if(win->alive)
01492         {
01493             for(int i=0;i<MAXCOLORS;i++)
01494             {
01495                 Str2RGB(colornames[i],&r,&g,&b);
01496                 colors[i]=win->AllocShortRGBColor(r,g,b);
01497             }
01498             
01499             Str2RGB(atomcolor[0],&r,&g,&b);
01500             colors[MAXCOLORS+0]=win->AllocShortRGBColor(r,g,b);
01501             Str2RGB(bondcolor,&r,&g,&b);
01502             colors[MAXCOLORS+1]=win->AllocShortRGBColor(r,g,b);
01503             Str2RGB(highlightcolor,&r,&g,&b);
01504             colors[MAXCOLORS+2]=win->AllocShortRGBColor(r,g,b);
01505             Str2RGB(fixatomcolor,&r,&g,&b);
01506             colors[MAXCOLORS+3]=win->AllocShortRGBColor(r,g,b);
01507             Str2RGB(backgroundcolor,&r,&g,&b);
01508             win->bgcolor=win->AllocShortRGBColor(r,g,b);
01509 
01510             for(int i=0;i<MAXSPECIES;i++)
01511             {
01512                 Str2RGB(atomcolor[i],&r,&g,&b);
01513                 colors[MAXCOLORS+5+i]=win->AllocShortRGBColor(r,g,b);
01514             }
01515         }
01516         else
01517         {
01518             WARNING("No window to allocate color for!");
01519         }
01520     }
01521     else
01522     {
01523         WARNING("No window to allocate color for!");
01524     }
01525 }
01526 void IsingFrame::alloccolorsX()
01527 {
01528     if(win!=NULL)
01529     {
01530         if(win->alive)
01531         {
01532             colors[0]=win->AllocNamedColor(atomcolor[0]);
01533             colors[1]=win->AllocNamedColor(bondcolor);
01534             colors[2]=win->AllocNamedColor(highlightcolor);
01535             colors[3]=win->AllocNamedColor(fixatomcolor);
01536             win->bgcolor=win->AllocNamedColor(backgroundcolor);
01537         }
01538         else
01539         {
01540             WARNING("No window to allocate color for!");
01541         }
01542     }
01543     else
01544     {
01545         WARNING("No window to allocate color for!");
01546     }
01547 }
01548 
01549 void IsingFrame::rotate()
01550 {
01551     if(win!=NULL)
01552     {
01553         if(win->alive)
01554         {
01555             win->horizontalRot(rotateangles[0]*M_PI/180);
01556             win->verticalRot(rotateangles[1]*M_PI/180);
01557             win->spinRot(rotateangles[2]*M_PI/180);
01558             if(rotateangles[3]!=0) win->zoom(rotateangles[3]);
01559             if(rotateangles[4]!=0) win->project(rotateangles[4]);
01560         }
01561         else
01562         {
01563             WARNING("No window to rotate for!");
01564         }
01565     }
01566     else
01567     {
01568         WARNING("No window to rotate for!");
01569     }
01570 }
01571 
01572 void IsingFrame::saverot()
01573 {
01574     if(win!=NULL)
01575     {
01576         if(win->alive)
01577         {
01578             //win->saveRot(); win->saveScale();
01579             win->saveView();
01580         }
01581         else
01582         {
01583             WARNING("No window to rotate for!");
01584         }
01585     }
01586     else
01587     {
01588         WARNING("No window to rotate for!");
01589     }
01590 }
01591 
01592 
01593 //#include "colormap.h"
01594   
01595 void IsingFrame::plot()
01596 {
01597     int i, j, k, spin;
01598     unsigned long c;
01599     
01600     if(win==NULL) return;
01601     if(!(win->alive)) return;
01602     
01603     win->Lock();
01604     win->Clear();
01605     
01606     /* draw atom s*/
01607     for(i=0;i<_NX;i++)
01608         for(j=0;j<_NY;j++)
01609             for(k=0;k<_NZ;k++)
01610             {
01611                 spin = _S[i*_NY*_NZ+j*_NZ+k];
01612 
01613                 if((_NZ>1)&&(spin<0)) continue;
01614                 
01615                 if(spin<0) c=colors[0]; else c=colors[1];
01616 
01617                 win->DrawPoint((2.0*i)/_NX-1,1-(2.0*j)/_NY,(2.0*k)/_NZ-1,
01618                                atomradius[0]/_NX,c);
01619             }
01620     
01621     win->Unlock();
01622     win->Refresh();
01623 
01624 }
01625 
01626 
01627 
01628 
01629 /***********************************************************/
01630 /*  Main Program Begins */
01631 
01632 #ifdef _TEST
01633 
01634 class IsingFrame sim;
01635 
01636 /* The main program is defined here */
01637 #include "main.cpp"
01638 
01639 #endif//_TEST