HomePatch.C File Reference

#include "time.h"
#include <math.h>
#include "charm++.h"
#include "qd.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 "Broadcasts.h"
#include "LdbCoordinator.h"
#include "ReductionMgr.h"
#include "Sync.h"
#include "Random.h"
#include "Priorities.h"
#include "ComputeNonbondedUtil.h"
#include "ComputeGBIS.inl"
#include "SortAtoms.h"
#include "ComputeQM.h"
#include "ComputeQMMgr.decl.h"
#include "Debug.h"
#include <vector>
#include <algorithm>
#include "ComputeNonbondedMICKernel.h"

Go to the source code of this file.

Defines

#define TINY   1.0e-20;
#define MAXHGS   10
#define MIN_DEBUG_LEVEL   2
#define MASS_EPSILON   (1.0e-35)
#define FIX_FOR_WATER

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)
static 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)


Define Documentation

#define FIX_FOR_WATER

Definition at line 1411 of file HomePatch.C.

#define MASS_EPSILON   (1.0e-35)

Definition at line 71 of file HomePatch.C.

#define MAXHGS   10

Definition at line 50 of file HomePatch.C.

#define MIN_DEBUG_LEVEL   2

Definition at line 51 of file HomePatch.C.

#define TINY   1.0e-20;

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 49 of file HomePatch.C.

Referenced by ludcmp().


Typedef Documentation

typedef BigReal HGArrayBigReal[MAXHGS]

Definition at line 60 of file HomePatch.C.

typedef int HGArrayInt[MAXHGS]

Definition at line 59 of file HomePatch.C.

typedef zVector HGArrayVector[MAXHGS]

Definition at line 61 of file HomePatch.C.

typedef BigReal HGMatrixBigReal[MAXHGS][MAXHGS]

Definition at line 62 of file HomePatch.C.

typedef zVector HGMatrixVector[MAXHGS][MAXHGS]

Definition at line 63 of file HomePatch.C.


Function Documentation

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 
)

Definition at line 4344 of file HomePatch.C.

References endi(), G_q(), iINFO(), iout, j, lubksb(), ludcmp(), and CompAtom::position.

Referenced by HomePatch::mollyAverage().

