NAMD
Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | Friends | List of all members
GridforceFullMainGrid Class Reference

#include <GridForceGrid.h>

Inheritance diagram for GridforceFullMainGrid:
GridforceGrid GridforceFullBaseGrid

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 GridforceGridnew_grid (int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
 
static void pack_grid (GridforceGrid *grid, MOStream *msg)
 
static GridforceGridunpack_grid (int gridnum, MIStream *msg)
 

Detailed Description

Definition at line 191 of file GridForceGrid.h.

Constructor & Destructor Documentation

GridforceFullMainGrid::GridforceFullMainGrid ( int  gridnum)
explicit
GridforceFullMainGrid::~GridforceFullMainGrid ( )
virtual

Definition at line 363 of file GridForceGrid.C.

References subgrids_flat.

364 {
365  delete[] subgrids_flat;
366 }
GridforceFullSubGrid ** 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().

408 {
409  DebugM(4, "buildSubgridsFlat() called, totalGrids-1 = " << totalGrids-1 << "\n" << endi);
410  delete[] subgrids_flat;
412  for (int i = 0; i < numSubgrids; i++) {
413  DebugM(3, "adding to subgridsFlat\n" << endi);
415  DebugM(3, "success!\n" << endi);
416  }
417  for (int i = 0; i < totalGrids-1; i++) {
418  DebugM(4, "subgrids_flat[" << i << "]->numSubgrids = " << subgrids_flat[i]->numSubgrids << "\n" << endi);
419  }
420  for (int i = 0; i < numSubgrids; i++) {
421  DebugM(4, "subgrids[" << i << "]->numSubgrids = " << subgrids[i]->numSubgrids << "\n" << endi);
422  }
423 }
#define DebugM(x, y)
Definition: Debug.h:59
GridforceFullSubGrid ** subgrids
GridforceFullSubGrid ** subgrids_flat
void addToSubgridsFlat(void)
infostream & endi(infostream &s)
Definition: InfoStream.C:38
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.

783 {
784  for (int i0 = 0; i0 < 8; i0++) {
785  int inds2[3];
786  int zero_derivs = FALSE;
787 
788  float voff = 0.0;
789  int bit = 1; // bit = 2^i1 in the below loop
790  for (int i1 = 0; i1 < 3; i1++) {
791  inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
792 
793  // Deal with voltage offsets
794  if (cont[i1] && inds[i1] == (k[i1]-1) && inds2[i1] == 0) {
795  voff += offset[i1];
796  DebugM(3, "offset[" << i1 << "] = " << offset[i1] << "\n" << endi);
797  }
798 
799  bit <<= 1; // i.e. multiply by 2
800  }
801 
802  DebugM(1, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
803 
804  // NOTE: leaving everything in terms of unit cell coordinates for now,
805  // eventually will multiply by inv tensor when applying the force
806 
807  // First set variables 'dk_{hi,lo}' (glob notation). The 'hi'
808  // ('lo') variable in a given dimension is the number added (subtracted)
809  // to go up (down) one grid point in that dimension; both are normally
810  // just the corresponding 'dk[i]'. However, if we are sitting on a
811  // boundary and we are using a continuous grid, then we want to map the
812  // next point off the grid back around to the other side. e.g. point
813  // (k[0], i1, k) maps to point (0, i1, k), which would be
814  // accomplished by changing 'dk1_hi' to -(k[0]-1)*dk1.
815 
816  int d_hi[3] = {1, 1, 1};
817  int d_lo[3] = {1, 1, 1};
818  float voffs[3];
819  float dscales[3] = {0.5, 0.5, 0.5};
820  for (int i1 = 0; i1 < 3; i1++) {
821  if (inds2[i1] == 0) {
822  if (cont[i1]) {
823  d_lo[i1] = -(k[i1]-1);
824  voffs[i1] = offset[i1];
825  dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
826  }
827  else zero_derivs = TRUE;
828  }
829  else if (inds2[i1] == k[i1]-1) {
830  if (cont[i1]) {
831  d_hi[i1] = -(k[i1]-1);
832  voffs[i1] = offset[i1];
833  dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
834  }
835  else zero_derivs = TRUE;
836  }
837  else {
838  voffs[i1] = 0.0;
839  }
840  }
841 
842 // DebugM(2, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << "\n" << endi);
843 // DebugM(2, "zero_derivs = " << zero_derivs << "\n" << endi);
844 // DebugM(2, "d_hi = " << d_hi[0] << " " << d_hi[1] << " " << d_hi[2] << "\n" << endi);
845 // DebugM(2, "d_lo = " << d_lo[0] << " " << d_lo[1] << " " << d_lo[2] << "\n" << endi);
846  DebugM(1, "dscales = " << dscales[0] << " " << dscales[1] << " " << dscales[2] << "\n" << endi);
847  DebugM(1, "voffs = " << voffs[0] << " " << voffs[1] << " " << voffs[2] << "\n" << endi);
848 
849  // V
850  b[i0] = get_grid(inds2[0],inds2[1],inds2[2]) + voff;
851 
852  if (zero_derivs) {
853  DebugM(2, "zero_derivs\n" << endi);
854  b[8+i0] = 0.0;
855  b[16+i0] = 0.0;
856  b[24+i0] = 0.0;
857  b[32+i0] = 0.0;
858  b[40+i0] = 0.0;
859  b[48+i0] = 0.0;
860  b[56+i0] = 0.0;
861  } else {
862  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
863  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
864  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
865  b[32+i0] = dscales[0] * dscales[1]
866  * (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]) -
867  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
868  b[40+i0] = dscales[0] * dscales[2]
869  * (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]) -
870  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
871  b[48+i0] = dscales[1] * dscales[2]
872  * (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]) -
873  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
874 
875  b[56+i0] = dscales[0] * dscales[1] * dscales[2] // d3V/dxdydz
876  * (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]) -
877  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]) +
878  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]) +
879  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]));
880  }
881 
882  DebugM(1, "V = " << b[i0] << "\n");
883 
884  DebugM(1, "dV/dx = " << b[8+i0] << "\n");
885  DebugM(1, "dV/dy = " << b[16+i0] << "\n");
886  DebugM(1, "dV/dz = " << b[24+i0] << "\n");
887 
888  DebugM(1, "d2V/dxdy = " << b[32+i0] << "\n");
889  DebugM(1, "d2V/dxdz = " << b[40+i0] << "\n");
890  DebugM(1, "d2V/dydz = " << b[48+i0] << "\n");
891 
892  DebugM(1, "d3V/dxdydz = " << b[56+i0] << "\n" << endi);
893  }
894 }
float get_grid(int i0, int i1, int i2) const
#define DebugM(x, y)
Definition: Debug.h:59
#define FALSE
Definition: common.h:118
double get_grid_d(int i0, int i1, int i2) const
infostream & endi(infostream &s)
Definition: InfoStream.C:38
#define TRUE
Definition: common.h:119
int GridforceFullMainGrid::compute_VdV ( Position  pos,
float &  V,
Vector dV 
) const
inline

