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

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

Referenced by HomePatch::mollyAverage().

04490                                                                                                                                                                                                                                                                          {
04491   //  input:  n = length of hydrogen group to be averaged (shaked)
04492   //          q[n] = vector array of original positions
04493   //          m = number of constraints
04494   //          imass[n] = inverse mass for each atom
04495   //          length2[m] = square of reference bond length for each constraint
04496   //          ial[m] = atom a in each constraint 
04497   //          ibl[m] = atom b in each constraint 
04498   //          refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
04499   //          tolf = function error tolerance for Newton's iteration
04500   //          ntrial = max number of Newton's iterations
04501   //  output: lambda[m] = double array of lagrange multipliers (used by mollify)
04502   //          qtilde[n] = vector array of averaged (shaked) positions
04503 
04504   int k,k1,i,j;
04505   BigReal errx,errf,d,tolx;
04506 
04507   HGArrayInt indx;
04508   HGArrayBigReal p;
04509   HGArrayBigReal fvec;
04510   HGMatrixBigReal fjac;
04511   HGArrayVector avgab;
04512   HGMatrixVector grhs;
04513   HGMatrixVector auxrhs;
04514   HGMatrixVector glhs;
04515 
04516   //  iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
04517   tolx=tolf; 
04518   
04519   // initialize lambda, globalGrhs
04520 
04521   for (i=0;i<m;i++) {
04522     lambda[i]=0.0;
04523   }
04524 
04525   // define grhs, auxrhs for all iterations
04526   // grhs= g_x(q)
04527   //
04528   G_q(refab,grhs,n,m,ial,ibl);
04529   for (k=1;k<=ntrial;k++) {
04530     //    usrfun(qtilde,q0,lambda,fvec,fjac,n,water); 
04531     HGArrayBigReal gij;
04532     // this used to be main loop of usrfun
04533     // compute qtilde given q0, lambda, IMASSes
04534     {
04535       BigReal multiplier;
04536       HGArrayVector tmp;
04537       for (i=0;i<m;i++) {
04538         multiplier = lambda[i];
04539         // auxrhs = M^{-1}grhs^{T}
04540         for (j=0;j<n;j++) {
04541           auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
04542         }
04543       }
04544       for (j=0;j<n;j++) {
04545         //      tmp[j]=0.0;      
04546         for (i=0;i<m;i++) {
04547           tmp[j]+=auxrhs[i][j];
04548         }
04549       }
04550  
04551       for (j=0;j<n;j++) {
04552         qtilde[j].position=q[j]+tmp[j];
04553       }
04554       //      delete [] tmp;
04555     }
04556   
04557     for ( i = 0; i < m; i++ ) {
04558       avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04559     }
04560 
04561     //  iout<<iINFO << "Calling Jac" << std::endl<<endi;
04562     //  Jac(qtilde, q0, fjac,n,water);
04563     {
04564       //  Vector glhs[3*n+3];
04565 
04566       HGMatrixVector grhs2;
04567 
04568       G_q(avgab,glhs,n,m,ial,ibl);
04569 #ifdef DEBUG0
04570       iout<<iINFO << "G_q:" << std::endl<<endi;
04571       for (i=0;i<m;i++) {
04572         iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
04573       }
04574 #endif
04575       //      G_q(refab,grhs2,m,ial,ibl);
04576       // update with the masses
04577       for (j=0; j<n; j++) { // number of atoms
04578         for (i=0; i<m;i++) { // number of constraints
04579           grhs2[i][j] = grhs[i][j]*imass[j];
04580         }
04581       }
04582 
04583       // G_q(qtilde) * M^-1 G_q'(q0) =
04584       // G_q(qtilde) * grhs'
04585       for (i=0;i<m;i++) { // number of constraints
04586         for (j=0;j<m;j++) { // number of constraints
04587           fjac[i][j] = 0; 
04588           for (k1=0;k1<n;k1++) {
04589             fjac[i][j] += glhs[i][k1]*grhs2[j][k1]; 
04590           }
04591         }
04592       }
04593 #ifdef DEBUG0  
04594       iout<<iINFO << "glhs" <<endi;
04595       for(i=0;i<9;i++) {
04596         iout<<iINFO << glhs[i] << ","<<endi;
04597       }
04598       iout<<iINFO << std::endl<<endi;
04599       for(i=0;i<9;i++) {
04600         iout<<iINFO << grhs2[i] << ","<<endi;
04601       }
04602       iout<<iINFO << std::endl<<endi;
04603 #endif
04604       //      delete[] grhs2;
04605     }
04606     // end of Jac calculation
04607 #ifdef DEBUG0
04608     iout<<iINFO << "Jac" << std::endl<<endi;
04609     for (i=0;i<m;i++) 
04610       for (j=0;j<m;j++)
04611         iout<<iINFO << fjac[i][j] << " "<<endi;
04612     iout<< std::endl<<endi;
04613 #endif
04614     // calculate constraints in gij for n constraints this being a water
04615     //  G(qtilde, gij, n, water);
04616     for (i=0;i<m;i++) {
04617       gij[i]=avgab[i]*avgab[i]-length2[i];
04618     }
04619 #ifdef DEBUG0
04620     iout<<iINFO << "G" << std::endl<<endi;
04621     iout<<iINFO << "( "<<endi;
04622     for(i=0;i<m-1;i++) {
04623       iout<<iINFO << gij[i] << ", "<<endi;
04624     }
04625     iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
04626 #endif
04627     // fill the return vector
04628     for(i=0;i<m;i++) {
04629       fvec[i] = gij[i];
04630     }
04631     // free up the constraints
04632     //    delete[] gij;
04633     // continue Newton's iteration    
04634     errf=0.0;
04635     for (i=0;i<m;i++) errf += fabs(fvec[i]);
04636 #ifdef DEBUG0
04637     iout<<iINFO << "errf: " << errf << std::endl<<endi;
04638 #endif
04639     if (errf <= tolf) {
04640       break;
04641     }
04642     for (i=0;i<m;i++) p[i] = -fvec[i];
04643     //    iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
04644     ludcmp(fjac,m,indx,&d);
04645     lubksb(fjac,m,indx,p);
04646 
04647     errx=0.0;
04648     for (i=0;i<m;i++) {
04649       errx += fabs(p[i]);
04650     }
04651     for (i=0;i<m;i++)  
04652       lambda[i] += p[i];
04653 
04654 #ifdef DEBUG0
04655     iout<<iINFO << "lambda:" << lambda[0] 
04656          << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04657     iout<<iINFO << "errx: " << errx << std::endl<<endi;
04658 #endif
04659     if (errx <= tolx) break;
04660 #ifdef DEBUG0
04661     iout<<iINFO << "Qtilde:" << std::endl<<endi;
04662     iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi; 
04663 #endif
04664   }
04665 #ifdef DEBUG
04666   iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
04667 #endif
04668 
04669   return k; // 
04670 }

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

