MD++: My Documents/Codes/MD++.svn/trunk/src/meam-marian.cpp Source File

00001 /*
00002   meam-marian.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Mon Jan  1 21:28:21 2007
00005 
00006   FUNCTION  : MD++ with MEAM potential (using MDCASK F90 Codes from Jaime Marian)
00007 */
00008 
00009 #include "meam-marian.h"
00010 
00011 
00012 void MEAMFrame::Alloc()
00013 {
00014     MDPARALLELFrame::Alloc();
00015 }    
00016 
00017 void MEAMFrame::potential()
00018 {
00019     MEAM();
00020 }
00021 
00022 void MEAMFrame::readMEAM()
00023 {
00024     meamdefine_();
00025     _RLIST=meamread_.rcutmeam;
00026     _SKIN=0;
00027 }
00028 
00029 void MEAMFrame::MEAM()
00030 {
00031     Vector3 h;
00032     int i;
00033     double eps;
00034     
00035     /*refreshneighborlist();*/
00036     h=_H.height();
00037     linklist_.nlcx=(int)floor(h[0]/((_RLIST)*1.05));if(linklist_.nlcx==0)linklist_.nlcx=1;
00038     linklist_.nlcy=(int)floor(h[1]/((_RLIST)*1.05));if(linklist_.nlcy==0)linklist_.nlcy=1;
00039     linklist_.nlcz=(int)floor(h[2]/((_RLIST)*1.05));if(linklist_.nlcz==0)linklist_.nlcz=1;
00040     linklist_.nlc=linklist_.nlcx*linklist_.nlcy*linklist_.nlcz;
00041     
00042     //INFO_Printf("nlcx=%d nlcy=%d nlcz=%d\n",linklist_.nlcx,linklist_.nlcy,linklist_.nlcz);
00043     
00044     np_.nm=_NP;
00045     for(i=0;i<_NP;i++)
00046     {
00047         adim_.x0[i]=_SR[i].x;
00048         adim_.y0[i]=_SR[i].y;
00049         adim_.z0[i]=_SR[i].z;
00050     }
00051     box_.Lx=_H[0][0];
00052     box_.Ly=_H[1][1];
00053     box_.Lz=_H[2][2];
00054 
00055     eps=1e-8;
00056     if ( (fabs(_H[0][1])>eps)||(fabs(_H[0][2])>eps)||
00057          (fabs(_H[1][0])>eps)||(fabs(_H[1][2])>eps)||
00058          (fabs(_H[2][0])>eps)||(fabs(_H[2][1])>eps) )
00059     {
00060         WARNING(" Only rectangular supercell is allowed!");
00061     }
00062     
00063     //INFO_Printf("np_.nm=%d\n",np_.nm);
00064     forces_();
00065 
00066     for(i=0;i<_NP;i++)
00067     {
00068         _F[i].set(area6_.fx[i],area6_.fy[i],area6_.fz[i]);
00069         _EPOT_IND[i] = ppote_.atpe[i]+ppote_.atpe3b[i];        
00070     }
00071 
00072     _EPOT=ppote_.pe+ppote_.petrip;
00073     //INFO_Printf("p2=%20.12e p3=%20.12e\n",ppote_.pe,ppote_.petrip);
00074 }
00075 
00076 void MEAMFrame::initvars()
00077 {
00078     MDPARALLELFrame::initvars();
00079 }
00080 
00081 void MEAMFrame::initparser()
00082 {
00083     MDPARALLELFrame::initparser();
00084     bindvar("meamfile",meamdatafile_.meamfile,STRING);
00085 }
00086 
00087 int MEAMFrame::exec(char *name)
00088 {
00089     if(MDPARALLELFrame::exec(name)==0) return 0;
00090     bindcommand(name,"readMEAM",readMEAM());
00091     return -1;
00092 }
00093    
00094 #ifdef _TEST
00095 
00096 /* Main Program Begins */
00097 class MEAMFrame sim;
00098 
00099 /* The main program is defined here */
00100 #include "main.cpp"
00101 
00102 #endif//_TEST
00103 
00104 
00105 
00106 
00107 
00108 
00109 
00110 
00111 
00112 
00113 
00114 
00115 
00116 
00117 
00118 
00119 
00120 
00121 
00122 
00123 
00124 
00125 
00126 
00127 
00128 
00129 
00130 
00131 
00132 
00133 
00134 
00135 
00136 
00137 
00138 
00139 
00140 
00141 
00142 
00143 
00144 
00145 
00146 
00147 
00148 
00149 
00150 
00151 
00152 
00153 
00154 
00155 
00156 
00157 
00158 
00159 
00160 
00161 
00162 
00163 /* old code */
00164 #if 0
00165 
00166 inline double interp(double func[],double deriv[],double dr,int ind,double qq)
00167 {
00168     double f, a, b, c, d, f1, p1, A1, A2, dr2, dr3, qq2, qq3;
00169 //    f = func[ind] + qq*deriv[ind];
00170     dr2=dr*dr; dr3=dr2*dr;
00171     qq2=qq*qq; qq3=qq2*qq;
00172     a = func[ind];
00173     b = deriv[ind];
00174     f1 = func[ind+1];
00175     p1 = deriv[ind+1];
00176     A1 = f1-a-b*dr;
00177     A2 = (p1-b)*dr;
00178     d = (A2-2*A1)/dr3;
00179     c = (3*A1-A2)/dr2;
00180     f=a+b*qq+c*qq2+d*qq3;
00181     return f;
00182 }
00183     
00184     
00185 inline double interp1(double func[],double deriv[],double dr,int ind,double qq)
00186 {
00187     double fp, a, b, c, d, f1, p1, A1, A2, dr2, dr3, qq2, qq3;
00188 //    fp = deriv[ind] + qq1/dr*(deriv[ind+1]-deriv[ind]);
00189     dr2=dr*dr; dr3=dr2*dr;
00190     qq2=qq*qq; qq3=qq2*qq;
00191     a = func[ind];
00192     b = deriv[ind];
00193     f1 = func[ind+1];
00194     p1 = deriv[ind+1];
00195     A1 = f1-a-b*dr;
00196     A2 = (p1-b)*dr;
00197     d = (A2-2*A1)/dr3;
00198     c = (3*A1-A2)/dr2;
00199     fp=b+2*c*qq+3*d*qq2;
00200     return fp;
00201 }
00202 #endif