Definition at line 216 of file GridForceGrid.h.

References GridforceFullBaseGrid::compute_VdV().

216 { return GridforceFullBaseGrid::compute_VdV(pos, V, dV); };
int compute_VdV(Position pos, float &V, Vector &dV) const
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.

719 {
720  // Creates a flat array of all grid values, including subgrids,
721  // and puts it in the value pointed to by the 'grids'
722  // argument. Returns the resulting array size. Caller is
723  // responsible for destroying the array via 'delete[]'
724 
725  DebugM(4, "get_all_gridvals called\n" << endi);
726 
727  long int sz = 0;
728  sz += size;
729  for (int i = 0; i < totalGrids-1; i++) {
730  sz += subgrids_flat[i]->size;
731  }
732  DebugM(4, "size = " << sz << "\n" << endi);
733 
734  float *grid_vals = new float[sz];
735  long int idx = 0;
736  for (long int i = 0; i < size; i++) {
737  grid_vals[idx++] = grid[i];
738  }
739  for (int j = 0; j < totalGrids-1; j++) {
740  for (long int i = 0; i < subgrids_flat[j]->size; i++) {
741  grid_vals[idx++] = subgrids_flat[j]->grid[i];
742  }
743  }
744  CmiAssert(idx == sz);
745 
746  *all_gridvals = grid_vals;
747 
748  DebugM(4, "get_all_gridvals finished\n" << endi);
749 
750  return sz;
751 }
#define DebugM(x, y)
Definition: Debug.h:59
GridforceFullSubGrid ** subgrids_flat
infostream & endi(infostream &s)
Definition: InfoStream.C:38
int GridforceFullMainGrid::get_border ( void  ) const
inlinevirtual

Implements GridforceFullBaseGrid.

Definition at line 214 of file GridForceGrid.h.

References border.

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

