NAMD
Macros | Typedefs | Functions
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.

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;

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

4493  {
4494  // input: n = length of hydrogen group to be averaged (shaked)
4495  // q[n] = vector array of original positions
4496  // m = number of constraints
4497  // imass[n] = inverse mass for each atom
4498  // length2[m] = square of reference bond length for each constraint
4499  // ial[m] = atom a in each constraint
4500  // ibl[m] = atom b in each constraint
4501  // refab[m] = vector of q_ial(i) - q_ibl(i) for each constraint
4502  // tolf = function error tolerance for Newton's iteration
4503  // ntrial = max number of Newton's iterations
4504  // output: lambda[m] = double array of lagrange multipliers (used by mollify)
4505  // qtilde[n] = vector array of averaged (shaked) positions
4506 
4507  int k,k1,i,j;
4508  BigReal errx,errf,d,tolx;
4509 
4510  HGArrayInt indx;
4511  HGArrayBigReal p;
4512  HGArrayBigReal fvec;
4513  HGMatrixBigReal fjac;
4514  HGArrayVector avgab;
4515  HGMatrixVector grhs;
4516  HGMatrixVector auxrhs;
4517  HGMatrixVector glhs;
4518 
4519  // iout <<iINFO << "average: n="<<n<<" m="<<m<<std::endl<<endi;
4520  tolx=tolf;
4521 
4522  // initialize lambda, globalGrhs
4523 
4524  for (i=0;i<m;i++) {
4525  lambda[i]=0.0;
4526  }
4527 
4528  // define grhs, auxrhs for all iterations
4529  // grhs= g_x(q)
4530  //
4531  G_q(refab,grhs,n,m,ial,ibl);
4532  for (k=1;k<=ntrial;k++) {
4533  // usrfun(qtilde,q0,lambda,fvec,fjac,n,water);
4534  HGArrayBigReal gij;
4535  // this used to be main loop of usrfun
4536  // compute qtilde given q0, lambda, IMASSes
4537  {
4538  BigReal multiplier;
4539  HGArrayVector tmp;
4540  for (i=0;i<m;i++) {
4541  multiplier = lambda[i];
4542  // auxrhs = M^{-1}grhs^{T}
4543  for (j=0;j<n;j++) {
4544  auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
4545  }
4546  }
4547  for (j=0;j<n;j++) {
4548  // tmp[j]=0.0;
4549  for (i=0;i<m;i++) {
4550  tmp[j]+=auxrhs[i][j];
4551  }
4552  }
4553 
4554  for (j=0;j<n;j++) {
4555  qtilde[j].position=q[j]+tmp[j];
4556  }
4557  // delete [] tmp;
4558  }
4559 
4560  for ( i = 0; i < m; i++ ) {
4561  avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
4562  }
4563 
4564  // iout<<iINFO << "Calling Jac" << std::endl<<endi;
4565  // Jac(qtilde, q0, fjac,n,water);
4566  {
4567  // Vector glhs[3*n+3];
4568 
4569  HGMatrixVector grhs2;
4570 
4571  G_q(avgab,glhs,n,m,ial,ibl);
4572 #ifdef DEBUG0
4573  iout<<iINFO << "G_q:" << std::endl<<endi;
4574  for (i=0;i<m;i++) {
4575  iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
4576  }
4577 #endif
4578  // G_q(refab,grhs2,m,ial,ibl);
4579  // update with the masses
4580  for (j=0; j<n; j++) { // number of atoms
4581  for (i=0; i<m;i++) { // number of constraints
4582  grhs2[i][j] = grhs[i][j]*imass[j];
4583  }
4584  }
4585 
4586  // G_q(qtilde) * M^-1 G_q'(q0) =
4587  // G_q(qtilde) * grhs'
4588  for (i=0;i<m;i++) { // number of constraints
4589  for (j=0;j<m;j++) { // number of constraints
4590  fjac[i][j] = 0;
4591  for (k1=0;k1<n;k1++) {
4592  fjac[i][j] += glhs[i][k1]*grhs2[j][k1];
4593  }
4594  }
4595  }
4596 #ifdef DEBUG0
4597  iout<<iINFO << "glhs" <<endi;
4598  for(i=0;i<9;i++) {
4599  iout<<iINFO << glhs[i] << ","<<endi;
4600  }
4601  iout<<iINFO << std::endl<<endi;
4602  for(i=0;i<9;i++) {
4603  iout<<iINFO << grhs2[i] << ","<<endi;
4604  }
4605  iout<<iINFO << std::endl<<endi;
4606 #endif
4607  // delete[] grhs2;
4608  }
4609  // end of Jac calculation
4610 #ifdef DEBUG0
4611  iout<<iINFO << "Jac" << std::endl<<endi;
4612  for (i=0;i<m;i++)
4613  for (j=0;j<m;j++)
4614  iout<<iINFO << fjac[i][j] << " "<<endi;
4615  iout<< std::endl<<endi;
4616 #endif
4617  // calculate constraints in gij for n constraints this being a water
4618  // G(qtilde, gij, n, water);
4619  for (i=0;i<m;i++) {
4620  gij[i]=avgab[i]*avgab[i]-length2[i];
4621  }
4622 #ifdef DEBUG0
4623  iout<<iINFO << "G" << std::endl<<endi;
4624  iout<<iINFO << "( "<<endi;
4625  for(i=0;i<m-1;i++) {
4626  iout<<iINFO << gij[i] << ", "<<endi;
4627  }
4628  iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
4629 #endif
4630  // fill the return vector
4631  for(i=0;i<m;i++) {
4632  fvec[i] = gij[i];
4633  }
4634  // free up the constraints
4635  // delete[] gij;
4636  // continue Newton's iteration
4637  errf=0.0;
4638  for (i=0;i<m;i++) errf += fabs(fvec[i]);
4639 #ifdef DEBUG0
4640  iout<<iINFO << "errf: " << errf << std::endl<<endi;
4641 #endif
4642  if (errf <= tolf) {
4643  break;
4644  }
4645  for (i=0;i<m;i++) p[i] = -fvec[i];
4646  // iout<<iINFO << "Doing dcmp in average " << std::endl<<endi;
4647  ludcmp(fjac,m,indx,&d);
4648  lubksb(fjac,m,indx,p);
4649 
4650  errx=0.0;
4651  for (i=0;i<m;i++) {
4652  errx += fabs(p[i]);
4653  }
4654  for (i=0;i<m;i++)
4655  lambda[i] += p[i];
4656 
4657 #ifdef DEBUG0
4658  iout<<iINFO << "lambda:" << lambda[0]
4659  << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
4660  iout<<iINFO << "errx: " << errx << std::endl<<endi;
4661 #endif
4662  if (errx <= tolx) break;
4663 #ifdef DEBUG0
4664  iout<<iINFO << "Qtilde:" << std::endl<<endi;
4665  iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi;
4666 #endif
4667  }
4668 #ifdef DEBUG
4669  iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
4670 #endif
4671 
4672  return k; //
4673 }
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Position position
Definition: NamdTypes.h:53
BigReal HGMatrixBigReal[MAXHGS][MAXHGS]
Definition: HomePatch.C:65
#define iout
Definition: InfoStream.h:87
BigReal HGArrayBigReal[MAXHGS]
Definition: HomePatch.C:63
zVector HGMatrixVector[MAXHGS][MAXHGS]
Definition: HomePatch.C:66
void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx, HGArrayBigReal &b)
Definition: HomePatch.C:4408
int HGArrayInt[MAXHGS]
Definition: HomePatch.C:62
void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
Definition: HomePatch.C:4431
void G_q(const HGArrayVector &refab, HGMatrixVector &gqij, const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl)
Definition: HomePatch.C:4481
infostream & endi(infostream &s)
Definition: InfoStream.C:38
zVector HGArrayVector[MAXHGS]
Definition: HomePatch.C:64
double BigReal
Definition: common.h:114
static int compDistance ( const void a,
const void b 
)
static

