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

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

Referenced by HomePatch::mollyAverage().

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

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

Referenced by average(), and mollify().

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

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

Definition at line 4408 of file HomePatch.C.

Referenced by average(), and mollify().

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

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

Definition at line 4431 of file HomePatch.C.

References NAMD_die(), and TINY.

Referenced by average(), and mollify().

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

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

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

Referenced by HomePatch::mollyMollify().

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

Generated on 16 Jun 2022 for NAMD by  doxygen 1.6.1