Position get_center(void) const
Definition: GridForceGrid.h:93
Bool GridforceFullMainGrid::get_checksize ( void  ) const
inlinevirtual

Implements GridforceGrid.

Definition at line 210 of file GridForceGrid.h.

References GridforceFullBaseGrid::get_checksize().

Bool get_checksize(void) const
Definition: GridForceGrid.h:98
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().

207 { return GridforceFullBaseGrid::get_e(); };
Tensor get_e(void) const
Definition: GridForceGrid.h:95
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().

208 { return GridforceFullBaseGrid::get_inv(); };
Tensor get_inv(void) const
Definition: GridForceGrid.h:96
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().

211 { return GridforceFullBaseGrid::get_k0(); };
int get_k0(void) const
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().

212 { return GridforceFullBaseGrid::get_k1(); };
int get_k1(void) const
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().

213 { return GridforceFullBaseGrid::get_k2(); };
int get_k2(void) const
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().

Position get_origin(void) const
Definition: GridForceGrid.h:94
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().

Vector get_scale(void) const
Definition: GridForceGrid.h:97
int GridforceFullMainGrid::get_total_grids ( void  ) const
inlinevirtual

Implements GridforceGrid.

Definition at line 218 of file GridForceGrid.h.

References totalGrids.

Referenced by GridforceLiteGrid::initialize().

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

