GridforceFullMainGrid Class Reference

#include <GridForceGrid.h>

Inheritance diagram for GridforceFullMainGrid:

Public Member Functions
	GridforceFullMainGrid (int gridnum)
virtual	~GridforceFullMainGrid ()
void	initialize (char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams, int border)
void	initialize (char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
void	reinitialize (SimParameters *simParams, MGridforceParams *mgridParams)
Position	get_center (void) const
Position	get_origin (void) const
Tensor	get_e (void) const
Tensor	get_inv (void) const
Vector	get_scale (void) const
Bool	get_checksize (void) const
int	get_k0 (void) const
int	get_k1 (void) const
int	get_k2 (void) const
int	get_border (void) const
int	compute_VdV (Position pos, float &V, Vector &dV) const
int	get_total_grids (void) const
void	set_scale (Vector s)
Public Member Functions inherited from GridforceGrid
virtual	~GridforceGrid ()
Position	wrap_position (const Position &pos, const Lattice &lattice)
bool	fits_lattice (const Lattice &lattice)
int	compute_VdV (Position pos, float &V, Vector &dV) const
GridforceGridType	get_grid_type (void)
Public Member Functions inherited from GridforceFullBaseGrid
	GridforceFullBaseGrid (void)
virtual	~GridforceFullBaseGrid ()
Position	get_center (void) const
Position	get_origin (void) const
Tensor	get_e (void) const
Tensor	get_inv (void) const
Vector	get_scale (void) const
Bool	get_checksize (void) const
float	get_grid (int i0, int i1, int i2) const
double	get_grid_d (int i0, int i1, int i2) const
void	set_grid (int i0, int i1, int i2, float V)
void	set_scale (Vector s)
int	compute_VdV (Position pos, float &V, Vector &dV) const
int	get_k0 (void) const
int	get_k1 (void) const
int	get_k2 (void) const

Protected Member Functions
void	pack (MOStream *msg) const
void	unpack (MIStream *msg)
long int	get_all_gridvals (float **all_gridvals) const
void	set_all_gridvals (float *all_gridvals, long int sz)
void	compute_b (float *b, int *inds, Vector gapscale) const
void	buildSubgridsFlat (void)
Protected Member Functions inherited from GridforceGrid
Position	get_corner (int idx)
	GridforceGrid ()
Protected Member Functions inherited from GridforceFullBaseGrid
void	readHeader (SimParameters *simParams, MGridforceParams *mgridParams)
long int	grid_index (int i0, int i1, int i2) const
int	get_inds (Position pos, int *inds, Vector &dg, Vector &gapscale) const
void	compute_a (float *a, float *b) const
float	compute_V (float *a, float *x, float *y, float *z) const
Vector	compute_dV (float *a, float *x, float *y, float *z) const
Vector	compute_d2V (float *a, float *x, float *y, float *z) const
float	compute_d3V (float *a, float *x, float *y, float *z) const
void	readSubgridHierarchy (FILE *poten, int &totalGrids)

Protected Attributes
char	filename [NAMD_FILENAME_BUFFER_SIZE]
int	totalGrids
GridforceFullSubGrid **	subgrids_flat
int	border
Protected Attributes inherited from GridforceGrid
GridforceGridType	type
int	mygridnum
Protected Attributes inherited from GridforceFullBaseGrid
FILE *	poten_fp
float *	grid
GridforceFullSubGrid **	subgrids
int	numSubgrids
int	generation
int	k [3]
int	k_nopad [3]
long int	size
long int	size_nopad
long int	dk [3]
long int	dk_nopad [3]
float	factor
Position	origin
Position	center
Tensor	e
Tensor	inv
double	p_sum [3]
double	n_sum [3]
double	pad_p [3]
double	pad_n [3]
Bool	cont [3]
float	offset [3]
float	gap [3]
float	gapinv [3]
Vector	scale
Bool	checksize

Static Protected Attributes
static const int	default_border = 1

Friends
class	GridforceFullBaseGrid
class	GridforceFullSubGrid

Additional Inherited Members
Public Types inherited from GridforceGrid
enum	GridforceGridType { GridforceGridTypeUndefined = 0, GridforceGridTypeFull, GridforceGridTypeLite }
Static Public Member Functions inherited from GridforceGrid
static GridforceGrid *	new_grid (int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
static void	pack_grid (GridforceGrid *grid, MOStream *msg)
static GridforceGrid *	unpack_grid (int gridnum, MIStream *msg)

Detailed Description

Definition at line 191 of file GridForceGrid.h.

Constructor & Destructor Documentation

GridforceFullMainGrid::GridforceFullMainGrid ( int  gridnum )

explicit

Definition at line 354 of file GridForceGrid.C.

References GridforceFullBaseGrid::generation, GridforceGrid::GridforceGridTypeFull, GridforceGrid::mygridnum, subgrids_flat, and GridforceGrid::type.

{
    mygridnum = gridnum;
    generation = 0;
    subgrids_flat = NULL;
    type = GridforceGridTypeFull;
}

GridforceFullMainGrid::~GridforceFullMainGrid ( )

virtual

Definition at line 363 of file GridForceGrid.C.

References subgrids_flat.

{
    delete[] subgrids_flat;
}

Member Function Documentation

void GridforceFullMainGrid::buildSubgridsFlat ( void  )

protected

Definition at line 407 of file GridForceGrid.C.

References GridforceFullSubGrid::addToSubgridsFlat(), DebugM, endi(), GridforceFullBaseGrid::numSubgrids, GridforceFullBaseGrid::subgrids, subgrids_flat, and totalGrids.

Referenced by initialize(), and unpack().

{
    DebugM(4, "buildSubgridsFlat() called, totalGrids-1 = " << totalGrids-1 << "\n" << endi);
    delete[] subgrids_flat;
    subgrids_flat = new GridforceFullSubGrid *[totalGrids-1];
    for (int i = 0; i < numSubgrids; i++) {
        DebugM(3, "adding to subgridsFlat\n" << endi);
        subgrids[i]->addToSubgridsFlat();
        DebugM(3, "success!\n" << endi);
    }
    for (int i = 0; i < totalGrids-1; i++) {
        DebugM(4, "subgrids_flat[" << i << "]->numSubgrids = " << subgrids_flat[i]->numSubgrids << "\n" << endi);
    }
    for (int i = 0; i < numSubgrids; i++) {
        DebugM(4, "subgrids[" << i << "]->numSubgrids = " << subgrids[i]->numSubgrids << "\n" << endi);
    }
}

void GridforceFullMainGrid::compute_b ( float *  b, int *  inds, Vector  gapscale  ) const

protectedvirtual

Implements GridforceFullBaseGrid.

Definition at line 782 of file GridForceGrid.C.

References GridforceFullBaseGrid::cont, DebugM, endi(), FALSE, GridforceFullBaseGrid::gap, GridforceFullBaseGrid::get_grid(), GridforceFullBaseGrid::get_grid_d(), GridforceFullBaseGrid::k, GridforceFullBaseGrid::offset, and TRUE.

{
    for (int i0 = 0; i0 < 8; i0++) {
        int inds2[3];
        int zero_derivs = FALSE;
        
        float voff = 0.0;
        int bit = 1;    // bit = 2^i1 in the below loop
        for (int i1 = 0; i1 < 3; i1++) {
            inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
            
            // Deal with voltage offsets
            if (cont[i1] && inds[i1] == (k[i1]-1) && inds2[i1] == 0) {
                voff += offset[i1];
                DebugM(3, "offset[" << i1 << "] = " << offset[i1] << "\n" << endi);
            }
            
            bit <<= 1;  // i.e. multiply by 2
        }
        
        DebugM(1, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
        
        // NOTE: leaving everything in terms of unit cell coordinates for now,
        // eventually will multiply by inv tensor when applying the force
        
        // First set variables 'dk_{hi,lo}' (glob notation). The 'hi'
        // ('lo') variable in a given dimension is the number added (subtracted)
        // to go up (down) one grid point in that dimension; both are normally
        // just the corresponding 'dk[i]'. However, if we are sitting on a
        // boundary and we are using a continuous grid, then we want to map the
        // next point off the grid back around to the other side. e.g. point
        // (k[0], i1, k) maps to point (0, i1, k), which would be
        // accomplished by changing 'dk1_hi' to -(k[0]-1)*dk1.
        
        int d_hi[3] = {1, 1, 1};
        int d_lo[3] = {1, 1, 1};
        float voffs[3];
        float dscales[3] = {0.5, 0.5, 0.5};
        for (int i1 = 0; i1 < 3; i1++) {
            if (inds2[i1] == 0) {
                if (cont[i1]) {
                    d_lo[i1] = -(k[i1]-1);
                    voffs[i1] = offset[i1];
                    dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
                }
                else zero_derivs = TRUE;
            }
            else if (inds2[i1] == k[i1]-1) {
                if (cont[i1]) {
                    d_hi[i1] = -(k[i1]-1);
                    voffs[i1] = offset[i1];
                    dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
                }
                else zero_derivs = TRUE;
            }
            else {
                voffs[i1] = 0.0;
            }
        }
        
//      DebugM(2, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << "\n" << endi);
//      DebugM(2, "zero_derivs = " << zero_derivs << "\n" << endi);
//      DebugM(2, "d_hi = " << d_hi[0] << " " << d_hi[1] << " " << d_hi[2] << "\n" << endi);
//      DebugM(2, "d_lo = " << d_lo[0] << " " << d_lo[1] << " " << d_lo[2] << "\n" << endi);
        DebugM(1, "dscales = " << dscales[0] << " " << dscales[1] << " " << dscales[2] << "\n" << endi);
        DebugM(1, "voffs = " << voffs[0] << " " << voffs[1] << " " << voffs[2] << "\n" << endi);
        
        // V
        b[i0] = get_grid(inds2[0],inds2[1],inds2[2]) + voff;
        
        if (zero_derivs) {
            DebugM(2, "zero_derivs\n" << endi);
            b[8+i0] = 0.0;
            b[16+i0] = 0.0;
            b[24+i0] = 0.0;
            b[32+i0] = 0.0;
            b[40+i0] = 0.0;
            b[48+i0] = 0.0;
            b[56+i0] = 0.0;
        } else {
            b[8+i0]  = dscales[0] * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]) + voffs[0]);       //  dV/dx
            b[16+i0] = dscales[1] * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]) + voffs[1]);       //  dV/dy
            b[24+i0] = dscales[2] * (get_grid_d(inds2[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1],inds2[2]-d_lo[2]) + voffs[2]);       //  dV/dz
            b[32+i0] = dscales[0] * dscales[1]
                * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]) -
                   get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]));    //  d2V/dxdy
            b[40+i0] = dscales[0] * dscales[2]
                * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]+d_hi[2]) -
                   get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]-d_lo[2]));    //  d2V/dxdz
            b[48+i0] = dscales[1] * dscales[2]
                * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) -
                   get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));    //  d2V/dydz
        
            b[56+i0] = dscales[0] * dscales[1] * dscales[2]                                     // d3V/dxdydz
                * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) -
                   get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) +
                   get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) +
                   get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));
        }
        
        DebugM(1, "V = " << b[i0] << "\n");
        
        DebugM(1, "dV/dx = " << b[8+i0] << "\n");
        DebugM(1, "dV/dy = " << b[16+i0] << "\n");
        DebugM(1, "dV/dz = " << b[24+i0] << "\n");
        
        DebugM(1, "d2V/dxdy = " << b[32+i0] << "\n");
        DebugM(1, "d2V/dxdz = " << b[40+i0] << "\n");
        DebugM(1, "d2V/dydz = " << b[48+i0] << "\n");
        
        DebugM(1, "d3V/dxdydz = " << b[56+i0] << "\n" << endi);
    }
}

