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 "NamdEventsProfiling.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

Redistribute forces from lone pair site onto its host atoms.

Atoms are "labeled" i, j, k, l, where atom i is the lone pair. Positions of atoms are ri, rj, rk, rl. Forces of atoms are fi, fj, fk, fl; updated forces are returned. Accumulate updates to the virial.

The forces on the atoms are updated so that:

  • the force fi on the lone pair site is 0
  • the net force fi+fj+fk+fl is conserved
  • the net torque cross(ri,fi)+cross(rj,fj)+cross(rk,fk)+cross(rl,fl) is conserved

If "midpt" is true (nonzero), then use the midpoint of rk and rl (e.g. rk and rl are the hydrogen atoms for water). Otherwise use planes defined by ri,rj,rk and rj,rk,rl.

Having "midpt" set true corresponds in CHARMM to having a negative distance parameter in Lphost.

Use FIX_FOR_WATER below to fix the case that occurs when the lone pair site for water lies on the perpendicular bisector of rk and rl, making the cross(r,s) zero.

Definition at line 1371 of file HomePatch.C.

#define MASS_EPSILON   (1.0e-35)

Definition at line 74 of file HomePatch.C.

#define MAXHGS   10

Definition at line 52 of file HomePatch.C.

#define MIN_DEBUG_LEVEL   2

Definition at line 53 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 51 of file HomePatch.C.

Referenced by ludcmp().


Typedef Documentation

typedef BigReal HGArrayBigReal[MAXHGS]

Definition at line 63 of file HomePatch.C.

typedef int HGArrayInt[MAXHGS]

Definition at line 62 of file HomePatch.C.

typedef zVector HGArrayVector[MAXHGS]

Definition at line 64 of file HomePatch.C.

typedef BigReal HGMatrixBigReal[MAXHGS][MAXHGS]

Definition at line 65 of file HomePatch.C.

typedef zVector HGMatrixVector[MAXHGS][MAXHGS]

Definition at line 66 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 4402 of file HomePatch.C.

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

Referenced by HomePatch::mollyAverage().