427 {
428  if (brd >= 0) {
429  border = brd;
430  } else {
432  }
433 
434  // FROM init1
435  //FILE *poten = Fopen(potfilename, "r");
436  poten_fp = Fopen(potfilename, "r");
437  if (!poten_fp) {
438  NAMD_die("Problem reading grid force potential file");
439  }
440 
441  // save file name so that grid can be re-read via Tcl
442  strcpy(filename, potfilename);
443 
444  // Read special comment fields and create subgrid objects
445  totalGrids = 1;
446  char line[256];
447  Bool flag = FALSE;
448  numSubgrids = 0;
449  float version;
450  long int poten_offset;
451  do {
452  poten_offset = ftell(poten_fp);
453  fgets(line, 256, poten_fp); // Read comment lines
454  //flag = sscanf(line, "# maingrid subgrids count %d\n", &numSubgrids);
455  flag = sscanf(line, "# namdnugrid version %f\n", &version);
456  } while (line[0] == '#' && !flag);
457 
458  if (flag) {
459  if (version != 1.0) {
460  NAMD_die("Unsupported version of non-uniform grid file format!");
461  }
462  fscanf(poten_fp, "# namdnugrid maingrid subgrids count %d\n", &numSubgrids);
465  } else {
466  fseek(poten_fp, poten_offset, SEEK_SET);
467  }
468 
469  // Read header
470  readHeader(simParams, mgridParams);
471 
472  factor = 1.0;
473  if (mgridParams->gridforceVolts)
474  {
475  factor /= 0.0434; // convert V -> kcal/mol*e
476  }
477  scale = mgridParams->gridforceScale;
478  checksize = mgridParams->gridforceCheckSize;
479 
480  // Allocate storage for potential and read it
481  float *grid_nopad = new float[size_nopad];
482 
483  float tmp2;
484  for (long int count = 0; count < size_nopad; count++) {
485  int err = fscanf(poten_fp, "%f", &tmp2);
486  if (err == EOF || err == 0) {
487  NAMD_die("Grid force potential file incorrectly formatted");
488  }
489  grid_nopad[count] = tmp2 * factor; // temporary, so just store flat
490  }
491  fscanf(poten_fp, "\n");
492 
493  // Shortcuts for accessing 1-D array with four indices
494  dk_nopad[0] = k_nopad[1] * k_nopad[2];
495  dk_nopad[1] = k_nopad[2];
496  dk_nopad[2] = 1;
497 
498  Vector Kvec[3];
499  Kvec[0] = e * Position(k_nopad[0]-1, 0, 0);
500  Kvec[1] = e * Position(0, k_nopad[1]-1, 0);
501  Kvec[2] = e * Position(0, 0, k_nopad[2]-1);
502  Vector Avec[3];
503  Avec[0] = simParams->lattice.a();
504  Avec[1] = simParams->lattice.b();
505  Avec[2] = simParams->lattice.c();
506 
507  // Decide whether we're wrapping
508  for (int i0 = 0; i0 < 3; i0++) {
509  if (mgridParams->gridforceCont[i0])
510  {
511  Bool found = FALSE;
512  for (int i1 = 0; i1 < 3; i1++) {
513  if (cross(Avec[i0].unit(), Kvec[i1].unit()).length() < 1e-4) {
514  found = TRUE;
515  cont[i1] = TRUE;
516  offset[i1] = mgridParams->gridforceVOffset[i0] * factor;
517  // want in grid-point units (normal = 1)
518  gap[i1] = (inv * (Avec[i0] - Kvec[i1])).length();
519  gapinv[i1] = 1.0/gap[i1];
520 
521  if (gap[i1] < 0) {
522  NAMD_die("Gridforce Grid overlap!");
523  }
524 
525  DebugM(4, "cont[" << i1 << "] = " << cont[i1] << "\n");
526  DebugM(4, "gap[" << i1 << "] = " << gap[i1] << "\n");
527  DebugM(4, "gapinv[" << i1 << "] = " << gapinv[i1] << "\n" << endi);
528  }
529  }
530 
531  if (!found) {
532  NAMD_die("No Gridforce unit vector found parallel to requested continuous grid direction!");
533  }
534  } else {
535  // check for grid overlap in non-wrapping dimensions
536  // occurs below
537  }
538  }
539 
540  // Figure out size of true grid (padded on non-periodic sides)
541  Vector delta = 0;
542  for (int i = 0; i < 3; i++) {
543  if (cont[i]) {
544  k[i] = k_nopad[i];
545  } else {
546  k[i] = k_nopad[i] + 2*border;
547  delta[i] -= border;
548  }
549  }
550  DebugM(4, "delta = " << e * delta << " (" << delta << ")\n" << endi);
551  origin += e * delta;
552 
553  // Check for grid overlap
554  if (!fits_lattice(simParams->lattice)) {
555  char errmsg[512];
556  if (checksize) {
557  sprintf(errmsg, "Warning: Periodic cell basis too small for Gridforce grid %d. Set gridforcechecksize off in configuration file to ignore.\n", mygridnum);
558  NAMD_die(errmsg);
559  }
560  }
561 
562  size = k[0] * k[1] * k[2];
563  dk[0] = k[1] * k[2];
564  dk[1] = k[2];
565  dk[2] = 1;
566 
567  DebugM(3, "size = " << size << ", size_nopad = " << size_nopad << "\n" << endi);
568 
569  delete[] grid;
570  grid = new float[size];
571 
572  n_sum[0] = n_sum[1] = n_sum[2] = 0;
573  p_sum[0] = p_sum[1] = p_sum[2] = 0;
574  for (int i0 = 0; i0 < k_nopad[0]; i0++) {
575  for (int i1 = 0; i1 < k_nopad[1]; i1++) {
576  for (int i2 = 0; i2 < k_nopad[2]; i2++) {
577  // Edges are special cases -- take force there to be
578  // zero for smooth transition across potential
579  // boundary
580 
581  long int ind_nopad = i0*dk_nopad[0] + i1*dk_nopad[1] + i2*dk_nopad[2];
582  int j0 = (cont[0]) ? i0 : i0 + border;
583  int j1 = (cont[1]) ? i1 : i1 + border;
584  int j2 = (cont[2]) ? i2 : i2 + border;
585  long int ind = j0*dk[0] + j1*dk[1] + j2*dk[2];
586 
587  if (i0 == 0) n_sum[0] += grid_nopad[ind_nopad];
588  else if (i0 == k_nopad[0]-1) p_sum[0] += grid_nopad[ind_nopad];
589  if (i1 == 0) n_sum[1] += grid_nopad[ind_nopad];
590  else if (i1 == k_nopad[1]-1) p_sum[1] += grid_nopad[ind_nopad];
591  if (i2 == 0) n_sum[2] += grid_nopad[ind_nopad];
592  else if (i2 == k_nopad[2]-1) p_sum[2] += grid_nopad[ind_nopad];
593 
594  //grid[ind] = grid_nopad[ind_nopad];
595  set_grid(j0, j1, j2, grid_nopad[ind_nopad]);
596  }
597  }
598  }
599 
600  const BigReal modThresh = 1.0;
601 
602  BigReal n_avg[3], p_avg[3];
603  int i0;
604  for (int i0 = 0; i0 < 3; i0++) {
605  int i1 = (i0 + 1) % 3;
606  int i2 = (i0 + 2) % 3;
607  n_avg[i0] = n_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
608  p_avg[i0] = p_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
609 
610  if (cont[i0] && fabs(offset[i0] - (p_avg[i0]-n_avg[i0])) > modThresh)
611  {
612  iout << iWARN << "GRID FORCE POTENTIAL DIFFERENCE IN K" << i0
613  << " DIRECTION IS "
614  << offset[i0] - (p_avg[i0]-n_avg[i0])
615  << " KCAL/MOL*E\n" << endi;
616  }
617  }
618 
619  Bool twoPadVals = (cont[0] + cont[1] + cont[2] == 2);
620  double padVal = 0.0;
621  long int weight = 0;
622  if (!twoPadVals) {
623  // Determine pad value (must average)
624  if (!cont[0]) {
625  padVal += p_sum[0] + n_sum[0];
626  weight += 2 * k_nopad[1] * k_nopad[2];
627  }
628  if (!cont[1]) {
629  padVal += p_sum[1] + n_sum[1];
630  weight += 2 * k_nopad[0] * k_nopad[2];
631  }
632  if (!cont[2]) {
633  padVal += p_sum[2] + n_sum[2];
634  weight += 2 * k_nopad[0] * k_nopad[1];
635  }
636  padVal /= weight;
637  }
638 
639  for (int i = 0; i < 3; i++) {
640  pad_n[i] = (cont[i]) ? 0.0 : (twoPadVals) ? n_avg[i] : padVal;
641  pad_p[i] = (cont[i]) ? 0.0 : (twoPadVals) ? p_avg[i] : padVal;
642  DebugM(4, "pad_n[" << i << "] = " << pad_n[i] << "\n");
643  DebugM(4, "pad_p[" << i << "] = " << pad_p[i] << "\n" << endi);
644  }
645 
646  if (cont[0] && cont[1] && cont[2]) {
647  // Nothing to do
648  return;
649  }
650 
651  // Now fill in rest of new grid
652  for (int i0 = 0; i0 < k[0]; i0++) {
653  for (int i1 = 0; i1 < k[1]; i1++) {
654  for (int i2 = 0; i2 < k[2]; i2++) {
655  if ( (cont[0] || (i0 >= border && i0 < k[0]-border))
656  && (cont[1] || (i1 >= border && i1 < k[1]-border))
657  && (cont[2] || i2 == border) )
658  {
659  i2 += k_nopad[2]-1;
660  continue;
661  }
662 
663  long int ind = i0*dk[0] + i1*dk[1] + i2*dk[2];
664 
665  Position pos = e * Position(i0, i1, i2);
666  int var[3] = {i0, i1, i2};
667 
668  for (int dir = 0; dir < 3; dir++) {
669  if (cont[dir])
670  continue;
671 
672  if (var[dir] < border) {
673  //grid[ind] = pad_n[dir];
674  set_grid(i0, i1, i2, pad_n[dir]);
675  } else if (var[dir] >= k[dir]-border) {
676  //grid[ind] = pad_p[dir];
677  set_grid(i0, i1, i2, pad_p[dir]);
678  }
679  }
680 
681 // DebugM(2, "grid[" << ind << "; " << i0 << ", " << i1
682 // << ", " << i2 << "] = " << get_grid(ind)
683 // << "\n" << endi);
684  }
685  }
686  }
687 
688  for (int i0 = 0; i0 < k[0]; i0++) {
689  for (int i1 = 0; i1 < k[1]; i1++) {
690  for (int i2 = 0; i2 < k[2]; i2++) {
691  DebugM(1, "grid[" << i0 << ", " << i1 << ", " << i2 << "] = " << get_grid(i0,i1,i2) << "\n" << endi);
692  }
693  }
694  }
695 
696  // Clean up
697  DebugM(3, "clean up\n" << endi);
698  delete[] grid_nopad;
699 
700  // Call initialize for each subgrid
701  for (int i = 0; i < numSubgrids; i++) {
702  subgrids[i]->poten_fp = poten_fp;
703  subgrids[i]->initialize(simParams, mgridParams);
704  }
705 
706  // close file pointer
707  fclose(poten_fp);
708 }
float get_grid(int i0, int i1, int i2) const
Definition: Vector.h:64
#define DebugM(x, y)
Definition: Debug.h:59
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
#define FALSE
Definition: common.h:118
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
#define iout
Definition: InfoStream.h:87
char filename[NAMD_FILENAME_BUFFER_SIZE]
void readHeader(SimParameters *simParams, MGridforceParams *mgridParams)
GridforceFullSubGrid ** subgrids
bool fits_lattice(const Lattice &lattice)
Definition: GridForceGrid.C:84
void set_grid(int i0, int i1, int i2, float V)
int Bool
Definition: common.h:133
FILE * Fopen(const char *filename, const char *mode)
Definition: common.C:265
Vector Position
Definition: NamdTypes.h:18
void readSubgridHierarchy(FILE *poten, int &totalGrids)
void NAMD_die(const char *err_msg)
Definition: common.C:83
void buildSubgridsFlat(void)
Vector b() const
Definition: Lattice.h:253
static const int default_border
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void initialize(SimParameters *simParams, MGridforceParams *mgridParams)
Vector a() const
Definition: Lattice.h:252
#define TRUE
Definition: common.h:119
Vector c() const
Definition: Lattice.h:254
double BigReal
Definition: common.h:114
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().

