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

00001 /*
00002   aluminum.cpp  translated from aluminum.f
00003   
00004   by Wei Cai  caiwei@mit.edu, Maurice de Koning maurice@mit.edu
00005   Last Modified : Mon Feb  4 09:58:29 2002
00006 
00007   FUNCTION  : MD++ implementation of Al Glue potential
00008               (Ercolessi-Adams Glue potential)
00009 */
00010 
00011 void seval(int n,double u,double *x, double *y,
00012            double *b, double *c, double *d,
00013            double *f, double *df, double *d2f);
00014 
00015 void v2(double arg, double *func, double *dfunc, double *d2func)
00016 {
00017 #define nv2 17
00018     double xv2[nv2]=
00019     {
00020          .202111069753385e+01  ,
00021          .227374953472558e+01  ,
00022          .252638837191732e+01  ,
00023          .277902720910905e+01  ,
00024          .303166604630078e+01  ,
00025          .328430488349251e+01  ,
00026          .353694372068424e+01  ,
00027          .378958255787597e+01  ,
00028          .404222139506771e+01  ,
00029          .429486023225944e+01  ,
00030          .454749906945117e+01  ,
00031          .480013790664290e+01  ,
00032          .505277674383463e+01  ,
00033          .530541558102636e+01  ,
00034          .555805441821810e+01  ,
00035          .555807968210182e+01  ,
00036          .555810494598553e+01  ,
00037     };
00038     double yv2[nv2]=
00039     {
00040          .196016472197158e+01  ,
00041          .682724240745344e+00  ,
00042          .147370824539188e+00  ,
00043         -.188188235860390e-01  ,
00044         -.576011902692490e-01  ,
00045         -.519846499644276e-01  ,
00046         -.376352484845919e-01  ,
00047         -.373737879689433e-01  ,
00048         -.531351030124350e-01  ,
00049         -.632864983555742e-01  ,
00050         -.548103623840369e-01  ,
00051         -.372889232343935e-01  ,
00052         -.188876517630154e-01  ,
00053         -.585239362533525e-02  ,
00054          .000000000000000e+00  ,
00055          .000000000000000e+00  ,
00056          .000000000000000e+00  ,
00057     };
00058     double bv2[nv2]=
00059     {
00060         -.702739315585347e+01  ,
00061         -.333140549270729e+01  ,
00062         -.117329394261502e+01  ,
00063         -.306003283486901e+00  ,
00064         -.366656699104026e-01  ,
00065          .588330899204400e-01  ,
00066          .384220572312032e-01  ,
00067         -.390223173707191e-01  ,
00068         -.663882722510521e-01  ,
00069         -.312918894386669e-02  ,
00070          .590118945294245e-01  ,
00071          .757939459148246e-01  ,
00072          .643822548468606e-01  ,
00073          .399750987463792e-01  ,
00074          .177103852679117e-05  ,
00075         -.590423369301474e-06  ,
00076          .590654950414731e-06  ,
00077     };
00078     double cv2[nv2]=
00079     {
00080          .877545959718548e+01  ,
00081          .585407125495837e+01  ,
00082          .268820820643116e+01  ,
00083          .744718689404422e+00  ,
00084          .321378734769888e+00  ,
00085          .566263292669091e-01  ,
00086         -.137417679148505e+00  ,
00087         -.169124163201523e+00  ,
00088          .608037039066423e-01  ,
00089          .189589640245655e+00  ,
00090          .563784150384640e-01  ,
00091          .100486298765028e-01  ,
00092         -.552186092621482e-01  ,
00093         -.413902746758285e-01  ,
00094         -.116832934994489e+00  ,
00095          .233610871054729e-01  ,
00096          .233885865725971e-01  ,
00097     };
00098     double dv2[nv2]=
00099     {
00100         -.385449887634130e+01  ,
00101         -.417706040200591e+01  ,
00102         -.256425277368288e+01  ,
00103         -.558557503589276e+00  ,
00104         -.349316054551627e+00  ,
00105         -.256022933201611e+00  ,
00106         -.418337423301704e-01  ,
00107          .303368330939646e+00  ,
00108          .169921006301015e+00  ,
00109         -.175759761362548e+00  ,
00110         -.611278214082881e-01  ,
00111         -.861140219824535e-01  ,
00112          .182451950513387e-01  ,
00113         -.995395392057973e-01  ,
00114          .184972909229936e+04  ,
00115          .362829766922787e+00  ,
00116          .362829766922787e+00  ,
00117     };
00118     double argmax;
00119 
00120     argmax=.555805441821810e+01;
00121     
00122     if(arg>=argmax)
00123     {
00124         *func = *dfunc = *d2func = 0;
00125     }
00126     else
00127     {
00128         seval(nv2,arg,xv2,yv2,bv2,cv2,dv2,func,dfunc,d2func);
00129     }
00130 }
00131 
00132 
00133 
00134 
00135 void rh(double arg, double *func, double *dfunc, double *d2func)
00136 {
00137 #define nrh 17
00138     double xrh[nrh]=
00139     {
00140          .202111069753385e+01  ,
00141          .227374953472558e+01  ,
00142          .252638837191732e+01  ,
00143          .277902720910905e+01  ,
00144          .303166604630078e+01  ,
00145          .328430488349251e+01  ,
00146          .353694372068424e+01  ,
00147          .378958255787597e+01  ,
00148          .404222139506771e+01  ,
00149          .429486023225944e+01  ,
00150          .454749906945117e+01  ,
00151          .480013790664290e+01  ,
00152          .505277674383463e+01  ,
00153          .530541558102636e+01  ,
00154          .555805441821810e+01  ,
00155          .555807968210182e+01  ,
00156          .555810494598553e+01  ,
00157     };
00158     double yrh[nrh]=
00159     {
00160          .865674623712589e-01  ,
00161          .925214702944478e-01  ,
00162          .862003123832002e-01  ,
00163          .762736292751052e-01  ,
00164          .606481841271735e-01  ,
00165          .466030959588197e-01  ,
00166          .338740138848363e-01  ,
00167          .232572661705343e-01  ,
00168          .109046405489829e-01  ,
00169          .524910605677597e-02  ,
00170          .391702419142291e-02  ,
00171          .308277776293383e-02  ,
00172          .250214745349505e-02  ,
00173          .147220513798186e-02  ,
00174          .000000000000000e+00  ,
00175          .000000000000000e+00  ,
00176          .000000000000000e+00  ,
00177     };
00178     double brh[nrh]=
00179     {
00180          .608555214104682e-01  ,
00181         -.800158928716306e-02  ,
00182         -.332089451111092e-01  ,
00183         -.521001991705069e-01  ,
00184         -.618130637429111e-01  ,
00185         -.529750064268036e-01  ,
00186         -.442210477548108e-01  ,
00187         -.473645664984640e-01  ,
00188         -.390741582571631e-01  ,
00189         -.101795580610560e-01  ,
00190         -.318316981110289e-02  ,
00191         -.281217210746153e-02  ,
00192         -.236932031483360e-02  ,
00193         -.683554708271547e-02  ,
00194         -.638718204858808e-06  ,
00195          .212925486831149e-06  ,
00196         -.212983742465787e-06  ,
00197     };
00198     double crh[nrh]=
00199     {
00200         -.170233687052940e+00  ,
00201         -.102317878901959e+00  ,
00202          .254162872544396e-02  ,
00203         -.773173610292656e-01  ,
00204          .388717099948882e-01  ,
00205         -.388873819867093e-02  ,
00206          .385388290924526e-01  ,
00207         -.509815666327127e-01  ,
00208          .837968231208082e-01  ,
00209          .305743500420042e-01  ,
00210         -.288110886134041e-02  ,
00211          .434959924771674e-02  ,
00212         -.259669459714693e-02  ,
00213         -.150816117849093e-01  ,
00214          .421356801161513e-01  ,
00215         -.842575249165724e-02  ,
00216         -.843267014952237e-02  ,
00217     };
00218     double drh[nrh]=
00219     {
00220          .896085612514625e-01  ,
00221          .138352319847830e+00  ,
00222         -.105366473134009e+00  ,
00223          .153300619856764e+00  ,
00224         -.564184148788224e-01  ,
00225          .559792096400504e-01  ,
00226         -.118113795329664e+00  ,
00227          .177827488509794e+00  ,
00228         -.702220789044304e-01  ,
00229         -.441413511810337e-01  ,
00230          .954024354744484e-02  ,
00231         -.916498550800407e-02  ,
00232         -.164726813535368e-01  ,
00233          .754928689733184e-01  ,
00234         -.667110847110954e+03  ,
00235         -.912720300911022e-01  ,
00236         -.912720300911022e-01  ,
00237     };
00238     double argmax;
00239 
00240     argmax=.555805441821810e+01;
00241     
00242     if(arg>=argmax)
00243     {
00244         *func = *dfunc = *d2func = 0;
00245     }
00246     else
00247     {
00248         seval(nrh,arg,xrh,yrh,brh,crh,drh,func,dfunc,d2func);
00249     }
00250 }
00251 
00252 
00253 
00254 
00255 
00256 void uu(double arg, double *func, double *dfunc, double *d2func)
00257 {
00258 #define nuu 13
00259     double xuu[nuu]=
00260     {
00261          .000000000000000e+00  ,
00262          .100000000000000e+00  ,
00263          .200000000000000e+00  ,
00264          .300000000000000e+00  ,
00265          .400000000000000e+00  ,
00266          .500000000000000e+00  ,
00267          .600000000000000e+00  ,
00268          .700000000000000e+00  ,
00269          .800000000000000e+00  ,
00270          .900000000000000e+00  ,
00271          .100000000000000e+01  ,
00272          .110000000000000e+01  ,
00273          .120000000000000e+01  ,
00274     };
00275     double yuu[nuu]=
00276     {
00277          .000000000000000e+00  ,
00278         -.113953324143752e+01  ,
00279         -.145709859805864e+01  ,
00280         -.174913308002738e+01  ,
00281         -.202960322136630e+01  ,
00282         -.225202324967546e+01  ,
00283         -.242723053979436e+01  ,
00284         -.255171976467357e+01  ,
00285         -.260521638832322e+01  ,
00286         -.264397894381693e+01  ,
00287         -.265707884842034e+01  ,
00288         -.264564149400021e+01  ,
00289         -.260870604452106e+01  ,
00290     };
00291     double buu[nuu]=
00292     {
00293         -.183757286015853e+02  ,
00294         -.574233124410516e+01  ,
00295         -.236790436375322e+01  ,
00296         -.307404645857774e+01  ,
00297         -.251104850116555e+01  ,
00298         -.196846462620234e+01  ,
00299         -.154391254686695e+01  ,
00300         -.846780636273251e+00  ,
00301         -.408540363905760e+00  ,
00302         -.286833282404628e+00  ,
00303         -.309389414590161e-06  ,
00304          .236958014464143e+00  ,
00305          .503352368511243e+00  ,
00306     };
00307     double cuu[nuu]=
00308     {
00309          .830779120415016e+02  ,
00310          .432560615333001e+02  ,
00311         -.951179272978074e+01  ,
00312          .245037178153561e+01  ,
00313          .317960779258630e+01  ,
00314          .224623095704576e+01  ,
00315          .199928983630817e+01  ,
00316          .497202926962879e+01  ,
00317         -.589626545953876e+00  ,
00318          .180669736096520e+01  ,
00319          .106163236918694e+01  ,
00320          .130795086934864e+01  ,
00321          .135599267112235e+01  ,
00322     };
00323     double duu[nuu]=
00324     {
00325         -.132739501694005e+03  ,
00326         -.175892847543603e+03  ,
00327          .398738817043878e+02  ,
00328          .243078670350231e+01  ,
00329         -.311125611846847e+01  ,
00330         -.823137069125319e+00  ,
00331          .990913144440207e+01  ,
00332         -.185388527186089e+02  ,
00333          .798774635639692e+01  ,
00334         -.248354997259420e+01  ,
00335          .821061667205675e+00  ,
00336          .160139339245701e+00  ,
00337          .160139339245701e+00  ,
00338     };        
00339     double argmin;
00340 
00341     argmin=0;
00342     
00343     if(arg<=argmin)
00344     {
00345         *func = *dfunc = *d2func = 0;
00346     }
00347     else
00348     {
00349         seval(nuu,arg,xuu,yuu,buu,cuu,duu,func,dfunc,d2func);
00350     }
00351 
00352 }
00353 
00354 
00355 void seval(int n,double u,double *x, double *y,
00356            double *b, double *c, double *d,
00357            double *f, double *df, double *d2f)
00358 {
00359     int i, j, k;
00360     double dx;
00361 
00362     i=0;
00363     if ( i>=n ) i = 0;    
00364     if ( u<x[i] ) goto line10;
00365     if ( u<=x[i+1] ) goto line30;
00366  line10:
00367     i = 0;
00368     j = n;
00369  line20:
00370     k = (i+j)/2;
00371     if ( u<x[k] ) j = k;
00372     if ( u>=x[k] ) i = k;
00373     if ( j>(i+1) ) goto line20;
00374  line30:
00375     dx = u - x[i];
00376     *f = y[i] + dx*(b[i] + dx*(c[i] + dx*d[i]));
00377     *df = b[i] + dx*(2.0*c[i] + 3.0*dx*d[i]);
00378     *d2f = 2.0*c[i] + 6.0*dx*d[i];
00379 }