HomePatch.C File Reference

#include "NamdTypes.h"
#include "common.h"
#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 "MigrationCUDAKernel.h"
#include "ComputeQM.h"
#include "ComputeQMMgr.decl.h"
#include "NamdEventsProfiling.h"
#include "Debug.h"
#include <vector>
#include <algorithm>

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)
void	MSHAKE_CUDA (int *, const int size, const RattleParam *rattleParam, BigReal *refx, BigReal *refy, BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, const BigReal tol2, const int maxiter, bool &done, bool &consFailure)
static int	compDistance (const void *a, const void *b)
void	PatchDataSOA_initialize (PatchDataSOA *p)
size_t	PatchDataSOA_set_size (PatchDataSOA *p, int natoms, int padding)
void	PatchDataSOA_set_buffer (PatchDataSOA *p, void *mybuffer)
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

FIX_FOR_WATER

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

MASS_EPSILON

#define MASS_EPSILON   (1.0e-35)

Definition at line 80 of file HomePatch.C.

MAXHGS

#define MAXHGS   10

Definition at line 55 of file HomePatch.C.

MIN_DEBUG_LEVEL

#define MIN_DEBUG_LEVEL   2

Definition at line 56 of file HomePatch.C.

TINY

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

Referenced by ludcmp().

Typedef Documentation

HGArrayBigReal

typedef BigReal HGArrayBigReal[MAXHGS]

Definition at line 66 of file HomePatch.C.

HGArrayInt

typedef int HGArrayInt[MAXHGS]

Definition at line 65 of file HomePatch.C.

HGArrayVector

typedef zVector HGArrayVector[MAXHGS]

Definition at line 67 of file HomePatch.C.

HGMatrixBigReal

typedef BigReal HGMatrixBigReal[MAXHGS][MAXHGS]

Definition at line 68 of file HomePatch.C.

HGMatrixVector

typedef zVector HGMatrixVector[MAXHGS][MAXHGS]

Definition at line 69 of file HomePatch.C.

Function Documentation

average()

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

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

Referenced by PressureProfileReduction::getData(), HomePatch::mollyAverage(), PressureProfileReduction::PressureProfileReduction(), and PressureProfileReduction::~PressureProfileReduction().

                                                                                                                                                                                                                                                                         {
  //  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; // 
}

compDistance()

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

static

Definition at line 497 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;
}

G_q()

void G_q ( const HGArrayVector &  refab, HGMatrixVector &  gqij, const int  n, const int  m, const HGArrayInt &  ial, const HGArrayInt &  ibl  )

inline

Definition at line 6458 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]];
  }
};

lubksb()

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

inline

Definition at line 6385 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];
        }
}

ludcmp()

void ludcmp ( HGMatrixBigReal &  a, int  n, HGArrayInt &  indx, BigReal *  d  )

inline

Definition at line 6408 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;
                }
        }
}

mollify()

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

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

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];
  }

}

MSHAKE_CUDA()

void MSHAKE_CUDA	(	int *	,
		const int	size,
		const RattleParam *	rattleParam,
		BigReal *	refx,
		BigReal *	refy,
		BigReal *	refz,
		BigReal *	posx,
		BigReal *	posy,
		BigReal *	posz,
		const BigReal	tol2,
		const int	maxiter,
		bool &	done,
		bool &	consFailure
	)

PatchDataSOA_initialize()

void PatchDataSOA_initialize

(

PatchDataSOA *

p

)

Initialize fundamental parts of PatchDataSOA.

Definition at line 2366 of file HomePatch.C.

References PatchDataSOA::buffer, PatchDataSOA::maxAtoms, PatchDataSOA::numAtoms, and PatchDataSOA::numBytes.

                                              {
  p->buffer = NULL;
  p->numBytes = 0;
  p->numAtoms = 0;
  p->maxAtoms = 0;
}

PatchDataSOA_set_buffer()

void PatchDataSOA_set_buffer	(	PatchDataSOA *	p,
		void *	mybuffer
	)

Caller provides the buffer space to be used. Set all pointers relative to this provided address, pointing into the buffer space.

Definition at line 2397 of file HomePatch.C.

