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

00001 // Last Modified : Fri Sep  5 18:03:42 2008
00002 
00003 #ifndef _ISING_H
00004 #define _ISING_H
00005 
00006 #ifdef _USETCL
00007 #include "tcl.h"
00008 #endif
00009 
00010 #include <math.h>
00011 #include <stdio.h>
00012 #include <stdlib.h>
00013 #include "general.h"
00014 #include "organizer.h"
00015 #include "display.h"
00016 
00017 /* Physical constants */
00018 #define EV   1.6021772E-19     /* (J)    */
00019 #define BOLZ 1.380658E-23      /* (J/K)  */
00020 #define KB   0.8617336E-4      /* (eV/K) */
00021 #ifndef max
00022 #define max(a,b) (((a)>(b))?(a):(b))
00023 #endif
00024 #ifndef min
00025 #define min(a,b) (((a)<(b))?(a):(b))
00026 #endif
00027 
00028 #define MAXSPECIES 10          /* maximum number of species */
00029 #define MAXCOLORS  10          /* maximum number of colors in display */
00030 
00031 /* Random number generator for platforms not supporting drand48 (e.g. cygwin) */
00032 #ifdef _NODRAND48
00033 #define drand48 (1.0/RAND_MAX)*rand
00034 #endif
00035 
00036 #define Realloc(p,t,s) {p=(t *)realloc(p,sizeof(t)*s);}
00037 
00038 
00039 /* Input/output files */
00040 class CNFile : public AUXFile /* configuration(cn) file */
00041 {
00042 public:
00043     CNFile(int t):AUXFile(t){};
00044     virtual char * describe();
00045 private:
00046     virtual int writeblock(void *p);
00047     virtual int readblock(void *p);
00048 };
00049 
00050 class PropFile : public AUXFile/* Output property file */
00051 {
00052 public:
00053     PropFile(int t):AUXFile(t){};
00054     virtual char * describe();
00055     virtual int writeentry(void *p);
00056 };
00057 
00058 class Path
00059 {
00060     enum { blocksize = 10000 };
00061 public:
00062     int maxlen;
00063     int len;
00064     int _NX, _NY, _NZ, Ntot;
00065     int *_S0, S0tot;
00066     int *ix, *iy, *iz, *newS, *newStot;
00067 
00068     Path(int nc=0,int nx=100,int ny=100,int nz=1):
00069         maxlen(0),len(nc),_NX(nx),_NY(ny),_NZ(nz),
00070         _S0(0),ix(0),iy(0),iz(0),newS(0),newStot(0) {};
00071     
00072     ~Path(){ Free(); };
00073 
00074     void Alloc(int nc,int nx,int ny,int nz);
00075     void Free();
00076 
00077     void Append(int i,int j,int k,int s);
00078     void SetS0(int *s);
00079 
00080 };
00081 
00082 /* MD simulation framework */
00083 class IsingFrame : public Organizer
00084 {
00085 public:
00086     int _NX, _NY, _NZ; /* number of spins in X, Y, Z direction */
00087                        /* _NZ = 1 for 2D problem */
00088 
00089     int *_S, *_S0;     /* spin array */
00090     int Stot;          /* sum of all spin */
00091     int Ntot;          /* total number of spins */
00092 
00093     double kBT, J, H;  /* parameters for the Hamiltonian */
00094     double E;          /* current energy (call Energy() to update */
00095     
00096     /* simulation setting */
00097     int totalsteps, curstep, continue_curstep, step0;
00098 
00099     double MC_spin, MC_accept, MC_accept_tot, MC_accept_ratio;
00100     int _RANDSEED;
00101     int _NSAMPLE, _NSUCCESS;  /* committor information for current spin configuration */
00102 
00103     /* paths */
00104     class Path pathA, pathB, pathC;
00105     
00106     /* Forward Flux Sampling (FFS) parameters */
00107     int N_lgst_cluster;/* total number of spins in the largest cluster */ 
00108     int N_lgst_avg;
00109     int N_cluster_temp;
00110     int *DFS_mark;     /* DFS algorithm marker */
00111     int *_Lam_array, *_Lam_check; /* lambda arrays for FFS */
00112     
00113     int     saveFFScn;  /* configuration save if saveFFScn=1 */
00114     int     FFSoption;  /* 0 for q/T, >1 for P(i+1/i) */
00115     int     FFSpruning; /* 0 for non, 1 for pruning */
00116     double  FFScn_weight; /* default = 1, pruning factor */
00117     double  FFS_Pp;     /* Pruning probability */
00118     int     FFS0_check; /* check for each FFS run */
00119     int     lambda_A;   /* region A if lambda < lambda_A */
00120     int     lambda_B;   /* region B if lambda > lambda_B */ 
00121     int     FFScurstep; /* step marker for time-limited FFS */ 
00122     int     FFSautoend; /* autoend mark */ 
00123    
00124     /* Visualization */
00125     YWindow *win;                   /* Display window object */
00126 
00127     int win_width, win_height;      /* size of the window */
00128     double rotateangles[5];         /* rotation angles and scaling parameter */
00129     int plotfreq;                   /* frequency of re-plotting */
00130     double atomradius[MAXSPECIES];  /* size of atoms */
00131     double bondradius, bondlength;  /* thickness and length of bonds */
00132 
00133     char atomcolor[MAXSPECIES][30]; /* color of atom species */
00134     char bondcolor[30];             /* color of bonds */    
00135     char fixatomcolor[30];          /* color of fixed atoms, fixed[i]=1 */
00136     char highlightcolor[30];        /* color of highlighted atoms */
00137     char backgroundcolor[30];       /* background color of window */
00138     char colornames[MAXCOLORS][30]; /* names of allocated colors */
00139     unsigned colors[MAXCOLORS+15];  /* value of allocated colors */
00140 
00141     /* File input and output */
00142     double input[20000];                    /* general input variable in script file */
00143     class CNFile initcn,intercn,FFScn, finalcn;    /* configuration file objects */
00144     class PropFile pf;                      /* property file objects */
00145     char incnfile[200], finalcnfile[200];   /* configuration file name */
00146     char intercnfile[200];                  /* configuration file name */
00147     char FFScnfile[200];                    /* cn file for FFS */
00148     char outpropfile[200];                  /* property file name */
00149     char myname[200];
00150     char pathfile[200];                     /* path file */
00151     char command[1000]; int ncom;           /* shell command */
00152     int  savecn, savecnfreq;                /* frequency of saving cn files */
00153     int  saveprop, savepropfreq;            /* frequency of saving prop files */
00154     int  printfreq, calcryfreq;             /* frequency of printing on scren */
00155     int  filecounter;                       /* number of current file */
00156     int  FFSfilecounter;
00157     int  zipfiles;
00158     
00159 
00160     /* Constructor (set initial values of member variables) */
00161     IsingFrame(): _NX(100), _NY(100), _NZ(1), _S(0),
00162                
00163                   /* simulation steps */
00164                   totalsteps(100), curstep(0), continue_curstep(0), step0(0),
00165                   
00166                   /* Monte Carlo simulation */
00167                   MC_spin(0),MC_accept(0),MC_accept_tot(0), MC_accept_ratio(0),
00168                   _RANDSEED(12345),
00169                   _NSAMPLE(100),_NSUCCESS(0),
00170 
00171                   /* Path sampling */
00172                   pathA(0,_NX,_NY,_NZ),
00173                   pathB(0,_NX,_NY,_NZ),
00174                   pathC(0,_NX,_NY,_NZ),
00175                   
00176                   /* FFS parameters */  
00177                    N_lgst_cluster(0), N_lgst_avg(0), N_cluster_temp(0),
00178                   DFS_mark(0), _Lam_array(0), _Lam_check(0), saveFFScn(0),
00179                   FFSoption(0), FFSpruning(0), FFScn_weight(1), FFS_Pp(0),
00180                   FFS0_check(1), lambda_A(0),lambda_B(0),FFScurstep(0),
00181                   FFSautoend(0),
00182                   
00183                   /* Visualization */
00184                   win(0),win_width(350),win_height(350),
00185                   plotfreq(1),bondradius(0.1),bondlength(0),
00186                   
00187                   /* File input and output */
00188                   initcn(AUXFile::BLOCK),intercn(AUXFile::SERIES),
00189                   FFScn(AUXFile::SERIES), finalcn(AUXFile::BLOCK),pf(AUXFile::ENTRIES),
00190                   ncom(0),savecn(0),savecnfreq(100),
00191                   saveprop(0),savepropfreq(100),
00192                   printfreq(100), calcryfreq(100), filecounter(1),FFSfilecounter(1),zipfiles(0)
00193     { 
00194         /* Input/output control */
00195         input[0]=0;
00196         
00197     };
00198 
00199     virtual ~IsingFrame() {delete(win);}
00200 
00201     /* Parser */    
00202     void initparser();
00203     int exec(char *name);            
00204     void initvars();
00205     void runcommand();        /* run shell command */    
00206 
00207 #ifdef _USETCL    
00208     virtual int Tcl_AppInit(Tcl_Interp *interp);
00209 #endif
00210 
00211     void Alloc();
00212     void initspin();               /* initialize _S */
00213     void initpath(int pn);         /* initialize path C */
00214     void copypath(int pdes, int psrc); /* copy path psrc to pdes */
00215     void cutpath(int pn, int istart, int newlen);
00216     void reversepath(int pn);
00217     void mergepath(int phead,int ptail);
00218     
00219     /* Monte Carlo simulation */
00220     void MCrun();              /* run the simulation */
00221     double Energy();           /* compute energy of current spin config */
00222    
00223     /* Path sampling */
00224     void SampleMCpath(double smin, double smax, int nmax, int pn);
00225     void WalkonChain(double smin,double smax,int nmax);
00226     int ComputeSuccessRate(int nsample,double smin,double smax,int nmax);
00227     void AnalyzeConfigs(int nstart,int nend,int nsample,double smin, double smax,int nmax);
00228     void copyPathtoS(int pn, int nstep); /* copy a state along path to S */
00229     void copyStoS0();
00230     void copyS0toS();
00231     
00232     /* Forward Flux Sampling (FFS) implementation */
00233     int FFSprocess();
00234     void calcrystalorder();
00235     void allocDFS();
00236     void assign_Lam(); 
00237     void DFS(int);
00238     
00239     /* File input and output */
00240     virtual void saveintercn(int);   /* save intermediate cn files */
00241     int readcn();
00242     int setfilecounter();
00243     int setFFSfilecounter();
00244     int openintercnfile();
00245     int openFFScnfile();
00246     int openpropfile();
00247     int writefinalcnfile(int zip=1,bool bg=true);
00248     int writeintercnfile(int zip=1,bool bg=false);
00249     int writeFFScnfile(int zip=1,bool bg=false); 
00250     int writepath(int pn);                 /* write pathC into a file */
00251     int readpath(int pn);                  /* read  pathC from a file */
00252     
00253     /* Interface with Fortran code (Relax) */
00254     void writefortraninifile(char *fname);
00255     void fortranrelax();
00256 
00257     /* Visualization */
00258     void winplot();        
00259     void winplot(int);
00260     void openwin();
00261     void plot();
00262     int  openwindow(int w,int h,const char *n);
00263     void closewindow();
00264     void wintogglepause();
00265     void alloccolors();
00266     void alloccolorsX();
00267     void rotate();
00268     void saverot();
00269     
00270 };    
00271 
00272 /* TCL driver program */
00273 #ifdef _USETCL
00274 
00275 int MD_Cmd(ClientData clientData,Tcl_Interp *interp,int argc, const char *argv[]);
00276 int MD_SetVar(ClientData clientData,Tcl_Interp *interp,int argc, const char *argv[]);
00277 int MD_GetVar(ClientData clientData,Tcl_Interp *interp,int argc, const char *argv[]);
00278 int Tcl_Main_Wrapper(int argc, char *argv[]);
00279 
00280 #define Tcl_Parser_Init(sim)                                         \
00281 int Tcl_AppInit(Tcl_Interp *interp)                                  \
00282 {                                                                    \
00283     return sim.Tcl_AppInit(interp);                                  \
00284 }                                                                    \
00285                                                                      \
00286 int MD_Cmd(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[])\
00287 {                                                                    \
00288     FILE *istream, *ostream; int mypipe[2];                          \
00289                                                                      \
00290     if (pipe (mypipe))                                               \
00291     {                                                                \
00292         fprintf (stderr, "Pipe failed.\n");                          \
00293         return EXIT_FAILURE;                                         \
00294     }                                                                \
00295     istream = fdopen(mypipe[0], "r");                                \
00296     ostream = fdopen(mypipe[1], "w");                                \
00297     for(int i=1;i<argc;i++)                                          \
00298         fprintf(ostream,"%s ",(char *)argv[i]);                      \
00299     fprintf(ostream,"\n");                                           \
00300     fclose(ostream);                                                 \
00301     sim.parse_line(istream);                                         \
00302     fclose(istream);                                                 \
00303     return TCL_OK;                                                   \
00304 }                                                                    \
00305                                                                      \
00306 int MD_SetVar(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[])\
00307 {                                                                    \
00308     int i, curn;                                                     \
00309     curn = sim.identify(argv[1]);                                    \
00310     if(curn<0)                                                       \
00311     {                                                                \
00312         WARNING("Unrecognized variable "<<argv[1]);                  \
00313         return TCL_OK;                                               \
00314     }                                                                \
00315     for(i=2;i<argc;i++)                                              \
00316     {                                                                \
00317         switch(sim.vartype[curn])                                    \
00318         {                                                            \
00319         case(INT): sscanf(argv[i],"%d",(int *)sim.varptr[curn]+i-2+sim.shift); break; \
00320         case(LONG): sscanf(argv[i],"%ld",(long *)sim.varptr[curn]+i-2+sim.shift); break;\
00321         case(DOUBLE): sscanf(argv[i],"%lf",(double *)sim.varptr[curn]+i-2+sim.shift); break;\
00322         case(STRING): strcpy((char *)sim.varptr[curn]+i-2+sim.shift,argv[i]); break;   \
00323         default: FATAL("unknown vartype ("<<sim.vartype[curn]<<")"); \
00324         }                                                            \
00325     }                                                                \
00326     return TCL_OK;                                                   \
00327 }                                                                    \
00328                                                                      \
00329 int MD_GetVar(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[])\
00330 {                                                                    \
00331     int i, curn, offset;  char buf[1000];                            \
00332     curn = sim.identify(argv[1]);                                    \
00333     if(curn<0)                                                       \
00334     {                                                                \
00335         WARNING("Unrecognized variable "<<argv[1]);                  \
00336         return TCL_OK;                                               \
00337     }                                                                \
00338     for(i=1;i<argc;i++)                                              \
00339     {                                                                \
00340         if(i==1) offset=0;                                           \
00341         else sscanf(argv[i],"%d",&offset);                           \
00342         if(i>2) sprintf(buf," ");                                    \
00343         switch(sim.vartype[curn])                                    \
00344         {                                                            \
00345         case(INT): sprintf(buf,"%d",*((int *)sim.varptr[curn]+offset+sim.shift)); break; \
00346         case(LONG): sprintf(buf,"%ld",*((long *)sim.varptr[curn]+offset+sim.shift)); break;\
00347         case(DOUBLE): sprintf(buf,"%.16g",*((double *)sim.varptr[curn]+offset+sim.shift)); break;\
00348         case(STRING): sprintf(buf,"%s",(char *)sim.varptr[curn]+offset+sim.shift); break;\
00349         default: FATAL("unknown vartype ("<<sim.vartype[curn]<<")"); \
00350         }                                                            \
00351     }                                                                \
00352     Tcl_SetResult(interp, buf, TCL_VOLATILE);                        \
00353     return TCL_OK;                                                   \
00354 }                                                                    
00355 
00356 #endif
00357 
00358 #endif // _ISING_H
00359