int GridforceFullMainGrid::compute_VdV ( Position  pos, float &  V, Vector &  dV  ) const

inline

Definition at line 216 of file GridForceGrid.h.

References GridforceFullBaseGrid::compute_VdV().

{ return GridforceFullBaseGrid::compute_VdV(pos, V, dV); };

long int GridforceFullMainGrid::get_all_gridvals ( float **  all_gridvals ) const

protectedvirtual

Implements GridforceGrid.

Definition at line 718 of file GridForceGrid.C.

References DebugM, endi(), GridforceFullBaseGrid::grid, GridforceFullBaseGrid::size, subgrids_flat, and totalGrids.

{
    // Creates a flat array of all grid values, including subgrids,
    // and puts it in the value pointed to by the 'grids'
    // argument. Returns the resulting array size. Caller is
    // responsible for destroying the array via 'delete[]'
    
    DebugM(4, "get_all_gridvals called\n" << endi);
    
    long int sz = 0;
    sz += size;
    for (int i = 0; i < totalGrids-1; i++) {
        sz += subgrids_flat[i]->size;
    }
    DebugM(4, "size = " << sz << "\n" << endi);
    
    float *grid_vals = new float[sz];
    long int idx = 0;
    for (long int i = 0; i < size; i++) {
        grid_vals[idx++] = grid[i];
    }
    for (int j = 0; j < totalGrids-1; j++) {
        for (long int i = 0; i < subgrids_flat[j]->size; i++) {
            grid_vals[idx++] = subgrids_flat[j]->grid[i];
        }
    }
    CmiAssert(idx == sz);
    
    *all_gridvals = grid_vals;
    
    DebugM(4, "get_all_gridvals finished\n" << endi);
    
    return sz;
}