Definition at line 456 of file HomePatch.C.

Referenced by HomePatch::buildSpanningTree().

457 {
458  int d1 = abs(*(int *)a - CkMyPe());
459  int d2 = abs(*(int *)b - CkMyPe());
460  if (d1 < d2)
461  return -1;
462  else if (d1 == d2)
463  return 0;
464  else
465  return 1;
466 }
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().

4482  {
4483  int i;
4484  // step through the rows of the matrix
4485  for(i=0;i<m;i++) {
4486  gqij[i][ial[i]]=2.0*refab[i];
4487  gqij[i][ibl[i]]=-gqij[i][ial[i]];
4488  }
4489 };
void lubksb ( HGMatrixBigReal a,
int  n,
HGArrayInt indx,
HGArrayBigReal b 
)
inline

Definition at line 4408 of file HomePatch.C.

Referenced by average(), and mollify().

4410 {
4411  int i,ii=-1,ip,j;
4412  double sum;
4413 
4414  for (i=0;i<n;i++) {
4415  ip=indx[i];
4416  sum=b[ip];
4417  b[ip]=b[i];
4418  if (ii >= 0)
4419  for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
4420  else if (sum) ii=i;
4421  b[i]=sum;
4422  }
4423  for (i=n-1;i>=0;i--) {
4424  sum=b[i];
4425  for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
4426  b[i]=sum/a[i][i];
4427  }
4428 }
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().