04402                                                                                                                                                                                                                                                                          {
04403   //  input:  n = length of hydrogen group to be averaged (shaked)
04404   //          q[n] = vector array of original positions
04405   //          m = number of constraints
04406   //          imass[n] = inverse mass for each atom
04407   //          length2[m] = square of reference bond length for each constraint
04408   //          ial[m] = atom a in each constraint 
04409   //          ibl[m] = atom b in each constraint 
04410   //          refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
04411   //          tolf = function error tolerance for Newton's iteration
04412   //          ntrial = max number of Newton's iterations
04413   //  output: lambda[m] = double array of lagrange multipliers (used by mollify)
04414   //          qtilde[n] = vector array of averaged (shaked) positions
04415 
04416   int k,k1,i,j;
04417   BigReal errx,errf,d,tolx;
04418 
04419   HGArrayInt indx;
04420   HGArrayBigReal p;
04421   HGArrayBigReal fvec;
04422   HGMatrixBigReal fjac;
04423   HGArrayVector avgab;
04424   HGMatrixVector grhs;
04425   HGMatrixVector auxrhs;
04426   HGMatrixVector glhs;
04427 
04428   //  iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
04429   tolx=tolf; 
04430   
04431   // initialize lambda, globalGrhs
04432 
04433   for (i=0;i<m;i++) {
04434     lambda[i]=0.0;
04435   }
04436 
04437   // define grhs, auxrhs for all iterations
04438   // grhs= g_x(q)
04439   //
04440   G_q(refab,grhs,n,m,ial,ibl);
04441   for (k=1;k<=ntrial;k++) {
04442     //    usrfun(qtilde,q0,lambda,fvec,fjac,n,water); 
04443     HGArrayBigReal gij;
04444     // this used to be main loop of usrfun
04445     // compute qtilde given q0, lambda, IMASSes
04446     {
04447       BigReal multiplier;
04448       HGArrayVector tmp;
04449       for (i=0;i<m;i++) {
04450         multiplier = lambda[i];
04451         // auxrhs = M^{-1}grhs^{T}
04452         for (j=0;j<n;j++) {
04453           auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
04454         }
04455       }
04456       for (j=0;j<n;j++) {
04457         //      tmp[j]=0.0;      
04458         for (i=0;i<m;i++) {
04459           tmp[j]+=auxrhs[i][j];
04460         }
04461       }
04462  
04463       for (j=0;j<n;j++) {
04464         qtilde[j].position=q[j]+tmp[j];
04465       }
04466       //      delete [] tmp;
04467     }
04468   
04469     for ( i = 0; i < m; i++ ) {
04470       avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04471     }
04472 
04473     //  iout<<iINFO << "Calling Jac" << std::endl<<endi;
04474     //  Jac(qtilde, q0, fjac,n,water);
04475     {
04476       //  Vector glhs[3*n+3];
04477 
04478       HGMatrixVector grhs2;
04479 
04480       G_q(avgab,glhs,n,m,ial,ibl);
04481 #ifdef DEBUG0
04482       iout<<iINFO << "G_q:" << std::endl<<endi;
04483       for (i=0;i<m;i++) {
04484         iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
04485       }
04486 #endif
04487       //      G_q(refab,grhs2,m,ial,ibl);
04488       // update with the masses
04489       for (j=0; j<n; j++) { // number of atoms
04490         for (i=0; i<m;i++) { // number of constraints
04491           grhs2[i][j] = grhs[i][j]*imass[j];
04492         }
04493       }
04494 
04495       // G_q(qtilde) * M^-1 G_q'(q0) =
04496       // G_q(qtilde) * grhs'
04497       for (i=0;i<m;i++) { // number of constraints
04498         for (j=0;j<m;j++) { // number of constraints
04499           fjac[i][j] = 0; 
04500           for (k1=0;k1<n;k1++) {
04501             fjac[i][j] += glhs[i][k1]*grhs2[j][k1]; 
04502           }
04503         }
04504       }
04505 #ifdef DEBUG0  
04506       iout<<iINFO << "glhs" <<endi;
04507       for(i=0;i<9;i++) {
04508         iout<<iINFO << glhs[i] << ","<<endi;
04509       }
04510       iout<<iINFO << std::endl<<endi;
04511       for(i=0;i<9;i++) {
04512         iout<<iINFO << grhs2[i] << ","<<endi;
04513       }
04514       iout<<iINFO << std::endl<<endi;
04515 #endif
04516       //      delete[] grhs2;
04517     }
04518     // end of Jac calculation
04519 #ifdef DEBUG0
04520     iout<<iINFO << "Jac" << std::endl<<endi;
04521     for (i=0;i<m;i++) 
04522       for (j=0;j<m;j++)
04523         iout<<iINFO << fjac[i][j] << " "<<endi;
04524     iout<< std::endl<<endi;
04525 #endif
04526     // calculate constraints in gij for n constraints this being a water
04527     //  G(qtilde, gij, n, water);
04528     for (i=0;i<m;i++) {
04529       gij[i]=avgab[i]*avgab[i]-length2[i];
04530     }
04531 #ifdef DEBUG0
04532     iout<<iINFO << "G" << std::endl<<endi;
04533     iout<<iINFO << "( "<<endi;
04534     for(i=0;i<m-1;i++) {
04535       iout<<iINFO << gij[i] << ", "<<endi;
04536     }
04537     iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
04538 #endif
04539     // fill the return vector
04540     for(i=0;i<m;i++) {
04541       fvec[i] = gij[i];
04542     }
04543     // free up the constraints
04544     //    delete[] gij;
04545     // continue Newton's iteration    
04546     errf=0.0;
04547     for (i=0;i<m;i++) errf += fabs(fvec[i]);
04548 #ifdef DEBUG0
04549     iout<<iINFO << "errf: " << errf << std::endl<<endi;
04550 #endif
04551     if (errf <= tolf) {
04552       break;
04553     }
04554     for (i=0;i<m;i++) p[i] = -fvec[i];
04555     //    iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
04556     ludcmp(fjac,m,indx,&d);
04557     lubksb(fjac,m,indx,p);
04558 
04559     errx=0.0;
04560     for (i=0;i<m;i++) {
04561       errx += fabs(p[i]);
04562     }
04563     for (i=0;i<m;i++)  
04564       lambda[i] += p[i];
04565 
04566 #ifdef DEBUG0
04567     iout<<iINFO << "lambda:" << lambda[0] 
04568          << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04569     iout<<iINFO << "errx: " << errx << std::endl<<endi;
04570 #endif
04571     if (errx <= tolx) break;
04572 #ifdef DEBUG0
04573     iout<<iINFO << "Qtilde:" << std::endl<<endi;
04574     iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi; 
04575 #endif
04576   }
04577 #ifdef DEBUG
04578   iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04579 #endif
04580 
04581   return k; // 
04582 }

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

