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 }