4432 {
4433 
4434  int i,imax,j,k;
4435  double big,dum,sum,temp;
4436  HGArrayBigReal vv;
4437  *d=1.0;
4438  for (i=0;i<n;i++) {
4439  big=0.0;
4440  for (j=0;j<n;j++)
4441  if ((temp=fabs(a[i][j])) > big) big=temp;
4442  if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
4443  vv[i]=1.0/big;
4444  }
4445  for (j=0;j<n;j++) {
4446  for (i=0;i<j;i++) {
4447  sum=a[i][j];
4448  for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
4449  a[i][j]=sum;
4450  }
4451  big=0.0;
4452  for (i=j;i<n;i++) {
4453  sum=a[i][j];
4454  for (k=0;k<j;k++)
4455  sum -= a[i][k]*a[k][j];
4456  a[i][j]=sum;
4457  if ( (dum=vv[i]*fabs(sum)) >= big) {
4458  big=dum;
4459  imax=i;
4460  }
4461  }
4462  if (j != imax) {
4463  for (k=0;k<n;k++) {
4464  dum=a[imax][k];
4465  a[imax][k]=a[j][k];
4466  a[j][k]=dum;
4467  }
4468  *d = -(*d);
4469  vv[imax]=vv[j];
4470  }
4471  indx[j]=imax;
4472  if (a[j][j] == 0.0) a[j][j]=TINY;
4473  if (j != n-1) {
4474  dum=1.0/(a[j][j]);
4475  for (i=j+1;i<n;i++) a[i][j] *= dum;
4476  }
4477  }
4478 }
#define TINY
Definition: HomePatch.C:51
BigReal HGArrayBigReal[MAXHGS]
Definition: HomePatch.C:63
void NAMD_die(const char *err_msg)
Definition: common.C:83
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().

4675  {
4676  int i,j,k;
4677  BigReal d;
4678  HGMatrixBigReal fjac;
4679  Vector zero(0.0,0.0,0.0);
4680 
4681  HGArrayVector tmpforce;
4682  HGArrayVector tmpforce2;
4683  HGArrayVector y;
4684  HGMatrixVector grhs;
4685  HGMatrixVector glhs;
4686  HGArrayBigReal aux;
4687  HGArrayInt indx;
4688 
4689  for(i=0;i<n;i++) {
4690  tmpforce[i]=imass[i]*force[i];
4691  }
4692 
4693  HGMatrixVector grhs2;
4694  HGArrayVector avgab;
4695 
4696  for ( i = 0; i < m; i++ ) {
4697  avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
4698  }
4699 
4700  G_q(avgab,glhs,n,m,ial,ibl);
4701  G_q(refab,grhs,n,m,ial,ibl);
4702  // update with the masses
4703  for (j=0; j<n; j++) { // number of atoms
4704  for (i=0; i<m;i++) { // number of constraints
4705  grhs2[i][j] = grhs[i][j]*imass[j];
4706  }
4707  }
4708 
4709  // G_q(qtilde) * M^-1 G_q'(q0) =
4710  // G_q(qtilde) * grhs'
4711  for (i=0;i<m;i++) { // number of constraints
4712  for (j=0;j<m;j++) { // number of constraints
4713  fjac[j][i] = 0;
4714  for (k=0;k<n;k++) {
4715  fjac[j][i] += glhs[i][k]*grhs2[j][k];
4716  }
4717  }
4718  }
4719 
4720  // aux=gqij*tmpforce
4721  // globalGrhs::computeGlobalGrhs(q0,n,water);
4722  // G_q(refab,grhs,m,ial,ibl);
4723  for(i=0;i<m;i++) {
4724  aux[i]=0.0;
4725  for(j=0;j<n;j++) {
4726  aux[i]+=grhs[i][j]*tmpforce[j];
4727  }
4728  }
4729 
4730  ludcmp(fjac,m,indx,&d);
4731  lubksb(fjac,m,indx,aux);
4732 
4733  for(j=0;j<n;j++) {
4734  y[j] = zero;
4735  for(i=0;i<m;i++) {
4736  y[j] += aux[i]*glhs[i][j];
4737  }
4738  }
4739  for(i=0;i<n;i++) {
4740  y[i]=force[i]-y[i];
4741  }
4742 
4743  // gqq12*y
4744  for(i=0;i<n;i++) {
4745  tmpforce2[i]=imass[i]*y[i];
4746  }
4747 
4748  // here we assume that tmpforce is initialized to zero.
4749  for (i=0;i<n;i++) {
4750  tmpforce[i]=zero;
4751  }
4752 
4753  for (j=0;j<m;j++) {
4754  Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
4755  tmpforce[ial[j]] += tmpf;
4756  tmpforce[ibl[j]] -= tmpf;
4757  }
4758  // c-ji the other bug for 2 constraint water was this line (2-4-99)
4759  // for(i=0;i<m;i++) {
4760  for(i=0;i<n;i++) {
4761  force[i]=tmpforce[i]+y[i];
4762  }
4763 
4764 }
Definition: Vector.h:64
Position position
Definition: NamdTypes.h:53
BigReal HGMatrixBigReal[MAXHGS][MAXHGS]
Definition: HomePatch.C:65
BigReal HGArrayBigReal[MAXHGS]
Definition: HomePatch.C:63
zVector HGMatrixVector[MAXHGS][MAXHGS]
Definition: HomePatch.C:66
void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx, HGArrayBigReal &b)
Definition: HomePatch.C:4408
int HGArrayInt[MAXHGS]
Definition: HomePatch.C:62
void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
Definition: HomePatch.C:4431
void G_q(const HGArrayVector &refab, HGMatrixVector &gqij, const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl)
Definition: HomePatch.C:4481
gridSize y
zVector HGArrayVector[MAXHGS]
Definition: HomePatch.C:64
double BigReal
Definition: common.h:114