200  {
201  initialize(potfilename, simParams, mgridParams, default_border);
202  }
void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams, int border)
static const int 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.

370 {
371  DebugM(4, "Packing maingrid\n" << endi);
372 
373 // msg->put(3*sizeof(float), (char*)pad_p);
374 // msg->put(3*sizeof(float), (char*)pad_n);
375  msg->put(totalGrids);
376  msg->put(mygridnum);
377  msg->put(129*sizeof(char), (char*)filename);
378 
379  DebugM(3, "calling GridforceFullBaseGrid::pack\n" << endi);
380 
382 }
#define DebugM(x, y)
Definition: Debug.h:59
char filename[NAMD_FILENAME_BUFFER_SIZE]
virtual void pack(MOStream *msg) const
MOStream * put(char data)
Definition: MStream.h:112
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void GridforceFullMainGrid::reinitialize ( SimParameters simParams,
MGridforceParams mgridParams 
)
virtual

Implements GridforceGrid.

Definition at line 711 of file GridForceGrid.C.

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

712 {
713  DebugM(4, "reinitializing grid\n" << endi);
714  initialize(filename, simParams, mgridParams);
715 }
void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams, int border)
#define DebugM(x, y)
Definition: Debug.h:59
char filename[NAMD_FILENAME_BUFFER_SIZE]
infostream & endi(infostream &s)
Definition: InfoStream.C:38
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.

