#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"

Macros
#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)

Macro Definition 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;

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().

                                                                                                                                                                                                                                                                          {
   //  input:  n = length of hydrogen group to be averaged (shaked)
   //          q[n] = vector array of original positions
   //          m = number of constraints
   //          imass[n] = inverse mass for each atom
   //          length2[m] = square of reference bond length for each constraint
   //          ial[m] = atom a in each constraint 
   //          ibl[m] = atom b in each constraint 
   //          refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
   //          tolf = function error tolerance for Newton's iteration
   //          ntrial = max number of Newton's iterations
   //  output: lambda[m] = double array of lagrange multipliers (used by mollify)
   //          qtilde[n] = vector array of averaged (shaked) positions
 
   int k,k1,i,j;
   BigReal errx,errf,d,tolx;
 
   HGArrayInt indx;
   HGArrayBigReal p;
   HGArrayBigReal fvec;
   HGMatrixBigReal fjac;
   HGArrayVector avgab;
   HGMatrixVector grhs;
   HGMatrixVector auxrhs;
   HGMatrixVector glhs;
 
   //  iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
   tolx=tolf; 
   
   // initialize lambda, globalGrhs
 
   for (i=0;i<m;i++) {
     lambda[i]=0.0;
   }
 
   // define grhs, auxrhs for all iterations
   // grhs= g_x(q)
   //
   G_q(refab,grhs,n,m,ial,ibl);
   for (k=1;k<=ntrial;k++) {
     //    usrfun(qtilde,q0,lambda,fvec,fjac,n,water); 
     HGArrayBigReal gij;
     // this used to be main loop of usrfun
     // compute qtilde given q0, lambda, IMASSes
     {
       BigReal multiplier;
       HGArrayVector tmp;
       for (i=0;i<m;i++) {
         multiplier = lambda[i];
         // auxrhs = M^{-1}grhs^{T}
         for (j=0;j<n;j++) {
           auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
         }
       }
       for (j=0;j<n;j++) {
         //      tmp[j]=0.0;      
         for (i=0;i<m;i++) {
           tmp[j]+=auxrhs[i][j];
         }
       }
  
       for (j=0;j<n;j++) {
         qtilde[j].position=q[j]+tmp[j];
       }
       //      delete [] tmp;
     }
   
     for ( i = 0; i < m; i++ ) {
       avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
     }
 
     //  iout<<iINFO << "Calling Jac" << std::endl<<endi;
     //  Jac(qtilde, q0, fjac,n,water);
     {
       //  Vector glhs[3*n+3];
 
       HGMatrixVector grhs2;
 
       G_q(avgab,glhs,n,m,ial,ibl);
 #ifdef DEBUG0
       iout<<iINFO << "G_q:" << std::endl<<endi;
       for (i=0;i<m;i++) {
         iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
       }
 #endif
       //      G_q(refab,grhs2,m,ial,ibl);
       // update with the masses
       for (j=0; j<n; j++) { // number of atoms
         for (i=0; i<m;i++) { // number of constraints
           grhs2[i][j] = grhs[i][j]*imass[j];
         }
       }
 
       // G_q(qtilde) * M^-1 G_q'(q0) =
       // G_q(qtilde) * grhs'
       for (i=0;i<m;i++) { // number of constraints
         for (j=0;j<m;j++) { // number of constraints
           fjac[i][j] = 0; 
           for (k1=0;k1<n;k1++) {
             fjac[i][j] += glhs[i][k1]*grhs2[j][k1]; 
           }
         }
       }
 #ifdef DEBUG0  
       iout<<iINFO << "glhs" <<endi;
       for(i=0;i<9;i++) {
         iout<<iINFO << glhs[i] << ","<<endi;
       }
       iout<<iINFO << std::endl<<endi;
       for(i=0;i<9;i++) {
         iout<<iINFO << grhs2[i] << ","<<endi;
       }
       iout<<iINFO << std::endl<<endi;
 #endif
       //      delete[] grhs2;
     }
     // end of Jac calculation
 #ifdef DEBUG0
     iout<<iINFO << "Jac" << std::endl<<endi;
     for (i=0;i<m;i++) 
       for (j=0;j<m;j++)
         iout<<iINFO << fjac[i][j] << " "<<endi;
     iout<< std::endl<<endi;
 #endif
     // calculate constraints in gij for n constraints this being a water
     //  G(qtilde, gij, n, water);
     for (i=0;i<m;i++) {
       gij[i]=avgab[i]*avgab[i]-length2[i];
     }
 #ifdef DEBUG0
     iout<<iINFO << "G" << std::endl<<endi;
     iout<<iINFO << "( "<<endi;
     for(i=0;i<m-1;i++) {
       iout<<iINFO << gij[i] << ", "<<endi;
     }
     iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
 #endif
     // fill the return vector
     for(i=0;i<m;i++) {
       fvec[i] = gij[i];
     }
     // free up the constraints
     //    delete[] gij;
     // continue Newton's iteration    
     errf=0.0;
     for (i=0;i<m;i++) errf += fabs(fvec[i]);
 #ifdef DEBUG0
     iout<<iINFO << "errf: " << errf << std::endl<<endi;
 #endif
     if (errf <= tolf) {
       break;
     }
     for (i=0;i<m;i++) p[i] = -fvec[i];
     //    iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
     ludcmp(fjac,m,indx,&d);
     lubksb(fjac,m,indx,p);
 
     errx=0.0;
     for (i=0;i<m;i++) {
       errx += fabs(p[i]);
     }
     for (i=0;i<m;i++)  
       lambda[i] += p[i];
 
 #ifdef DEBUG0
     iout<<iINFO << "lambda:" << lambda[0] 
          << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
     iout<<iINFO << "errx: " << errx << std::endl<<endi;
 #endif
     if (errx <= tolx) break;
 #ifdef DEBUG0
     iout<<iINFO << "Qtilde:" << std::endl<<endi;
     iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi; 
 #endif
   }
 #ifdef DEBUG
   iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
 #endif
 
   return k; // 
 }

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