04344                                                                                                                                                                                                                                                                          {
04345   //  input:  n = length of hydrogen group to be averaged (shaked)
04346   //          q[n] = vector array of original positions
04347   //          m = number of constraints
04348   //          imass[n] = inverse mass for each atom
04349   //          length2[m] = square of reference bond length for each constraint
04350   //          ial[m] = atom a in each constraint 
04351   //          ibl[m] = atom b in each constraint 
04352   //          refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
04353   //          tolf = function error tolerance for Newton's iteration
04354   //          ntrial = max number of Newton's iterations
04355   //  output: lambda[m] = double array of lagrange multipliers (used by mollify)
04356   //          qtilde[n] = vector array of averaged (shaked) positions
04357 
04358   int k,k1,i,j;
04359   BigReal errx,errf,d,tolx;
04360 
04361   HGArrayInt indx;
04362   HGArrayBigReal p;
04363   HGArrayBigReal fvec;
04364   HGMatrixBigReal fjac;
04365   HGArrayVector avgab;
04366   HGMatrixVector grhs;
04367   HGMatrixVector auxrhs;
04368   HGMatrixVector glhs;
04369 
04370   //  iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
04371   tolx=tolf; 
04372   
04373   // initialize lambda, globalGrhs
04374 
04375   for (i=0;i<m;i++) {
04376     lambda[i]=0.0;
04377   }
04378 
04379   // define grhs, auxrhs for all iterations
04380   // grhs= g_x(q)
04381   //
04382   G_q(refab,grhs,n,m,ial,ibl);
04383   for (k=1;k<=ntrial;k++) {
04384     //    usrfun(qtilde,q0,lambda,fvec,fjac,n,water); 
04385     HGArrayBigReal gij;
04386     // this used to be main loop of usrfun
04387     // compute qtilde given q0, lambda, IMASSes
04388     {
04389       BigReal multiplier;
04390       HGArrayVector tmp;
04391       for (i=0;i<m;i++) {
04392         multiplier = lambda[i];
04393         // auxrhs = M^{-1}grhs^{T}
04394         for (j=0;j<n;j++) {
04395           auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
04396         }
04397       }
04398       for (j=0;j<n;j++) {
04399         //      tmp[j]=0.0;      
04400         for (i=0;i<m;i++) {
04401           tmp[j]+=auxrhs[i][j];
04402         }
04403       }
04404  
04405       for (j=0;j<n;j++) {
04406         qtilde[j].position=q[j]+tmp[j];
04407       }
04408       //      delete [] tmp;
04409     }
04410   
04411     for ( i = 0; i < m; i++ ) {
04412       avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04413     }
04414 
04415     //  iout<<iINFO << "Calling Jac" << std::endl<<endi;
04416     //  Jac(qtilde, q0, fjac,n,water);
04417     {
04418       //  Vector glhs[3*n+3];
04419 
04420       HGMatrixVector grhs2;
04421 
04422       G_q(avgab,glhs,n,m,ial,ibl);
04423 #ifdef DEBUG0
04424       iout<<iINFO << "G_q:" << std::endl<<endi;
04425       for (i=0;i<m;i++) {
04426         iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
04427       }
04428 #endif
04429       //      G_q(refab,grhs2,m,ial,ibl);
04430       // update with the masses
04431       for (j=0; j<n; j++) { // number of atoms
04432         for (i=0; i<m;i++) { // number of constraints
04433           grhs2[i][j] = grhs[i][j]*imass[j];
04434         }
04435       }
04436 
04437       // G_q(qtilde) * M^-1 G_q'(q0) =
04438       // G_q(qtilde) * grhs'
04439       for (i=0;i<m;i++) { // number of constraints
04440         for (j=0;j<m;j++) { // number of constraints
04441           fjac[i][j] = 0; 
04442           for (k1=0;k1<n;k1++) {
04443             fjac[i][j] += glhs[i][k1]*grhs2[j][k1]; 
04444           }
04445         }
04446       }
04447 #ifdef DEBUG0  
04448       iout<<iINFO << "glhs" <<endi;
04449       for(i=0;i<9;i++) {
04450         iout<<iINFO << glhs[i] << ","<<endi;
04451       }
04452       iout<<iINFO << std::endl<<endi;
04453       for(i=0;i<9;i++) {
04454         iout<<iINFO << grhs2[i] << ","<<endi;
04455       }
04456       iout<<iINFO << std::endl<<endi;
04457 #endif
04458       //      delete[] grhs2;
04459     }
04460     // end of Jac calculation
04461 #ifdef DEBUG0
04462     iout<<iINFO << "Jac" << std::endl<<endi;
04463     for (i=0;i<m;i++) 
04464       for (j=0;j<m;j++)
04465         iout<<iINFO << fjac[i][j] << " "<<endi;
04466     iout<< std::endl<<endi;
04467 #endif
04468     // calculate constraints in gij for n constraints this being a water
04469     //  G(qtilde, gij, n, water);
04470     for (i=0;i<m;i++) {
04471       gij[i]=avgab[i]*avgab[i]-length2[i];
04472     }
04473 #ifdef DEBUG0
04474     iout<<iINFO << "G" << std::endl<<endi;
04475     iout<<iINFO << "( "<<endi;
04476     for(i=0;i<m-1;i++) {
04477       iout<<iINFO << gij[i] << ", "<<endi;
04478     }
04479     iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
04480 #endif
04481     // fill the return vector
04482     for(i=0;i<m;i++) {
04483       fvec[i] = gij[i];
04484     }
04485     // free up the constraints
04486     //    delete[] gij;
04487     // continue Newton's iteration    
04488     errf=0.0;
04489     for (i=0;i<m;i++) errf += fabs(fvec[i]);
04490 #ifdef DEBUG0
04491     iout<<iINFO << "errf: " << errf << std::endl<<endi;
04492 #endif
04493     if (errf <= tolf) {
04494       break;
04495     }
04496     for (i=0;i<m;i++) p[i] = -fvec[i];
04497     //    iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
04498     ludcmp(fjac,m,indx,&d);
04499     lubksb(fjac,m,indx,p);
04500 
04501     errx=0.0;
04502     for (i=0;i<m;i++) {
04503       errx += fabs(p[i]);
04504     }
04505     for (i=0;i<m;i++)  
04506       lambda[i] += p[i];
04507 
04508 #ifdef DEBUG0
04509     iout<<iINFO << "lambda:" << lambda[0] 
04510          << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04511     iout<<iINFO << "errx: " << errx << std::endl<<endi;
04512 #endif
04513     if (errx <= tolx) break;
04514 #ifdef DEBUG0
04515     iout<<iINFO << "Qtilde:" << std::endl<<endi;
04516     iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi; 
04517 #endif
04518   }
04519 #ifdef DEBUG
04520   iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04521 #endif
04522 
04523   return k; // 
04524 }