References PatchDataSOA::atomFixed, PatchDataSOA::buffer, PatchDataSOA::charge, PatchDataSOA::exclId, PatchDataSOA::f_global_x, PatchDataSOA::f_global_y, PatchDataSOA::f_global_z, PatchDataSOA::f_nbond_x, PatchDataSOA::f_nbond_y, PatchDataSOA::f_nbond_z, PatchDataSOA::f_normal_x, PatchDataSOA::f_normal_y, PatchDataSOA::f_normal_z, PatchDataSOA::f_saved_nbond_x, PatchDataSOA::f_saved_nbond_y, PatchDataSOA::f_saved_nbond_z, PatchDataSOA::f_saved_slow_x, PatchDataSOA::f_saved_slow_y, PatchDataSOA::f_saved_slow_z, PatchDataSOA::f_slow_x, PatchDataSOA::f_slow_y, PatchDataSOA::f_slow_z, PatchDataSOA::fixedPosition_x, PatchDataSOA::fixedPosition_y, PatchDataSOA::fixedPosition_z, PatchDataSOA::gaussrand_x, PatchDataSOA::gaussrand_y, PatchDataSOA::gaussrand_z, PatchDataSOA::groupFixed, PatchDataSOA::hydrogenGroupSize, PatchDataSOA::id, PatchDataSOA::isWater, PatchDataSOA::langevinParam, PatchDataSOA::langScalRandBBK2, PatchDataSOA::langScalVelBBK2, PatchDataSOA::mass, PatchDataSOA::maxAtoms, PatchDataSOA::migrationGroupSize, PatchDataSOA::nonbondedGroupSize, Node::Object(), PatchDataSOA::partition, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, PatchDataSOA::posNew_x, PatchDataSOA::posNew_y, PatchDataSOA::posNew_z, PatchDataSOA::recipMass, PatchDataSOA::rigidBondLength, PatchDataSOA::sigId, Node::simParameters, simParams, PatchDataSOA::sortOrder, PatchDataSOA::status, PatchDataSOA::transform_i, PatchDataSOA::transform_j, PatchDataSOA::transform_k, PatchDataSOA::unsortOrder, PatchDataSOA::vdwType, PatchDataSOA::vel_x, PatchDataSOA::vel_y, PatchDataSOA::vel_z, PatchDataSOA::velNew_x, PatchDataSOA::velNew_y, and PatchDataSOA::velNew_z.

      {
  const SimParameters *simParams = Node::Object()->simParameters;
  p->buffer = (unsigned char *) mybuffer;
  unsigned char *t = p->buffer;
  p->pos_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->pos_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->pos_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->charge = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->vdwType = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->partition = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->nonbondedGroupSize = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->hydrogenGroupSize = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->isWater = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->sortOrder = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->unsortOrder = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->id = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->exclId = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->sigId = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->atomFixed = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->groupFixed = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->vel_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->vel_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->vel_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->fixedPosition_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->fixedPosition_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->fixedPosition_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->recipMass = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->mass = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->langevinParam = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->status = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->transform_i = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->transform_j = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->transform_k = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->migrationGroupSize = (int *) t;
  t += p->maxAtoms * sizeof(int);
  p->rigidBondLength = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->langScalVelBBK2 = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->langScalRandBBK2 = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->gaussrand_x = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->gaussrand_y = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->gaussrand_z = (float *) t;
  t += p->maxAtoms * sizeof(float);
  p->f_normal_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_normal_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_normal_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_nbond_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_nbond_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_nbond_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_slow_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_slow_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_slow_z = (double *) t;
  t += p->maxAtoms * sizeof(float);
  if (simParams->colvarsOn || simParams->tclForcesOn || (simParams->IMDon && ! (simParams->IMDignore || simParams->IMDignoreForces))){
    p->f_global_x = (double *) t;
    t += p->maxAtoms * sizeof(double);
    p->f_global_y = (double *) t;
    t += p->maxAtoms * sizeof(double);
    p->f_global_z = (double *) t;
    t += p->maxAtoms * sizeof(double);
  }
  p->f_saved_nbond_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_saved_nbond_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_saved_nbond_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_saved_slow_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_saved_slow_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->f_saved_slow_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->velNew_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->velNew_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->velNew_z = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->posNew_x = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->posNew_y = (double *) t;
  t += p->maxAtoms * sizeof(double);
  p->posNew_z = (double *) t;
}

PatchDataSOA_set_size()

size_t PatchDataSOA_set_size	(	PatchDataSOA *	p,
		int	natoms,
		int	pad = MAXFACTOR
	)

Set sizes based on number of atoms and padding factor. Sets internally numAtoms, maxAtoms, and numBytes. Returns numBytes.

Definition at line 2373 of file HomePatch.C.

References PatchDataSOA::maxAtoms, PatchDataSOA::numAtoms, PatchDataSOA::numBytes, Node::Object(), Node::simParameters, and simParams.

      {
  if (natoms > p->maxAtoms) {
    // maxAtoms extends numAtoms to next multiple of padding
    p->maxAtoms = ((natoms + padding - 1) / padding) * padding;
    // set numBytes to allow each array to have maxAtoms space
    // count up total number of double, float, and int arrays
    // XXX TODO : these magic numbers are bad for maintability
    const SimParameters *simParams = Node::Object()->simParameters;
    // XXX TODO: For Colvars, tclForces, IMD or any GlobalMaster
    // client that requires the global forces, we need extra 3 buffers
    // to save the x, y and z components of the forces, and therefore
    // the number of double fields should be increased by 3 (from 31 to 34).
    int numdoubles = (simParams->colvarsOn || simParams->tclForcesOn || simParams->IMDon) ? 34: 31;
    p->numBytes = p->maxAtoms *
      (numdoubles * sizeof(double) + 9*sizeof(float) + 17*sizeof(int));
  }
  p->numAtoms = natoms;
  return p->numBytes;
}