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 1302 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 4235 of file HomePatch.C.

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

Referenced by HomePatch::mollyAverage().

04235                                                                                                                                                                                                                                                                          {
04236   //  input:  n = length of hydrogen group to be averaged (shaked)
04237   //          q[n] = vector array of original positions
04238   //          m = number of constraints
04239   //          imass[n] = inverse mass for each atom
04240   //          length2[m] = square of reference bond length for each constraint
04241   //          ial[m] = atom a in each constraint 
04242   //          ibl[m] = atom b in each constraint 
04243   //          refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
04244   //          tolf = function error tolerance for Newton's iteration
04245   //          ntrial = max number of Newton's iterations
04246   //  output: lambda[m] = double array of lagrange multipliers (used by mollify)
04247   //          qtilde[n] = vector array of averaged (shaked) positions
04248 
04249   int k,k1,i,j;
04250   BigReal errx,errf,d,tolx;
04251 
04252   HGArrayInt indx;
04253   HGArrayBigReal p;
04254   HGArrayBigReal fvec;
04255   HGMatrixBigReal fjac;
04256   HGArrayVector avgab;
04257   HGMatrixVector grhs;
04258   HGMatrixVector auxrhs;
04259   HGMatrixVector glhs;
04260 
04261   //  iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
04262   tolx=tolf; 
04263   
04264   // initialize lambda, globalGrhs
04265 
04266   for (i=0;i<m;i++) {
04267     lambda[i]=0.0;
04268   }
04269 
04270   // define grhs, auxrhs for all iterations
04271   // grhs= g_x(q)
04272   //
04273   G_q(refab,grhs,n,m,ial,ibl);
04274   for (k=1;k<=ntrial;k++) {
04275     //    usrfun(qtilde,q0,lambda,fvec,fjac,n,water); 
04276     HGArrayBigReal gij;
04277     // this used to be main loop of usrfun
04278     // compute qtilde given q0, lambda, IMASSes
04279     {
04280       BigReal multiplier;
04281       HGArrayVector tmp;
04282       for (i=0;i<m;i++) {
04283         multiplier = lambda[i];
04284         // auxrhs = M^{-1}grhs^{T}
04285         for (j=0;j<n;j++) {
04286           auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
04287         }
04288       }
04289       for (j=0;j<n;j++) {
04290         //      tmp[j]=0.0;      
04291         for (i=0;i<m;i++) {
04292           tmp[j]+=auxrhs[i][j];
04293         }
04294       }
04295  
04296       for (j=0;j<n;j++) {
04297         qtilde[j].position=q[j]+tmp[j];
04298       }
04299       //      delete [] tmp;
04300     }
04301   
04302     for ( i = 0; i < m; i++ ) {
04303       avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04304     }
04305 
04306     //  iout<<iINFO << "Calling Jac" << std::endl<<endi;
04307     //  Jac(qtilde, q0, fjac,n,water);
04308     {
04309       //  Vector glhs[3*n+3];
04310 
04311       HGMatrixVector grhs2;
04312 
04313       G_q(avgab,glhs,n,m,ial,ibl);
04314 #ifdef DEBUG0
04315       iout<<iINFO << "G_q:" << std::endl<<endi;
04316       for (i=0;i<m;i++) {
04317         iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
04318       }
04319 #endif
04320       //      G_q(refab,grhs2,m,ial,ibl);
04321       // update with the masses
04322       for (j=0; j<n; j++) { // number of atoms
04323         for (i=0; i<m;i++) { // number of constraints
04324           grhs2[i][j] = grhs[i][j]*imass[j];
04325         }
04326       }
04327 
04328       // G_q(qtilde) * M^-1 G_q'(q0) =
04329       // G_q(qtilde) * grhs'
04330       for (i=0;i<m;i++) { // number of constraints
04331         for (j=0;j<m;j++) { // number of constraints
04332           fjac[i][j] = 0; 
04333           for (k1=0;k1<n;k1++) {
04334             fjac[i][j] += glhs[i][k1]*grhs2[j][k1]; 
04335           }
04336         }
04337       }
04338 #ifdef DEBUG0  
04339       iout<<iINFO << "glhs" <<endi;
04340       for(i=0;i<9;i++) {
04341         iout<<iINFO << glhs[i] << ","<<endi;
04342       }
04343       iout<<iINFO << std::endl<<endi;
04344       for(i=0;i<9;i++) {
04345         iout<<iINFO << grhs2[i] << ","<<endi;
04346       }
04347       iout<<iINFO << std::endl<<endi;
04348 #endif
04349       //      delete[] grhs2;
04350     }
04351     // end of Jac calculation
04352 #ifdef DEBUG0
04353     iout<<iINFO << "Jac" << std::endl<<endi;
04354     for (i=0;i<m;i++) 
04355       for (j=0;j<m;j++)
04356         iout<<iINFO << fjac[i][j] << " "<<endi;
04357     iout<< std::endl<<endi;
04358 #endif
04359     // calculate constraints in gij for n constraints this being a water
04360     //  G(qtilde, gij, n, water);
04361     for (i=0;i<m;i++) {
04362       gij[i]=avgab[i]*avgab[i]-length2[i];
04363     }
04364 #ifdef DEBUG0
04365     iout<<iINFO << "G" << std::endl<<endi;
04366     iout<<iINFO << "( "<<endi;
04367     for(i=0;i<m-1;i++) {
04368       iout<<iINFO << gij[i] << ", "<<endi;
04369     }
04370     iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
04371 #endif
04372     // fill the return vector
04373     for(i=0;i<m;i++) {
04374       fvec[i] = gij[i];
04375     }
04376     // free up the constraints
04377     //    delete[] gij;
04378     // continue Newton's iteration    
04379     errf=0.0;
04380     for (i=0;i<m;i++) errf += fabs(fvec[i]);
04381 #ifdef DEBUG0
04382     iout<<iINFO << "errf: " << errf << std::endl<<endi;
04383 #endif
04384     if (errf <= tolf) {
04385       break;
04386     }
04387     for (i=0;i<m;i++) p[i] = -fvec[i];
04388     //    iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
04389     ludcmp(fjac,m,indx,&d);
04390     lubksb(fjac,m,indx,p);
04391 
04392     errx=0.0;
04393     for (i=0;i<m;i++) {
04394       errx += fabs(p[i]);
04395     }
04396     for (i=0;i<m;i++)  
04397       lambda[i] += p[i];
04398 
04399 #ifdef DEBUG0
04400     iout<<iINFO << "lambda:" << lambda[0] 
04401          << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04402     iout<<iINFO << "errx: " << errx << std::endl<<endi;
04403 #endif
04404     if (errx <= tolx) break;
04405 #ifdef DEBUG0
04406     iout<<iINFO << "Qtilde:" << std::endl<<endi;
04407     iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi; 
04408 #endif
04409   }
04410 #ifdef DEBUG
04411   iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04412 #endif
04413 
04414   return k; // 
04415 }

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