static int compDistance ( const void *  a,
const void *  b 
) [static]

Definition at line 548 of file HomePatch.C.

Referenced by HomePatch::buildSpanningTree().

00549 {
00550   int d1 = abs(*(int *)a - CkMyPe());
00551   int d2 = abs(*(int *)b - CkMyPe());
00552   if (d1 < d2) 
00553     return -1;
00554   else if (d1 == d2) 
00555     return 0;
00556   else 
00557     return 1;
00558 }

void G_q ( const HGArrayVector refab,
HGMatrixVector gqij,
const int  n,
const int  m,
const HGArrayInt ial,
const HGArrayInt ibl 
) [inline]

Definition at line 4332 of file HomePatch.C.

Referenced by average(), and mollify().

04333                                                                              {
04334   int i; 
04335   // step through the rows of the matrix
04336   for(i=0;i<m;i++) {
04337     gqij[i][ial[i]]=2.0*refab[i];
04338     gqij[i][ibl[i]]=-gqij[i][ial[i]];
04339   }
04340 };

void lubksb ( HGMatrixBigReal a,
int  n,
HGArrayInt indx,
HGArrayBigReal b 
) [inline]

Definition at line 4259 of file HomePatch.C.

References j.

Referenced by average(), and mollify().

04261 {
04262         int i,ii=-1,ip,j;
04263         double sum;
04264 
04265         for (i=0;i<n;i++) {
04266                 ip=indx[i];
04267                 sum=b[ip];
04268                 b[ip]=b[i];
04269                 if (ii >= 0)
04270                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
04271                 else if (sum) ii=i;
04272                 b[i]=sum;
04273         }
04274         for (i=n-1;i>=0;i--) {
04275                 sum=b[i];
04276                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
04277                 b[i]=sum/a[i][i];
04278         }
04279 }

void ludcmp ( HGMatrixBigReal a,
int  n,
HGArrayInt indx,
BigReal d 
) [inline]

Definition at line 4282 of file HomePatch.C.

References j, NAMD_die(), and TINY.

Referenced by average(), and mollify().

04283 {
04284 
04285         int i,imax,j,k;
04286         double big,dum,sum,temp;
04287         HGArrayBigReal vv;
04288         *d=1.0;
04289         for (i=0;i<n;i++) {
04290                 big=0.0;
04291                 for (j=0;j<n;j++)
04292                         if ((temp=fabs(a[i][j])) > big) big=temp;
04293                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
04294                 vv[i]=1.0/big;
04295         }
04296         for (j=0;j<n;j++) {
04297                 for (i=0;i<j;i++) {
04298                         sum=a[i][j];
04299                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
04300                         a[i][j]=sum;
04301                 }
04302                 big=0.0;
04303                 for (i=j;i<n;i++) {
04304                         sum=a[i][j];
04305                         for (k=0;k<j;k++)
04306                                 sum -= a[i][k]*a[k][j];
04307                         a[i][j]=sum;
04308                         if ( (dum=vv[i]*fabs(sum)) >= big) {
04309                                 big=dum;
04310                                 imax=i;
04311                         }
04312                 }
04313                 if (j != imax) {
04314                         for (k=0;k<n;k++) {
04315                                 dum=a[imax][k];
04316                                 a[imax][k]=a[j][k];
04317                                 a[j][k]=dum;
04318                         }
04319                         *d = -(*d);
04320                         vv[imax]=vv[j];
04321                 }
04322                 indx[j]=imax;
04323                 if (a[j][j] == 0.0) a[j][j]=TINY;
04324                 if (j != n-1) {
04325                         dum=1.0/(a[j][j]);
04326                         for (i=j+1;i<n;i++) a[i][j] *= dum;
04327                 }
04328         }
04329 }

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 
)