755 {
756  DebugM(4, "set_all_gridvals called\n" << endi);
757 
758  long int sz_calc = 0;
759  sz_calc += size;
760  for (int i = 0; i < totalGrids-1; i++) {
761  sz_calc += subgrids_flat[i]->size;
762  }
763  CmiAssert(sz == sz_calc);
764 
765  long int idx = 0;
766  for (long int i = 0; i < size; i++) {
767  DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
768  grid[i] = all_gridvals[idx++];
769  }
770  for (int j = 0; j < totalGrids-1; j++) {
771  for (long int i = 0; i < subgrids_flat[j]->size; i++) {
772  DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
773  subgrids_flat[j]->grid[i] = all_gridvals[idx++];
774  }
775  }
776  CmiAssert(idx == sz);
777 
778  DebugM(4, "set_all_gridvals finished\n" << endi);
779 }
#define DebugM(x, y)
Definition: Debug.h:59
GridforceFullSubGrid ** subgrids_flat
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void GridforceFullMainGrid::set_scale ( Vector  s)
inlinevirtual

Implements GridforceGrid.

Definition at line 219 of file GridForceGrid.h.

References GridforceFullBaseGrid::scale.

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

386 {
387  DebugM(4, "Unpacking maingrid\n" << endi);
388 
389 // msg->get(3*sizeof(float), (char*)pad_p);
390 // msg->get(3*sizeof(float), (char*)pad_n);
391  msg->get(totalGrids);
392  msg->get(mygridnum);
393  msg->get(129*sizeof(char), (char*)filename);
394 
396 
397  DebugM(4, "size = " << size << "\n");
398  DebugM(4, "numSubgrids = " << numSubgrids << "\n");
399  DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
400  DebugM(4, "generation = " << generation << "\n" << endi);
401  DebugM(4, "filename = " << filename << "\n" << endi);
402 
404 }
#define DebugM(x, y)
Definition: Debug.h:59
MIStream * get(char &data)
Definition: MStream.h:29
char filename[NAMD_FILENAME_BUFFER_SIZE]
virtual void unpack(MIStream *msg)
void buildSubgridsFlat(void)
infostream & endi(infostream &s)
Definition: InfoStream.C:38

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
int GridforceFullMainGrid::totalGrids
protected

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