static

Definition at line 456 of file HomePatch.C.

Referenced by HomePatch::buildSpanningTree().

 {
   int d1 = abs(*(int *)a - CkMyPe());
   int d2 = abs(*(int *)b - CkMyPe());
   if (d1 < d2) 
     return -1;
   else if (d1 == d2) 
     return 0;
   else 
     return 1;
 }

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().

                                                                              {
   int i; 
   // step through the rows of the matrix
   for(i=0;i<m;i++) {
     gqij[i][ial[i]]=2.0*refab[i];
     gqij[i][ibl[i]]=-gqij[i][ial[i]];
   }
 };

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

inline

Definition at line 4408 of file HomePatch.C.

Referenced by average(), and mollify().

 {
         int i,ii=-1,ip,j;
         double sum;
 
         for (i=0;i<n;i++) {
                 ip=indx[i];
                 sum=b[ip];
                 b[ip]=b[i];
                 if (ii >= 0)
                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
                 else if (sum) ii=i;
                 b[i]=sum;
         }
         for (i=n-1;i>=0;i--) {
                 sum=b[i];
                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
                 b[i]=sum/a[i][i];
         }
 }

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().

 {
 
         int i,imax,j,k;
         double big,dum,sum,temp;
         HGArrayBigReal vv;
         *d=1.0;
         for (i=0;i<n;i++) {
                 big=0.0;
                 for (j=0;j<n;j++)
                         if ((temp=fabs(a[i][j])) > big) big=temp;
                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
                 vv[i]=1.0/big;
         }
         for (j=0;j<n;j++) {
                 for (i=0;i<j;i++) {
                         sum=a[i][j];
                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
                         a[i][j]=sum;
                 }
                 big=0.0;
                 for (i=j;i<n;i++) {
                         sum=a[i][j];
                         for (k=0;k<j;k++)
                                 sum -= a[i][k]*a[k][j];
                         a[i][j]=sum;
                         if ( (dum=vv[i]*fabs(sum)) >= big) {
                                 big=dum;
                                 imax=i;
                         }
                 }
                 if (j != imax) {
                         for (k=0;k<n;k++) {
                                 dum=a[imax][k];
                                 a[imax][k]=a[j][k];
                                 a[j][k]=dum;
                         }
                         *d = -(*d);
                         vv[imax]=vv[j];
                 }
                 indx[j]=imax;
                 if (a[j][j] == 0.0) a[j][j]=TINY;
                 if (j != n-1) {
                         dum=1.0/(a[j][j]);
                         for (i=j+1;i<n;i++) a[i][j] *= dum;
                 }
         }
 }

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(), CompAtom::position, and y.