Definition at line 456 of file HomePatch.C.

Referenced by HomePatch::buildSpanningTree().

00457 {
00458   int d1 = abs(*(int *)a - CkMyPe());
00459   int d2 = abs(*(int *)b - CkMyPe());
00460   if (d1 < d2) 
00461     return -1;
00462   else if (d1 == d2) 
00463     return 0;
00464   else 
00465     return 1;
00466 }

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

Definition at line 4390 of file HomePatch.C.

Referenced by average(), and mollify().

04391                                                                              {
04392   int i; 
04393   // step through the rows of the matrix
04394   for(i=0;i<m;i++) {
04395     gqij[i][ial[i]]=2.0*refab[i];
04396     gqij[i][ibl[i]]=-gqij[i][ial[i]];
04397   }
04398 };

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

Definition at line 4317 of file HomePatch.C.

Referenced by average(), and mollify().

04319 {
04320         int i,ii=-1,ip,j;
04321         double sum;
04322 
04323         for (i=0;i<n;i++) {
04324                 ip=indx[i];
04325                 sum=b[ip];
04326                 b[ip]=b[i];
04327                 if (ii >= 0)
04328                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
04329                 else if (sum) ii=i;
04330                 b[i]=sum;
04331         }
04332         for (i=n-1;i>=0;i--) {
04333                 sum=b[i];
04334                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
04335                 b[i]=sum/a[i][i];
04336         }
04337 }

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

Definition at line 4340 of file HomePatch.C.

References NAMD_die(), and TINY.

Referenced by average(), and mollify().

04341 {
04342 
04343         int i,imax,j,k;
04344         double big,dum,sum,temp;
04345         HGArrayBigReal vv;
04346         *d=1.0;
04347         for (i=0;i<n;i++) {
04348                 big=0.0;
04349                 for (j=0;j<n;j++)
04350                         if ((temp=fabs(a[i][j])) > big) big=temp;
04351                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
04352                 vv[i]=1.0/big;
04353         }
04354         for (j=0;j<n;j++) {
04355                 for (i=0;i<j;i++) {
04356                         sum=a[i][j];
04357                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
04358                         a[i][j]=sum;
04359                 }
04360                 big=0.0;
04361                 for (i=j;i<n;i++) {
04362                         sum=a[i][j];
04363                         for (k=0;k<j;k++)
04364                                 sum -= a[i][k]*a[k][j];
04365                         a[i][j]=sum;
04366                         if ( (dum=vv[i]*fabs(sum)) >= big) {
04367                                 big=dum;
04368                                 imax=i;
04369                         }
04370                 }
04371                 if (j != imax) {
04372                         for (k=0;k<n;k++) {
04373                                 dum=a[imax][k];
04374                                 a[imax][k]=a[j][k];
04375                                 a[j][k]=dum;
04376                         }
04377                         *d = -(*d);
04378                         vv[imax]=vv[j];
04379                 }
04380                 indx[j]=imax;
04381                 if (a[j][j] == 0.0) a[j][j]=TINY;
04382                 if (j != n-1) {
04383                         dum=1.0/(a[j][j]);
04384                         for (i=j+1;i<n;i++) a[i][j] *= dum;
04385                 }
04386         }
04387 }

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

