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

00001 // Last Modified : Fri Apr  4 11:10:35 2008
00002 
00003 #include "stringmethod.h"
00004 
00005 void MDFrame::runstringmethod()
00006 {/*
00007   * The string method
00008   * J. of Chemical Physics, vol.126, 164103 (2007)
00009   * Apr. 01 2008, Keonowook Kang
00010   * Must constrain whole atoms.
00011   */
00012     int i, j, k, m, n, ipt, size;
00013     Vector3 ds, ds0, dr, dt, dR, **_RcStore;
00014     Matrix33 hinv;
00015     double s, y_csp, dRnorm, dRnorm2, Xi0, XiN, s_j, hi, w, wbar, yi, yim1;
00016     double *Ec, *Lc;
00017     void *c;
00018     FILE *fp;
00019 
00020     INFO("The String Method");
00021 
00022     if(_SR1 == NULL || _SR2 == NULL)
00023     {
00024         if(_SR1 == NULL)
00025             ERROR("config1 is not set!");
00026         else
00027             ERROR("config2 is not set!");
00028         return;
00029     }
00030 
00031     if(_Rc==NULL || _Fc==NULL)
00032     {
00033         AllocChain();
00034         caldscom12();
00035         initRchain();
00036     }
00037     
00038     n=constrainatoms[0]; 
00039     Ec=(double *)malloc(sizeof(double)*(_CHAINLENGTH+1));
00040     Lc=(double *)malloc(sizeof(double)*(_CHAINLENGTH+1));
00041     size=sizeof(double *)*(3*n) + sizeof(double)*(3*n)*(_CHAINLENGTH+1) + 10;
00042 
00043     size=sizeof(Vector3 *)*(_CHAINLENGTH+1)+sizeof(Vector3)*(_CHAINLENGTH+1)*n + 10;
00044     c=malloc(size); memset(c,0,size);
00045     _RcStore=(Vector3 **) c;   _RcStore[0]=(Vector3 *)(_RcStore+(_CHAINLENGTH+1));
00046     for(i=1;i<=_CHAINLENGTH;i++) _RcStore[i]=_RcStore[i-1]+n;
00047 
00048     fp=fopen("stringeng.out","w");
00049 
00050     /* compute current constrain value */
00051     constrain_dist=0; s=0;
00052     for(i=0;i<n;i++)
00053     {
00054         ipt=constrainatoms[i+1];
00055         ds0=(_SR2[ipt]-dscom12)-_SR1[ipt];
00056         constrain_dist+=ds0.norm2();
00057     }
00058     INFO_Printf("string method: constrain_dist=%e\n",constrain_dist);
00059 
00060     if(nebspec[0]==1)
00061     {/* fix constrained atoms for relaxing surrounding atoms */
00062         for(i=0;i<n;i++)
00063         {
00064             ipt=constrainatoms[i+1];
00065             fixed[ipt]=1;
00066         }
00067     }
00068     
00069     LaVectorDouble DD(_CHAINLENGTH+1);
00070     LaVectorDouble DL(_CHAINLENGTH);
00071     LaVectorDouble DU(_CHAINLENGTH);
00072     LaVectorDouble b(_CHAINLENGTH+1), X(_CHAINLENGTH+1);
00073     LaGenMatDouble sigma(_CHAINLENGTH+1,3*n);
00074     LaTridiagFactDouble Tf;
00075     
00076 #define get_Xi0(y0,y1,y2,y3) protect(\
00077         Xi0 = DL(0)*( y0/(_nebinterp[0]-_nebinterp[1])/(_nebinterp[0]-_nebinterp[2])/(_nebinterp[0]-_nebinterp[3]) \
00078                     + y1/(_nebinterp[1]-_nebinterp[0])/(_nebinterp[1]-_nebinterp[2])/(_nebinterp[1]-_nebinterp[3]) \
00079                     + y2/(_nebinterp[2]-_nebinterp[0])/(_nebinterp[2]-_nebinterp[1])/(_nebinterp[2]-_nebinterp[3]) \
00080                     + y3/(_nebinterp[3]-_nebinterp[0])/(_nebinterp[3]-_nebinterp[1])/(_nebinterp[3]-_nebinterp[2])); \
00081         )
00082 #define get_XiN(y_nm3,y_nm2,y_nm1,y_n) protect(\
00083         XiN = DU(_CHAINLENGTH-1)*( y_nm3/(_nebinterp[_CHAINLENGTH-3]-_nebinterp[_CHAINLENGTH-2])/(_nebinterp[_CHAINLENGTH-3]-_nebinterp[_CHAINLENGTH-1])/(_nebinterp[_CHAINLENGTH-3]-_nebinterp[_CHAINLENGTH]) \
00084                     + y_nm2/(_nebinterp[_CHAINLENGTH-2]-_nebinterp[_CHAINLENGTH-3])/(_nebinterp[_CHAINLENGTH-2]-_nebinterp[_CHAINLENGTH-1])/(_nebinterp[_CHAINLENGTH-2]-_nebinterp[_CHAINLENGTH]) \
00085                     + y_nm1/(_nebinterp[_CHAINLENGTH-1]-_nebinterp[_CHAINLENGTH-3])/(_nebinterp[_CHAINLENGTH-1]-_nebinterp[_CHAINLENGTH-2])/(_nebinterp[_CHAINLENGTH-1]-_nebinterp[_CHAINLENGTH]) \
00086                     + y_n/(_nebinterp[_CHAINLENGTH]-_nebinterp[_CHAINLENGTH-3])/(_nebinterp[_CHAINLENGTH]-_nebinterp[_CHAINLENGTH-2])/(_nebinterp[_CHAINLENGTH]-_nebinterp[_CHAINLENGTH-1])); \
00087         )
00088 #define get_cubicspline(h_i,w,wbar,sigma_i,sigma_im1,y_i,y_im1) protect(\
00089         y_csp = w*y_i + wbar*y_im1 + SQR(h_i)*((pow(w,3)-w)*sigma_i + (pow(wbar,3)-wbar)*sigma_im1); \
00090         )
00091                
00092     step0 = curstep;
00093     for(curstep=step0;curstep<(step0 + totalsteps);curstep++)    
00094     {
00095         /* get forces */
00096         for(j=0;j<=_CHAINLENGTH;j++)
00097         {
00098             s=_nebinterp[j];
00099 
00100             interpCN(s); /* interpolate all atoms for a given s */
00101             /* constrained atoms at chain value */
00102             copyRchaintoCN(j);
00103 
00104             if(nebspec[0]==1)
00105             {
00106                 INFO("Relaxation for the point "<<j<<" in the energy hill");
00107                 relax(); 
00108             }
00109 
00110             call_potential();
00111             Ec[j]=_EPOT;
00112 
00113             for(i=0;i<n;i++)
00114             {
00115                 ipt=constrainatoms[i+1];
00116                 _Fc[j][i]=_F[ipt];
00117             }
00118         }
00119 
00120         if(curstep%printfreq==0) 
00121         {
00122             INFO_Printf("curstep = %d\n",curstep);
00123             for(j=0;j<=_CHAINLENGTH;j++)
00124                 INFO_Printf("%20.12e %25.15e %25.15e\n",_nebinterp[j],Ec[j]-Ec[0]);
00125             fprintf(fp,"%d ",curstep);
00126             for(j=0;j<=_CHAINLENGTH;j++)
00127                 fprintf(fp,"%20.12e %25.15e ",_nebinterp[j], Ec[j]-Ec[0]);
00128             fprintf(fp,"\n");
00129             fflush(fp);            
00130         }
00131         
00132         /* move along force - This step is more like steepest descent relaxation. */
00133         for(j=0;j<=_CHAINLENGTH;j++)
00134            for(i=0;i<n;i++)
00135            {
00136                 _Rc[j][i]+=_Fc[j][i]*_TIMESTEP;
00137            }
00138 
00139         copyRchaintoCN(0); setconfig1();
00140         copyRchaintoCN(_CHAINLENGTH); setconfig2();
00141         caldscom12();
00142         
00143         /* compute current constrain value */
00144         constrain_dist=0;
00145         for(i=0;i<n;i++)
00146         {
00147             ipt=constrainatoms[i+1];
00148             ds0=(_SR2[ipt]-dscom12)-_SR1[ipt];
00149             constrain_dist+=ds0.norm2();
00150         }
00151         INFO_Printf("string method: constrain_dist=%e\n",constrain_dist);
00152         
00153         /* calculate current reaction coordinates s */
00154         Lc[0]=0;
00155         for(j=1;j<=_CHAINLENGTH;j++)
00156         {
00157             dRnorm2 = 0;
00158             for(i=0;i<n;i++)
00159             {
00160                 dR = _Rc[j][i]-_Rc[j-1][i];
00161                 dRnorm2 += dR.norm2();
00162             }
00163             dRnorm = sqrt(dRnorm2);
00164             Lc[j]=Lc[j-1]+dRnorm;
00165         }
00166         _nebinterp[0]=0;
00167         for(j=1;j<=_CHAINLENGTH;j++)
00168             _nebinterp[j] = Lc[j]/Lc[_CHAINLENGTH];
00169 
00170         /* compute DD, DL and DU which are tri-diagonals of Matrix T */
00171         DD(0)=-1; DU(0)=1;
00172         for(j=1;j<_CHAINLENGTH;j++)
00173         {
00174             DL(j-1)=_nebinterp[j]-_nebinterp[j-1];
00175             DU(j)=_nebinterp[j+1]-_nebinterp[j];
00176             DD(j)=2*(DL(j-1)+DU(j));
00177         }
00178         DL(_CHAINLENGTH-1)=-1; DD(_CHAINLENGTH)=1;
00179         LaTridiagMatDouble T(DD,DL,DU);
00180         LaTridiagMatFactorize (T, Tf);
00181         
00182         for(k=0;k<n;k++)
00183             for(i=0;i<3;i++)
00184             {/* compute the vector b, T*sigma=b */
00185                 get_Xi0(_Rc[0][k][i],_Rc[1][k][i],_Rc[2][k][i],_Rc[3][k][i]);
00186                 b(0) = Xi0;
00187                 for(j=1;j<_CHAINLENGTH;j++)
00188                     b(j)=(_Rc[j+1][k][i]-_Rc[j][k][i])/DU(j) - (_Rc[j][k][i]-_Rc[j-1][k][i])/DL(j-1);
00189                 get_XiN(_Rc[_CHAINLENGTH-3][k][i],_Rc[_CHAINLENGTH-2][k][i],_Rc[_CHAINLENGTH-1][k][i],_Rc[_CHAINLENGTH][k][i]);
00190                 b(_CHAINLENGTH)=XiN;
00191             /* store sigma */
00192                 LaLinearSolve(Tf,X,b);
00193                 sigma(LaIndex(),3*k+i).inject(X);
00194             }
00195 
00196         for(k=0;k<n;k++)
00197             for(j=0;j<_CHAINLENGTH;j++)
00198                 _RcStore[j][k] = _Rc[j][k];
00199         
00200         /* redistribut string evenly using spline interpolation
00201            - parametrization by equal arc length */
00202         for(j=1;j<_CHAINLENGTH;j++)
00203         {
00204             s_j = j*1.0/_CHAINLENGTH;
00205             for(i=j;i<_CHAINLENGTH;i++)
00206                 if(_nebinterp[i]>s_j)
00207                 {
00208                     hi = _nebinterp[i]-_nebinterp[i-1];
00209                     w = (s_j - _nebinterp[i-1])/hi; wbar = 1-w;
00210                     for(m=0;m<n;m++)
00211                         for(k=0;k<3;k++)
00212                         {
00213                             yi = _RcStore[i][m][k]; yim1 = _RcStore[i-1][m][k];
00214                             get_cubicspline(hi,w,wbar,sigma(3*m+k,i),sigma(3*m+k,i-1),yi,yim1);
00215                             _Rc[j][m][k]=y_csp;
00216                         }
00217                     
00218                     break;
00219                 }
00220         }
00221 
00222         for(j=1;j<_CHAINLENGTH;j++)
00223             _nebinterp[j]=j*1.0/_CHAINLENGTH;
00224     }
00225 
00226     INFO_Printf("curstep = %d\n",curstep);
00227     for(j=0;j<=_CHAINLENGTH;j++)
00228         INFO_Printf("%20.12e %25.15e %25.15e\n",_nebinterp[j],Ec[j]-Ec[0]);
00229     fprintf(fp,"%d ",curstep);
00230     for(j=0;j<=_CHAINLENGTH;j++)
00231         fprintf(fp,"%20.12e %25.15e ",_nebinterp[j], Ec[j]-Ec[0]);
00232     fprintf(fp,"\n");
00233     fflush(fp); 
00234 
00235     free(Ec);free(Lc);free(_RcStore);
00236     fclose(fp);    
00237 }
00238 
00239