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.vecAddSync(Res[0]); 315 Oforce.vecAddSync(Res[1]); 316 H2force.vecAddSync(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 #pragma omp parallel for schedule(static) 1148 for(i=0; i < max; i++) { 1149 potengLoop2(i); 1150 } 1151 } 1152 1153 void ensemble::potengLoop2(int idx){ 1154 int i; 1155 1156 1157 int max = numMol; 1158 for(i=idx+1; i < max; i++) { 1159 double tmp[3]; 1160 interPoteng2Aux(&pad[idx],&pad[i],tmp); 1161 POT.vecAdd(tmp); 1162 } 1163 } 1164 1165 void ensemble::INTER_POTENG() { 1166 loadData(); 1167 potengLoop1(); 1168 } 1169 1170 1171 // --------------------------------------------------------------------------- 1172 void ensemble::interPoteng2Aux(skratch_pad *p1, skratch_pad *p2, double *res){ 1173 int KC,K; 1174 double CL[3][15], RS[14], RL[14], PotR, PotF, S[2]; 1175 double Cut2, Cutoff, Ref1, Ref2; 1176 1177 PotF = 0.0; 1178 PotR = 0.0; 1179 S[0] = parms->getBOXH(); // BoxH 1180 S[1] = parms->getBOXL(); // BoxL 1181 Cut2 = parms->getCUT2(); 1182 Cutoff = parms->getCUTOFF(); 1183 Ref1 = parms->getREF1(); 1184 Ref2 = parms->getREF2(); 1185 1186 CSHIFT2(p1,p2,CL,S); 1187 KC=0; 1188 for(K = 0; K < 9; K++) { 1189 RS[K]=CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K]; 1190 if(RS[K] > Cut2) 1191 KC=KC+1; 1192 } 1193 1194 if (KC != 9) { 1195 for(K = 0; K < 9; K++) { 1196 if(RS[K] <= Cut2){ 1197 RL[K] = sqrt(RS[K]); 1198 } else { 1199 RL[K] = Cutoff; 1200 RS[K] = Cut2; 1201 } 1202 } 1203 1204 PotR += (-QQ2/RL[1]-QQ2/RL[2]-QQ2/RL[3]-QQ2/RL[4] 1205 +QQ /RL[5]+QQ /RL[6]+QQ /RL[7]+QQ /RL[8]+QQ4/RL[0]); 1206 PotF += (-Ref2*RS[0]-Ref1*((RS[5]+RS[6]+RS[7]+RS[8])*0.5 1207 -RS[1]-RS[2]-RS[3]-RS[4])); 1208 1209 if (KC <= 0) { 1210 for (K = 9; K < 14; K++) 1211 RL[K]=sqrt(CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K]); 1212 PotR += A1 * exp(-B1*RL[9]) +A2*(exp(-B2*RL[5])+ 1213 exp(-B2*RL[6])+exp(-B2*RL[7])+exp(-B2*RL[ 8])) 1214 +A3*(exp(-B3*RL[10])+exp(-B3*RL[11])+exp(-B3*RL[12]) 1215 +exp(-B3*RL[13]))-A4*(exp(-B4*RL[10])+exp(-B4*RL[11]) 1216 +exp(-B4*RL[12])+exp(-B4*RL[13])); 1217 } 1218 } 1219 res[0] = 0.0; 1220 res[1] = PotR; 1221 res[2] = PotF; 1222 } 1223 1224 1225 // ---------------------------------------------------------------------- 1226 void ensemble::INTRA_POTENG(){ 1227 int i; 1228 for(i=0; i < numMol; i++){ 1229 molecule[i].intra_poteng(&POT); 1230 } 1231 } 1232 1233 // ---------------------------------------------------------------------- 1234 void ensemble::POTENG(){ 1235 POT.vecClr(); 1236 INTRA_POTENG(); 1237 INTER_POTENG(); 1238 POT.vecScale(parms->getFPOT()); 1239 } 1240 1241 // ---------------------------------------------------------------------- 1242 // KINETI 1243 // ---------------------------------------------------------------------- 1244 void ensemble::clearTKIN(){ 1245 TKIN = 0.0; 1246 } 1247 1248 void ensemble::updateTKIN(){ 1249 TKIN += KIN.vecNorm1(); 1250 } 1251 1252 void ensemble::KINETI(){ 1253 int i; 1254 KIN.vecClr(); 1255 for(i=0; i < numMol; i++){ 1256 molecule[i].kineti(&KIN); 1257 } 1258 updateTKIN(); 1259 } 1260 1261 // ---------------------------------------------------------------------- 1262 // VIR 1263 // ---------------------------------------------------------------------- 1264 void ensemble::computeVIR(){ 1265 int i; 1266 1267 for(i=0; i < numMol; i++){ 1268 molecule[i].vir(&VIR); 1269 } 1270 } 1271 1272 // ---------------------------------------------------------------------- 1273 void ensemble::printENERGY(int iter){ 1274 double XVIR, TENN, XTT, AVGT; 1275 double loc_pot[3]; 1276 1277 POT.vecStore(loc_pot); 1278 XVIR = TVIR*parms->getFPOT()*0.50/((double)TTMV); 1279 AVGT = TKIN*parms->getFKIN()*(parms->getTEMP())*2.00/(3.00*((double)TTMV)); 1280 TENN = KIN.vecNorm1() * parms->getFKIN(); 1281 XTT = tvecNorm1(loc_pot)+TENN; 1282 1283 if ((iter % parms->getNPRINT()) == 0){ 1284 cout << " "<<iter; 1285 cout << " "<<TENN; 1286 cout << " "<<loc_pot[0]; 1287 cout << " "<<loc_pot[1]; 1288 cout << " "<<loc_pot[2]; 1289 cout << " "<<XTT; 1290 cout << " "<<AVGT; 1291 cout << " "<<XVIR<< "\n";; 1292 } 1293 } 1294 1295 // ---------------------------------------------------------------------------- 1296 void ensemble::PREDIC(){ 1297 int i, ord; 1298 double *coeffs; 1299 1300 for(i=0; i < numMol; i++){ 1301 ord = parms->getNORDER(); 1302 coeffs = parms->getTLC(); 1303 molecule[i].predic(ord,coeffs); 1304 } 1305 } 1306 1307 // ---------------------------------------------------------------------------- 1308 void ensemble::CORREC(){ 1309 int i, ord; 1310 double *coeffs; 1311 1312 for(i=0; i < numMol; i++){ 1313 ord = parms->getNORDER(); 1314 coeffs = parms->getPCC(); 1315 molecule[i].correc(ord,coeffs); 1316 } 1317 } 1318 1319 // ---------------------------------------------------------------------------- 1320 // MDMAIN --- MOLECULAR DYNAMICS LOOP 1321 // ---------------------------------------------------------------------------- 1322 void ensemble::MDMAIN(){ 1323 1324 cout << "RESTART "<< parms->getIRST(); 1325 cout << " AFTER "<< parms->getELPST() << " SECONDS\n"; 1326 clearTVIR(); 1327 clearTKIN(); 1328 1329 if(parms->getNSAVE() > 0){ 1330 cout << "COLLECTING X AND V DATA AT EVERY %4d TIME STEPS\n"; 1331 cout << parms->getNSAVE(); 1332 } 1333 cout << "INTERMEDIATE RESULTS (ENERGY EXPRESSED IN UNITS OF KT ATOM) \n"; 1334 cout << " TIME KINETIC INTRA POT INTER POT REACT POT "; 1335 cout << "TOTAL \n<TEMPERATURE> <VIR>\n"; 1336 1337 stepsystem(); 1338 1339 } 1340 1341 // --- MAIN LOOP ------------------------------------------------------- 1342 void ensemble::stepsystem(){ 1343 int i, n, start, stop; 1344 int start_serial, stop_serial, total_serial; 1345 int ticks; 1346 double dticks; 1347 1348 total_serial = 0; 1349 get_time(&start); 1350 1351 n = parms->getNSTEP(); 1352 for(i=1;i <= n; i++) { 1353 1354 TTMV += 1.00; 1355 #if defined(DEBUG_OPT) 1356 cout << " >>> Step: " << i << "\n"; 1357 #endif 1358 1359 get_time(&start_serial); 1360 PREDIC(); 1361 #if defined(DEBUG_OPT) 1362 cout << " >>> Predic: done.\n"; 1363 #endif 1364 INTRAF(); 1365 VIR.writeval(0.0); 1366 computeVIR(); 1367 #if defined(DEBUG_OPT) 1368 cout << " >>> Intraf: done.\n"; 1369 #endif 1370 get_time(&stop_serial); 1371 total_serial += (stop_serial - start_serial); 1372 INTERF(FORCES); 1373 SCALEFORCES(FORCES); 1374 1375 #if defined(DEBUG_OPT) 1376 cout << " >>> Interf: done.\n"; 1377 #endif 1378 get_time(&start_serial); 1379 CORREC(); 1380 #if defined(DEBUG_OPT) 1381 cout << " >>> Correc: done.\n"; 1382 #endif 1383 BNDRY(parms->getBOXL()); 1384 #if defined(DEBUG_OPT) 1385 cout << " >>> Bondary: done.\n"; 1386 #endif 1387 KINETI(); 1388 #if defined(DEBUG_OPT) 1389 cout << " >>> Kinetic: done.\n"; 1390 #endif 1391 updateTVIR(); 1392 1393 get_time(&stop_serial); 1394 total_serial += (stop_serial - start_serial); 1395 1396 if(( (i % parms->getNPRINT()) == 0) || 1397 ((parms->getNSAVE() > 0) && ((i % parms->getNSAVE()) == 0))) { 1398 POTENG(); 1399 printENERGY(i); 1400 } 1401 1402 } 1403 get_time(&stop); 1404 get_ticks_time(&ticks); 1405 dticks = ticks; 1406 1407 cout << "execution time = " << (stop-start) / dticks << "\n"; 1408 } 1409 1410 void ensemble::CSHIFT2(skratch_pad *p1, skratch_pad *p2, double L[3][15], double *S){ 1411 int i; 1412 1413 // --- XDIR --- 1414 double vm1,vm2,h1pos1,h1pos2,h2pos1,h2pos2,opos1,opos2; 1415 vm1 = p1->VM.vecField(0); 1416 vm2 = p2->VM.vecField(0); 1417 L[0][0] = vm1-vm2; 1418 h1pos2 = p2->H1pos.vecField(0); 1419 L[0][1] = vm1-h1pos2; 1420 h2pos2 = p2->H2pos.vecField(0); 1421 L[0][2] = vm1-h2pos2; 1422 h1pos1 = p1->H1pos.vecField(0); 1423 L[0][3] = h1pos1-vm2; 1424 h2pos1 = p1->H2pos.vecField(0); 1425 L[0][4] = h2pos1-vm2; 1426 L[0][5] = h1pos1-h1pos2; 1427 L[0][6] = h1pos1-h2pos2; 1428 L[0][7] = h2pos1-h1pos2; 1429 L[0][8] = h2pos1-h2pos2; 1430 opos1 = p1->Opos.vecField(0); 1431 opos2 = p2->Opos.vecField(0); 1432 L[0][9] = opos1-opos2; 1433 L[0][10] = opos1-h1pos2; 1434 L[0][11] = opos1-h2pos2; 1435 L[0][12] = h1pos1-opos2; 1436 L[0][13] = h2pos1-opos2; 1437 for (i = 0; i < 14; i++) { 1438 if (fabs(L[0][i]) > S[0]) 1439 L[0][i] -= (sign(S[1],L[0][i])); 1440 } 1441 1442 // --- YDIR --- 1443 vm1 = p1->VM.vecField(1); 1444 vm2 = p2->VM.vecField(1); 1445 h1pos1 = p1->H1pos.vecField(1); 1446 L[1][0] = vm1-vm2; 1447 h1pos2 = p2->H1pos.vecField(1); 1448 L[1][1] = vm1-h1pos2; 1449 h2pos2 = p2->H2pos.vecField(1); 1450 L[1][2] = vm1-h2pos2; 1451 L[1][3] = h1pos1-vm2; 1452 h2pos1 = p1->H2pos.vecField(1); 1453 L[1][4] = h2pos1 - vm2; 1454 L[1][5] = h1pos1-h1pos2; 1455 L[1][6] = h1pos1-h2pos2; 1456 L[1][7] = h2pos1-h1pos2; 1457 L[1][8] = h2pos1-h2pos2; 1458 opos1 = p1->Opos.vecField(1); 1459 opos2 = p2->Opos.vecField(1); 1460 L[1][9] = opos1-opos2; 1461 L[1][10] = opos1-h1pos2; 1462 L[1][11] = opos1-h2pos2; 1463 L[1][12] = h1pos1-opos2; 1464 L[1][13] = h2pos1-opos2; 1465 for (i = 0; i < 14; i++) { 1466 if (fabs(L[1][i]) > S[0]) 1467 L[1][i] -= (sign(S[1],L[1][i])); 1468 } 1469 1470 // --- ZDIR --- 1471 vm1 = p1->VM.vecField(2); 1472 vm2 = p2->VM.vecField(2); 1473 L[2][0] = vm1-vm2; 1474 h1pos2 = p2->H1pos.vecField(2); 1475 L[2][1] = vm1-h1pos2; 1476 h2pos2 = p2->H2pos.vecField(2); 1477 L[2][2] = vm1-h2pos2; 1478 h1pos1 = p1->H1pos.vecField(2); 1479 L[2][3] = h1pos1-vm2; 1480 h2pos1 = p1->H2pos.vecField(2); 1481 L[2][4] = h2pos1-vm2; 1482 L[2][5] = h1pos1-h1pos2; 1483 L[2][6] = h1pos1-h2pos2; 1484 L[2][7] = h2pos1-h1pos2; 1485 L[2][8] = h2pos1-h2pos2; 1486 opos1 = p1->Opos.vecField(2); 1487 opos2 = p2->Opos.vecField(2); 1488 L[2][9] = opos1-opos2; 1489 L[2][10] = opos1-h1pos2; 1490 L[2][11] = opos1-h2pos2; 1491 L[2][12] = h1pos1-opos2; 1492 L[2][13] = h2pos1-opos2; 1493 for (i = 0; i < 14; i++) { 1494 if (fabs(L[2][i]) > S[0]) 1495 L[2][i] -= (sign(S[1],L[2][i])); 1496 } 1497 } 1498 1499 // ----------------------------------------------------------------------- 1500 // INTRAF: 1501 // ----------------------------------------------------------------------- 1502 void ensemble::INTRAF(){ 1503 int i; 1504 for (i = 0; i < numMol; i++){ 1505 molecule[i].intraf(); 1506 } 1507 } 1508 1509 // ------------------------------------------------------------------------ 1510 // interf2 - operates on the skratch pad areas 1511 // ------------------------------------------------------------------------ 1512 double ensemble::interf2_aux(skratch_pad *p1, skratch_pad *p2, double Res1[3][3], double Res2[3][3]){ 1513 1514 int KC, K; 1515 double CL[3][15], RS[15], FF[15], RL[15], GG[15]; 1516 double G110, G23, G45, TT1, TT, FTEMP; 1517 double gCUT2, gREF1, gREF2, gREF4; 1518 double loc_vir; 1519 double S[2]; 1520 1521 S[0] = parms->getBOXH(); 1522 S[1] = parms->getBOXL(); 1523 gCUT2 = parms->getCUT2(); 1524 gREF1 = parms->getREF1(); 1525 gREF2 = parms->getREF2(); 1526 gREF4 = parms->getREF4(); 1527 loc_vir = 0.0; 1528 1529 CSHIFT2(p1,p2,CL,S); 1530 KC=0; 1531 for (K = 0; K < 9; K++) { 1532 RS[K]=CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K]; 1533 if (RS[K] > gCUT2) 1534 KC++; 1535 } 1536 1537 if(KC != 9) { 1538 for (K = 0; K < 14; K++) 1539 FF[K]=0.0; 1540 if(RS[0] < gCUT2) { 1541 FF[0]=QQ4/(RS[0]*sqrt(RS[0]))+gREF4; 1542 loc_vir = loc_vir + FF[0]*RS[0]; 1543 } 1544 for (K = 1; K < 5; K++) { 1545 if(RS[K] < gCUT2) { 1546 FF[K]= -QQ2/(RS[K]*sqrt(RS[K]))-gREF2; 1547 loc_vir = loc_vir + FF[K]*RS[K]; 1548 } 1549 if(RS[K+4] <= gCUT2) { 1550 RL[K+4]=sqrt(RS[K+4]); 1551 FF[K+4]=QQ/(RS[K+4]*RL[K+4])+gREF1; 1552 loc_vir = loc_vir + FF[K+4]*RS[K+4]; 1553 } 1554 } 1555 1556 if(KC == 0) { 1557 RS[9]=CL[0][9]*CL[0][9]+CL[1][9]*CL[1][9]+CL[2][9]*CL[2][9]; 1558 RL[9]=sqrt(RS[9]); 1559 FF[9]=AB1*exp(-B1*RL[9])/RL[9]; 1560 loc_vir = loc_vir + FF[9]*RS[9]; 1561 for (K = 10; K < 14; K++) { 1562 FTEMP=AB2*exp(-B2*RL[K-5])/RL[K-5]; 1563 FF[K-5]=FF[K-5]+FTEMP; 1564 loc_vir= loc_vir+FTEMP*RS[K-5]; 1565 RS[K]=CL[0][K]*CL[0][K]+CL[1][K]*CL[1][K]+CL[2][K]*CL[2][K]; 1566 RL[K]=sqrt(RS[K]); 1567 FF[K]=(AB3*exp(-B3*RL[K])-AB4*exp(-B4*RL[K]))/RL[K]; 1568 loc_vir = loc_vir + FF[K]*RS[K]; 1569 } 1570 } 1571 1572 /* CALCULATE X-COMPONENT FORCES */ 1573 for (K = 0; K < 14; K++) 1574 GG[K+1]=FF[K]*CL[0][K]; 1575 1576 G110=GG[10]+GG[1]*C1; 1577 G23=GG[2]+GG[3]; 1578 G45=GG[4]+GG[5]; 1579 Res1[1][0] = G110+GG[11]+GG[12]+C1*G23; 1580 Res2[1][0] = -G110-GG[13]-GG[14]-C1*G45; 1581 TT1=GG[1]*C2; 1582 TT=G23*C2+TT1; 1583 Res1[0][0] = GG[6]+GG[7]+GG[13]+TT+GG[4]; 1584 // Res2[0][0] = GG[8]+GG[9]+GG[14]+TT+GG[5]; 1585 Res1[2][0] = GG[8]+GG[9]+GG[14]+TT+GG[5]; 1586 TT=G45*C2+TT1; 1587 // Res1[2][0] = -GG[6]-GG[8]-GG[11]-TT-GG[2]; 1588 Res2[0][0] = -GG[6]-GG[8]-GG[11]-TT-GG[2]; 1589 Res2[2][0] = -GG[7]-GG[9]-GG[12]-TT-GG[3]; 1590 1591 /* CALCULATE Y-COMPONENT FORCES */ 1592 for(K = 0; K < 14; K++) 1593 GG[K+1]=FF[K]*CL[1][K]; 1594 1595 G110=GG[10]+GG[1]*C1; 1596 G23=GG[2]+GG[3]; 1597 G45=GG[4]+GG[5]; 1598 Res1[1][1] = G110+GG[11]+GG[12]+C1*G23; 1599 Res2[1][1] = -G110-GG[13]-GG[14]-C1*G45; 1600 TT1=GG[1]*C2; 1601 TT=G23*C2+TT1; 1602 Res1[0][1] = GG[6]+GG[7]+GG[13]+TT+GG[4]; 1603 // Res2[0][1] = GG[8]+GG[9]+GG[14]+TT+GG[5]; 1604 Res1[2][1] = GG[8]+GG[9]+GG[14]+TT+GG[5]; 1605 TT=G45*C2+TT1; 1606 // Res1[2][1] = -GG[6]-GG[8]-GG[11]-TT-GG[2]; 1607 Res2[0][1] = -GG[6]-GG[8]-GG[11]-TT-GG[2]; 1608 Res2[2][1] = -GG[7]-GG[9]-GG[12]-TT-GG[3]; 1609 1610 /* CALCULATE Z-COMPONENT FORCES */ 1611 for(K = 0; K < 14; K++) 1612 GG[K+1]=FF[K]*CL[2][K]; 1613 1614 G110=GG[10]+GG[1]*C1; 1615 G23=GG[2]+GG[3]; 1616 G45=GG[4]+GG[5]; 1617 Res1[1][2] = G110+GG[11]+GG[12]+C1*G23; 1618 Res2[1][2] = -G110-GG[13]-GG[14]-C1*G45; 1619 TT1=GG[1]*C2; 1620 TT=G23*C2+TT1; 1621 Res1[0][2] = GG[6]+GG[7]+GG[13]+TT+GG[4]; 1622 // Res2[0][2] = GG[8]+GG[9]+GG[14]+TT+GG[5]; 1623 Res1[2][2] = GG[8]+GG[9]+GG[14]+TT+GG[5]; 1624 TT=G45*C2+TT1; 1625 // Res1[2][2] = -GG[6]-GG[8]-GG[11]-TT-GG[2]; 1626 Res2[0][2] = -GG[6]-GG[8]-GG[11]-TT-GG[2]; 1627 Res2[2][2] = -GG[7]-GG[9]-GG[12]-TT-GG[3]; 1628 1629 /* Update DESTination of calculations */ 1630 // p1->update_forces(Res1); 1631 // p2->update_forces(Res2); 1632 } else { 1633 tvecClr(Res1[0]); tvecClr(Res1[1]); tvecClr(Res1[2]); 1634 tvecClr(Res2[0]); tvecClr(Res2[1]); tvecClr(Res2[2]); 1635 } 1636 return loc_vir; 1637 } 1638 1639 1640 void ensemble::interf2(skratch_pad *p1, skratch_pad *p2){ 1641 double incr, Res1[3][3], Res2[3][3]; 1642 1643 incr = interf2_aux(p1,p2,Res1,Res2); 1644 p1->update_forces(Res1); 1645 p2->update_forces(Res2); 1646 VIR.addvalRepl(incr); 1647 } 1648 1649 void ensemble::loadData(){ 1650 int i; 1651 h2o *mol; 1652 skratch_pad *p; 1653 1654 for(i=0; i < numMol; i++){ 1655 mol = getMolecule(i); 1656 p = getPad(i); 1657 p->read_data(mol); 1658 } 1659 } 1660 1661 void ensemble::storeData(int dest){ 1662 int i; 1663 h2o *p1; 1664 skratch_pad *p2; 1665 1666 for(i=0; i < numMol; i++){ 1667 p1 = getMolecule(i); 1668 p2 = getPad(i); 1669 p1->updateFields(dest,p2); 1670 } 1671 } 1672 1673 void ensemble::interfLoop2(int idx){ 1674 int i; 1675 skratch_pad *p1, *p2; 1676 int max = numMol; 1677 1678 #pragma omp parallel for private (p2, p1) schedule(static) 1679 for(i = idx+1; i < max; i++){ 1680 p1 = getPad(idx); 1681 p2 = getPad(i); 1682 interf2(p1,p2); 1683 } 1684 } 1685 1686 void ensemble::interfLoop1(){ 1687 int i; 1688 int max = numMol-1; 1689 for(i = 0; i < max; i++) { 1690 interfLoop2(i); 1691 } 1692 } 1693 1694 1695 void ensemble::INTERF(int DEST){ 1696 loadData(); 1697 interfLoop1(); 1698 storeData(DEST); 1699 } 1700 1701 // -------------------------------------------------------------------- 1702 void ensemble::SCALEFORCES(int Dest){ 1703 int i; 1704 double HM, OM; 1705 1706 for(i=0; i < numMol; i++){ 1707 HM = parms->getFHM(); 1708 OM = parms->getFOM(); 1709 molecule[i].scaleMomenta(Dest,HM,OM); 1710 } 1711 } 1712 1713 // -------------------------------------------------------------------- 1714 void ensemble::DUMP(int iter){ 1715 int i; 1716 1717 printf(" >> STEP: %d\n",iter); 1718 printf(" NMOL: %d\n",numMol); 1719 printf(" TVIR: %6.5f\n",TVIR); 1720 printf(" VIR: %6.5f\n",VIR.readval()); 1721 printf(" FKIN: %6.5f\n",parms->getFKIN()); 1722 printf(" TKIN: %6.5f\n",TKIN); 1723 printf(" KIN: "); KIN.vecPrint(); printf("\n"); 1724 printf(" FPOT: %6.5f\n",parms->getFPOT()); 1725 printf(" POT: "); POT.vecPrint(); printf("\n"); 1726 printf("\n"); 1727 for (i = 0; i < numMol; i++){ 1728 molecule[i].dump(); 1729 printf("\n"); 1730 } 1731 } 1732 1733 1734 // ---------------------------------------------------------------------------- 1735 // MAIN 1736 // ---------------------------------------------------------------------------- 1737 int main(int argc, char **argv) { 1738 1739 int n, start_time, stop_time; 1740 int ticks; 1741 double dticks; 1742 1743 const char * filename, * randomfilename; 1744 1745 if (argc != 3) { 1746 printf("Usage: %s input_param_file random_file\n", argv[0]); 1747 return(1); 1748 } 1749 filename = argv[1]; 1750 randomfilename = argv[2]; 1751 1752 cout << " >>> Program Started\n"; 1753 cout.flush(); 1754 get_ticks_time(&ticks); 1755 dticks = ticks; 1756 1757 get_time(&start_time); 1758 1759 parms = new simparm; 1760 parms->loadParms(filename); 1761 n = parms->getNMOL(); 1762 liquid = new ensemble(n); 1763 1764 liquid->INITIA(randomfilename); 1765 liquid->INTRAF(); 1766 liquid->computeVIR(); 1767 liquid->INTERF(ACC); 1768 liquid->SCALEFORCES(ACC); 1769 1770 parms->setNFMC(-1); 1771 parms->resetStat(); 1772 liquid->MDMAIN(); 1773 get_time(&stop_time); 1774 1775 get_ticks_time(&ticks); 1776 dticks = ticks; 1777 cout << "TOTAL CPU USED = " << ((stop_time-start_time)/dticks) << " SECONDS\n"; 1778 cout << "\nELAPSED CPU TIME IN SECONDS: " << ((stop_time-start_time)/dticks) << "\n"; 1779 cout << "\nMFLOP RATE: "<< 3432.550/((stop_time-start_time)/dticks) << "\n"; 1780 cout << "\nTotal Time = " << ((stop_time-start_time)/dticks) << "\n"; 1781 cout.flush(); 1782 }