References G_q(), lubksb(), ludcmp(), and CompAtom::position.

Referenced by HomePatch::mollyMollify().

04584                                                                                                                                                                                                                                 {
04585   int i,j,k;
04586   BigReal d;
04587   HGMatrixBigReal fjac;
04588   Vector zero(0.0,0.0,0.0);
04589   
04590   HGArrayVector tmpforce;
04591   HGArrayVector tmpforce2;
04592   HGArrayVector y;
04593   HGMatrixVector grhs;
04594   HGMatrixVector glhs;
04595   HGArrayBigReal aux;
04596   HGArrayInt indx;
04597 
04598   for(i=0;i<n;i++) {
04599     tmpforce[i]=imass[i]*force[i];
04600   }
04601 
04602   HGMatrixVector grhs2;
04603   HGArrayVector avgab;
04604 
04605   for ( i = 0; i < m; i++ ) {
04606         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04607   }
04608 
04609   G_q(avgab,glhs,n,m,ial,ibl);
04610   G_q(refab,grhs,n,m,ial,ibl);
04611   // update with the masses
04612   for (j=0; j<n; j++) { // number of atoms
04613         for (i=0; i<m;i++) { // number of constraints
04614           grhs2[i][j] = grhs[i][j]*imass[j];
04615         }
04616   }
04617 
04618   // G_q(qtilde) * M^-1 G_q'(q0) =
04619   // G_q(qtilde) * grhs'
04620   for (i=0;i<m;i++) { // number of constraints
04621         for (j=0;j<m;j++) { // number of constraints
04622           fjac[j][i] = 0; 
04623           for (k=0;k<n;k++) {
04624             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
04625           }
04626         }
04627   }
04628 
04629   // aux=gqij*tmpforce
04630   //  globalGrhs::computeGlobalGrhs(q0,n,water);
04631   //  G_q(refab,grhs,m,ial,ibl);
04632   for(i=0;i<m;i++) {
04633     aux[i]=0.0;
04634     for(j=0;j<n;j++) {
04635       aux[i]+=grhs[i][j]*tmpforce[j];
04636     }
04637   }
04638 
04639   ludcmp(fjac,m,indx,&d);
04640   lubksb(fjac,m,indx,aux);
04641 
04642   for(j=0;j<n;j++) {
04643     y[j] = zero;
04644     for(i=0;i<m;i++) {
04645       y[j] += aux[i]*glhs[i][j];
04646     }
04647   }
04648   for(i=0;i<n;i++) {
04649     y[i]=force[i]-y[i];
04650   }
04651     
04652   // gqq12*y
04653   for(i=0;i<n;i++) {
04654     tmpforce2[i]=imass[i]*y[i];
04655   }
04656 
04657   // here we assume that tmpforce is initialized to zero.
04658   for (i=0;i<n;i++) {
04659     tmpforce[i]=zero;
04660   }
04661   
04662   for (j=0;j<m;j++) {
04663     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
04664     tmpforce[ial[j]] += tmpf;
04665     tmpforce[ibl[j]] -= tmpf;
04666   }
04667   // c-ji the other bug for 2 constraint water was this line (2-4-99)
04668   //  for(i=0;i<m;i++) {
04669   for(i=0;i<n;i++) {
04670     force[i]=tmpforce[i]+y[i];
04671   }
04672 
04673 }


Generated on 18 Sep 2019 for NAMD by  doxygen 1.6.1