int GridforceFullMainGrid::get_border ( void  ) const

inlinevirtual

Implements GridforceFullBaseGrid.

Definition at line 214 of file GridForceGrid.h.

References border.

{ return border; }

Position GridforceFullMainGrid::get_center ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 205 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_center().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_center(); };

Bool GridforceFullMainGrid::get_checksize ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 210 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_checksize().

{ return GridforceFullBaseGrid::get_checksize(); };

Tensor GridforceFullMainGrid::get_e ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 207 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_e().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_e(); };

Tensor GridforceFullMainGrid::get_inv ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 208 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_inv().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_inv(); };

int GridforceFullMainGrid::get_k0 ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 211 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_k0().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_k0(); };

int GridforceFullMainGrid::get_k1 ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 212 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_k1().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_k1(); };

int GridforceFullMainGrid::get_k2 ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 213 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_k2().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_k2(); };

Position GridforceFullMainGrid::get_origin ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 206 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_origin().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_origin(); };

Vector GridforceFullMainGrid::get_scale ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 209 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_scale().

Referenced by GridforceLiteGrid::initialize().

{ return GridforceFullBaseGrid::get_scale(); };

int GridforceFullMainGrid::get_total_grids ( void  ) const

inlinevirtual

Implements GridforceGrid.

Definition at line 218 of file GridForceGrid.h.