Referenced by average(), and mollify().

04479                                                                              {
04480   int i; 
04481   // step through the rows of the matrix
04482   for(i=0;i<m;i++) {
04483     gqij[i][ial[i]]=2.0*refab[i];
04484     gqij[i][ibl[i]]=-gqij[i][ial[i]];
04485   }
04486 };

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

Definition at line 4405 of file HomePatch.C.

Referenced by average(), and mollify().

04407 {
04408         int i,ii=-1,ip,j;
04409         double sum;
04410 
04411         for (i=0;i<n;i++) {
04412                 ip=indx[i];
04413                 sum=b[ip];
04414                 b[ip]=b[i];
04415                 if (ii >= 0)
04416                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
04417                 else if (sum) ii=i;
04418                 b[i]=sum;
04419         }
04420         for (i=n-1;i>=0;i--) {
04421                 sum=b[i];
04422                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
04423                 b[i]=sum/a[i][i];
04424         }
04425 }

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

Definition at line 4428 of file HomePatch.C.

References NAMD_die(), and TINY.

Referenced by average(), and mollify().

04429 {
04430 
04431         int i,imax,j,k;
04432         double big,dum,sum,temp;
04433         HGArrayBigReal vv;
04434         *d=1.0;
04435         for (i=0;i<n;i++) {
04436                 big=0.0;
04437                 for (j=0;j<n;j++)
04438                         if ((temp=fabs(a[i][j])) > big) big=temp;
04439                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
04440                 vv[i]=1.0/big;
04441         }
04442         for (j=0;j<n;j++) {
04443                 for (i=0;i<j;i++) {
04444                         sum=a[i][j];
04445                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
04446                         a[i][j]=sum;
04447                 }
04448                 big=0.0;
04449                 for (i=j;i<n;i++) {
04450                         sum=a[i][j];
04451                         for (k=0;k<j;k++)
04452                                 sum -= a[i][k]*a[k][j];
04453                         a[i][j]=sum;
04454                         if ( (dum=vv[i]*fabs(sum)) >= big) {
04455                                 big=dum;
04456                                 imax=i;
04457                         }
04458                 }
04459                 if (j != imax) {
04460                         for (k=0;k<n;k++) {
04461                                 dum=a[imax][k];
04462                                 a[imax][k]=a[j][k];
04463                                 a[j][k]=dum;
04464                         }
04465                         *d = -(*d);
04466                         vv[imax]=vv[j];
04467                 }
04468                 indx[j]=imax;
04469                 if (a[j][j] == 0.0) a[j][j]=TINY;
04470                 if (j != n-1) {
04471                         dum=1.0/(a[j][j]);
04472                         for (i=j+1;i<n;i++) a[i][j] *= dum;
04473                 }
04474         }
04475 }

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

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

