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

00001 /*
00002   meam-lenosky.cpp
00003   by Wei Cai  caiwei@mit.edu
00004   Last Modified : Mon Jan  1 21:27:48 2007
00005 
00006   FUNCTION  : MD++ with Lenosky potential for Si
00007 
00008   Note: The Lenosky potential is for Si only
00009         Only rectangular supercell is allowed
00010         Virial stress is not computed
00011 */
00012 
00013 #include "meam-lenosky.h"
00014 
00015 void MEAMFrame::potential()
00016 {
00017     MEAM();
00018 }
00019 
00020 void MEAMFrame::MEAM()
00021 {
00022     double alat[3], coord, ener_var, coord_var, eps;
00023 
00024     alat[0]=_H[0][0];
00025     alat[1]=_H[1][1];
00026     alat[2]=_H[2][2];
00027 
00028     SHtoR();
00029 
00030     eps=1e-8;
00031     if ( (fabs(_H[0][1])>eps)||(fabs(_H[0][2])>eps)||
00032          (fabs(_H[1][0])>eps)||(fabs(_H[1][2])>eps)||
00033          (fabs(_H[2][0])>eps)||(fabs(_H[2][1])>eps) )
00034     {
00035         WARNING(" Only rectangular supercell is allowed!");
00036     }
00037     lenosky_(&_NP, alat, (double *)_R, (double *)_F, _EPOT_IND, &_EPOT,
00038              &coord, &ener_var, &coord_var, &count);
00039 }
00040 
00041 void MEAMFrame::Alloc()
00042 {
00043     MDPARALLELFrame::Alloc();
00044 }    
00045 
00046 void MEAMFrame::initvars()
00047 {
00048     MDPARALLELFrame::initvars();
00049 }
00050 
00051 void MEAMFrame::initparser()
00052 {
00053     MDPARALLELFrame::initparser();
00054 }
00055 
00056 int MEAMFrame::exec(char *name)
00057 {
00058     if(MDPARALLELFrame::exec(name)==0) return 0;
00059     return -1;
00060 }
00061    
00062 #ifdef _TEST
00063 
00064 /* Main Program Begins */
00065 class MEAMFrame sim;
00066 
00067 /* The main program is defined here */
00068 #include "main.cpp"
00069 
00070 #endif//_TEST
00071 
00072 
00073 
00074 
00075 
00076 
00077 
00078 
00079 
00080 
00081 
00082 
00083 
00084 
00085 
00086 
00087 
00088 
00089 
00090 
00091 
00092 
00093 
00094 
00095 
00096 
00097 
00098 
00099 
00100 
00101 
00102 
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 /* old code */
00134 #if 0
00135 
00136 inline double interp(double func[],double deriv[],double dr,int ind,double qq)
00137 {
00138     double f, a, b, c, d, f1, p1, A1, A2, dr2, dr3, qq2, qq3;
00139 //    f = func[ind] + qq*deriv[ind];
00140     dr2=dr*dr; dr3=dr2*dr;
00141     qq2=qq*qq; qq3=qq2*qq;
00142     a = func[ind];
00143     b = deriv[ind];
00144     f1 = func[ind+1];
00145     p1 = deriv[ind+1];
00146     A1 = f1-a-b*dr;
00147     A2 = (p1-b)*dr;
00148     d = (A2-2*A1)/dr3;
00149     c = (3*A1-A2)/dr2;
00150     f=a+b*qq+c*qq2+d*qq3;
00151     return f;
00152 }
00153     
00154     
00155 inline double interp1(double func[],double deriv[],double dr,int ind,double qq)
00156 {
00157     double fp, a, b, c, d, f1, p1, A1, A2, dr2, dr3, qq2, qq3;
00158 //    fp = deriv[ind] + qq1/dr*(deriv[ind+1]-deriv[ind]);
00159     dr2=dr*dr; dr3=dr2*dr;
00160     qq2=qq*qq; qq3=qq2*qq;
00161     a = func[ind];
00162     b = deriv[ind];
00163     f1 = func[ind+1];
00164     p1 = deriv[ind+1];
00165     A1 = f1-a-b*dr;
00166     A2 = (p1-b)*dr;
00167     d = (A2-2*A1)/dr3;
00168     c = (3*A1-A2)/dr2;
00169     fp=b+2*c*qq+3*d*qq2;
00170     return fp;
00171 }
00172 #endif