#include <math.h>#include "charm++.h"#include "SimParameters.h"#include "HomePatch.h"#include "AtomMap.h"#include "Node.h"#include "PatchMap.inl"#include "main.h"#include "ProxyMgr.decl.h"#include "ProxyMgr.h"#include "Migration.h"#include "Molecule.h"#include "PatchMgr.h"#include "Sequencer.h"#include "LdbCoordinator.h"#include "Settle.h"#include "ReductionMgr.h"#include "Sync.h"#include "Random.h"#include "Priorities.h"#include "Debug.h"Go to the source code of this file.
Defines | |
| #define | TINY 1.0e-20; |
| #define | MAXHGS 10 |
| #define | MIN_DEBUG_LEVEL 2 |
Typedefs | |
| typedef int | HGArrayInt [MAXHGS] |
| typedef BigReal | HGArrayBigReal [MAXHGS] |
| typedef zVector | HGArrayVector [MAXHGS] |
| typedef BigReal | HGMatrixBigReal [MAXHGS][MAXHGS] |
| typedef zVector | HGMatrixVector [MAXHGS][MAXHGS] |
Functions | |
| int | average (CompAtom *qtilde, const HGArrayVector &q, BigReal *lambda, const int n, const int m, const HGArrayBigReal &imass, const HGArrayBigReal &length2, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab, const BigReal tolf, const int ntrial) |
| void | mollify (CompAtom *qtilde, const HGArrayVector &q0, const BigReal *lambda, HGArrayVector &force, const int n, const int m, const HGArrayBigReal &imass, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab) |
| int | compDistance (const void *a, const void *b) |
| void | lubksb (HGMatrixBigReal &a, int n, HGArrayInt &indx, HGArrayBigReal &b) |
| void | ludcmp (HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d) |
| void | G_q (const HGArrayVector &refab, HGMatrixVector &gqij, const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl) |
|
|
Definition at line 40 of file HomePatch.C. |
|
|
Definition at line 41 of file HomePatch.C. |
|
|
Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved. Definition at line 39 of file HomePatch.C. |
|
|
Definition at line 46 of file HomePatch.C. Referenced by average(), ludcmp(), mollify(), HomePatch::mollyAverage(), and HomePatch::mollyMollify(). |
|
|
Definition at line 45 of file HomePatch.C. Referenced by average(), mollify(), HomePatch::mollyAverage(), and HomePatch::mollyMollify(). |
|
|
Definition at line 47 of file HomePatch.C. Referenced by average(), mollify(), HomePatch::mollyAverage(), and HomePatch::mollyMollify(). |
|
|
Definition at line 48 of file HomePatch.C. |
|
|
Definition at line 49 of file HomePatch.C. |
|
||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2498 of file HomePatch.C. References BigReal, G_q(), HGArrayBigReal, HGArrayInt, HGArrayVector, HGMatrixBigReal, HGMatrixVector, iINFO(), iout, lubksb(), ludcmp(), and CompAtom::position. Referenced by PressureProfileReduction::getData(), HomePatch::mollyAverage(), and PressureProfileReduction::PressureProfileReduction(). 02498 {
02499 // input: n = length of hydrogen group to be averaged (shaked)
02500 // q[n] = vector array of original positions
02501 // m = number of constraints
02502 // imass[n] = inverse mass for each atom
02503 // length2[m] = square of reference bond length for each constraint
02504 // ial[m] = atom a in each constraint
02505 // ibl[m] = atom b in each constraint
02506 // refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
02507 // tolf = function error tolerance for Newton's iteration
02508 // ntrial = max number of Newton's iterations
02509 // output: lambda[m] = double array of lagrange multipliers (used by mollify)
02510 // qtilde[n] = vector array of averaged (shaked) positions
02511
02512 int k,k1,i,j;
02513 BigReal errx,errf,d,tolx;
02514
02515 HGArrayInt indx;
02516 HGArrayBigReal p;
02517 HGArrayBigReal fvec;
02518 HGMatrixBigReal fjac;
02519 HGArrayVector avgab;
02520 HGMatrixVector grhs;
02521 HGMatrixVector auxrhs;
02522 HGMatrixVector glhs;
02523
02524 // iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
02525 tolx=tolf;
02526
02527 // initialize lambda, globalGrhs
02528
02529 for (i=0;i<m;i++) {
02530 lambda[i]=0.0;
02531 }
02532
02533 // define grhs, auxrhs for all iterations
02534 // grhs= g_x(q)
02535 //
02536 G_q(refab,grhs,n,m,ial,ibl);
02537 for (k=1;k<=ntrial;k++) {
02538 // usrfun(qtilde,q0,lambda,fvec,fjac,n,water);
02539 HGArrayBigReal gij;
02540 // this used to be main loop of usrfun
02541 // compute qtilde given q0, lambda, IMASSes
02542 {
02543 BigReal multiplier;
02544 HGArrayVector tmp;
02545 for (i=0;i<m;i++) {
02546 multiplier = lambda[i];
02547 // auxrhs = M^{-1}grhs^{T}
02548 for (j=0;j<n;j++) {
02549 auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
02550 }
02551 }
02552 for (j=0;j<n;j++) {
02553 // tmp[j]=0.0;
02554 for (i=0;i<m;i++) {
02555 tmp[j]+=auxrhs[i][j];
02556 }
02557 }
02558
02559 for (j=0;j<n;j++) {
02560 qtilde[j].position=q[j]+tmp[j];
02561 }
02562 // delete [] tmp;
02563 }
02564
02565 for ( i = 0; i < m; i++ ) {
02566 avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
02567 }
02568
02569 // iout<<iINFO << "Calling Jac" << std::endl<<endi;
02570 // Jac(qtilde, q0, fjac,n,water);
02571 {
02572 // Vector glhs[3*n+3];
02573
02574 HGMatrixVector grhs2;
02575
02576 G_q(avgab,glhs,n,m,ial,ibl);
02577 #ifdef DEBUG0
02578 iout<<iINFO << "G_q:" << std::endl<<endi;
02579 for (i=0;i<m;i++) {
02580 iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
02581 }
02582 #endif
02583 // G_q(refab,grhs2,m,ial,ibl);
02584 // update with the masses
02585 for (j=0; j<n; j++) { // number of atoms
02586 for (i=0; i<m;i++) { // number of constraints
02587 grhs2[i][j] = grhs[i][j]*imass[j];
02588 }
02589 }
02590
02591 // G_q(qtilde) * M^-1 G_q'(q0) =
02592 // G_q(qtilde) * grhs'
02593 for (i=0;i<m;i++) { // number of constraints
02594 for (j=0;j<m;j++) { // number of constraints
02595 fjac[i][j] = 0;
02596 for (k1=0;k1<n;k1++) {
02597 fjac[i][j] += glhs[i][k1]*grhs2[j][k1];
02598 }
02599 }
02600 }
02601 #ifdef DEBUG0
02602 iout<<iINFO << "glhs" <<endi;
02603 for(i=0;i<9;i++) {
02604 iout<<iINFO << glhs[i] << ","<<endi;
02605 }
02606 iout<<iINFO << std::endl<<endi;
02607 for(i=0;i<9;i++) {
02608 iout<<iINFO << grhs2[i] << ","<<endi;
02609 }
02610 iout<<iINFO << std::endl<<endi;
02611 #endif
02612 // delete[] grhs2;
02613 }
02614 // end of Jac calculation
02615 #ifdef DEBUG0
02616 iout<<iINFO << "Jac" << std::endl<<endi;
02617 for (i=0;i<m;i++)
02618 for (j=0;j<m;j++)
02619 iout<<iINFO << fjac[i][j] << " "<<endi;
02620 iout<< std::endl<<endi;
02621 #endif
02622 // calculate constraints in gij for n constraints this being a water
02623 // G(qtilde, gij, n, water);
02624 for (i=0;i<m;i++) {
02625 gij[i]=avgab[i]*avgab[i]-length2[i];
02626 }
02627 #ifdef DEBUG0
02628 iout<<iINFO << "G" << std::endl<<endi;
02629 iout<<iINFO << "( "<<endi;
02630 for(i=0;i<m-1;i++) {
02631 iout<<iINFO << gij[i] << ", "<<endi;
02632 }
02633 iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
02634 #endif
02635 // fill the return vector
02636 for(i=0;i<m;i++) {
02637 fvec[i] = gij[i];
02638 }
02639 // free up the constraints
02640 // delete[] gij;
02641 // continue Newton's iteration
02642 errf=0.0;
02643 for (i=0;i<m;i++) errf += fabs(fvec[i]);
02644 #ifdef DEBUG0
02645 iout<<iINFO << "errf: " << errf << std::endl<<endi;
02646 #endif
02647 if (errf <= tolf) {
02648 break;
02649 }
02650 for (i=0;i<m;i++) p[i] = -fvec[i];
02651 // iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
02652 ludcmp(fjac,m,indx,&d);
02653 lubksb(fjac,m,indx,p);
02654
02655 errx=0.0;
02656 for (i=0;i<m;i++) {
02657 errx += fabs(p[i]);
02658 }
02659 for (i=0;i<m;i++)
02660 lambda[i] += p[i];
02661
02662 #ifdef DEBUG0
02663 iout<<iINFO << "lambda:" << lambda[0]
02664 << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
02665 iout<<iINFO << "errx: " << errx << std::endl<<endi;
02666 #endif
02667 if (errx <= tolx) break;
02668 #ifdef DEBUG0
02669 iout<<iINFO << "Qtilde:" << std::endl<<endi;
02670 iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi;
02671 #endif
02672 }
02673 #ifdef DEBUG
02674 iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
02675 #endif
02676
02677 return k; //
02678 }
|
|
||||||||||||
|
Definition at line 432 of file HomePatch.C. Referenced by HomePatch::buildSpanningTree(). 00433 {
00434 int d1 = abs(*(int *)a - CkMyPe());
00435 int d2 = abs(*(int *)b - CkMyPe());
00436 if (d1 < d2)
00437 return -1;
00438 else if (d1 == d2)
00439 return 0;
00440 else
00441 return 1;
00442 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 2486 of file HomePatch.C. Referenced by average(), and mollify(). 02487 {
02488 int i;
02489 // step through the rows of the matrix
02490 for(i=0;i<m;i++) {
02491 gqij[i][ial[i]]=2.0*refab[i];
02492 gqij[i][ibl[i]]=-gqij[i][ial[i]];
02493 }
02494 };
|
|
||||||||||||||||||||
|
Definition at line 2413 of file HomePatch.C. Referenced by average(), and mollify(). 02415 {
02416 int i,ii=-1,ip,j;
02417 double sum;
02418
02419 for (i=0;i<n;i++) {
02420 ip=indx[i];
02421 sum=b[ip];
02422 b[ip]=b[i];
02423 if (ii >= 0)
02424 for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
02425 else if (sum) ii=i;
02426 b[i]=sum;
02427 }
02428 for (i=n-1;i>=0;i--) {
02429 sum=b[i];
02430 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
02431 b[i]=sum/a[i][i];
02432 }
02433 }
|
|
||||||||||||||||||||
|
Definition at line 2436 of file HomePatch.C. References HGArrayBigReal, and NAMD_die(). Referenced by average(), and mollify(). 02437 {
02438
02439 int i,imax,j,k;
02440 double big,dum,sum,temp;
02441 HGArrayBigReal vv;
02442 *d=1.0;
02443 for (i=0;i<n;i++) {
02444 big=0.0;
02445 for (j=0;j<n;j++)
02446 if ((temp=fabs(a[i][j])) > big) big=temp;
02447 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
02448 vv[i]=1.0/big;
02449 }
02450 for (j=0;j<n;j++) {
02451 for (i=0;i<j;i++) {
02452 sum=a[i][j];
02453 for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
02454 a[i][j]=sum;
02455 }
02456 big=0.0;
02457 for (i=j;i<n;i++) {
02458 sum=a[i][j];
02459 for (k=0;k<j;k++)
02460 sum -= a[i][k]*a[k][j];
02461 a[i][j]=sum;
02462 if ( (dum=vv[i]*fabs(sum)) >= big) {
02463 big=dum;
02464 imax=i;
02465 }
02466 }
02467 if (j != imax) {
02468 for (k=0;k<n;k++) {
02469 dum=a[imax][k];
02470 a[imax][k]=a[j][k];
02471 a[j][k]=dum;
02472 }
02473 *d = -(*d);
02474 vv[imax]=vv[j];
02475 }
02476 indx[j]=imax;
02477 if (a[j][j] == 0.0) a[j][j]=TINY;
02478 if (j != n-1) {
02479 dum=1.0/(a[j][j]);
02480 for (i=j+1;i<n;i++) a[i][j] *= dum;
02481 }
02482 }
02483 }
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2680 of file HomePatch.C. References BigReal, G_q(), HGArrayBigReal, HGArrayInt, HGArrayVector, HGMatrixBigReal, HGMatrixVector, lubksb(), ludcmp(), and CompAtom::position. Referenced by HomePatch::mollyMollify(). 02680 {
02681 int i,j,k;
02682 BigReal d;
02683 HGMatrixBigReal fjac;
02684 Vector zero(0.0,0.0,0.0);
02685
02686 HGArrayVector tmpforce;
02687 HGArrayVector tmpforce2;
02688 HGArrayVector y;
02689 HGMatrixVector grhs;
02690 HGMatrixVector glhs;
02691 HGArrayBigReal aux;
02692 HGArrayInt indx;
02693
02694 for(i=0;i<n;i++) {
02695 tmpforce[i]=imass[i]*force[i];
02696 }
02697
02698 HGMatrixVector grhs2;
02699 HGArrayVector avgab;
02700
02701 for ( i = 0; i < m; i++ ) {
02702 avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
02703 }
02704
02705 G_q(avgab,glhs,n,m,ial,ibl);
02706 G_q(refab,grhs,n,m,ial,ibl);
02707 // update with the masses
02708 for (j=0; j<n; j++) { // number of atoms
02709 for (i=0; i<m;i++) { // number of constraints
02710 grhs2[i][j] = grhs[i][j]*imass[j];
02711 }
02712 }
02713
02714 // G_q(qtilde) * M^-1 G_q'(q0) =
02715 // G_q(qtilde) * grhs'
02716 for (i=0;i<m;i++) { // number of constraints
02717 for (j=0;j<m;j++) { // number of constraints
02718 fjac[j][i] = 0;
02719 for (k=0;k<n;k++) {
02720 fjac[j][i] += glhs[i][k]*grhs2[j][k];
02721 }
02722 }
02723 }
02724
02725 // aux=gqij*tmpforce
02726 // globalGrhs::computeGlobalGrhs(q0,n,water);
02727 // G_q(refab,grhs,m,ial,ibl);
02728 for(i=0;i<m;i++) {
02729 aux[i]=0.0;
02730 for(j=0;j<n;j++) {
02731 aux[i]+=grhs[i][j]*tmpforce[j];
02732 }
02733 }
02734
02735 ludcmp(fjac,m,indx,&d);
02736 lubksb(fjac,m,indx,aux);
02737
02738 for(j=0;j<n;j++) {
02739 y[j] = zero;
02740 for(i=0;i<m;i++) {
02741 y[j] += aux[i]*glhs[i][j];
02742 }
02743 }
02744 for(i=0;i<n;i++) {
02745 y[i]=force[i]-y[i];
02746 }
02747
02748 // gqq12*y
02749 for(i=0;i<n;i++) {
02750 tmpforce2[i]=imass[i]*y[i];
02751 }
02752
02753 // here we assume that tmpforce is initialized to zero.
02754 for (i=0;i<n;i++) {
02755 tmpforce[i]=zero;
02756 }
02757
02758 for (j=0;j<m;j++) {
02759 Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
02760 tmpforce[ial[j]] += tmpf;
02761 tmpforce[ibl[j]] -= tmpf;
02762 }
02763 // c-ji the other bug for 2 constraint water was this line (2-4-99)
02764 // for(i=0;i<m;i++) {
02765 for(i=0;i<n;i++) {
02766 force[i]=tmpforce[i]+y[i];
02767 }
02768
02769 }
|
1.3.9.1