Definition at line 4526 of file HomePatch.C.

References G_q(), j, lubksb(), ludcmp(), and y.

Referenced by HomePatch::mollyMollify().

04526                                                                                                                                                                                                                                 {
04527   int i,j,k;
04528   BigReal d;
04529   HGMatrixBigReal fjac;
04530   Vector zero(0.0,0.0,0.0);
04531   
04532   HGArrayVector tmpforce;
04533   HGArrayVector tmpforce2;
04534   HGArrayVector y;
04535   HGMatrixVector grhs;
04536   HGMatrixVector glhs;
04537   HGArrayBigReal aux;
04538   HGArrayInt indx;
04539 
04540   for(i=0;i<n;i++) {
04541     tmpforce[i]=imass[i]*force[i];
04542   }
04543 
04544   HGMatrixVector grhs2;
04545   HGArrayVector avgab;
04546 
04547   for ( i = 0; i < m; i++ ) {
04548         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04549   }
04550 
04551   G_q(avgab,glhs,n,m,ial,ibl);
04552   G_q(refab,grhs,n,m,ial,ibl);
04553   // update with the masses
04554   for (j=0; j<n; j++) { // number of atoms
04555         for (i=0; i<m;i++) { // number of constraints
04556           grhs2[i][j] = grhs[i][j]*imass[j];
04557         }
04558   }
04559 
04560   // G_q(qtilde) * M^-1 G_q'(q0) =
04561   // G_q(qtilde) * grhs'
04562   for (i=0;i<m;i++) { // number of constraints
04563         for (j=0;j<m;j++) { // number of constraints
04564           fjac[j][i] = 0; 
04565           for (k=0;k<n;k++) {
04566             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
04567           }
04568         }
04569   }
04570 
04571   // aux=gqij*tmpforce
04572   //  globalGrhs::computeGlobalGrhs(q0,n,water);
04573   //  G_q(refab,grhs,m,ial,ibl);
04574   for(i=0;i<m;i++) {
04575     aux[i]=0.0;
04576     for(j=0;j<n;j++) {
04577       aux[i]+=grhs[i][j]*tmpforce[j];
04578     }
04579   }
04580 
04581   ludcmp(fjac,m,indx,&d);
04582   lubksb(fjac,m,indx,aux);
04583 
04584   for(j=0;j<n;j++) {
04585     y[j] = zero;
04586     for(i=0;i<m;i++) {
04587       y[j] += aux[i]*glhs[i][j];
04588     }
04589   }
04590   for(i=0;i<n;i++) {
04591     y[i]=force[i]-y[i];
04592   }
04593     
04594   // gqq12*y
04595   for(i=0;i<n;i++) {
04596     tmpforce2[i]=imass[i]*y[i];
04597   }
04598 
04599   // here we assume that tmpforce is initialized to zero.
04600   for (i=0;i<n;i++) {
04601     tmpforce[i]=zero;
04602   }
04603   
04604   for (j=0;j<m;j++) {
04605     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
04606     tmpforce[ial[j]] += tmpf;
04607     tmpforce[ibl[j]] -= tmpf;
04608   }
04609   // c-ji the other bug for 2 constraint water was this line (2-4-99)
04610   //  for(i=0;i<m;i++) {
04611   for(i=0;i<n;i++) {
04612     force[i]=tmpforce[i]+y[i];
04613   }
04614 
04615 }


Generated on Mon Nov 20 01:17:15 2017 for NAMD by  doxygen 1.4.7