References totalGrids.

Referenced by GridforceLiteGrid::initialize().

{ return totalGrids; }    

void GridforceFullMainGrid::initialize	(	char *	potfilename,
		SimParameters *	simParams,
		MGridforceParams *	mgridParams,
		int	border
	)

Definition at line 426 of file GridForceGrid.C.

References Lattice::a(), Lattice::b(), border, buildSubgridsFlat(), Lattice::c(), GridforceFullBaseGrid::checksize, GridforceFullBaseGrid::cont, cross(), DebugM, default_border, GridforceFullBaseGrid::dk, GridforceFullBaseGrid::dk_nopad, GridforceFullBaseGrid::e, endi(), GridforceFullBaseGrid::factor, FALSE, filename, GridforceGrid::fits_lattice(), Fopen(), GridforceFullBaseGrid::gap, GridforceFullBaseGrid::gapinv, GridforceFullBaseGrid::get_grid(), GridforceFullBaseGrid::grid, MGridforceParams::gridforceCheckSize, MGridforceParams::gridforceCont, MGridforceParams::gridforceScale, MGridforceParams::gridforceVOffset, MGridforceParams::gridforceVolts, GridforceFullSubGrid::initialize(), GridforceFullBaseGrid::inv, iout, iWARN(), GridforceFullBaseGrid::k, GridforceFullBaseGrid::k_nopad, SimParameters::lattice, GridforceGrid::mygridnum, GridforceFullBaseGrid::n_sum, NAMD_die(), GridforceFullBaseGrid::numSubgrids, GridforceFullBaseGrid::offset, GridforceFullBaseGrid::origin, GridforceFullBaseGrid::p_sum, GridforceFullBaseGrid::pad_n, GridforceFullBaseGrid::pad_p, GridforceFullBaseGrid::poten_fp, GridforceFullBaseGrid::readHeader(), GridforceFullBaseGrid::readSubgridHierarchy(), GridforceFullBaseGrid::scale, GridforceFullBaseGrid::set_grid(), GridforceFullBaseGrid::size, GridforceFullBaseGrid::size_nopad, GridforceFullBaseGrid::subgrids, totalGrids, and TRUE.