Definition at line 439 of file HomePatch.C.

Referenced by HomePatch::buildSpanningTree().

00440 {
00441   int d1 = abs(*(int *)a - CkMyPe());
00442   int d2 = abs(*(int *)b - CkMyPe());
00443   if (d1 < d2) 
00444     return -1;
00445   else if (d1 == d2) 
00446     return 0;
00447   else 
00448     return 1;
00449 }

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

Definition at line 4223 of file HomePatch.C.

Referenced by average(), and mollify().

04224                                                                              {
04225   int i; 
04226   // step through the rows of the matrix
04227   for(i=0;i<m;i++) {
04228     gqij[i][ial[i]]=2.0*refab[i];
04229     gqij[i][ibl[i]]=-gqij[i][ial[i]];
04230   }
04231 };

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

Definition at line 4150 of file HomePatch.C.

References j.

Referenced by average(), and mollify().

04152 {
04153         int i,ii=-1,ip,j;
04154         double sum;
04155 
04156         for (i=0;i<n;i++) {
04157                 ip=indx[i];
04158                 sum=b[ip];
04159                 b[ip]=b[i];
04160                 if (ii >= 0)
04161                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
04162                 else if (sum) ii=i;
04163                 b[i]=sum;
04164         }
04165         for (i=n-1;i>=0;i--) {
04166                 sum=b[i];
04167                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
04168                 b[i]=sum/a[i][i];
04169         }
04170 }

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

Definition at line 4173 of file HomePatch.C.

References j, NAMD_die(), and TINY.

Referenced by average(), and mollify().