Referenced by HomePatch::mollyMollify().

                                                                                                                                                                                                                                 {
   int i,j,k;
   BigReal d;
   HGMatrixBigReal fjac;
   Vector zero(0.0,0.0,0.0);
   
   HGArrayVector tmpforce;
   HGArrayVector tmpforce2;
   HGArrayVector y;
   HGMatrixVector grhs;
   HGMatrixVector glhs;
   HGArrayBigReal aux;
   HGArrayInt indx;
 
   for(i=0;i<n;i++) {
     tmpforce[i]=imass[i]*force[i];
   }
 
   HGMatrixVector grhs2;
   HGArrayVector avgab;
 
   for ( i = 0; i < m; i++ ) {
         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
   }
 
   G_q(avgab,glhs,n,m,ial,ibl);
   G_q(refab,grhs,n,m,ial,ibl);
   // update with the masses
   for (j=0; j<n; j++) { // number of atoms
         for (i=0; i<m;i++) { // number of constraints
           grhs2[i][j] = grhs[i][j]*imass[j];
         }
   }
 
   // G_q(qtilde) * M^-1 G_q'(q0) =
   // G_q(qtilde) * grhs'
   for (i=0;i<m;i++) { // number of constraints
         for (j=0;j<m;j++) { // number of constraints
           fjac[j][i] = 0; 
           for (k=0;k<n;k++) {
             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
           }
         }
   }
 
   // aux=gqij*tmpforce
   //  globalGrhs::computeGlobalGrhs(q0,n,water);
   //  G_q(refab,grhs,m,ial,ibl);
   for(i=0;i<m;i++) {
     aux[i]=0.0;
     for(j=0;j<n;j++) {
       aux[i]+=grhs[i][j]*tmpforce[j];
     }
   }
 
   ludcmp(fjac,m,indx,&d);
   lubksb(fjac,m,indx,aux);
 
   for(j=0;j<n;j++) {
     y[j] = zero;
     for(i=0;i<m;i++) {
       y[j] += aux[i]*glhs[i][j];
     }
   }
   for(i=0;i<n;i++) {
     y[i]=force[i]-y[i];
   }
     
   // gqq12*y
   for(i=0;i<n;i++) {
     tmpforce2[i]=imass[i]*y[i];
   }
 
   // here we assume that tmpforce is initialized to zero.
   for (i=0;i<n;i++) {
     tmpforce[i]=zero;
   }
   
   for (j=0;j<m;j++) {
     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
     tmpforce[ial[j]] += tmpf;
     tmpforce[ibl[j]] -= tmpf;
   }
   // c-ji the other bug for 2 constraint water was this line (2-4-99)
   //  for(i=0;i<m;i++) {
   for(i=0;i<n;i++) {
     force[i]=tmpforce[i]+y[i];
   }
 
 }

Macros

Typedefs

Functions

Macro Definition Documentation

Typedef Documentation

Function Documentation