Referenced by initialize(), GridforceLiteGrid::initialize(), and reinitialize().

{
    if (brd >= 0) {
        border = brd;
    } else {
        border = default_border;
    }
    
    // FROM init1
    //FILE *poten = Fopen(potfilename, "r");
    poten_fp = Fopen(potfilename, "r");
    if (!poten_fp) {
        NAMD_die("Problem reading grid force potential file");
    }
    
    // save file name so that grid can be re-read via Tcl
    strcpy(filename, potfilename);
    
    // Read special comment fields and create subgrid objects
    totalGrids = 1;
    char line[256];
    Bool flag = FALSE;
    numSubgrids = 0;
    float version;
    long int poten_offset;
    do {
        poten_offset = ftell(poten_fp);
        fgets(line, 256, poten_fp);     // Read comment lines
        //flag = sscanf(line, "# maingrid subgrids count %d\n", &numSubgrids);
        flag = sscanf(line, "# namdnugrid version %f\n", &version);
    } while (line[0] == '#' && !flag);
    
    if (flag) {
        if (version != 1.0) {
            NAMD_die("Unsupported version of non-uniform grid file format!");
        }
        fscanf(poten_fp, "# namdnugrid maingrid subgrids count %d\n", &numSubgrids);
        readSubgridHierarchy(poten_fp, totalGrids);
        buildSubgridsFlat();
    } else {
        fseek(poten_fp, poten_offset, SEEK_SET);
    }
    
    // Read header
    readHeader(simParams, mgridParams);
    
    factor = 1.0;
    if (mgridParams->gridforceVolts)
    {
        factor /= 0.0434;  // convert V -> kcal/mol*e
    }
    scale = mgridParams->gridforceScale;
    checksize = mgridParams->gridforceCheckSize;

    // Allocate storage for potential and read it
    float *grid_nopad = new float[size_nopad];
    
    float tmp2;
    for (long int count = 0; count < size_nopad; count++) {
        int err = fscanf(poten_fp, "%f", &tmp2);
        if (err == EOF || err == 0) {
            NAMD_die("Grid force potential file incorrectly formatted");
        }
        grid_nopad[count] = tmp2 * factor;      // temporary, so just store flat
    }
    fscanf(poten_fp, "\n");
    
    // Shortcuts for accessing 1-D array with four indices
    dk_nopad[0] = k_nopad[1] * k_nopad[2];
    dk_nopad[1] = k_nopad[2];
    dk_nopad[2] = 1;
    
    Vector Kvec[3];
    Kvec[0] = e * Position(k_nopad[0]-1, 0, 0);
    Kvec[1] = e * Position(0, k_nopad[1]-1, 0);
    Kvec[2] = e * Position(0, 0, k_nopad[2]-1);
    Vector Avec[3];
    Avec[0] = simParams->lattice.a();
    Avec[1] = simParams->lattice.b();
    Avec[2] = simParams->lattice.c();
    
    // Decide whether we're wrapping
    for (int i0 = 0; i0 < 3; i0++) {
        if (mgridParams->gridforceCont[i0])
        {
            Bool found = FALSE;
            for (int i1 = 0; i1 < 3; i1++) {
                if (cross(Avec[i0].unit(), Kvec[i1].unit()).length() < 1e-4) {
                    found = TRUE;
                    cont[i1] = TRUE;
                    offset[i1] = mgridParams->gridforceVOffset[i0] * factor;
                    // want in grid-point units (normal = 1)
                    gap[i1] = (inv * (Avec[i0] - Kvec[i1])).length();   
                    gapinv[i1] = 1.0/gap[i1];
                    
                    if (gap[i1] < 0) {
                        NAMD_die("Gridforce Grid overlap!");
                    }
                    
                    DebugM(4, "cont[" << i1 << "] = " << cont[i1] << "\n");
                    DebugM(4, "gap[" << i1 << "] = " << gap[i1] << "\n");
                    DebugM(4, "gapinv[" << i1 << "] = " << gapinv[i1] << "\n" << endi);
                }
            }
            
            if (!found) {
                NAMD_die("No Gridforce unit vector found parallel to requested continuous grid direction!");
            }
        } else {
            // check for grid overlap in non-wrapping dimensions
            // occurs below
        }
    }
    
    // Figure out size of true grid (padded on non-periodic sides)
    Vector delta = 0;
    for (int i = 0; i < 3; i++) {
        if (cont[i]) {
            k[i] = k_nopad[i];
        } else {
            k[i] = k_nopad[i] + 2*border;
            delta[i] -= border;
        }
    }
    DebugM(4, "delta = " << e * delta << " (" << delta << ")\n" << endi);
    origin += e * delta;
    
    // Check for grid overlap
    if (!fits_lattice(simParams->lattice)) {
      char errmsg[512];
      if (checksize) {
        sprintf(errmsg, "Warning: Periodic cell basis too small for Gridforce grid %d.  Set gridforcechecksize off in configuration file to ignore.\n", mygridnum);
        NAMD_die(errmsg);      
      }
    }
    
    size = k[0] * k[1] * k[2];
    dk[0] = k[1] * k[2];
    dk[1] = k[2];
    dk[2] = 1;
    
    DebugM(3, "size = " << size << ", size_nopad = " << size_nopad << "\n" << endi);
    
    delete[] grid;
    grid = new float[size];
    
    n_sum[0] = n_sum[1] = n_sum[2] = 0;
    p_sum[0] = p_sum[1] = p_sum[2] = 0;
    for (int i0 = 0; i0 < k_nopad[0]; i0++) {
        for (int i1 = 0; i1 < k_nopad[1]; i1++) {
            for (int i2 = 0; i2 < k_nopad[2]; i2++) {
                // Edges are special cases -- take force there to be
                // zero for smooth transition across potential
                // boundary
                
                long int ind_nopad = i0*dk_nopad[0] + i1*dk_nopad[1] + i2*dk_nopad[2];
                int j0 = (cont[0]) ? i0 : i0 + border;
                int j1 = (cont[1]) ? i1 : i1 + border;
                int j2 = (cont[2]) ? i2 : i2 + border;
                long int ind = j0*dk[0] + j1*dk[1] + j2*dk[2];
                
                if (i0 == 0)                    n_sum[0] += grid_nopad[ind_nopad];
                else if (i0 == k_nopad[0]-1)    p_sum[0] += grid_nopad[ind_nopad];
                if (i1 == 0)                    n_sum[1] += grid_nopad[ind_nopad];
                else if (i1 == k_nopad[1]-1)    p_sum[1] += grid_nopad[ind_nopad];
                if (i2 == 0)                    n_sum[2] += grid_nopad[ind_nopad];
                else if (i2 == k_nopad[2]-1)    p_sum[2] += grid_nopad[ind_nopad];
                
                //grid[ind] = grid_nopad[ind_nopad];
                set_grid(j0, j1, j2, grid_nopad[ind_nopad]);
            }
        }
    }
    
    const BigReal modThresh = 1.0;
    
    BigReal n_avg[3], p_avg[3];
    int i0;
    for (int i0 = 0; i0 < 3; i0++) {
        int i1 = (i0 + 1) % 3;
        int i2 = (i0 + 2) % 3;
        n_avg[i0] = n_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
        p_avg[i0] = p_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
        
        if (cont[i0] && fabs(offset[i0] - (p_avg[i0]-n_avg[i0])) > modThresh) 
        {
            iout << iWARN << "GRID FORCE POTENTIAL DIFFERENCE IN K" << i0
                 << " DIRECTION IS " 
                 << offset[i0] - (p_avg[i0]-n_avg[i0]) 
                 << " KCAL/MOL*E\n" << endi;
        }
    }
    
    Bool twoPadVals = (cont[0] + cont[1] + cont[2] == 2);
    double padVal = 0.0;
    long int weight = 0;
    if (!twoPadVals) {
        // Determine pad value (must average)
        if (!cont[0]) {
            padVal += p_sum[0] + n_sum[0];
            weight += 2 * k_nopad[1] * k_nopad[2];
        }
        if (!cont[1]) {
            padVal += p_sum[1] + n_sum[1];
            weight += 2 * k_nopad[0] * k_nopad[2];
        }
        if (!cont[2]) {
            padVal += p_sum[2] + n_sum[2];
            weight += 2 * k_nopad[0] * k_nopad[1];
        }
        padVal /= weight;
    }
    
    for (int i = 0; i < 3; i++) {
        pad_n[i] = (cont[i]) ? 0.0 : (twoPadVals) ? n_avg[i] : padVal;
        pad_p[i] = (cont[i]) ? 0.0 : (twoPadVals) ? p_avg[i] : padVal;
        DebugM(4, "pad_n[" << i << "] = " << pad_n[i] << "\n");
        DebugM(4, "pad_p[" << i << "] = " << pad_p[i] << "\n" << endi);
    }
    
    if (cont[0] && cont[1] && cont[2]) {
        // Nothing to do
        return;
    }
    
    // Now fill in rest of new grid
    for (int i0 = 0; i0 < k[0]; i0++) {
        for (int i1 = 0; i1 < k[1]; i1++) {
            for (int i2 = 0; i2 < k[2]; i2++) {
                if ( (cont[0] || (i0 >= border && i0 < k[0]-border)) 
                     && (cont[1] || (i1 >= border && i1 < k[1]-border)) 
                     && (cont[2] || i2 == border) )
                {
                    i2 += k_nopad[2]-1;
                    continue;
                }
                
                long int ind = i0*dk[0] + i1*dk[1] + i2*dk[2];

                Position pos = e * Position(i0, i1, i2);
                int var[3] = {i0, i1, i2};
                
                for (int dir = 0; dir < 3; dir++) {
                    if (cont[dir]) 
                        continue;
                    
                    if (var[dir] < border) {
                        //grid[ind] = pad_n[dir];
                        set_grid(i0, i1, i2, pad_n[dir]);
                    } else if (var[dir] >= k[dir]-border) {
                        //grid[ind] = pad_p[dir];
                        set_grid(i0, i1, i2, pad_p[dir]);
                    }
                }
                
//              DebugM(2, "grid[" << ind << "; " << i0 << ", " << i1
//                     << ", " << i2 << "] = " << get_grid(ind)
//                     << "\n" << endi);
            }
        }
    }
    
    for (int i0 = 0; i0 < k[0]; i0++) {
        for (int i1 = 0; i1 < k[1]; i1++) {
            for (int i2 = 0; i2 < k[2]; i2++) {
                DebugM(1, "grid[" << i0 << ", " << i1 << ", " << i2 << "] = " << get_grid(i0,i1,i2) << "\n" << endi);
            }
        }
    }
    
    // Clean up
    DebugM(3, "clean up\n" << endi);
    delete[] grid_nopad;
    
    // Call initialize for each subgrid
    for (int i = 0; i < numSubgrids; i++) {
        subgrids[i]->poten_fp = poten_fp;
        subgrids[i]->initialize(simParams, mgridParams);
    }
    
    // close file pointer
    fclose(poten_fp);
}

