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

00001 /*
00002   ffs.cpp
00003   by Seunghwa Ryu, shryu@stanford.edu
00004   Last Modified : Mon Sep  8 17:16:38 2008
00005 
00006   FUNCTION  :  Foward Flux Sampling for rare transition events
00007 
00008   To Do:
00009       1.  calcrystalorder() needs to be improved
00010 
00011   Note:
00012       This file is included by md.cpp.  This file cannot be compiled by itself.
00013 */
00014 
00015 /********************************************************************/
00016 /* Forward Flux Sampling (FFS) function (Seunghwa Ryu)              */
00017 /********************************************************************/
00018 void MDFrame::calcrystalorder()
00019 {
00020     int i, j, jpt, m, M, Nb,Nc_i;
00021     Vector3 sij, rij;
00022     double rrij, xij, yij, zij, r_l;
00023     double metric[30];
00024 
00025     int N_cluster, Size_cluster[_NP];    
00026     
00027     M = 2*L_for_QLM + 2;
00028     N_solid_P = 0;
00029     N_skin_P = 0;
00030     N_lgst_cluster = 0;
00031     N_lgst_skin = 0;
00032     N_cluster = 0;
00033 
00034 //      printf ("debug\n");
00035     if (L_for_QLM==3)
00036     {
00037         metric[0]=5.0;
00038         metric[1]=metric[0];
00039         metric[2]=30.0;
00040         metric[3]=metric[2];
00041         metric[4]=3.0;
00042         metric[5]=metric[4];
00043         metric[6]=2.0;
00044     }
00045      else if (L_for_QLM==6)
00046     {
00047         metric[0] = 231.0;
00048         metric[1] = metric[0];
00049         metric[2] = 2772.0;
00050         metric[3] = metric[2];
00051         metric[4] = 126.0;
00052         metric[5] = metric[4];
00053         metric[6] = 420.0;
00054         metric[7] = metric[6];
00055         metric[8] = 105.0;
00056         metric[9] = metric[8];
00057         metric[10]= 168.0;
00058         metric[11]= metric[10];
00059         metric[12]= 2.0;
00060     }
00061     else
00062     {
00063         ERROR("calcrystalorder(): L = "<<L_for_QLM<<" is not implemented");
00064         return;
00065     }
00066 //      printf("debug1 \n");
00067     /* initialize QLM and QL to zero */
00068     for(i=0 ; i < _NP ; i++)
00069     {
00070         QLM_mark[i]=0;
00071         QLM_solid[i]=0;
00072         Size_cluster[i]=0;
00073         _QL[i]=0;
00074         for(m=0;m<M;m++)
00075             _QLM[i*M+m]=0;
00076     }
00077 //      printf("debug1.5\n");
00078     for(i=0;i<_NP;i++)
00079     {
00080         Nb=0;
00081         for(j=0;j<nn[i];j++)
00082         {
00083             jpt=nindex[i][j];
00084             if(jpt==i) continue;
00085             sij=_SR[jpt]-_SR[i]; sij.subint();
00086             rij=_H*sij; rrij=rij.norm();
00087 
00088             if(rrij > Rc_for_QLM) continue;
00089             xij=rij.x;
00090             yij=rij.y;
00091             zij=rij.z;
00092             r_l=pow(rrij,L_for_QLM);
00093 
00094             /* rij.x, rij.y, rij.z
00095              * rij[0], rij[1], rij[2]
00096              */
00097             Nb++;
00098             if (L_for_QLM == 3)
00099             {   
00100                  _QLM[i*M+ 0 ] += ( pow(xij,3.0) - 3.0*xij*pow(yij,2.0) )/r_l;
00101                  _QLM[i*M+ 1 ] += ( pow(yij,3.0) - 3.0*yij*pow(xij,2.0) )/r_l;
00102                  _QLM[i*M+ 2 ] += ( pow(xij,2.0) - pow(yij,2.0) )*zij/r_l;
00103                  _QLM[i*M+ 3 ] += -2.0*xij*yij*zij/r_l;
00104                  _QLM[i*M+ 4 ] += (4.0*zij*zij-xij*xij-yij*yij)*xij/r_l;
00105                  _QLM[i*M+ 5 ] += -(4.0*zij*zij-xij*xij-yij*yij)*yij/r_l;
00106                  _QLM[i*M+ 6 ] += zij*(2.0*zij*zij-3.0*xij*xij-3.0*yij*yij)/r_l;            }
00107             else if (L_for_QLM ==6)
00108             {
00109                 _QLM[i*M+ 0 ] += (pow(xij,6)-15*pow(xij,4)*pow(yij,2)+15*pow(xij,2)*pow(yij,4)-pow(yij,6))/r_l;
00110                 _QLM[i*M+ 1 ] += 2*xij*yij*(3*pow(xij,4)-10*pow(xij,2)*pow(yij,2)+3*pow(yij,4))/r_l;
00111                 _QLM[i*M+ 2 ] += xij*zij*(pow(xij,4)-10*pow(xij,2)*pow(yij,2)+5*pow(yij,4))/r_l;
00112                 _QLM[i*M+ 3 ] += yij*zij*(5*pow(xij,4)-10*pow(xij,2)*pow(yij,2)+pow(yij,4))/r_l;
00113                 _QLM[i*M+ 4 ] += -(-10*pow(xij,4)*pow(zij,2)+pow(xij,6)-5*pow(xij,4)*pow(yij,2)+60*pow(xij*yij*zij,2)-5*pow(xij,2)*pow(yij,4)-10*pow(yij,4)*pow(zij,2)+pow(yij,6))/r_l;
00114                 _QLM[i*M+ 5 ] += -4*(-10*pow(zij,2)+xij*xij+yij*yij)*xij*yij*(xij*xij-yij*yij)/r_l;
00115                 _QLM[i*M+ 6 ] += -xij*zij*(-8*pow(xij,2)*pow(zij,2)+3*pow(xij,4)-6*pow(xij,2)*pow(yij,2)+24*pow(yij,2)*pow(zij,2)-9*pow(yij,4))/r_l;
00116                 _QLM[i*M+ 7 ] += -zij*(-8*pow(zij,2)+3*xij*xij+3*yij*yij)*yij*(3*xij*xij-yij*yij)/r_l;
00117                 _QLM[i*M+ 8 ] += (16*pow(xij,2)*pow(zij,4)-16*pow(xij,4)*pow(zij,2)+pow(xij,6)+pow(xij,4)*pow(yij,2)-pow(xij,2)*pow(yij,4)-16*pow(yij,2)*pow(zij,4)+16*pow(yij,4)*pow(zij,2)-pow(yij,6))/r_l;
00118                 _QLM[i*M+ 9 ] += 2*(16*pow(zij,4)-16*pow(xij,2)*pow(zij,2)-16*pow(yij,2)*pow(zij,2)+pow(xij,4)+2*pow(xij*yij,2)+pow(yij,4))*xij*yij/r_l;
00119                 _QLM[i*M+ 10] += xij*zij*(8*pow(zij,4)-20*pow(xij*zij,2)-20*pow(yij*zij,2)+5*pow(xij,4)+10*pow(xij*yij,2)+5*pow(yij,4))/r_l;
00120                 _QLM[i*M+ 11] += yij*zij*(8*pow(zij,4)-20*pow(xij*zij,2)-20*pow(yij*zij,2)+5*pow(xij,4)+10*pow(xij*yij,2)+5*pow(yij,4))/r_l;
00121                 _QLM[i*M+ 12] += -(-16*pow(zij,6)+120*pow(xij,2)*pow(zij,4)+120*pow(yij,2)*pow(zij,4)-90*pow(xij,4)*pow(zij,2)-180*pow(xij*yij*zij,2)-90*pow(yij,4)*pow(zij,2)+5*pow(xij,6)+15*pow(xij,4)*pow(yij,2)+15*pow(xij,2)*pow(yij,4)+5*pow(yij,6))/r_l;
00122             }
00123         }
00124         /* the last entry is the magnitude of the vector */
00125         for (m=0;m<M-1;m++)
00126             _QLM[i*M+M-1] += metric[m]*_QLM[i*M+m]*_QLM[i*M+m];
00127         _QLM[i*M+M-1]=sqrt(_QLM[i*M+M-1]);
00128 //        printf("%d %f\n",Nb,_QLM[i*M+M-1]);
00129     }
00130 
00131 //      printf("debug2\n");
00132     for(i=0;i<_NP;i++)
00133     {
00134         Nb=0;
00135         for(j=0;j<nn[i];j++)
00136         {
00137             jpt=nindex[i][j];
00138             if(jpt==i) continue;
00139             sij=_SR[jpt]-_SR[i]; sij.subint();
00140             rij=_H*sij; rrij=rij.norm();  
00141             if(rrij > Rc_for_QLM) continue;
00142             Nb++; 
00143             for(m=0;m<M-1;m++)
00144             {
00145                 _QL[i] += _QLM[i*M+m] * _QLM[jpt*M+m] * metric[m] / _QLM[jpt*M+M-1];
00146             }
00147         }
00148         if (Nb>0)
00149         {
00150         _QL[i]/= fabs((_QLM[i*M+M-1]*Nb));
00151                 if (L_for_QLM==6)
00152                         _QL[i]=-_QL[i];
00153         }
00154         else
00155         _QL[i] = 2; 
00156 //        printf("%d %f\n",Nb,_QL[i]);
00157 
00158         /* put results in TOPOL array for visualization */
00159         _TOPOL[i] = _QL[i];
00160         if (_QL[i] <= QLM_cutoff )
00161         {
00162                 N_solid_P++;
00163                 QLM_solid[i]=1;
00164         }
00165 
00166 //      _TOPOL[i] = _QLM[i*M+M-1];      
00167     }
00168     for (i=0;i<_NP;i++)
00169     {
00170         if ( QLM_mark[i]==0 && QLM_solid[i]==1 )
00171         {
00172             N_cluster++;
00173             N_cluster_temp = 0;
00174             if (N_cluster == 1)
00175                 Nc_i=i;
00176             if (wSKIN>0) DFS1(i,0,0,0);
00177             else DFS(i,0,0);
00178             Size_cluster[N_cluster]=N_cluster_temp;
00179             if (N_cluster_temp > N_lgst_cluster)
00180             {
00181                 N_lgst_cluster = N_cluster_temp;
00182                 Nc_i=i;
00183             }
00184 //          printf("Size_cluster = %d\n", N_cluster_temp);
00185 //          printf("N_largest_cluster = %d\n", N_lgst_cluster);     
00186         }
00187     }
00188     if ( N_lgst_index != 0 )
00189     {
00190 //        N_lgst_cluster = 0;
00191         for (i=0;i<_NP;i++)
00192             QLM_mark[i]=0;
00193         if (wSKIN >0 )
00194             DFS1(Nc_i,1,N_lgst_index,N_lgst_skin_index);
00195         else
00196             DFS(Nc_i,1,N_lgst_index);
00197     }
00198 //    printf("N_solid_P = %d\n",N_solid_P);
00199 //    printf("N_skin_P = %d\n",N_skin_P);
00200 //    printf("N_largest_cluster = %d\n",N_lgst_cluster);
00201 }
00202 
00203 void MDFrame::DFS(int i, int On_Off, double index )
00204 {
00205     int j, jpt;
00206     double rrij;
00207     Vector3 sij, rij;
00208 
00209     if (QLM_mark[i]==1 || QLM_solid[i]==0)
00210         return;
00211     else
00212     {
00213 //      printf("QLM_mark=%d , QLM_solid=%d, i = %d \n",QLM_mark[i],QLM_solid[i],i);
00214         QLM_mark[i]=1;
00215         if ( On_Off > 0 ) 
00216         {
00217             _QL[i]=index;
00218             _TOPOL[i]=_QL[i];
00219         }
00220         N_cluster_temp++;
00221         for (j=0;j<nn[i];j++)
00222         {
00223             jpt = nindex[i][j];
00224             if ( jpt == i || QLM_solid[jpt] == 0 || QLM_mark[jpt] == 1)
00225                 continue;
00226             sij=_SR[jpt]-_SR[i]; sij.subint();
00227             rij=_H*sij; rrij=rij.norm();
00228            
00229             if (rrij > Rc_for_QLM) continue;
00230             DFS(jpt, On_Off, index);
00231         }
00232     }
00233 }
00234 
00235 void MDFrame::DFS1(int i, int On_Off, double index, double index1 )
00236 {
00237     int j, jpt;
00238     double rrij;
00239     Vector3 sij, rij;
00240 
00241     if (QLM_mark[i]==1 && QLM_solid[i]==0)
00242         return;
00243     else if (QLM_mark[i]==0 && QLM_solid[i]==0)
00244     {
00245         QLM_mark[i]=1;
00246         if (On_Off > 0 )
00247         {
00248             _QL[i]=index1;
00249             _TOPOL[i]=_QL[i];
00250             N_lgst_skin++;
00251         }
00252 //        _QL[i]=-1.5;
00253 //        _TOPOL[i]=_QL[i];
00254         N_skin_P++;
00255         N_cluster_temp++;
00256         return; 
00257     }
00258     else
00259     {
00260 //      printf("QLM_mark=%d , QLM_solid=%d, i = %d \n",QLM_mark[i],QLM_solid[i],i);
00261         QLM_mark[i]=1;
00262         if ( On_Off > 0 )
00263         {
00264             _QL[i]=index;
00265             _TOPOL[i]=_QL[i];
00266         }
00267         N_cluster_temp++;
00268         for (j=0;j<nn[i];j++)
00269         {
00270             jpt = nindex[i][j];
00271             if ( jpt == i || QLM_mark[jpt] == 1)
00272                 continue;
00273             sij=_SR[jpt]-_SR[i]; sij.subint();
00274             rij=_H*sij; rrij=rij.norm();
00275            
00276             if (rrij > Rc_for_QLM) continue;
00277             DFS1(jpt, On_Off, index, index1);
00278         }
00279     }
00280 }
00281 
00282 void MDFrame::allocQLM()
00283 {
00284     int size;
00285 
00286     if (L_for_QLM <= 0)
00287     {
00288         ERROR("L_for_QLM = "<<L_for_QLM<<" skip allocQLM");
00289         return;
00290     }
00291     if (_NP <= 0)
00292     {
00293         ERROR("NP = "<<_NP<<" skip allocQLM");
00294         return;
00295     }
00296     size = (2*L_for_QLM + 2)*_NP;
00297 
00298     Realloc(_QLM,double,size);
00299    
00300     Realloc(QLM_solid,int,_NP);
00301 
00302     Realloc(QLM_mark,int,_NP);
00303 
00304     Realloc(_QL,double,_NP);
00305 
00306 }
00307 
00308 void MDFrame::assign_Lam()
00309 {
00310     int i, n;
00311     n = (int)input[0];
00312     Realloc(_Lam_array,int,n);
00313     Realloc(_Lam_check,int,n);
00314     for (i=1;i<=n;i++)
00315     {
00316         _Lam_array[i-1]=(int)input[i];
00317         _Lam_check[i-1]=0;
00318     }
00319     FFS0_check=1;
00320 }
00321 
00322 int MDFrame::FFSprocess()
00323 {
00324     int i;
00325     if (FFSoption==0) 
00326     {
00327         if (N_lgst_cluster >= _Lam_array[0] && FFS0_check == 0 )
00328         {
00329             FFScn.write(this,zipfiles,true);
00330             FFS0_check=1;
00331         }
00332         if (FFS0_check==1 && N_lgst_cluster <= FFS_i )
00333             FFS0_check=0;
00334         return 0;
00335     }
00336     else
00337     {
00338         if (N_lgst_cluster >= _Lam_array[FFSoption])
00339         {
00340             FFS0_check = 1;
00341             FFScn.write(this,zipfiles,true);
00342             return 1;
00343         }
00344         
00345         if ( FFSpruning > 0 )
00346         {
00347             for (i=FFSoption;i>1;i--)
00348             {
00349                 if ( N_lgst_cluster<=_Lam_array[i-2] && _Lam_check[i-2]==0 )
00350                 {
00351                     if (drand48()>FFS_Pp)
00352                     {
00353                         FFScn_weight=FFScn_weight/(1-FFS_Pp);
00354                         _Lam_check[i-2]=1;
00355                     }
00356                     else
00357                     {
00358                         FFS0_check = -1;
00359                     }
00360                 }
00361                 if ( N_lgst_cluster <= FFS_i ) FFS0_check=-1;
00362             }
00363             if ( N_lgst_cluster <= FFS_i ) FFS0_check=-1;
00364         }
00365         else
00366         {
00367             if ( FFS_i_minus > step0 + totalsteps || N_lgst_cluster <= FFS_i )
00368                 FFS0_check = -1;
00369         }
00370         
00371         if (FFS0_check<0) return 1;
00372         else return 0;
00373     }
00374 }
00375