04174 {
04175 
04176         int i,imax,j,k;
04177         double big,dum,sum,temp;
04178         HGArrayBigReal vv;
04179         *d=1.0;
04180         for (i=0;i<n;i++) {
04181                 big=0.0;
04182                 for (j=0;j<n;j++)
04183                         if ((temp=fabs(a[i][j])) > big) big=temp;
04184                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
04185                 vv[i]=1.0/big;
04186         }
04187         for (j=0;j<n;j++) {
04188                 for (i=0;i<j;i++) {
04189                         sum=a[i][j];
04190                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
04191                         a[i][j]=sum;
04192                 }
04193                 big=0.0;
04194                 for (i=j;i<n;i++) {
04195                         sum=a[i][j];
04196                         for (k=0;k<j;k++)
04197                                 sum -= a[i][k]*a[k][j];
04198                         a[i][j]=sum;
04199                         if ( (dum=vv[i]*fabs(sum)) >= big) {
04200                                 big=dum;
04201                                 imax=i;
04202                         }
04203                 }
04204                 if (j != imax) {
04205                         for (k=0;k<n;k++) {
04206                                 dum=a[imax][k];
04207                                 a[imax][k]=a[j][k];
04208                                 a[j][k]=dum;
04209                         }
04210                         *d = -(*d);
04211                         vv[imax]=vv[j];
04212                 }
04213                 indx[j]=imax;
04214                 if (a[j][j] == 0.0) a[j][j]=TINY;
04215                 if (j != n-1) {
04216                         dum=1.0/(a[j][j]);
04217                         for (i=j+1;i<n;i++) a[i][j] *= dum;
04218                 }
04219         }
04220 }

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 4417 of file HomePatch.C.

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

Referenced by HomePatch::mollyMollify().

04417                                                                                                                                                                                                                                 {
04418   int i,j,k;
04419   BigReal d;
04420   HGMatrixBigReal fjac;
04421   Vector zero(0.0,0.0,0.0);
04422   
04423   HGArrayVector tmpforce;
04424   HGArrayVector tmpforce2;
04425   HGArrayVector y;
04426   HGMatrixVector grhs;
04427   HGMatrixVector glhs;
04428   HGArrayBigReal aux;
04429   HGArrayInt indx;
04430 
04431   for(i=0;i<n;i++) {
04432     tmpforce[i]=imass[i]*force[i];
04433   }
04434 
04435   HGMatrixVector grhs2;
04436   HGArrayVector avgab;
04437 
04438   for ( i = 0; i < m; i++ ) {
04439         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04440   }
04441 
04442   G_q(avgab,glhs,n,m,ial,ibl);
04443   G_q(refab,grhs,n,m,ial,ibl);
04444   // update with the masses
04445   for (j=0; j<n; j++) { // number of atoms
04446         for (i=0; i<m;i++) { // number of constraints
04447           grhs2[i][j] = grhs[i][j]*imass[j];
04448         }
04449   }
04450 
04451   // G_q(qtilde) * M^-1 G_q'(q0) =
04452   // G_q(qtilde) * grhs'
04453   for (i=0;i<m;i++) { // number of constraints
04454         for (j=0;j<m;j++) { // number of constraints
04455           fjac[j][i] = 0; 
04456           for (k=0;k<n;k++) {
04457             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
04458           }
04459         }
04460   }
04461 
04462   // aux=gqij*tmpforce
04463   //  globalGrhs::computeGlobalGrhs(q0,n,water);
04464   //  G_q(refab,grhs,m,ial,ibl);
04465   for(i=0;i<m;i++) {
04466     aux[i]=0.0;
04467     for(j=0;j<n;j++) {
04468       aux[i]+=grhs[i][j]*tmpforce[j];
04469     }
04470   }
04471 
04472   ludcmp(fjac,m,indx,&d);
04473   lubksb(fjac,m,indx,aux);
04474 
04475   for(j=0;j<n;j++) {
04476     y[j] = zero;
04477     for(i=0;i<m;i++) {
04478       y[j] += aux[i]*glhs[i][j];
04479     }
04480   }
04481   for(i=0;i<n;i++) {
04482     y[i]=force[i]-y[i];
04483   }
04484     
04485   // gqq12*y
04486   for(i=0;i<n;i++) {
04487     tmpforce2[i]=imass[i]*y[i];
04488   }
04489 
04490   // here we assume that tmpforce is initialized to zero.
04491   for (i=0;i<n;i++) {
04492     tmpforce[i]=zero;
04493   }
04494   
04495   for (j=0;j<m;j++) {
04496     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
04497     tmpforce[ial[j]] += tmpf;
04498     tmpforce[ibl[j]] -= tmpf;
04499   }
04500   // c-ji the other bug for 2 constraint water was this line (2-4-99)
04501   //  for(i=0;i<m;i++) {
04502   for(i=0;i<n;i++) {
04503     force[i]=tmpforce[i]+y[i];
04504   }
04505 
04506 }


Generated on Fri Jul 20 01:17:17 2018 for NAMD by  doxygen 1.4.7