void GridforceFullMainGrid::initialize ( char *  potfilename, SimParameters *  simParams, MGridforceParams *  mgridParams  )

inlinevirtual

Implements GridforceGrid.

Definition at line 200 of file GridForceGrid.h.

References default_border, and initialize().

                                                                                                       {
        initialize(potfilename, simParams, mgridParams, default_border);
    }

void GridforceFullMainGrid::pack ( MOStream *  msg ) const

protectedvirtual

Implements GridforceGrid.

Definition at line 369 of file GridForceGrid.C.

References DebugM, endi(), filename, GridforceGrid::mygridnum, GridforceFullBaseGrid::pack(), MOStream::put(), and totalGrids.

{
    DebugM(4, "Packing maingrid\n" << endi);
    
//     msg->put(3*sizeof(float), (char*)pad_p);
//     msg->put(3*sizeof(float), (char*)pad_n);
    msg->put(totalGrids);
    msg->put(mygridnum);
    msg->put(129*sizeof(char), (char*)filename);
    
    DebugM(3, "calling GridforceFullBaseGrid::pack\n" << endi);
    
    GridforceFullBaseGrid::pack(msg);
}

void GridforceFullMainGrid::reinitialize ( SimParameters *  simParams, MGridforceParams *  mgridParams  )

virtual

