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

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

Referenced by average(), and mollify().

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

Definition at line 4405 of file HomePatch.C.

Referenced by average(), and mollify().

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

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

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

Referenced by HomePatch::mollyMollify().

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