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

00001 /*
00002   rods.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Wed May  2 19:03:18 2007
00005 
00006   FUNCTION  :  MD simulation package of Lennard-Jones (LJ) potential
00007                LJ parameters depend on species
00008                Atoms can also form covalent bond with neighbors
00009 
00010   Test cases: scripts/work/ar/lj2-md.script, ljbond.tcl
00011 */
00012 
00013 #include "rods.h"
00014 
00015 void RODSFrame::initparser()
00016 {
00017     MDPARALLELFrame::initparser();
00018     bindvar("bond_k",&BOND_K,DOUBLE);
00019     bindvar("repul",&REPUL,DOUBLE);
00020 }
00021 
00022 int RODSFrame::exec(char *name)
00023 {
00024     if(MDPARALLELFrame::exec(name)==0) return 0;
00025     bindcommand(name,"init_structure",init_structure());
00026     bindcommand(name,"final_structure",final_structure());
00027     bindcommand(name,"connect_nodes",connect_nodes());
00028     bindcommand(name,"zero_com_rotation",zero_com_rotation());
00029     return -1;
00030 }
00031 
00032 void RODSFrame::Alloc()
00033 {
00034     int size;
00035     MDPARALLELFrame::Alloc();
00036 
00037     size = _NP*allocmultiple;    
00038     Realloc(_Fext,Vector3,size);
00039     bindvar("Fext",_Fext,DOUBLE);
00040 }
00041 
00042 void RODSFrame::Alloc_Bonds()
00043 {
00044     int size;
00045     size = NUM_BONDS;
00046     Realloc(BOND_R0, double,size);
00047     Realloc(BOND_INDEX, int,size*2);
00048 
00049     memset(BOND_R0, 0,sizeof(double)*size);
00050     memset(BOND_INDEX,  0,sizeof(int)*size*2);
00051 }
00052 
00053 
00054 void RODSFrame::initvars()
00055 {
00056     _RLIST=1.1;
00057     _SKIN =1.0;
00058     MDPARALLELFrame::initvars();
00059 }
00060 
00061 void RODSFrame::init_structure()
00062 {
00063     _NP = 24;
00064     
00065     Alloc();
00066 
00067     /* Assign positions of the nodes */
00068     _R[ 0].set(-1,  -1, 0);
00069     _R[ 1].set( 0,  -1, 0);
00070     _R[ 2].set( 1,  -1, 0);
00071     _R[ 3].set(-1,-0.6, 0);
00072     _R[ 4].set( 0,-0.4, 0);
00073     _R[ 5].set( 1,-0.6, 0);
00074     _R[ 6].set(-1,   0, 0);
00075     _R[ 7].set( 0,   0, 0);
00076     _R[ 8].set( 1,   0, 0);
00077     _R[ 9].set(-1, 0.6, 0);
00078     _R[10].set( 0, 0.4, 0);
00079     _R[11].set( 1, 0.6, 0);
00080     _R[12].set(-1,   1, 0);
00081     _R[13].set( 0,   1, 0);
00082     _R[14].set( 1,   1, 0);
00083 
00084     _R[15].set(-1,-0.6, 0);
00085     _R[16].set( 0,-0.4, 0);
00086     _R[17].set( 1,-0.6, 0);
00087     _R[18].set(-1,   0, 0);
00088     _R[19].set( 0,   0, 0);
00089     _R[20].set( 1,   0, 0);
00090     _R[21].set(-1, 0.6, 0);
00091     _R[22].set( 0, 0.4, 0);
00092     _R[23].set( 1, 0.6, 0);
00093     
00094     _H.clear(); _H[0][0]=4; _H[1][1]=4; _H[2][2]=4;
00095 
00096     RHtoS();
00097 
00098     connect_nodes();
00099 }
00100 
00101 void RODSFrame::final_structure()
00102 {
00103     _NP = 24;
00104     
00105      /* Assign positions of the nodes */
00106     _R[ 0].set(-1,  -0.6, 0);
00107     _R[ 1].set( 0,  -0.4, 0);
00108     _R[ 2].set( 1,  -0.6, 0);
00109     _R[ 3].set(-1,-0.6, 0.4);
00110     _R[ 4].set( 0,-0.4, 0.6);
00111     _R[ 5].set( 1,-0.6, 0.4);
00112     _R[ 6].set(-1,   0, 0.4);
00113     _R[ 7].set( 0,   0, 0.6);
00114     _R[ 8].set( 1,   0, 0.4);
00115     _R[ 9].set(-1, 0.6, 0.4);
00116     _R[10].set( 0, 0.4, 0.6);
00117     _R[11].set( 1, 0.6, 0.4);
00118     _R[12].set(-1,   0.6, 0);
00119     _R[13].set( 0,   0.4, 0);
00120     _R[14].set( 1,   0.6, 0);
00121 
00122     _R[15].set(-1,-0.6, -0.4);
00123     _R[16].set( 0,-0.4, -0.6);
00124     _R[17].set( 1,-0.6, -0.4);
00125     _R[18].set(-1,   0, -0.4);
00126     _R[19].set( 0,   0, -0.6);
00127     _R[20].set( 1,   0, -0.4);
00128     _R[21].set(-1, 0.6, -0.4);
00129     _R[22].set( 0, 0.4, -0.6);
00130     _R[23].set( 1, 0.6, -0.4);
00131     
00132     _H.clear(); _H[0][0]=4; _H[1][1]=4; _H[2][2]=4;
00133 
00134     RHtoS();
00135 }
00136 
00137 void RODSFrame::connect_nodes()
00138 {
00139     int i;
00140     Vector3 dr;
00141     int bond_index_data[112] = { 0,1,
00142                                  1,2,
00143                                  0,3,
00144                                  0,4,
00145                                  1,4,
00146                                  2,4,
00147                                  2,5,
00148                                  3,4,
00149                                  4,5,
00150                                  3,6,
00151                                  3,7,
00152                                  4,7,
00153                                  5,7,
00154                                  5,8,
00155                                  6,7,
00156                                  7,8,
00157                                  6,9,
00158                                  7,9,
00159                                  7,10,
00160                                  7,11,
00161                                  8,11,
00162                                  9,10,
00163                                  10,11,
00164                                  9,12,
00165                                  10,12,
00166                                  10,13,
00167                                  10,14,
00168                                  11,14,
00169                                  12,13,
00170                                  13,14,
00171                                  0,15,
00172                                  0,16,
00173                                  1,16,
00174                                  2,16,
00175                                  2,17,
00176                                  15,16,
00177                                  16,17,
00178                                  15,18,
00179                                  15,19,
00180                                  16,19,
00181                                  17,19,
00182                                  17,20,
00183                                  18,19,
00184                                  19,20,
00185                                  18,21,
00186                                  19,21,
00187                                  19,22,
00188                                  19,23,
00189                                  20,23,
00190                                  21,22,
00191                                  22,23,
00192                                  21,12,
00193                                  22,12,
00194                                  22,13,
00195                                  22,14,
00196                                  23,14 };
00197     
00198     NUM_BONDS = 56;
00199     Alloc_Bonds();
00200     
00201     /* Assign bond indices */
00202     for(i=0;i<NUM_BONDS;i++)
00203     {
00204         BOND_INDEX[i*2]   = bond_index_data[i*2];
00205         BOND_INDEX[i*2+1] = bond_index_data[i*2+1];
00206     }
00207     
00208     /* Compute distances of the bonds */
00209     SHtoR();
00210     for(i=0;i<NUM_BONDS;i++)
00211     {
00212         dr = _R[BOND_INDEX[i*2]]-_R[BOND_INDEX[i*2+1]];
00213         BOND_R0[i] = dr.norm();
00214         INFO_Printf("BOND_R0[%d]=%20.10e\n",i,BOND_R0[i]);
00215     }
00216 
00217 }
00218 
00219 void RODSFrame::calcprop()
00220 {
00221     MDPARALLELFrame::calcprop();
00222 }
00223 
00224 void RODSFrame::rods_potential()
00225 {
00226     int i,ipt,jpt;
00227     double U,r,r2;
00228     Vector3 rij, fij;
00229     
00230     DUMP(HIG"Lennard Jones"NOR);
00231         
00232     refreshneighborlist();
00233     
00234     _EPOT=0;
00235 
00236     for(i=0;i<_NP;i++)
00237     {_F[i].clear(); _EPOT_IND[i]=0;}
00238     _VIRIAL.clear();
00239 
00240     SHtoR();
00241     for(i=0;i<NUM_BONDS;i++)
00242     {
00243         ipt = BOND_INDEX[i*2];
00244         jpt = BOND_INDEX[i*2+1];
00245         rij = _R[jpt] - _R[ipt];
00246 
00247         r2=rij.norm2();
00248         r=sqrt(r2);
00249                 
00250         U = 0.5* BOND_K * SQR(r-BOND_R0[i]);
00251         fij = rij * (-1.0*BOND_K*(1-BOND_R0[i]/r));
00252 
00253         _F[ipt]-=fij;
00254         _F[jpt]+=fij;
00255         _EPOT_IND[ipt]+=U*0.5;
00256         _EPOT_IND[jpt]+=U*0.5;
00257         _EPOT+=U;
00258         _VIRIAL.addnvv(1.,fij,rij);
00259     }
00260 
00261     /* Adding repulsive force between all atom pairs */
00262     for(ipt=0;ipt<_NP;ipt++)
00263         for(jpt=ipt+1;jpt<_NP;jpt++)
00264         {
00265             rij = _R[jpt] - _R[ipt];
00266 
00267             r2=rij.norm2();
00268             
00269             U = -0.5*REPUL* r2;
00270             fij = rij*REPUL;
00271 
00272             _F[ipt]-=fij;
00273             _F[jpt]+=fij;
00274             _EPOT_IND[ipt]+=U*0.5;
00275             _EPOT_IND[jpt]+=U*0.5;
00276             _EPOT+=U;
00277             _VIRIAL.addnvv(1.,fij,rij);
00278         }
00279     
00280     /* external force */
00281     for(i=0;i<_NP;i++)
00282     {
00283         if(!fixed[i])
00284         {
00285             _F[i] += _Fext[i];
00286             _EPOT -= dot(_Fext[i],_R[i]);
00287         }
00288     }
00289 
00290     /* fixed atoms */
00291     for(i=0;i<_NP;i++)
00292     {
00293         if(fixed[i]) _F[i].clear();
00294     }
00295     
00296 }
00297 
00298 
00299 void RODSFrame::potential()
00300 {
00301     rods_potential();
00302 }
00303 
00304 void RODSFrame::zero_com_rotation()
00305 {
00306     int i, n;
00307     double Jx, Jy, Jz;
00308     Matrix33 hinv;
00309     
00310     for(i=0;i<_NP;i++)
00311     {
00312         _R[i]=_H*_SR[i];
00313         _VR[i]=_H*_VSR[i];
00314     }
00315 
00316     /* computing COM angular momentum */
00317     Jz = 0;  n = 0;
00318     for(i=0;i<_NP;i++)
00319         if(!fixed[i])
00320         {
00321             Jz += _R[i].x*_VR[i].y - _R[i].y*_VR[i].x;
00322             n++;
00323         }
00324     Jz /= n;
00325 
00326     Jx = 0;  n = 0;
00327     for(i=0;i<_NP;i++)
00328         if(!fixed[i])
00329         {
00330             Jx += _R[i].y*_VR[i].z - _R[i].z*_VR[i].y;
00331             n++;
00332         }
00333     Jx /= n;
00334 
00335     Jy = 0;  n = 0;
00336     for(i=0;i<_NP;i++)
00337         if(!fixed[i])
00338         {
00339             Jy += _R[i].z*_VR[i].x - _R[i].x*_VR[i].z;
00340             n++;
00341         }
00342     Jy /= n;
00343     
00344     for(i=0;i<_NP;i++)
00345         if(!fixed[i])
00346         {
00347             _VR[i].y -= Jz*_R[i].x/(_R[i].x*_R[i].x+_R[i].y*_R[i].y);
00348             _VR[i].x += Jz*_R[i].y/(_R[i].x*_R[i].x+_R[i].y*_R[i].y);
00349         }
00350     
00351 #if 0 /* not working yet, need to debug */
00352     for(i=0;i<_NP;i++)
00353         if(!fixed[i])
00354         {
00355             _VR[i].z -= Jx*_R[i].y/(_R[i].y*_R[i].y+_R[i].z*_R[i].z);
00356             _VR[i].y += Jx*_R[i].z/(_R[i].y*_R[i].y+_R[i].z*_R[i].z);
00357         }
00358 
00359     for(i=0;i<_NP;i++)
00360         if(!fixed[i])
00361         {
00362             _VR[i].x -= Jy*_R[i].z/(_R[i].x*_R[i].x+_R[i].z*_R[i].z);
00363             _VR[i].z += Jy*_R[i].x/(_R[i].x*_R[i].x+_R[i].z*_R[i].z);
00364         }
00365 #endif
00366     
00367     
00368     hinv = _H.inv();
00369     for(i=0;i<_NP;i++)
00370     {
00371         if(!fixed[i])
00372         {
00373             _SR[i]=hinv*_R[i];
00374             _VSR[i]=hinv*_VR[i];
00375         }
00376     }
00377 }
00378 
00379 void RODSFrame::plot()
00380 {
00381     int i,ipt,jpt;
00382     double L;
00383     char s[100];
00384     double x1,y1,z1,x2,y2,z2,dx,dy,dz,dr;
00385     unsigned ce;
00386     Vector3 s1, s2, vr1, vr2, sri, srj, ri, rj;
00387 
00388     if(win==NULL) return;
00389     if(!(win->alive)) return;
00390     
00391     L=max(_H[0][0],_H[1][1]);
00392     L=max(L,_H[2][2])*.5;
00393 
00394     switch(plot_color_axis){
00395      case(0): color_ind=_EPOT_IND; break;
00396      case(1): color_ind=_EPOT_IND; break;
00397      case(2): color_ind=_TOPOL; break;
00398      case(3): color_ind=_TOPOL; break;
00399      default: ERROR("plot() unknown coloraxis "<<plot_color_axis); return;
00400     }
00401     
00402     SHtoR();
00403     win->Lock();
00404     win->Clear();
00405     
00406     /* draw atom s*/
00407     for(i=0;i<_NP;i++)
00408     {
00409         sri=_SR[i];
00410         if(plot_map_pbc==1) sri.subint();
00411         ri = _H*sri;
00412                 
00413         if (!plot_atom_info) s[0]=0;
00414         else
00415         {   /* when an atom is clicked */
00416             if(plot_atom_info==1)    /* print scaled coordinates of the atom */
00417                 sprintf(s,"%d,%6.4f,%6.4f,%6.4f",
00418                         i,sri.x,sri.y,sri.z);
00419             else if(plot_atom_info==2) /* print real coordinates of the atom */
00420                 sprintf(s,"%d,%6.4f,%6.4f,%6.4f",
00421                         i,ri.x,ri.y,ri.z);
00422             else if(plot_atom_info==3) /* print color_ind, fixed, species info of the atom */
00423                 sprintf(s,"%d,%8.6e %d %d",
00424                         i,color_ind[i], fixed[i], species[i]);
00425             else if(plot_atom_info==4) /* print force on the atom */
00426                 sprintf(s,"%d,%6.4f,%6.4f,%6.4f",
00427                         i,_F[i].x,_F[i].y,_F[i].z);
00428             else
00429                 s[0]=0;
00430         }
00431         /* otherwise color atoms according to their species */        
00432         ce=colors[MAXCOLORS+5+species[i]];
00433         if(plot_atom_info==5) /* print atom species and color code if clicked */
00434             sprintf(s,"%d,%d,%x",i,species[i],ce);
00435         
00436         win->DrawPoint(ri.x/L,ri.y/L,ri.z/L,
00437                        atomradius[species[i]]/L,ce,s,2);
00438     }
00439         
00440     /* Draw Bonds */
00441     for(i=0;i<NUM_BONDS;i++)
00442     {
00443         ipt = BOND_INDEX[i*2];
00444         jpt = BOND_INDEX[i*2+1];
00445         ri = _R[ipt];
00446         rj = _R[jpt];
00447 
00448         x1=ri.x/L;y1=ri.y/L;z1=ri.z/L;
00449         x2=rj.x/L;y2=rj.y/L;z2=rj.z/L;
00450         dx=x2-x1;dy=y2-y1;dz=z2-z1;dr=sqrt(dx*dx+dy*dy+dz*dz);
00451         dx/=dr;dy/=dr;dz/=dr;
00452         win->DrawLine(x1+dx*atomradius[species[ipt]]/L,
00453                       y1+dy*atomradius[species[ipt]]/L,
00454                       z1+dz*atomradius[species[ipt]]/L,
00455                       x2-dx*atomradius[species[jpt]]/L,
00456                       y2-dy*atomradius[species[jpt]]/L,
00457                       z2-dz*atomradius[species[jpt]]/L,
00458                       colors[MAXCOLORS+1],bondradius/L,1);
00459     }
00460 
00461     /* draw frame */
00462 #define drawsline(a,b,c,d,e,f,g,h,i) s1.set(a,b,c); s2.set(d,e,f);\
00463         vr1=_H*s1; vr2=_H*s2; vr1/=2*L; vr2/=2*L;\
00464         win->DrawLine(vr1.x,vr1.y,vr1.z,vr2.x,vr2.y,vr2.z,g,h,i);
00465 #if 0    
00466     drawsline(-1,-1,-1,-1,-1, 1,colors[MAXCOLORS+2],0,0);
00467     drawsline(-1,-1, 1,-1, 1, 1,colors[MAXCOLORS+2],0,0);
00468     drawsline(-1, 1, 1,-1, 1,-1,colors[MAXCOLORS+2],0,0);
00469     drawsline(-1, 1,-1,-1,-1,-1,colors[MAXCOLORS+2],0,0);
00470     drawsline( 1,-1,-1, 1,-1, 1,colors[MAXCOLORS+2],0,0);
00471     drawsline( 1,-1, 1, 1, 1, 1,colors[MAXCOLORS+2],0,0);
00472     drawsline( 1, 1, 1, 1, 1,-1,colors[MAXCOLORS+2],0,0);
00473     drawsline( 1, 1,-1, 1,-1,-1,colors[MAXCOLORS+2],0,0);
00474     drawsline(-1,-1,-1, 1,-1,-1,colors[MAXCOLORS+2],0,0);
00475     drawsline(-1,-1, 1, 1,-1, 1,colors[MAXCOLORS+2],0,0);
00476     drawsline(-1, 1, 1, 1, 1, 1,colors[MAXCOLORS+2],0,0);
00477     drawsline(-1, 1,-1, 1, 1,-1,colors[MAXCOLORS+2],0,0);
00478 #endif
00479     
00480     win->Unlock();
00481     win->Refresh();
00482 }
00483 
00484     
00485 
00486     
00487 #ifdef _TEST
00488 
00489 /* Main Program Begins */
00490 class RODSFrame sim;
00491 
00492 /* The main program is defined here */
00493 #include "main.cpp"
00494 
00495 #endif//_TEST
00496