Implements GridforceGrid.

Definition at line 711 of file GridForceGrid.C.

References DebugM, endi(), filename, and initialize().

{
    DebugM(4, "reinitializing grid\n" << endi);
    initialize(filename, simParams, mgridParams);
}

void GridforceFullMainGrid::set_all_gridvals ( float *  all_gridvals, long int  sz  )

protectedvirtual

Implements GridforceGrid.

Definition at line 754 of file GridForceGrid.C.

References DebugM, endi(), GridforceFullBaseGrid::grid, GridforceFullBaseGrid::size, subgrids_flat, and totalGrids.

{
    DebugM(4, "set_all_gridvals called\n" << endi);
    
    long int sz_calc = 0;
    sz_calc += size;
    for (int i = 0; i < totalGrids-1; i++) {
        sz_calc += subgrids_flat[i]->size;
    }
    CmiAssert(sz == sz_calc);
    
    long int idx = 0;
    for (long int i = 0; i < size; i++) {
        DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
        grid[i] = all_gridvals[idx++];
    }
    for (int j = 0; j < totalGrids-1; j++) {
        for (long int i = 0; i < subgrids_flat[j]->size; i++) {
            DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
            subgrids_flat[j]->grid[i] = all_gridvals[idx++];
        }
    }
    CmiAssert(idx == sz);

    DebugM(4, "set_all_gridvals finished\n" << endi);
}

