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

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

Referenced by HomePatch::mollyAverage().

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

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

Referenced by average(), and mollify().

04447                                                                              {
04448   int i; 
04449   // step through the rows of the matrix
04450   for(i=0;i<m;i++) {
04451     gqij[i][ial[i]]=2.0*refab[i];
04452     gqij[i][ibl[i]]=-gqij[i][ial[i]];
04453   }
04454 };

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

Definition at line 4373 of file HomePatch.C.

Referenced by average(), and mollify().

04375 {
04376         int i,ii=-1,ip,j;
04377         double sum;
04378 
04379         for (i=0;i<n;i++) {
04380                 ip=indx[i];
04381                 sum=b[ip];
04382                 b[ip]=b[i];
04383                 if (ii >= 0)
04384                         for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
04385                 else if (sum) ii=i;
04386                 b[i]=sum;
04387         }
04388         for (i=n-1;i>=0;i--) {
04389                 sum=b[i];
04390                 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
04391                 b[i]=sum/a[i][i];
04392         }
04393 }

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

Definition at line 4396 of file HomePatch.C.

References NAMD_die(), and TINY.

Referenced by average(), and mollify().

04397 {
04398 
04399         int i,imax,j,k;
04400         double big,dum,sum,temp;
04401         HGArrayBigReal vv;
04402         *d=1.0;
04403         for (i=0;i<n;i++) {
04404                 big=0.0;
04405                 for (j=0;j<n;j++)
04406                         if ((temp=fabs(a[i][j])) > big) big=temp;
04407                 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
04408                 vv[i]=1.0/big;
04409         }
04410         for (j=0;j<n;j++) {
04411                 for (i=0;i<j;i++) {
04412                         sum=a[i][j];
04413                         for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
04414                         a[i][j]=sum;
04415                 }
04416                 big=0.0;
04417                 for (i=j;i<n;i++) {
04418                         sum=a[i][j];
04419                         for (k=0;k<j;k++)
04420                                 sum -= a[i][k]*a[k][j];
04421                         a[i][j]=sum;
04422                         if ( (dum=vv[i]*fabs(sum)) >= big) {
04423                                 big=dum;
04424                                 imax=i;
04425                         }
04426                 }
04427                 if (j != imax) {
04428                         for (k=0;k<n;k++) {
04429                                 dum=a[imax][k];
04430                                 a[imax][k]=a[j][k];
04431                                 a[j][k]=dum;
04432                         }
04433                         *d = -(*d);
04434                         vv[imax]=vv[j];
04435                 }
04436                 indx[j]=imax;
04437                 if (a[j][j] == 0.0) a[j][j]=TINY;
04438                 if (j != n-1) {
04439                         dum=1.0/(a[j][j]);
04440                         for (i=j+1;i<n;i++) a[i][j] *= dum;
04441                 }
04442         }
04443 }

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

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

Referenced by HomePatch::mollyMollify().

04640                                                                                                                                                                                                                                 {
04641   int i,j,k;
04642   BigReal d;
04643   HGMatrixBigReal fjac;
04644   Vector zero(0.0,0.0,0.0);
04645   
04646   HGArrayVector tmpforce;
04647   HGArrayVector tmpforce2;
04648   HGArrayVector y;
04649   HGMatrixVector grhs;
04650   HGMatrixVector glhs;
04651   HGArrayBigReal aux;
04652   HGArrayInt indx;
04653 
04654   for(i=0;i<n;i++) {
04655     tmpforce[i]=imass[i]*force[i];
04656   }
04657 
04658   HGMatrixVector grhs2;
04659   HGArrayVector avgab;
04660 
04661   for ( i = 0; i < m; i++ ) {
04662         avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
04663   }
04664 
04665   G_q(avgab,glhs,n,m,ial,ibl);
04666   G_q(refab,grhs,n,m,ial,ibl);
04667   // update with the masses
04668   for (j=0; j<n; j++) { // number of atoms
04669         for (i=0; i<m;i++) { // number of constraints
04670           grhs2[i][j] = grhs[i][j]*imass[j];
04671         }
04672   }
04673 
04674   // G_q(qtilde) * M^-1 G_q'(q0) =
04675   // G_q(qtilde) * grhs'
04676   for (i=0;i<m;i++) { // number of constraints
04677         for (j=0;j<m;j++) { // number of constraints
04678           fjac[j][i] = 0; 
04679           for (k=0;k<n;k++) {
04680             fjac[j][i] += glhs[i][k]*grhs2[j][k]; 
04681           }
04682         }
04683   }
04684 
04685   // aux=gqij*tmpforce
04686   //  globalGrhs::computeGlobalGrhs(q0,n,water);
04687   //  G_q(refab,grhs,m,ial,ibl);
04688   for(i=0;i<m;i++) {
04689     aux[i]=0.0;
04690     for(j=0;j<n;j++) {
04691       aux[i]+=grhs[i][j]*tmpforce[j];
04692     }
04693   }
04694 
04695   ludcmp(fjac,m,indx,&d);
04696   lubksb(fjac,m,indx,aux);
04697 
04698   for(j=0;j<n;j++) {
04699     y[j] = zero;
04700     for(i=0;i<m;i++) {
04701       y[j] += aux[i]*glhs[i][j];
04702     }
04703   }
04704   for(i=0;i<n;i++) {
04705     y[i]=force[i]-y[i];
04706   }
04707     
04708   // gqq12*y
04709   for(i=0;i<n;i++) {
04710     tmpforce2[i]=imass[i]*y[i];
04711   }
04712 
04713   // here we assume that tmpforce is initialized to zero.
04714   for (i=0;i<n;i++) {
04715     tmpforce[i]=zero;
04716   }
04717   
04718   for (j=0;j<m;j++) {
04719     Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
04720     tmpforce[ial[j]] += tmpf;
04721     tmpforce[ibl[j]] -= tmpf;
04722   }
04723   // c-ji the other bug for 2 constraint water was this line (2-4-99)
04724   //  for(i=0;i<m;i++) {
04725   for(i=0;i<n;i++) {
04726     force[i]=tmpforce[i]+y[i];
04727   }
04728 
04729 }


Generated on 29 May 2020 for NAMD by  doxygen 1.6.1