Referenced by HomePatch::mollyMollify().

04672                                                                                                                                                                                                                                 {
04673   int i,j,k;
04674   BigReal d;
04675   HGMatrixBigReal fjac;
04676   Vector zero(0.0,0.0,0.0);
04677   
04678   HGArrayVector tmpforce;
04679   HGArrayVector tmpforce2;
04680   HGArrayVector y;
04681   HGMatrixVector grhs;
04682   HGMatrixVector glhs;
04683   HGArrayBigReal aux;
04684   HGArrayInt indx;
04685 
04686   for(i=0;i<n;i++) {
04687     tmpforce[i]=imass[i]*force[i];
04688   }
04689 
04690   HGMatrixVector grhs2;
04691   HGArrayVector avgab;
04692 
04693   for ( i = 0; i < m; i++ ) {
04694         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04695   }
04696 
04697   G_q(avgab,glhs,n,m,ial,ibl);
04698   G_q(refab,grhs,n,m,ial,ibl);
04699   // update with the masses
04700   for (j=0; j<n; j++) { // number of atoms
04701         for (i=0; i<m;i++) { // number of constraints
04702           grhs2[i][j] = grhs[i][j]*imass[j];
04703         }
04704   }
04705 
04706   // G_q(qtilde) * M^-1 G_q'(q0) =
04707   // G_q(qtilde) * grhs'
04708   for (i=0;i<m;i++) { // number of constraints
04709         for (j=0;j<m;j++) { // number of constraints
04710           fjac[j][i] = 0; 
04711           for (k=0;k<n;k++) {
04712             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
04713           }
04714         }
04715   }
04716 
04717   // aux=gqij*tmpforce
04718   //  globalGrhs::computeGlobalGrhs(q0,n,water);
04719   //  G_q(refab,grhs,m,ial,ibl);
04720   for(i=0;i<m;i++) {
04721     aux[i]=0.0;
04722     for(j=0;j<n;j++) {
04723       aux[i]+=grhs[i][j]*tmpforce[j];
04724     }
04725   }
04726 
04727   ludcmp(fjac,m,indx,&d);
04728   lubksb(fjac,m,indx,aux);
04729 
04730   for(j=0;j<n;j++) {
04731     y[j] = zero;
04732     for(i=0;i<m;i++) {
04733       y[j] += aux[i]*glhs[i][j];
04734     }
04735   }
04736   for(i=0;i<n;i++) {
04737     y[i]=force[i]-y[i];
04738   }
04739     
04740   // gqq12*y
04741   for(i=0;i<n;i++) {
04742     tmpforce2[i]=imass[i]*y[i];
04743   }
04744 
04745   // here we assume that tmpforce is initialized to zero.
04746   for (i=0;i<n;i++) {
04747     tmpforce[i]=zero;
04748   }
04749   
04750   for (j=0;j<m;j++) {
04751     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
04752     tmpforce[ial[j]] += tmpf;
04753     tmpforce[ibl[j]] -= tmpf;
04754   }
04755   // c-ji the other bug for 2 constraint water was this line (2-4-99)
04756   //  for(i=0;i<m;i++) {
04757   for(i=0;i<n;i++) {
04758     force[i]=tmpforce[i]+y[i];
04759   }
04760 
04761 }


Generated on 19 Sep 2020 for NAMD by  doxygen 1.6.1