void GridforceFullMainGrid::set_scale ( Vector  s )

inlinevirtual

Implements GridforceGrid.

Definition at line 219 of file GridForceGrid.h.

References GridforceFullBaseGrid::scale.

{ scale = s; }

void GridforceFullMainGrid::unpack ( MIStream *  msg )

protectedvirtual

Implements GridforceGrid.

Definition at line 385 of file GridForceGrid.C.

References buildSubgridsFlat(), DebugM, endi(), filename, GridforceFullBaseGrid::gapinv, GridforceFullBaseGrid::generation, MIStream::get(), GridforceGrid::mygridnum, GridforceFullBaseGrid::numSubgrids, GridforceFullBaseGrid::size, totalGrids, and GridforceFullBaseGrid::unpack().

{
    DebugM(4, "Unpacking maingrid\n" << endi);
    
//     msg->get(3*sizeof(float), (char*)pad_p);
//     msg->get(3*sizeof(float), (char*)pad_n);
    msg->get(totalGrids);
    msg->get(mygridnum);
    msg->get(129*sizeof(char), (char*)filename);
    
    GridforceFullBaseGrid::unpack(msg);
    
    DebugM(4, "size  = " << size << "\n");
    DebugM(4, "numSubgrids = " << numSubgrids << "\n");
    DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
    DebugM(4, "generation = " << generation << "\n" << endi);
    DebugM(4, "filename = " << filename << "\n" << endi);
    
    buildSubgridsFlat();
}

Friends And Related Function Documentation

friend class GridforceFullBaseGrid

friend

Definition at line 192 of file GridForceGrid.h.

friend class GridforceFullSubGrid

friend

Definition at line 193 of file GridForceGrid.h.

Member Data Documentation

int GridforceFullMainGrid::border

protected

Definition at line 238 of file GridForceGrid.h.

Referenced by get_border(), and initialize().

const int GridforceFullMainGrid::default_border = 1

staticprotected

Definition at line 237 of file GridForceGrid.h.

Referenced by initialize().

char GridforceFullMainGrid::filename[NAMD_FILENAME_BUFFER_SIZE]

protected

Definition at line 232 of file GridForceGrid.h.

Referenced by initialize(), pack(), reinitialize(), and unpack().

GridforceFullSubGrid** GridforceFullMainGrid::subgrids_flat

protected

Definition at line 234 of file GridForceGrid.h.

Referenced by GridforceFullSubGrid::addToSubgridsFlat(), buildSubgridsFlat(), get_all_gridvals(), GridforceFullMainGrid(), set_all_gridvals(), and ~GridforceFullMainGrid().

int GridforceFullMainGrid::totalGrids

protected

Definition at line 233 of file GridForceGrid.h.

Referenced by buildSubgridsFlat(), get_all_gridvals(), get_total_grids(), initialize(), pack(), set_all_gridvals(), and unpack().

---

The documentation for this class was generated from the following files:

• GridForceGrid.h
• GridForceGrid.C