1 #include <omp.h>
    2 /******************************************************************************/
    3 /*                                                                            */
    4 /*  Copyright (c) 1996 University of California at Santa Barbara              */
    5 /*                                                                            */
    6 /*  All rights reserved.                                                      */
    7 /*                                                                            */
    8 /*  Permission is given to use, copy, and modify this software for any        */
    9 /*  non-commercial purpose as long as this copyright notice is not removed    */
   10 /*  All other uses, including redistribution in whole or in part, are         */
   11 /*  forbidden without prior written permission.                               */
   12 /*                                                                            */
   13 /*  This software is provided with absolutely no warranty and no support.     */
   14 /*                                                                            */
   15 /*  For further information contact: pedro@cs.ucsb.edu or martin@cs.ucsb.edu  */
   16 /*                                                                            */
   17 /******************************************************************************/
   18 
   19 // -------- WATER_H -----------------------------------------------------
   20 
   21 
   22 using namespace std;
   23 
   24 #include <stdio.h>
   25 #include <math.h>
   26 #include <stdlib.h>
   27 #include <iostream>
   28 #include <fstream>
   29 
   30 // -------- CONSTANT_H ---------------------------------------------------
   31 
   32 const double ONE 	= (double) 1.00;
   33 const double TWO	= (double) 2.00;
   34 const double THREE	= (double) 3.00;
   35 const double FIVE	= (double) 5.00;
   36 
   37 const double ROH = 0.9572;
   38 const double ROHI = (1.00/ROH);
   39 const double ROHI2 = (ROHI*ROHI);
   40 const double ANGLE = 1.824218;
   41 const double OMAS = 15.99945;
   42 const double HMAS = 1.007825;
   43 const double WTMOL = (OMAS+2.00*HMAS);
   44 
   45 /*.....UNITS USeD TO SCALe VARIABLeS (IN C.G.S.) */
   46 
   47 const double UNITT = (1.0e-15);
   48 const double UNITL = (1.0e-8);
   49 const double UNITM = (1.6605655e-24);
   50 const double BOLTZ = (1.380662e-16);
   51 const double AVGNO = (6.0220045e23);
   52 
   53 /* .....FORCe CONSTANTS SCALeD(DIVIDeD) BY (UNITM/UNITT**2) */
   54 
   55 const double FC11  = (0.512596);
   56 const double FC33  = (0.048098);
   57 const double FC12 = (-0.0058230);
   58 const double FC13  = (0.016452);
   59 const double FC111 = (-0.57191);
   60 const double FC333 = (-0.007636);
   61 const double FC112 = (-0.001867);
   62 const double FC113 = (-0.002047);
   63 const double FC123 = (-0.03083);
   64 const double FC133 = (-0.0094245);
   65 const double FC1111  = (0.8431);
   66 const double FC3333 = (-0.00193);
   67 const double FC1112 = (-0.0030);
   68 const double FC1122  = (0.0036);
   69 const double FC1113 = (-0.012);
   70 const double FC1123  = (0.0060);
   71 const double FC1133 = (-0.0048);
   72 const double FC1233  = (0.0211);
   73 const double FC1333  = (0.006263);
   74 
   75 /*.....WATeR-WATeR INTeRACTION PARAMeTeRS */
   76 
   77 const double QQ = 0.07152158;
   78 const double A1 = 455.313100;
   79 const double B1 = 5.15271070;
   80 const double A2 = 0.27879839;
   81 const double B2 = 2.76084370;
   82 const double A3 = 0.60895706;
   83 const double B3 = 2.96189550;
   84 const double A4 = 0.11447336;
   85 const double B4 = 2.23326410;
   86 const double CM = 0.45682590;
   87 const double AB1 = (A1*B1);
   88 const double AB2 = (A2*B2);
   89 const double AB3 = (A3*B3);
   90 const double AB4 = (A4*B4);
   91 const double C1 = (1.00-CM);
   92 const double C2 = (0.50*CM);
   93 const double QQ2 = (2.00*QQ);
   94 const double QQ4 = (2.00*QQ2);
   95 
   96 // ------ EXTERN_H ---------
   97 
   98 
   99 extern "C" void get_time(int *v);
  100 extern "C" void get_ticks_time(int *v);
  101 
  102 inline double sign(double a, double b){
  103   if(b < 0) {
  104     if(a < 0)
  105       return a;
  106     return -a;
  107   } else {
  108     if(a < 0)
  109       return -a;
  110     return a;
  111   }
  112 }
  113 
  114 // ------ VECTOR_H ---------
  115 
  116 #define _OPENMP_
  117 #include "water.h"
  118 
  119 // ----------------------------------------------------------------------------
  120 
  121 inline void tvecClr(double *v) { int i; for(i=0; i < NDIM; i++) v[i] = 0.0; }
  122 
  123 inline void tvecAdd(double *v1, double *v2) { int i;
  124 for(i=0; i < NDIM; i++) v1[i] += v2[i];
  125 }
  126 
  127 inline void tvecSub(double *v1, double *v2) { int i;
  128 for(i=0; i < NDIM; i++) v1[i] -= v2[i];
  129 }
  130 
  131 inline void tvecLoad(double *v1, double *v2) { int i;
  132 for(i=0; i < NDIM; i++) v1[i] = v2[i];
  133 }
  134 
  135 inline void tvecStore(double *v1, double *v2) { int i;
  136 for(i=0; i < NDIM; i++) v2[i] = v1[i];
  137 }
  138 
  139 inline void tvecDiv(double *v1, double *v2) { int i;
  140 for(i=0; i < NDIM; i++) v1[i] /= v2[i];
  141 }
  142 
  143 inline void tvecProd(double *v1, double *v2) { int i;
  144 for(i=0; i < NDIM; i++) v1[i] *= v2[i];
  145 }
  146 
  147 inline void tvecScale(double *v, double s) { int i;
  148 for(i=0; i < NDIM; i++) v[i] *= s;
  149 }
  150 
  151 inline double tvecDotProd(double *v1,double *v2) {
  152   int i; double d = 0.0;
  153   for(i=0; i < NDIM; i++)
  154     d += (v1[i] * v2[i]);
  155   return d;
  156 }
  157 
  158 inline void tvecMin(double *v1,double  *v2){ int i;
  159 for(i=0; i < NDIM; i++)
  160   if(v1[i] > v2[i]) v1[i] = v2[i];
  161 }
  162 
  163 inline double tvecNorm1(double *v){ int i; double d=0.0;
  164 for(i=0; i < NDIM; i++) d += v[i];
  165 return d;
  166 }
  167 
  168 // ------ H2O_H ---------
  169 
  170 const int NDIR = 3;             /* XDIR, YDIR, ZDIR             */
  171 const int XDIR = 0;
  172 const int YDIR = 1;
  173 const int ZDIR = 2;
  174 
  175 const int NATOM = 3;		/* H1, Oxy, H2 			*/
  176 
  177 const int MAXODR = 7;		/* Predictor/corrector order	*/
  178 const int MAXODR2 = 9;		// (MAXODR+2);
  179 
  180 const int DISP   = 0;		/* Displacement (position)      */
  181 const int VEL    = 1;           /* Velocity                     */
  182 const int ACC    = 2;           /* Acceleration                 */
  183 const int DER_3  = 3;		/* Higher order derivatives ... */
  184 const int DER_4  = 4;
  185 const int DER_5  = 5;
  186 const int DER_6  = 6;
  187 const int FORCES = 7;
  188 
  189 // -----------------------------------------------------------------------------
  190 // ATOM
  191 // -----------------------------------------------------------------------------
  192 class atom {
  193   vector M[MAXODR2];
  194 public:
  195   atom() { }
  196   ~atom() { }
  197 vector *getM(int idx)		{ return &M[idx]; }
  198 //void setM(int idx, double *v)	{ M[idx].vecLoad(v); }
  199 //void addM(int idx, double *v)	{ M[idx].vecAdd(v); }
  200 //void clearM(int idx) 		{ M[idx].vecClr(); }
  201 void clearAtom() 		{
  202   int i;
  203   for(i=0; i < MAXODR2;i++)
  204    M[i].vecClr();
  205 }
  206 void predic(int n, double *c);
  207 void correc(int n, double *c);
  208 };
  209 
  210 // -----------------------------------------------------------------------------
  211 // WATER
  212 // -----------------------------------------------------------------------------
  213 class skratch_pad;
  214 class h2o {
  215 public:
  216   atom H1;
  217   atom O;
  218   atom H2;
  219   vector V;
  220   h2o() { }
  221   ~h2o() { }
  222 void  clear();
  223 
  224 //vector *getH1M(int idx)   { return &(H1.M[idx]); }
  225 //vector *getOM(int idx)    { return &(O.M[idx]); }
  226 //vector *getH2M(int idx)   { return &(H2.M[idx]); }
  227 
  228 void loadH1Pos(double *v) { H1.getM(DISP)->vecLoad(v); }
  229 void loadOPos(double *v)  {  O.getM(DISP)->vecLoad(v); }
  230 void loadH2Pos(double *v) { H2.getM(DISP)->vecLoad(v); }
  231 
  232 void storeH1Pos(double *v) { H1.getM(DISP)->vecStore(v); }
  233 void storeOPos(double *v)  {  O.getM(DISP)->vecStore(v); }
  234 void storeH2Pos(double *v) { H2.getM(DISP)->vecStore(v); }
  235 
  236 void loadH1Vel(double *v) { H1.getM(VEL)->vecLoad(v); }
  237 void loadOVel(double *v)  {  O.getM(VEL)->vecLoad(v); }
  238 void loadH2Vel(double *v) { H2.getM(VEL)->vecLoad(v); }
  239 
  240 void storeH1Vel(double *v) { H1.getM(VEL)->vecStore(v); }
  241 void storeOVel(double *v)  {  O.getM(VEL)->vecStore(v); }
  242 void storeH2Vel(double *v) { H2.getM(VEL)->vecStore(v); }
  243 
  244 //vector *getV()            { return &V; }
  245 void loadV(double *v)     { V.vecLoad(v); }
  246 void storeV(double *v)    { V.vecStore(v); }
  247 
  248 // --------------------------------------------------------------------
  249 void 	loadDirPos(int dir,double *v);
  250 void	storeDirVel(int dir,double *v);
  251 void	shiftAxis(int dir, double v);
  252 void	shift(double *v);
  253 // void	PosMin(vector *min);
  254 
  255 void	kineti(vector *v);
  256 void	bndry(double b);
  257 
  258 void	intra_poteng(vector *v);
  259 void	intraf();
  260 void    vir(acc_double *v);
  261 
  262 void	updateFields(int dest, skratch_pad *p);
  263 
  264 void	predic(int n, double *c);
  265 void	correc(int n, double *c);
  266 void	scaleMomenta(int Dest, double HM, double OM);
  267 void	dump();
  268 };
  269 
  270 // --------------------------------------------------------------------
  271 // PAD: to perform the accumulation of forces during INTERF
  272 // --------------------------------------------------------------------
  273 class skratch_pad {
  274 public:
  275   vector H1pos;	// input data
  276   vector  Opos;
  277   vector H2pos;
  278   vector    VM;
  279   vector H1force; // output data
  280   vector  Oforce;
  281   vector H2force;
  282   skratch_pad() { H1pos.vecClr(); Opos.vecClr(); H2pos.vecClr(); VM.vecClr();
  283 		H1force.vecClr(); Oforce.vecClr(); H2force.vecClr(); }
  284   ~skratch_pad() { }
  285   void storeVM(double *v)	{ VM.vecStore(v); }
  286   void storeH1pos(double *v) { H1pos.vecStore(v); }
  287   void storeOpos(double *v)  { Opos.vecStore(v); }
  288   void storeH2pos(double *v) { H2pos.vecStore(v); }
  289   void storeH1force(double *v) { H1force.vecStore(v); }
  290   void storeOforce(double *v)  { Oforce.vecStore(v); }
  291   void storeH2force(double *v) { H2force.vecStore(v); }
  292   void read_data(h2o *m);
  293   void update_forces(double Res[3][3]);
  294 };
  295 
  296 // ---- H2O.C -----------------------------------------------------
  297 void skratch_pad::read_data(h2o *m) {
  298   double tmp[3];
  299 
  300   m->storeH1Pos(tmp);
  301   H1pos.vecLoad(tmp);
  302   m->storeOPos(tmp);
  303   Opos.vecLoad(tmp);
  304   m->storeH2Pos(tmp);
  305   H2pos.vecLoad(tmp);
  306   m->storeV(tmp);
  307   VM.vecLoad(tmp);
  308   H1force.vecClr();
  309   Oforce.vecClr();
  310   H2force.vecClr();
  311 }
  312 
  313 void skratch_pad::update_forces(double Res[3][3]) {
  314   H1force.vecAdd(Res[0]);
  315   Oforce.vecAdd(Res[1]);
  316   H2force.vecAdd(Res[2]);
  317 }
  318 
  319 void h2o::clear(){
  320   H1.clearAtom();
  321    O.clearAtom();
  322   H2.clearAtom();
  323 }
  324 
  325 void h2o::predic(int n, double *c){
  326   H1.predic(n,c);
  327    O.predic(n,c);
  328   H2.predic(n,c);
  329 }
  330 
  331 void h2o::correc(int n, double *c){
  332   H1.correc(n,c);
  333    O.correc(n,c);
  334   H2.correc(n,c);
  335 }
  336 
  337 void h2o::scaleMomenta(int Dest, double HM, double OM){
  338   H1.getM(Dest)->vecScale(HM);
  339    O.getM(Dest)->vecScale(OM);
  340   H2.getM(Dest)->vecScale(HM);
  341 }
  342 
  343 void h2o::loadDirPos(int dir,double *v){
  344   H1.getM(DISP)->vecFieldSet(dir,v[0]);
  345    O.getM(DISP)->vecFieldSet(dir,v[1]);
  346   H2.getM(DISP)->vecFieldSet(dir,v[2]);
  347 }
  348 
  349 // ---------------------------------------------------------------------------
  350 void h2o::storeDirVel(int dir,double *v){
  351   v[0] = H1.getM(VEL)->vecField(dir);
  352   v[1] =  O.getM(VEL)->vecField(dir);
  353   v[2] = H2.getM(VEL)->vecField(dir);
  354 }
  355 // ---------------------------------------------------------------------------
  356 void h2o::shiftAxis(int dir, double v){
  357   H1.getM(DISP)->vecFieldAdd(dir,v);
  358    O.getM(DISP)->vecFieldAdd(dir,v);
  359   H2.getM(DISP)->vecFieldAdd(dir,v);
  360 }
  361 // ---------------------------------------------------------------------------
  362 void h2o::shift(double *v){
  363   H1.getM(DISP)->vecSub(v);
  364    O.getM(DISP)->vecSub(v);
  365   H2.getM(DISP)->vecSub(v);
  366 }
  367 // ---------------------------------------------------------------------------
  368 /*
  369 void h2o::PosMin(vector *min){
  370   min->vecMin(&(H1.M[DISP]));
  371   min->vecMin(&(O.M[DISP]));
  372   min->vecMin(&(H2.M[DISP]));
  373 }
  374 */
  375 // ---------------------------------------------------------------------------
  376 void h2o::kineti(vector *s){
  377   double T1[3], T2[3], v1[3], v2[3], v3[3];
  378   H1.getM(VEL)->vecStore(v1);
  379    O.getM(VEL)->vecStore(v2);
  380   H2.getM(VEL)->vecStore(v3);
  381   tvecLoad(T1,v1);
  382   tvecProd(T1,T1);
  383   tvecLoad(T2,v3);
  384   tvecProd(T2,T2);
  385   tvecAdd(T1,T2);
  386   tvecScale(T1,HMAS);
  387   tvecProd(v2,v2);
  388   tvecScale(v2,OMAS);
  389   tvecAdd(T1,v2);
  390   s->vecAdd(T1);
  391 }
  392 
  393 // ---------------------------------------------------------------------------
  394 void h2o::bndry(double b){
  395   int i;
  396   double t[3];
  397 
  398   O.getM(DISP)->vecStore(t);
  399   for(i = XDIR; i <= ZDIR; i++ ) {
  400     if(t[i] > b) {
  401       shiftAxis(i,-b);
  402     } else if (t[i] < 0.00) {
  403       shiftAxis(i,b);
  404     }
  405   }
  406 }
  407 
  408 // ---------------------------------------------------------------------------
  409 //  INTRA-MOLECULAR POTENTIAL ENERGY
  410 // ---------------------------------------------------------------------------
  411 void h2o::intra_poteng(vector *v){
  412   double LocPot, DTS, R1, R2, RX, COS, DT, DR1, DR2, DR1S, DR2S, DRP;
  413   double t1[3], t2[3], t3[3], t4[3], h1[3], Ox[3], h2[3], VM[3];
  414   double tmp[3];
  415 
  416   tvecClr(tmp);
  417   LocPot = 0.0;
  418 
  419   H1.getM(DISP)->vecStore(h1);
  420    O.getM(DISP)->vecStore(Ox);
  421   H2.getM(DISP)->vecStore(h2);
  422 
  423   // Compute VM = Ox.Pos * C1 + (H1.Pos + H2.Pos) *C2
  424   tvecLoad(t1,h1);
  425   tvecLoad(t2,Ox);
  426   tvecLoad(t3,h2);
  427   tvecScale(t2,C1);
  428   tvecAdd(t3,t1);
  429   tvecScale(t3,C2);
  430   tvecAdd(t2,t3);
  431   tvecLoad(VM,t2);
  432 
  433   // R1 = OxH1Potential();
  434   tvecLoad(t1,Ox);
  435   tvecSub(t1,h1);
  436   R1 = tvecDotProd(t1,t1);
  437 
  438   // R2 = OxH2Potential();
  439   tvecLoad(t1,Ox);
  440   tvecSub(t1,h2);
  441   R2 = tvecDotProd(t1,t1);
  442 
  443   // RX = H1xOxH2Potential();
  444   tvecLoad(t1,h1);
  445   tvecLoad(t2,Ox);
  446   tvecLoad(t3,h2);
  447   tvecLoad(t4,t2);
  448   tvecSub(t2,t1);
  449   tvecSub(t4,t3);
  450   RX = tvecDotProd(t2,t4);
  451 
  452   R1 = sqrt(R1);
  453   R2 = sqrt(R2);
  454   COS = RX/(R1*R2);
  455   DT = (acos(COS)-ANGLE)*ROH;
  456   DR1 = R1-ROH;
  457   DR2 = R2-ROH;
  458   DR1S = DR1*DR1;
  459   DR2S = DR2*DR2;
  460   DRP = DR1+DR2;
  461   DTS = DT*DT;
  462 
  463   LocPot += (FC11*(DR1S+DR2S)+FC33*DTS)*0.5+FC12*DR1*DR2+FC13*DRP*DT
  464      +(FC111*(DR1S*DR1+DR2S*DR2)+FC333*DTS*DT+FC112*DRP*DR1*DR2
  465      +FC113*(DR1S+DR2S)*DT+FC123*DR1*DR2*DT+FC133*DRP*DTS)*ROHI;
  466   LocPot += (FC1111*(DR1S*DR1S+DR2S*DR2S)+FC3333*DTS*DTS+
  467      FC1112*(DR1S+DR2S)*DR1*DR2+FC1122*DR1S*DR2S+
  468      FC1113*(DR1S*DR1+DR2S*DR2)*DT+FC1123*DRP*DR1*DR2*DT+
  469      FC1133*(DR1S+DR2S)*DTS+FC1233*DR1*DR2*DTS+
  470      FC1333*DRP*DTS*DT)*ROHI2;
  471 
  472   tmp[0] = LocPot;
  473   V.vecLoad(VM);
  474   v->vecAdd(tmp);
  475 }
  476 
  477 // ---------------------------------------------------------------------------
  478 void atom::predic(int norder, double *coeffs){
  479   int JIZ, JI, L, f;
  480   double S[3], T[3];
  481 
  482   JIZ=2;
  483   for (f = 0; f < norder; f++) {
  484     JI = JIZ;
  485     tvecClr(S);
  486     for (L = f; L < norder; L++) {
  487       M[L+1].vecStore(T);
  488       tvecScale(T,coeffs[JI]);
  489       tvecAdd(S,T);
  490       JI++;
  491     }
  492     M[f].vecAdd(S);
  493     JIZ += (norder+1);
  494   }
  495 }
  496 
  497 // ---------------------------------------------------------------------------
  498 void atom::correc(int norder, double *coeffs){
  499   int f;
  500   double S[3], T[3];
  501 
  502   M[FORCES].vecStore(S);
  503   M[ACC].vecStore(T);
  504   tvecSub(S,T);
  505   for(f = 0; f < (norder+1); f++){
  506     tvecLoad(T,S);
  507     tvecScale(T,coeffs[f]);
  508     M[f].vecAdd(T);
  509   }
  510 }
  511 
  512 // ------------------------------------------------------------------------
  513 // computeVIR:
  514 // ------------------------------------------------------------------------
  515 void h2o::vir(acc_double *v){
  516   double loc_vir;
  517   double tmp1[3], tmp2[3];
  518 
  519   loc_vir = 0.0;
  520   H1.getM(DISP)->vecStore(tmp1);
  521   H1.getM(FORCES)->vecStore(tmp2);
  522   loc_vir += tvecDotProd(tmp1,tmp2);
  523   O.getM(DISP)->vecStore(tmp1);
  524   O.getM(FORCES)->vecStore(tmp2);
  525   loc_vir += tvecDotProd(tmp1,tmp2);
  526   H2.getM(DISP)->vecStore(tmp1);
  527   H2.getM(FORCES)->vecStore(tmp2);
  528   loc_vir += tvecDotProd(tmp1,tmp2);
  529   v->addval(loc_vir);
  530 }
  531 
  532 
  533 // -----------------------------------------------------------------------
  534 // intraf: Sets V and FORCES fields in each ATOM
  535 // -----------------------------------------------------------------------
  536 void h2o::intraf(){
  537   double SUM, R1, R2, COS, SIN;
  538   double DT, DTS, DR1, DR1S, DR2, DR2S, R1S, R2S, F1, F2, F3;
  539   double vr1[3], vr2[3], dt1[3], dt3[3], dr11[3], dr23[3];
  540   double s[3], v1[3], v2[3], v3[3], h1[3], Ox[3], h2[3];
  541 
  542   SUM=0.0;
  543   R1=0.0;
  544   R2=0.0;
  545 
  546   tvecClr(s);
  547 
  548   H1.getM(DISP)->vecStore(h1);
  549    O.getM(DISP)->vecStore(Ox);
  550   H2.getM(DISP)->vecStore(h2);
  551 
  552   tvecLoad(v1,Ox);
  553   tvecScale(v1,C1);
  554   tvecLoad(v2,h1);
  555   tvecLoad(v3,h2);
  556   tvecAdd(v2,v3);
  557   tvecScale(v2,C2);
  558   tvecAdd(v1,v2);
  559   V.vecLoad(v1);
  560 
  561   tvecLoad(vr1,Ox);
  562   tvecLoad(v1,h1);
  563   tvecSub(vr1,v1);
  564   R1 = tvecDotProd(vr1,vr1);
  565   tvecLoad(vr2,Ox);
  566   tvecLoad(v2,h2);
  567   tvecSub(vr2,v2);
  568   R2 = tvecDotProd(vr2,vr2);
  569   SUM = tvecDotProd(vr1,vr2);
  570 
  571   R1=sqrt(R1);
  572   R2=sqrt(R2);
  573 
  574 /* CALCULATE COS(THETA), SIN(THETA), DELTA(R1), DELTA(R2), AND DELTA(THETA) */
  575   COS=SUM/(R1*R2);
  576   SIN=sqrt(ONE-COS*COS);
  577   DT=(acos(COS)-ANGLE)*ROH;
  578   DTS=DT*DT;
  579   DR1=R1-ROH;
  580   DR1S=DR1*DR1;
  581   DR2=R2-ROH;
  582   DR2S=DR2*DR2;
  583 
  584 /* CALCULATE DERIVATIVES OF R1/X1, R2/X3, THETA/X1, AND THETA/X3 */
  585   R1S=ROH/(R1*SIN);
  586   R2S=ROH/(R2*SIN);
  587   tvecLoad(dr11,vr1); tvecScale(dr11,(1.0/R1));
  588   tvecLoad(dr23,vr2); tvecScale(dr23,(1.0/R2));
  589   tvecLoad(dt1,dr11); tvecScale(dt1,COS);
  590   tvecSub(dt1,dr23); tvecScale(dt1,R1S);
  591   tvecLoad(dt3,dr23); tvecScale(dt3,COS);
  592   tvecSub(dt3,dr11); tvecScale(dt3,R2S);
  593 
  594 /* CALCULATE FORCES */
  595   F1=FC11*DR1+FC12*DR2+FC13*DT;
  596   F2=FC33*DT +FC13*(DR1+DR2);
  597   F3=FC11*DR2+FC12*DR1+FC13*DT;
  598   F1=F1+(3.0*FC111*DR1S+FC112*(2.0*DR1+DR2)*DR2
  599      +2.0*FC113*DR1*DT+FC123*DR2*DT+FC133*DTS)*ROHI;
  600   F2=F2+(3.0*FC333*DTS+FC113*(DR1S+DR2S)
  601      +FC123*DR1*DR2+2.0*FC133*(DR1+DR2)*DT)*ROHI;
  602   F3=F3+(3.0*FC111*DR2S+FC112*(2.0*DR2+DR1)*DR1
  603      +2.0*FC113*DR2*DT+FC123*DR1*DT+FC133*DTS)*ROHI;
  604   F1=F1+(4.0*FC1111*DR1S*DR1+FC1112*(3.0*DR1S+DR2S)
  605      *DR2+2.0*FC1122*DR1*DR2S+3.0*FC1113*DR1S*DT
  606      +FC1123*(2.0*DR1+DR2)*DR2*DT+(2.0*FC1133*DR1
  607      +FC1233*DR2+FC1333*DT)*DTS)*ROHI2;
  608   F2=F2+(4.0*FC3333*DTS*DT+FC1113*(DR1S*DR1+DR2S*DR2)
  609      +FC1123*(DR1+DR2)*DR1*DR2+2.0*FC1133*(DR1S+DR2S)
  610      *DT+2.0*FC1233*DR1*DR2*DT+3.0*FC1333*(DR1+DR2)*DTS)*ROHI2;
  611   F3=F3+(4.0*FC1111*DR2S*DR2+FC1112*(3.0*DR2S+DR1S)
  612      *DR1+2.0*FC1122*DR1S*DR2+3.0*FC1113*DR2S*DT
  613      +FC1123*(2.0*DR2+DR1)*DR1*DT+(2.0*FC1133*DR2
  614      +FC1233*DR1+FC1333*DT)*DTS)*ROHI2;
  615 
  616   tvecLoad(v1,dr11);
  617   tvecScale(v1,F1);
  618   tvecLoad(v2,dt1);
  619   tvecScale(v2,F2);
  620   tvecAdd(v1,v2);
  621   tvecLoad(v3,v1);
  622   H1.getM(FORCES)->vecLoad(v1);
  623 
  624   tvecLoad(v1,dr23);
  625   tvecScale(v1,F3);
  626   tvecLoad(v2,dt3);
  627   tvecScale(v2,F2);
  628   tvecAdd(v1,v2);
  629   H2.getM(FORCES)->vecLoad(v1);
  630 
  631   tvecAdd(v3,v1);
  632   tvecScale(v3,-1.00);
  633   O.getM(FORCES)->vecLoad(v3);
  634 }
  635 
  636 // --------------------------------------------------------------------
  637 void h2o::updateFields(int d, skratch_pad *p){
  638   double tmp[3];
  639    p->storeH1force(tmp);
  640    H1.getM(d)->vecAdd(tmp);
  641    p->storeOforce(tmp);
  642    O.getM(d)->vecAdd(tmp);
  643   p->storeH2force(tmp);
  644   H2.getM(d)->vecAdd(tmp);
  645 }
  646 
  647 // --------------------------------------------------------------------
  648 void h2o::dump(){
  649   int i;
  650   for(i=0; i <= FORCES; i++){
  651     H1.getM(i)->vecPrint();
  652     O.getM(i)->vecPrint();
  653     H2.getM(i)->vecPrint();
  654     printf("\n");
  655   }
  656   printf("\n V:");
  657   V.vecPrint();
  658   printf("\n");
  659 }
  660 
  661 // --- SIMPARM.H --------------------------------------------------------------
  662 class simparm {
  663   double 	TLC[100];		/* Taylor Series Coeffs		*/
  664   double	PCC[11];		/* Predictor/Corrector Coeffs	*/
  665 
  666   double        R1;
  667   double        ELPST;
  668 
  669   int	IRST;
  670   int	NVAR;
  671   int	NXYZ;
  672   int 	NXV;
  673   int 	IXF;
  674   int 	IYF;
  675   int 	IZF;
  676   int 	IMY;
  677   int 	IMZ;
  678 
  679   double TEMP;		/* Temperature			*/
  680   double RHO;		/* Density			*/
  681   double TSTEP;		/* Time step			*/
  682   double BOXL;		/* Box length			*/
  683   double BOXH;		/* Box height			*/
  684   double CUTOFF;	/* Cut-Off radius		*/
  685   double CUT2;		/* Square of Cut-Off radius	*/
  686 
  687   double FKIN;
  688   double FPOT;
  689 
  690   int NMOL;			/* Number of Molecules		*/
  691   int NORDER;			/* Integration order		*/
  692   int NATMO;			/* Number of Atoms		*/
  693   int NATMO3;
  694   int NMOL1;
  695   int NSTEP;			/* Number of time steps		*/
  696   int NSAVE;
  697   int NRST;
  698   int NPRINT;
  699   int NFMC;
  700   int I2;
  701 
  702   double FHM;
  703   double FOM;
  704   double REF1;
  705   double REF2;
  706   double REF4;
  707 
  708 public:
  709   simparm();
  710   ~simparm();
  711   void SYSCNS();
  712   void CNSTNT();
  713   void loadParms(const char *);
  714 
  715   double *getTLC()		{ return TLC; }
  716   double *getPCC()		{ return PCC; }
  717 
  718   double getR1()         { return R1; }
  719   double getELPST()      { return ELPST; }
  720 
  721   double getFKIN()      { return FKIN; }
  722   double getFPOT()      { return FPOT; }
  723 
  724   void setPS1(double *v) { R1 =v[0]; ELPST=v[1]; }
  725   void resetStat()       { ELPST=0.00; }
  726 
  727  int getIRST() 	{ return IRST; }
  728  int getNVAR() 	{ return NVAR; }
  729  int getNXYZ() 	{ return NXYZ; }
  730  int getNXV() 	{ return NXV; }
  731  int getIXF() 	{ return IXF; }
  732  int getIYF() 	{ return IYF; }
  733  int getIZF() 	{ return IZF; }
  734  int getIMY() 	{ return IMY; }
  735  int getIMZ() 	{ return IMZ; }
  736 
  737  void setPS2(int *v){ IRST=v[0]; NVAR=v[1]; NXYZ=v[2];
  738 NXV=v[3]; IXF=v[4]; IYF=v[5]; IZF=v[6]; IMY=v[7]; IMZ=v[8]; }
  739 
  740  double getTEMP()	{ return TEMP; }
  741  double getRHO() 	{ return RHO; }
  742  double getTSTEP()	{ return TSTEP; }
  743  double getBOXL()	{ return BOXL; }
  744  double getBOXH()	{ return BOXH; }
  745  double getCUTOFF()	{ return CUTOFF; }
  746  double getCUT2()  	{ return CUT2; }
  747 
  748  void setPS3(double *v) { TEMP =v[0]; RHO=v[1];
  749 TSTEP=v[2]; BOXL=v[3]; BOXH=v[4]; CUTOFF=v[5]; CUT2=v[6];}
  750 
  751  int getNMOL()		{ return  NMOL; }
  752  int getNORDER()	{ return  NORDER; }
  753  int getNATMO()		{ return  NATMO; }
  754  int getNATMO3()	{ return  NATMO3; }
  755  int getNMOL1()		{ return  NMOL1; }
  756  int getNSTEP()		{ return  NSTEP; }
  757  int getNSAVE()		{ return  NSAVE; }
  758  int getNRST()		{ return  NRST; }
  759  int getNPRINT()	{ return  NPRINT; }
  760  int getNFMC()		{ return  NFMC; }
  761  int getI2()		{ return  I2; }
  762 
  763  void setPS4(int *v){ NMOL=v[0]; NORDER=v[1]; NATMO=v[2];
  764 NATMO3=v[3]; NMOL1=v[4]; NSTEP=v[5]; NSAVE=v[6]; NRST=v[7]; NPRINT=v[8];
  765 NFMC=v[9]; I2=v[10]; }
  766  void setNFMC(int v) { NFMC = v; }
  767 
  768  double getFHM() { return FHM; }
  769  double getFOM() { return FOM; }
  770  double getREF1() { return REF1; }
  771  double getREF2() { return REF2; }
  772  double getREF4() { return REF4; }
  773 
  774  void   setPS5(double *v) {FHM=v[0];FOM=v[1];REF1=v[2]; REF2=v[3]; REF4=v[4]; }
  775 
  776  double computeFAC(){
  777    return BOLTZ*TEMP*NATMO/UNITM * pow((UNITT*TSTEP/UNITL),2.0);
  778  }
  779 
  780 // void	dump();
  781 
  782 };
  783 
  784 // --- SIMPARM.C --------------------------------------------------------------
  785 
  786 simparm::simparm(){
  787   int i;
  788 
  789   for(i=0; i < 100; i++)
  790     TLC[i] = 0.0;
  791   for(i=0; i < 11; i++)
  792     PCC[i] = 0.0;
  793   R1 = ELPST = 0.0;
  794   IRST = NVAR = NXYZ = NXV = IXF = IYF = IZF = IMY = IMZ = 0;
  795   TEMP = RHO = TSTEP = BOXL = BOXH = CUTOFF = CUT2 = 0.0;
  796   NMOL = NORDER = NATMO = NATMO3 = NMOL1 = NSTEP = NSAVE = 0;
  797   NRST = NPRINT = NFMC = I2 = 0;
  798   FHM = FOM = REF1 = REF2 = REF4 = 0.0;
  799 }
  800 
  801 void simparm::CNSTNT(){
  802   int NN,N1,K1,N;
  803   double TN, TK;
  804 
  805   N = NORDER+1;
  806   TLC[1] = 1.00;
  807   for (N1 = 2; N1<=N; N1++){
  808     NN = N1-1;
  809     TN = NN;
  810     TLC[N1] = 1.00;
  811     TK = 1.00;
  812     for (K1=2;K1<=N1;K1++) {
  813       TLC[(K1-1)*N+NN] = TLC[(K1-2)*N+NN+1]*TN/TK;
  814       NN = NN-1;
  815       TN = TN-ONE;
  816       TK = TK+ONE;
  817     }
  818   }
  819 
  820   PCC[2] = ONE;
  821   N1 = N-1;
  822   if((N1 == 1) || (N1 == 2)){
  823     cout << "***** ERROR: THE ORDER HAS TO BE GREATER THAN 2 ****\n";
  824   }
  825   if(N1 == 3){
  826     PCC[0] = ONE/6.00;
  827     PCC[1] = FIVE/6.00;
  828     PCC[3] = ONE/3.00;
  829   }
  830   if(N1 == 4){
  831     PCC[0] = (double) 19.00/120.00;
  832     PCC[1] = THREE/4.00;
  833     PCC[3] = ONE/2.00;
  834     PCC[4] = ONE/12.00;
  835   }
  836   if(N1 == 5){
  837     PCC[0] = THREE/20.00;
  838     PCC[1] = (double) 251.00/360.00;
  839     PCC[3] = (double) 11.00/18.00;
  840     PCC[4] = ONE/6.00;
  841     PCC[5] = ONE/60.00;
  842   }
  843   if(N1 ==  6){
  844 	PCC[0] = (double) 863.00/6048.00;
  845 	PCC[1] = (double) 665.00/1008.00;
  846 	PCC[3] = (double) 25.00/36.00;
  847 	PCC[4] = (double) 35.00/144.00;
  848 	PCC[5] = ONE/24.00;
  849 	PCC[6] = ONE/360.00;
  850   }
  851     if(N1 ==  7){
  852 	PCC[0] = (double) 275.00/2016.00;
  853 	PCC[1] = (double) 19087.00/30240.00;
  854 	PCC[3] = (double) 137.00/180.00;
  855 	PCC[4] = FIVE/16.00;
  856 	PCC[5] = (double) 17.00/240.00;
  857 	PCC[6] = ONE/120.00;
  858 	PCC[7] = ONE/2520.00;
  859     }
  860 }
  861 
  862 void simparm::SYSCNS() {
  863   TSTEP=TSTEP/UNITT;
  864   NATMO=NATOM*NMOL;
  865   NATMO3=NATMO*3;
  866   BOXL=pow((NMOL*WTMOL*UNITM/RHO),(1.00/3.00));
  867   BOXL=BOXL/UNITL;
  868   BOXH=BOXL*0.50;
  869   CUTOFF=max(BOXH,CUTOFF);
  870   REF1= -QQ/(CUTOFF*CUTOFF*CUTOFF);
  871   REF2=TWO*REF1;
  872   REF4=TWO*REF2;
  873   CUT2=CUTOFF*CUTOFF;
  874   FPOT= UNITM * pow((UNITL/UNITT),2.0) / (BOLTZ*TEMP*NATMO);
  875   FKIN=FPOT*0.50/(TSTEP*TSTEP);
  876   FHM=(TSTEP*TSTEP*0.50)/HMAS;
  877   FOM=(TSTEP*TSTEP*0.50)/OMAS;
  878   NMOL1=NMOL-1;
  879 }
  880 
  881 void simparm::loadParms(const char * file){
  882   FILE * infile;
  883 
  884   cout << "OUTPUT FOR PERFECT CLUB BENCHMARK: MDG Revision: 1.*  Author: kipp\n";
  885   NORDER = 5;
  886   TEMP = 298.0;
  887   RHO = 0.9980;
  888   TSTEP = 1.5e-16;
  889   CUTOFF = 0.0;
  890 
  891   infile = fopen(file, "r");
  892   if(infile == NULL){
  893     cout << " *** ERROR: Cannot open LWI5\n";
  894     exit(-1);
  895   }
  896   fscanf(infile,"%lf",&TSTEP);
  897   fscanf(infile,"%d", &NMOL);
  898   fscanf(infile,"%d", &NSTEP);
  899   fscanf(infile,"%d", &NORDER);
  900   fscanf(infile,"%d", &NSAVE);
  901   fscanf(infile,"%d", &NRST);
  902   fscanf(infile,"%d", &NPRINT);
  903   fscanf(infile,"%d", &NFMC);
  904 
  905 
  906   CNSTNT();
  907   cout << "\nTEMPERATURE                = "<< TEMP << "\n";
  908   cout << "DENSITY                    = "<< RHO << " G/C.C.\n";
  909   cout << "NUMBER OF MOLECULES        = "<< NMOL << "\n";
  910   cout << "TIME STEP                  = "<< TSTEP << " SEC\n";
  911   cout << "ORDER USED TO SOLVE F=MA   = "<< NORDER << "\n";
  912   cout << "NO. OF TIME STEPS          = "<< NSTEP <<"\n";
  913   cout << "FREQUENCY OF DATA SAVING   = "<< NSAVE <<"\n";
  914   cout << "FREQUENCY TO WRITE RST FILE= "<< NRST <<"\n";
  915   cout.flush();
  916   SYSCNS();
  917   cout << "SPHERICAL CUTOFF RADIUS    = "<< CUTOFF <<" ANGSTROM\n";
  918   IRST=0;
  919   fclose(infile);
  920 }
  921 
  922 // --- ENSEMBLE_H -------------------------
  923 
  924 class ensemble {
  925 public:
  926   int numMol;
  927   h2o *molecule;
  928   skratch_pad *pad;
  929   double TTMV;
  930   double TVIR;
  931   acc_double VIR;
  932   vector POT;
  933   double TKIN;
  934   vector KIN;
  935 h2o* getMolecule(int idx)  	{ return &(molecule[idx]); }
  936 skratch_pad* getPad(int idx)	{ return &(pad[idx]); }
  937 void INTER_POTENG();
  938 void INTRA_POTENG();
  939 void interfLoop1();
  940 void interfLoop2(int idx);
  941 void potengLoop1();
  942 void potengLoop2(int idx);
  943 double interf2_aux(skratch_pad *p1, skratch_pad *p2,
  944 double Res1[3][3], double Res2[3][3]);
  945 void interPoteng2Aux(skratch_pad *p1, skratch_pad *p2, double *res);
  946 void interf2(skratch_pad *p1, skratch_pad *p2);
  947   ensemble(int n);
  948   ~ensemble();
  949 
  950 void BNDRY(double size);
  951 void INITIA(const char * );
  952 void clearTVIR();
  953 void updateTVIR();
  954 void computeVIR();
  955 void POTENG();
  956 void KINETI();
  957 void clearTKIN();
  958 void updateTKIN();
  959 void INTRAF();
  960 void loadData();
  961 void storeData(int dest);
  962 void INTERF(int dest);
  963 void SCALEFORCES(int dest);
  964 void printENERGY(int iter);
  965 void PREDIC();
  966 void CORREC();
  967 void CSHIFT2(skratch_pad *p1, skratch_pad *p2,double CL[3][15],double *S);
  968 void MDMAIN();
  969 void stepsystem();
  970 void DUMP(int iter);
  971 };
  972 
  973 // -----------------------------------------------------------------------
  974 // GLOBALS
  975 // -----------------------------------------------------------------------
  976 simparm *parms;
  977 ensemble *liquid;
  978 
  979 ensemble::ensemble(int num){
  980 
  981   numMol = num;
  982   molecule = new h2o[numMol];
  983   pad = new skratch_pad[numMol];
  984 
  985   TTMV = 0.0;
  986   TVIR = 0.0;
  987   TKIN = 0.0;
  988   VIR.writeval(0.0);
  989   POT.vecClr();
  990   KIN.vecClr();
  991 }
  992 
  993 ensemble::~ensemble(){
  994   delete[] molecule;
  995   delete[] pad;
  996 }
  997 
  998 void ensemble::clearTVIR(){
  999   TVIR = 0.0;
 1000 }
 1001 
 1002 void ensemble::updateTVIR(){
 1003   TVIR -= VIR.readval();
 1004 }
 1005 
 1006 void ensemble::BNDRY(double size) {
 1007   int i;
 1008 
 1009   for(i = 0; i < numMol; i++){
 1010     molecule[i].bndry(size);
 1011   }
 1012 }
 1013 
 1014 void ensemble::INITIA(const char * file){
 1015   vector minimum;
 1016   double t1[3];
 1017   FILE *random_numbers;
 1018   double XMAS[3],XS,ZERO,WCOS,WSIN,XT[3],YT[3],ZT[3],SU[3],FAC,SUM[3],t[3],NS;
 1019   int i, j, k, mol;
 1020 
 1021   random_numbers = fopen(file,"r");
 1022   if (random_numbers == NULL) {
 1023     cout << " *** Error: Can't open file random.in (ensemble::INITIA)\n";
 1024     cout.flush();
 1025     exit(0);
 1026   }
 1027 
 1028      NS = pow((double) numMol, 1.0/3.0) - 0.00001;
 1029      XS = parms->getBOXL()/NS;
 1030      ZERO = XS * 0.50;
 1031      WCOS = ROH * cos(ANGLE * 0.5);
 1032      WSIN = ROH * sin(ANGLE * 0.5);
 1033 
 1034      cout << "***** NEW RUN STARTING FROM REGULAR LATTICE *****\n";
 1035      cout.flush();
 1036      XT[1] = ZERO;
 1037      mol = 0;
 1038      for(i=0; i < NS; i++) {
 1039        XT[0]=XT[1]+WCOS;
 1040        XT[2]=XT[0];
 1041        YT[1]=ZERO;
 1042        for(j=0; j < NS; j++) {
 1043          YT[0]=YT[1]+WSIN;
 1044          YT[2]=YT[1]-WSIN;
 1045          ZT[0]=ZT[1]=ZT[2]=ZERO;
 1046          for(k = 0; k < NS; k++) {
 1047 	   molecule[mol].loadDirPos(XDIR,XT);
 1048 	   molecule[mol].loadDirPos(YDIR,YT);
 1049 	   molecule[mol].loadDirPos(ZDIR,ZT);
 1050            mol++;
 1051 	   if (mol == numMol) goto exitLoop;
 1052            ZT[0] += XS; ZT[1] += XS; ZT[2] += XS;
 1053          }
 1054          YT[1] += XS;
 1055        }
 1056        XT[1] += XS;
 1057      }
 1058 
 1059  exitLoop: if(numMol != mol) {
 1060        cout << " *** Error: Lattice number of mol mismatch \n";
 1061        exit(-1);
 1062      }
 1063 
 1064   /* ASSIGN RANDOM MOMENTA */
 1065   fscanf(random_numbers,"%lf",&SU[0]);
 1066 
 1067   tvecClr(SUM);
 1068   for(i = 0; i < numMol; i++) {
 1069     fscanf(random_numbers,"%lf",&t[0]);
 1070     fscanf(random_numbers,"%lf",&t[1]);
 1071     fscanf(random_numbers,"%lf",&t[2]);
 1072 //    fscanf(random_numbers,"%lf%lf%lf",&t[0],&t[1],&t[2]);
 1073 
 1074     molecule[i].loadH1Vel(t);
 1075     tvecAdd(SUM,t);
 1076     fscanf(random_numbers,"%lf",&t[0]);
 1077     fscanf(random_numbers,"%lf",&t[1]);
 1078     fscanf(random_numbers,"%lf",&t[2]);
 1079 
 1080 //    fscanf(random_numbers,"%lf%lf%lf",&t[0],&t[1],&t[2]);
 1081     molecule[i].loadOVel(t);
 1082     tvecAdd(SUM,t);
 1083     fscanf(random_numbers,"%lf",&t[0]);
 1084     fscanf(random_numbers,"%lf",&t[1]);
 1085     fscanf(random_numbers,"%lf",&t[2]);
 1086 
 1087 //    fscanf(random_numbers,"%lf%lf%lf",&t[0],&t[1],&t[2]);
 1088     molecule[i].loadH2Vel(t);
 1089     tvecAdd(SUM,t);
 1090   }
 1091   /* .....FIND AVERAGE MOMENTA PER ATOM */
 1092   SUM[0] /= (NATOM*parms->getNMOL());
 1093   SUM[1] /= (NATOM*parms->getNMOL());
 1094   SUM[2] /= (NATOM*parms->getNMOL());
 1095 
 1096   /* FIND NORMALIZATION FACTOR SO THAT <K.E.>=KT/2 */
 1097 
 1098   tvecClr(SU);
 1099   for (i = 0; i < numMol; i++) {
 1100     molecule[i].storeDirVel(XDIR,t);
 1101     SU[0] += (pow((t[0]-SUM[0]),2.0) + pow((t[2]-SUM[0]),2.0))/HMAS
 1102 	+ pow((t[1]-SUM[0]),2.0)/OMAS;
 1103     molecule[i].storeDirVel(YDIR,t);
 1104     SU[1] += (pow((t[0]-SUM[1]),2.0) + pow((t[2]-SUM[1]),2.0))/HMAS
 1105 	+ pow((t[1]-SUM[1]),2.0)/OMAS;
 1106     molecule[i].storeDirVel(ZDIR,t);
 1107     SU[2] += (pow((t[0]-SUM[2]),2.0) + pow((t[2]-SUM[2]),2.0))/HMAS
 1108 	+ pow((t[1]-SUM[2]),2.0)/OMAS;
 1109   }
 1110   FAC=parms->computeFAC();
 1111   SU[0]=sqrt(FAC/SU[0]);
 1112   SU[1]=sqrt(FAC/SU[1]);
 1113   SU[2]=sqrt(FAC/SU[2]);
 1114   /* NORMALIZE INDIVIDUAL VELOCITIES SO THAT THERE IS NO BULK MOMENTA */
 1115 
 1116   XMAS[0] = 1.0/HMAS;
 1117   XMAS[1] = 1.0/OMAS;
 1118   XMAS[2] = 1.0/HMAS;
 1119 
 1120   for(i = 0; i < numMol; i++) {
 1121     molecule[i].storeH1Vel(t1);
 1122     tvecSub(t1,SUM);
 1123     tvecProd(t1,SU);
 1124     tvecScale(t1,XMAS[0]);
 1125     molecule[i].loadH1Vel(t1);
 1126 
 1127     molecule[i].storeOVel(t1);
 1128     tvecSub(t1,SUM);
 1129     tvecProd(t1,SU);
 1130     tvecScale(t1,XMAS[1]);
 1131     molecule[i].loadOVel(t1);
 1132 
 1133     molecule[i].storeH2Vel(t1);
 1134     tvecSub(t1,SUM);
 1135     tvecProd(t1,SU);
 1136     tvecScale(t1,XMAS[2]);
 1137     molecule[i].loadH2Vel(t1);
 1138   }
 1139 }
 1140 
 1141 // ----------------------------------------------------------------------
 1142 // POTENG
 1143 // ----------------------------------------------------------------------
 1144 void ensemble::potengLoop1(){
 1145   int i;
 1146   int max =  numMol-1;
 1147   for(i=0; i < max; i++) {
 1148     potengLoop2(i);
 1149   }
 1150 }
 1151 
 1152 void ensemble::potengLoop2(int idx){
 1153   int i;
 1154 
 1155 
 1156   int max =  numMol;
 1157   for(i=idx+1; i < max; i++) {
 1158      double tmp[3];
 1159     interPoteng2Aux(&pad[idx],&pad[i],tmp);
 1160     POT.vecAdd(tmp);
 1161   }
 1162 }
 1163 
 1164 void ensemble::INTER_POTENG() {
 1165   loadData();
 1166   potengLoop1();
 1167 }
 1168 
 1169 
 1170 // ---------------------------------------------------------------------------
 1171 void ensemble::interPoteng2Aux(skratch_pad *p1, skratch_pad *p2, double *res){
 1172   int KC,K;
 1173   double CL[3][15], RS[14], RL[14], PotR, PotF, S[2];
 1174   double Cut2, Cutoff, Ref1, Ref2;
 1175 
 1176   PotF = 0.0;
 1177   PotR = 0.0;
 1178   S[0] = parms->getBOXH();	// BoxH
 1179   S[1] = parms->getBOXL();	// BoxL
 1180   Cut2 = parms->getCUT2();
 1181   Cutoff = parms->getCUTOFF();
 1182   Ref1 = parms->getREF1();
 1183   Ref2 = parms->getREF2();
 1184 
 1185   CSHIFT2(p1,p2,CL,S);
 1186   KC=0;
 1187   for(K = 0; K < 9; K++) {
 1188     RS[K]=CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K];
 1189     if(RS[K] > Cut2)
 1190       KC=KC+1;
 1191   }
 1192 
 1193   if (KC != 9) {
 1194      for(K = 0; K < 9; K++) {
 1195       if(RS[K] <= Cut2){
 1196         RL[K] = sqrt(RS[K]);
 1197       } else {
 1198     	RL[K] = Cutoff;
 1199 	RS[K] = Cut2;
 1200       }
 1201     }
 1202 
 1203     PotR += (-QQ2/RL[1]-QQ2/RL[2]-QQ2/RL[3]-QQ2/RL[4]
 1204       +QQ /RL[5]+QQ /RL[6]+QQ /RL[7]+QQ /RL[8]+QQ4/RL[0]);
 1205     PotF += (-Ref2*RS[0]-Ref1*((RS[5]+RS[6]+RS[7]+RS[8])*0.5
 1206       -RS[1]-RS[2]-RS[3]-RS[4]));
 1207 
 1208     if (KC <= 0) {
 1209       for (K = 9; K <  14; K++)
 1210 	RL[K]=sqrt(CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K]);
 1211       PotR += A1 * exp(-B1*RL[9]) +A2*(exp(-B2*RL[5])+
 1212           exp(-B2*RL[6])+exp(-B2*RL[7])+exp(-B2*RL[ 8]))
 1213      	+A3*(exp(-B3*RL[10])+exp(-B3*RL[11])+exp(-B3*RL[12])
 1214 	+exp(-B3*RL[13]))-A4*(exp(-B4*RL[10])+exp(-B4*RL[11])
 1215 	+exp(-B4*RL[12])+exp(-B4*RL[13]));
 1216     }
 1217   }
 1218   res[0] = 0.0;
 1219   res[1] = PotR;
 1220   res[2] = PotF;
 1221 }
 1222 
 1223 
 1224 // ----------------------------------------------------------------------
 1225 void ensemble::INTRA_POTENG(){
 1226   int i;
 1227   for(i=0; i < numMol; i++){
 1228     molecule[i].intra_poteng(&POT);
 1229   }
 1230 }
 1231 
 1232 // ----------------------------------------------------------------------
 1233 void ensemble::POTENG(){
 1234   POT.vecClr();
 1235   INTRA_POTENG();
 1236   INTER_POTENG();
 1237   POT.vecScale(parms->getFPOT());
 1238 }
 1239 
 1240 // ----------------------------------------------------------------------
 1241 // KINETI
 1242 // ----------------------------------------------------------------------
 1243 void ensemble::clearTKIN(){
 1244   TKIN = 0.0;
 1245 }
 1246 
 1247 void ensemble::updateTKIN(){
 1248   TKIN += KIN.vecNorm1();
 1249 }
 1250 
 1251 void ensemble::KINETI(){
 1252   int i;
 1253   KIN.vecClr();
 1254   for(i=0; i < numMol; i++){
 1255     molecule[i].kineti(&KIN);
 1256   }
 1257   updateTKIN();
 1258 }
 1259 
 1260 // ----------------------------------------------------------------------
 1261 // VIR
 1262 // ----------------------------------------------------------------------
 1263 void ensemble::computeVIR(){
 1264   int i;
 1265 
 1266   for(i=0; i < numMol; i++){
 1267     molecule[i].vir(&VIR);
 1268   }
 1269 }
 1270 
 1271 // ----------------------------------------------------------------------
 1272 void ensemble::printENERGY(int iter){
 1273   double XVIR, TENN, XTT, AVGT;
 1274   double loc_pot[3];
 1275 
 1276   POT.vecStore(loc_pot);
 1277   XVIR = TVIR*parms->getFPOT()*0.50/((double)TTMV);
 1278   AVGT = TKIN*parms->getFKIN()*(parms->getTEMP())*2.00/(3.00*((double)TTMV));
 1279   TENN = KIN.vecNorm1() * parms->getFKIN();
 1280   XTT = tvecNorm1(loc_pot)+TENN;
 1281 
 1282   if ((iter % parms->getNPRINT()) == 0){
 1283     cout << "     "<<iter;
 1284     cout << " "<<TENN;
 1285     cout << " "<<loc_pot[0];
 1286     cout << " "<<loc_pot[1];
 1287     cout << " "<<loc_pot[2];
 1288     cout << " "<<XTT;
 1289     cout << " "<<AVGT;
 1290     cout << " "<<XVIR<< "\n";;
 1291   }
 1292 }
 1293 
 1294 // ----------------------------------------------------------------------------
 1295 void ensemble::PREDIC(){
 1296   int i, ord;
 1297   double *coeffs;
 1298 
 1299   for(i=0; i < numMol; i++){
 1300     ord = parms->getNORDER();
 1301     coeffs = parms->getTLC();
 1302     molecule[i].predic(ord,coeffs);
 1303   }
 1304 }
 1305 
 1306 // ----------------------------------------------------------------------------
 1307 void ensemble::CORREC(){
 1308   int i, ord;
 1309   double *coeffs;
 1310 
 1311   for(i=0; i < numMol; i++){
 1312     ord = parms->getNORDER();
 1313     coeffs = parms->getPCC();
 1314     molecule[i].correc(ord,coeffs);
 1315   }
 1316 }
 1317 
 1318 // ----------------------------------------------------------------------------
 1319 // MDMAIN --- MOLECULAR DYNAMICS LOOP
 1320 // ----------------------------------------------------------------------------
 1321 void ensemble::MDMAIN(){
 1322 
 1323   cout << "RESTART "<< parms->getIRST();
 1324   cout << " AFTER "<< parms->getELPST() << " SECONDS\n";
 1325   clearTVIR();
 1326   clearTKIN();
 1327 
 1328   if(parms->getNSAVE() > 0){
 1329     cout << "COLLECTING X AND V DATA AT EVERY %4d TIME STEPS\n";
 1330     cout << parms->getNSAVE();
 1331   }
 1332   cout << "INTERMEDIATE RESULTS (ENERGY EXPRESSED IN UNITS OF KT ATOM) \n";
 1333   cout << "  TIME       KINETIC   INTRA POT   INTER POT   REACT POT       ";
 1334   cout << "TOTAL  \n<TEMPERATURE>   <VIR>\n";
 1335 
 1336   stepsystem();
 1337 
 1338 }
 1339 
 1340 // --- MAIN LOOP -------------------------------------------------------
 1341 void ensemble::stepsystem(){
 1342   int i, n, start, stop;
 1343   int start_serial, stop_serial, total_serial;
 1344   int ticks;
 1345   double dticks;
 1346 
 1347   total_serial = 0;
 1348   get_time(&start);
 1349 
 1350   n = parms->getNSTEP();
 1351   for(i=1;i <= n; i++) {
 1352 
 1353     TTMV += 1.00;
 1354 #if defined(DEBUG_OPT)
 1355   cout << " >>> Step: " << i << "\n";
 1356 #endif
 1357 
 1358     get_time(&start_serial);
 1359     PREDIC();
 1360 #if defined(DEBUG_OPT)
 1361   cout << "    >>> Predic: done.\n";
 1362 #endif
 1363     INTRAF();
 1364     VIR.writeval(0.0);
 1365     computeVIR();
 1366 #if defined(DEBUG_OPT)
 1367   cout << "    >>> Intraf: done.\n";
 1368 #endif
 1369     get_time(&stop_serial);
 1370     total_serial += (stop_serial - start_serial);
 1371     INTERF(FORCES);
 1372     SCALEFORCES(FORCES);
 1373 
 1374 #if defined(DEBUG_OPT)
 1375   cout << "    >>> Interf: done.\n";
 1376 #endif
 1377     get_time(&start_serial);
 1378     CORREC();
 1379 #if defined(DEBUG_OPT)
 1380   cout << "    >>> Correc: done.\n";
 1381 #endif
 1382     BNDRY(parms->getBOXL());
 1383 #if defined(DEBUG_OPT)
 1384   cout << "    >>> Bondary: done.\n";
 1385 #endif
 1386     KINETI();
 1387 #if defined(DEBUG_OPT)
 1388   cout << "    >>> Kinetic: done.\n";
 1389 #endif
 1390     updateTVIR();
 1391 
 1392     get_time(&stop_serial);
 1393     total_serial += (stop_serial - start_serial);
 1394 
 1395     if(( (i % parms->getNPRINT()) == 0) ||
 1396 	((parms->getNSAVE() > 0) && ((i % parms->getNSAVE()) == 0))) {
 1397       POTENG();
 1398       printENERGY(i);
 1399     }
 1400 
 1401   }
 1402   get_time(&stop);
 1403   get_ticks_time(&ticks);
 1404   dticks = ticks;
 1405 
 1406   cout << "execution time = " << (stop-start) / dticks << "\n";
 1407 }
 1408 
 1409 void ensemble::CSHIFT2(skratch_pad *p1, skratch_pad *p2, double L[3][15], double *S){
 1410   int i;
 1411 
 1412   // ---  XDIR ---
 1413   double vm1,vm2,h1pos1,h1pos2,h2pos1,h2pos2,opos1,opos2;
 1414   vm1 = p1->VM.vecField(0);
 1415   vm2 = p2->VM.vecField(0);
 1416   L[0][0] = vm1-vm2;
 1417   h1pos2  = p2->H1pos.vecField(0);
 1418   L[0][1] = vm1-h1pos2;
 1419   h2pos2  = p2->H2pos.vecField(0);
 1420   L[0][2] = vm1-h2pos2;
 1421   h1pos1 = p1->H1pos.vecField(0);
 1422   L[0][3] = h1pos1-vm2;
 1423   h2pos1 = p1->H2pos.vecField(0);
 1424   L[0][4] = h2pos1-vm2;
 1425   L[0][5] = h1pos1-h1pos2;
 1426   L[0][6] = h1pos1-h2pos2;
 1427   L[0][7] = h2pos1-h1pos2;
 1428   L[0][8] = h2pos1-h2pos2;
 1429   opos1 = p1->Opos.vecField(0);
 1430   opos2 = p2->Opos.vecField(0);
 1431   L[0][9] = opos1-opos2;
 1432   L[0][10] = opos1-h1pos2;
 1433   L[0][11] = opos1-h2pos2;
 1434   L[0][12] = h1pos1-opos2;
 1435   L[0][13] = h2pos1-opos2;
 1436   for (i = 0; i < 14; i++) {
 1437     if (fabs(L[0][i]) > S[0])
 1438       L[0][i] -= (sign(S[1],L[0][i]));
 1439   }
 1440 
 1441   // ---  YDIR ---
 1442   vm1 = p1->VM.vecField(1);
 1443   vm2 = p2->VM.vecField(1);
 1444   h1pos1 = p1->H1pos.vecField(1);
 1445   L[1][0] = vm1-vm2;
 1446   h1pos2 = p2->H1pos.vecField(1);
 1447   L[1][1] = vm1-h1pos2;
 1448   h2pos2 = p2->H2pos.vecField(1);
 1449   L[1][2] = vm1-h2pos2;
 1450   L[1][3] = h1pos1-vm2;
 1451   h2pos1 = p1->H2pos.vecField(1);
 1452   L[1][4] = h2pos1 - vm2;
 1453   L[1][5] = h1pos1-h1pos2;
 1454   L[1][6] = h1pos1-h2pos2;
 1455   L[1][7] = h2pos1-h1pos2;
 1456   L[1][8] = h2pos1-h2pos2;
 1457   opos1 = p1->Opos.vecField(1);
 1458   opos2 = p2->Opos.vecField(1);
 1459   L[1][9] = opos1-opos2;
 1460   L[1][10] = opos1-h1pos2;
 1461   L[1][11] = opos1-h2pos2;
 1462   L[1][12] = h1pos1-opos2;
 1463   L[1][13] = h2pos1-opos2;
 1464   for (i = 0; i < 14; i++) {
 1465     if (fabs(L[1][i]) > S[0])
 1466       L[1][i] -= (sign(S[1],L[1][i]));
 1467   }
 1468 
 1469   // ---  ZDIR ---
 1470   vm1 = p1->VM.vecField(2);
 1471   vm2 = p2->VM.vecField(2);
 1472   L[2][0] = vm1-vm2;
 1473   h1pos2 = p2->H1pos.vecField(2);
 1474   L[2][1] = vm1-h1pos2;
 1475   h2pos2 = p2->H2pos.vecField(2);
 1476   L[2][2] = vm1-h2pos2;
 1477   h1pos1 = p1->H1pos.vecField(2);
 1478   L[2][3] = h1pos1-vm2;
 1479   h2pos1 = p1->H2pos.vecField(2);
 1480   L[2][4] = h2pos1-vm2;
 1481   L[2][5] = h1pos1-h1pos2;
 1482   L[2][6] = h1pos1-h2pos2;
 1483   L[2][7] = h2pos1-h1pos2;
 1484   L[2][8] = h2pos1-h2pos2;
 1485   opos1 = p1->Opos.vecField(2);
 1486   opos2 = p2->Opos.vecField(2);
 1487   L[2][9] = opos1-opos2;
 1488   L[2][10] = opos1-h1pos2;
 1489   L[2][11] = opos1-h2pos2;
 1490   L[2][12] = h1pos1-opos2;
 1491   L[2][13] = h2pos1-opos2;
 1492   for (i = 0; i < 14; i++) {
 1493     if (fabs(L[2][i]) > S[0])
 1494       L[2][i] -= (sign(S[1],L[2][i]));
 1495   }
 1496 }
 1497 
 1498 // -----------------------------------------------------------------------
 1499 // INTRAF:
 1500 // -----------------------------------------------------------------------
 1501 void ensemble::INTRAF(){
 1502   int i;
 1503   for (i = 0; i < numMol; i++){
 1504     molecule[i].intraf();
 1505   }
 1506 }
 1507 
 1508 // ------------------------------------------------------------------------
 1509 // interf2 - operates on the skratch pad areas
 1510 // ------------------------------------------------------------------------
 1511 double ensemble::interf2_aux(skratch_pad *p1, skratch_pad *p2, double Res1[3][3], double Res2[3][3]){
 1512 
 1513   int KC, K;
 1514   double CL[3][15], RS[15], FF[15], RL[15], GG[15];
 1515   double G110, G23, G45, TT1, TT, FTEMP;
 1516   double gCUT2, gREF1, gREF2, gREF4;
 1517   double loc_vir;
 1518   double S[2];
 1519 
 1520   S[0] = parms->getBOXH();
 1521   S[1] = parms->getBOXL();
 1522   gCUT2 = parms->getCUT2();
 1523   gREF1 = parms->getREF1();
 1524   gREF2 = parms->getREF2();
 1525   gREF4 = parms->getREF4();
 1526   loc_vir = 0.0;
 1527 
 1528   CSHIFT2(p1,p2,CL,S);
 1529   KC=0;
 1530   for (K = 0; K < 9; K++) {
 1531     RS[K]=CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K];
 1532     if (RS[K] > gCUT2)
 1533       KC++;
 1534   }
 1535 
 1536   if(KC != 9) {
 1537     for (K = 0; K < 14; K++)
 1538         FF[K]=0.0;
 1539     if(RS[0] < gCUT2) {
 1540       FF[0]=QQ4/(RS[0]*sqrt(RS[0]))+gREF4;
 1541       loc_vir = loc_vir + FF[0]*RS[0];
 1542     }
 1543     for (K = 1; K < 5; K++) {
 1544       if(RS[K] < gCUT2) {
 1545 	FF[K]= -QQ2/(RS[K]*sqrt(RS[K]))-gREF2;
 1546 	loc_vir = loc_vir + FF[K]*RS[K];
 1547       }
 1548       if(RS[K+4] <= gCUT2) {
 1549 	RL[K+4]=sqrt(RS[K+4]);
 1550  	FF[K+4]=QQ/(RS[K+4]*RL[K+4])+gREF1;
 1551 	loc_vir = loc_vir + FF[K+4]*RS[K+4];
 1552       }
 1553     }
 1554 
 1555     if(KC == 0) {
 1556       RS[9]=CL[0][9]*CL[0][9]+CL[1][9]*CL[1][9]+CL[2][9]*CL[2][9];
 1557       RL[9]=sqrt(RS[9]);
 1558       FF[9]=AB1*exp(-B1*RL[9])/RL[9];
 1559       loc_vir = loc_vir + FF[9]*RS[9];
 1560       for (K = 10; K < 14; K++) {
 1561 	FTEMP=AB2*exp(-B2*RL[K-5])/RL[K-5];
 1562 	FF[K-5]=FF[K-5]+FTEMP;
 1563 	loc_vir= loc_vir+FTEMP*RS[K-5];
 1564 	RS[K]=CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K];
 1565 	RL[K]=sqrt(RS[K]);
 1566 	FF[K]=(AB3*exp(-B3*RL[K])-AB4*exp(-B4*RL[K]))/RL[K];
 1567 	loc_vir = loc_vir + FF[K]*RS[K];
 1568       }
 1569     }
 1570 
 1571   /* CALCULATE X-COMPONENT FORCES */
 1572     for (K = 0; K < 14; K++)
 1573       GG[K+1]=FF[K]*CL[0][K];
 1574 
 1575     G110=GG[10]+GG[1]*C1;
 1576     G23=GG[2]+GG[3];
 1577     G45=GG[4]+GG[5];
 1578     Res1[1][0] = G110+GG[11]+GG[12]+C1*G23;
 1579     Res2[1][0] = -G110-GG[13]-GG[14]-C1*G45;
 1580     TT1=GG[1]*C2;
 1581     TT=G23*C2+TT1;
 1582     Res1[0][0] = GG[6]+GG[7]+GG[13]+TT+GG[4];
 1583 //    Res2[0][0] = GG[8]+GG[9]+GG[14]+TT+GG[5];
 1584     Res1[2][0] = GG[8]+GG[9]+GG[14]+TT+GG[5];
 1585     TT=G45*C2+TT1;
 1586 //    Res1[2][0] = -GG[6]-GG[8]-GG[11]-TT-GG[2];
 1587     Res2[0][0] = -GG[6]-GG[8]-GG[11]-TT-GG[2];
 1588     Res2[2][0] = -GG[7]-GG[9]-GG[12]-TT-GG[3];
 1589 
 1590   /* CALCULATE Y-COMPONENT FORCES */
 1591     for(K = 0; K < 14; K++)
 1592       GG[K+1]=FF[K]*CL[1][K];
 1593 
 1594     G110=GG[10]+GG[1]*C1;
 1595     G23=GG[2]+GG[3];
 1596     G45=GG[4]+GG[5];
 1597     Res1[1][1] = G110+GG[11]+GG[12]+C1*G23;
 1598     Res2[1][1] = -G110-GG[13]-GG[14]-C1*G45;
 1599     TT1=GG[1]*C2;
 1600     TT=G23*C2+TT1;
 1601     Res1[0][1] = GG[6]+GG[7]+GG[13]+TT+GG[4];
 1602 //    Res2[0][1] = GG[8]+GG[9]+GG[14]+TT+GG[5];
 1603     Res1[2][1] = GG[8]+GG[9]+GG[14]+TT+GG[5];
 1604     TT=G45*C2+TT1;
 1605 //    Res1[2][1] = -GG[6]-GG[8]-GG[11]-TT-GG[2];
 1606     Res2[0][1] = -GG[6]-GG[8]-GG[11]-TT-GG[2];
 1607     Res2[2][1] = -GG[7]-GG[9]-GG[12]-TT-GG[3];
 1608 
 1609   /* CALCULATE Z-COMPONENT FORCES */
 1610     for(K = 0; K < 14; K++)
 1611       GG[K+1]=FF[K]*CL[2][K];
 1612 
 1613     G110=GG[10]+GG[1]*C1;
 1614     G23=GG[2]+GG[3];
 1615     G45=GG[4]+GG[5];
 1616     Res1[1][2] = G110+GG[11]+GG[12]+C1*G23;
 1617     Res2[1][2] = -G110-GG[13]-GG[14]-C1*G45;
 1618     TT1=GG[1]*C2;
 1619     TT=G23*C2+TT1;
 1620     Res1[0][2] = GG[6]+GG[7]+GG[13]+TT+GG[4];
 1621 //    Res2[0][2] = GG[8]+GG[9]+GG[14]+TT+GG[5];
 1622     Res1[2][2] = GG[8]+GG[9]+GG[14]+TT+GG[5];
 1623     TT=G45*C2+TT1;
 1624 //    Res1[2][2] = -GG[6]-GG[8]-GG[11]-TT-GG[2];
 1625     Res2[0][2] = -GG[6]-GG[8]-GG[11]-TT-GG[2];
 1626     Res2[2][2] = -GG[7]-GG[9]-GG[12]-TT-GG[3];
 1627 
 1628     /* Update DESTination of calculations */
 1629 //    p1->update_forces(Res1);
 1630 //    p2->update_forces(Res2);
 1631   } else {
 1632     tvecClr(Res1[0]); tvecClr(Res1[1]); tvecClr(Res1[2]);
 1633     tvecClr(Res2[0]); tvecClr(Res2[1]); tvecClr(Res2[2]);
 1634   }
 1635   return loc_vir;
 1636 }
 1637 
 1638 
 1639 void ensemble::interf2(skratch_pad *p1, skratch_pad *p2){
 1640   double incr, Res1[3][3], Res2[3][3];
 1641 
 1642   incr = interf2_aux(p1,p2,Res1,Res2);
 1643   p1->update_forces(Res1);
 1644   p2->update_forces(Res2);
 1645   VIR.addvalSync(incr);
 1646 }
 1647 
 1648 void ensemble::loadData(){
 1649   int i;
 1650   h2o *mol;
 1651   skratch_pad *p;
 1652 
 1653   for(i=0; i < numMol; i++){
 1654     mol = getMolecule(i);
 1655     p = getPad(i);
 1656     p->read_data(mol);
 1657   }
 1658 }
 1659 
 1660 void ensemble::storeData(int dest){
 1661   int i;
 1662   h2o *p1;
 1663   skratch_pad *p2;
 1664 
 1665   for(i=0; i < numMol; i++){
 1666     p1 = getMolecule(i);
 1667     p2 = getPad(i);
 1668     p1->updateFields(dest,p2);
 1669   }
 1670 }
 1671 
 1672 void ensemble::interfLoop2(int idx){
 1673   int i;
 1674   skratch_pad *p1, *p2;
 1675   int max =  numMol;
 1676 
 1677   for(i = idx+1; i < max; i++){
 1678     p1 = getPad(idx);
 1679     p2 = getPad(i);
 1680     interf2(p1,p2);
 1681   }
 1682 }
 1683 
 1684 void ensemble::interfLoop1(){
 1685   int i;
 1686   int max =  numMol-1;
 1687    #pragma omp parallel for schedule(static)
 1688   for(i = 0; i < max; i++) {
 1689     interfLoop2(i);
 1690   }
 1691 }
 1692 
 1693 
 1694 void ensemble::INTERF(int DEST){
 1695   loadData();
 1696   interfLoop1();
 1697   storeData(DEST);
 1698 }
 1699 
 1700 // --------------------------------------------------------------------
 1701 void ensemble::SCALEFORCES(int Dest){
 1702   int i;
 1703   double HM, OM;
 1704 
 1705   for(i=0; i < numMol; i++){
 1706     HM = parms->getFHM();
 1707     OM = parms->getFOM();
 1708     molecule[i].scaleMomenta(Dest,HM,OM);
 1709   }
 1710 }
 1711 
 1712 // --------------------------------------------------------------------
 1713 void ensemble::DUMP(int iter){
 1714   int i;
 1715 
 1716   printf(" >> STEP: %d\n",iter);
 1717   printf(" NMOL: %d\n",numMol);
 1718   printf(" TVIR: %6.5f\n",TVIR);
 1719   printf("  VIR: %6.5f\n",VIR.readval());
 1720   printf(" FKIN: %6.5f\n",parms->getFKIN());
 1721   printf(" TKIN: %6.5f\n",TKIN);
 1722   printf("  KIN: "); KIN.vecPrint(); printf("\n");
 1723   printf(" FPOT: %6.5f\n",parms->getFPOT());
 1724   printf("  POT: "); POT.vecPrint(); printf("\n");
 1725   printf("\n");
 1726   for (i = 0; i < numMol; i++){
 1727     molecule[i].dump();
 1728     printf("\n");
 1729   }
 1730 }
 1731 
 1732 
 1733 // ----------------------------------------------------------------------------
 1734 // MAIN
 1735 // ----------------------------------------------------------------------------
 1736 int main(int argc, char **argv) {
 1737 
 1738   int n, start_time, stop_time;
 1739   int ticks;
 1740   double dticks;
 1741 
 1742   const char * filename, * randomfilename;
 1743 
 1744   if (argc != 3) {
 1745     printf("Usage: %s input_param_file random_file\n", argv[0]);
 1746     return(1);
 1747   }
 1748   filename = argv[1];
 1749   randomfilename = argv[2];
 1750 
 1751   cout << " >>> Program Started\n";
 1752   cout.flush();
 1753   get_ticks_time(&ticks);
 1754   dticks = ticks;
 1755 
 1756   get_time(&start_time);
 1757 
 1758   parms = new simparm;
 1759   parms->loadParms(filename);
 1760   n = parms->getNMOL();
 1761   liquid = new ensemble(n);
 1762 
 1763   liquid->INITIA(randomfilename);
 1764   liquid->INTRAF();
 1765   liquid->computeVIR();
 1766   liquid->INTERF(ACC);
 1767   liquid->SCALEFORCES(ACC);
 1768 
 1769   parms->setNFMC(-1);
 1770   parms->resetStat();
 1771   liquid->MDMAIN();
 1772   get_time(&stop_time);
 1773 
 1774   get_ticks_time(&ticks);
 1775   dticks = ticks;
 1776   cout << "TOTAL CPU USED = " << ((stop_time-start_time)/dticks) << " SECONDS\n";
 1777   cout << "\nELAPSED CPU TIME IN SECONDS: " << ((stop_time-start_time)/dticks) << "\n";
 1778   cout << "\nMFLOP RATE: "<< 3432.550/((stop_time-start_time)/dticks) << "\n";
 1779   cout << "\nTotal Time = " << ((stop_time-start_time)/dticks) << "\n";
 1780   cout.flush();
 1781 }