#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 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:
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 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.
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().
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.
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 }