NAMD
ComputeGridForce.C
Go to the documentation of this file.
1 
7 #include "ComputeGridForce.h"
8 #include "GridForceGrid.h"
9 #include "Node.h"
10 
11 #define MIN_DEBUG_LEVEL 2
12 //#define DEBUGM
13 #include "Debug.h"
14 
15 #include "GridForceGrid.inl"
16 #include "MGridforceParams.h"
17 
18 //#define GF_FORCE_OUTPUT
19 //#define GF_FORCE_OUTPUT_FREQ 100
20 #define GF_OVERLAPCHECK_FREQ 1000
21 
22 
24  : ComputeHomePatch(c,pid)
25 {
26 
28 
29 }
30 /* END OF FUNCTION ComputeGridForce */
31 
32 
34 {
35  delete reduction;
36 }
37 /* END OF FUNCTION ~ComputeGridForce */
38 
39 template <class T> void ComputeGridForce::do_calc(T *grid, int gridnum, FullAtom *p, int numAtoms, Molecule *mol, Force *forces, BigReal &energy, Force &extForce, Tensor &extVirial)
40 {
41  Real scale; // Scaling factor
42  Charge charge; // Charge
43  Vector dV;
44  float V;
45 
46  Vector gfScale = grid->get_scale();
47  DebugM(3, "doCalc()\n" << endi);
48 
49  // Loop through and check each atom
50  for (int i = 0; i < numAtoms; i++) {
51  if (mol->is_atom_gridforced(p[i].id, gridnum))
52  {
53  DebugM(1, "Atom " << p[i].id << " is gridforced\n" << endi);
54 
55  mol->get_gridfrc_params(scale, charge, p[i].id, gridnum);
56 
57  // Wrap coordinates using grid center
58  Position pos = grid->wrap_position(p[i].position, homePatch->lattice);
59  DebugM(1, "pos = " << pos << "\n" << endi);
60 
61  // Here's where the action happens
62  int err = grid->compute_VdV(pos, V, dV);
63 
64  if (err) {
65  DebugM(2, "V = 0\n" << endi);
66  DebugM(2, "dV = 0 0 0\n" << endi);
67  continue; // This means the current atom is outside the potential
68  }
69 
70  //Force force = scale * Tensor::diagonal(gfScale) * (-charge * dV);
71  Force force = -charge * scale * Vector(gfScale.x * dV.x, gfScale.y * dV.y, gfScale.z * dV.z);
72 
73 #ifdef DEBUGM
74  DebugM(2, "scale = " << scale << " gfScale = " << gfScale << " charge = " << charge << "\n" << endi);
75 
76  DebugM(2, "V = " << V << "\n" << endi);
77  DebugM(2, "dV = " << dV << "\n" << endi);
78  DebugM(2, "grid = " << gridnum << " force = " << force << " pos = " << pos << " V = " << V << " dV = " << dV << " step = " << homePatch->flags.step << " index = " << p[i].id << "\n" << endi);
79 
80  DebugM(1, "transform = " << (int)p[i].transform.i << " "
81  << (int)p[i].transform.j << " " << (int)p[i].transform.k << "\n" << endi);
82 
83  if (V != V) {
84  iout << iWARN << "V is NaN!\natomid = " << p[i].id << " loc = " << p[i].position << " V = " << V << "\n" << endi;
85  }
86 #endif
87 
88  forces[i] += force;
89  extForce += force;
90  Position vpos = homePatch->lattice.reverse_transform(p[i].position, p[i].transform);
91 
92  //energy -= force * (vpos - homePatch->lattice.origin());
93  if (gfScale.x == gfScale.y && gfScale.x == gfScale.z)
94  {
95  // only makes sense when scaling is isotropic
96  energy += scale * gfScale.x * (charge * V);
97 
98  // add something when we're off the grid? I'm thinking no
99  }
100  extVirial += outer(force,vpos);
101  }
102  }
103  DebugM(3, "doCalc() done\n" << endi);
104 }
105 
107 {
109  Molecule *mol = Node::Object()->molecule;
110 
111  Force *forces = r->f[Results::normal];
112  BigReal energy = 0;
113  Force extForce = 0.;
114  Tensor extVirial;
115 
116  int numAtoms = homePatch->getNumAtoms();
117 
118  if ( mol->numGridforceGrids < 1 ) NAMD_bug("No grids loaded in ComputeGridForce::doForce()");
119 
120  DebugM(3, "doForce()\n" << endi);
121 
122  for (int gridnum = 0; gridnum < mol->numGridforceGrids; gridnum++) {
123  GridforceGrid *grid = mol->get_gridfrc_grid(gridnum);
124 
125  const Vector gfScale = grid->get_scale();
126  if ((gfScale.x == 0.0) && (gfScale.y == 0.0) && (gfScale.z == 0.0)) {
127  DebugM(3, "Skipping grid index " << gridnum << "\n" << endi);
128  continue;
129  }
130 
132  // only check on node 0 and every GF_OVERLAPCHECK_FREQ steps
133  if (simParams->langevinPistonOn || simParams->berendsenPressureOn) {
134  // check for grid overlap if pressure control is on
135  // not needed without pressure control, since the check is also performed on startup
136  if (!grid->fits_lattice(homePatch->lattice)) {
137  char errmsg[512];
138  if (grid->get_checksize()) {
139  sprintf(errmsg, "Warning: Periodic cell basis too small for Gridforce grid %d. Set gridforcechecksize off in configuration file to ignore.\n", gridnum);
140  NAMD_die(errmsg);
141  }
142  }
143  }
144  }
145 
146  if (homePatch->flags.step % 100 == 1) {
147  Position center = grid->get_center();
148  DebugM(3, "center = " << center << "\n" << endi);
149  DebugM(3, "e = " << grid->get_e() << "\n" << endi);
150  }
151 
154  do_calc(g, gridnum, p, numAtoms, mol, forces, energy, extForce, extVirial);
155  } else if (grid->get_grid_type() == GridforceGrid::GridforceGridTypeLite) {
157  do_calc(g, gridnum, p, numAtoms, mol, forces, energy, extForce, extVirial);
158  }
159  }
161  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
162  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
163  reduction->submit();
164  DebugM(3, "doForce() done\n" << endi);
165 }
166 /* END OF FUNCTION force */
static Node * Object()
Definition: Node.h:86
Bool berendsenPressureOn
signed char j
Definition: NamdTypes.h:38
int ComputeID
Definition: NamdTypes.h:183
Lattice & lattice
Definition: Patch.h:126
Definition: Vector.h:64
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
SimParameters * simParameters
Definition: Node.h:178
ComputeGridForce(ComputeID c, PatchID pid)
static __thread float4 * forces
float Real
Definition: common.h:109
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
#define GF_OVERLAPCHECK_FREQ
int numGridforceGrids
Definition: Molecule.h:590
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:108
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
void doForce(FullAtom *p, Results *r)
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define iout
Definition: InfoStream.h:87
void get_gridfrc_params(Real &k, Charge &q, int atomnum, int gridnum) const
Definition: Molecule.h:1270
Flags flags
Definition: Patch.h:127
signed char k
Definition: NamdTypes.h:38
bool fits_lattice(const Lattice &lattice)
Definition: GridForceGrid.C:84
void NAMD_bug(const char *err_msg)
Definition: common.C:123
virtual Position get_center(void) const =0
Force * f[maxNumForces]
Definition: PatchTypes.h:67
virtual ~ComputeGridForce()
SubmitReduction * reduction
BigReal x
Definition: Vector.h:66
int PatchID
Definition: NamdTypes.h:182
GridforceGrid * get_gridfrc_grid(int gridnum) const
Definition: Molecule.h:1276
void NAMD_die(const char *err_msg)
Definition: common.C:83
virtual Bool get_checksize(void) const =0
Tensor outer(const Vector &v1, const Vector &v2)
Definition: Tensor.h:241
virtual Vector get_scale(void) const =0
#define simParams
Definition: Output.C:127
Position reverse_transform(Position data, const Transform &t) const
Definition: Lattice.h:138
Definition: Tensor.h:15
BigReal y
Definition: Vector.h:66
int getNumAtoms()
Definition: Patch.h:105
virtual Tensor get_e(void) const =0
#define ADD_VECTOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:27
void do_calc(T *grid, int gridnum, FullAtom *p, int numAtoms, Molecule *mol, Force *forces, BigReal &energy, Force &extForce, Tensor &extVirial)
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
GridforceGridType get_grid_type(void)
Definition: GridForceGrid.h:65
void submit(void)
Definition: ReductionMgr.h:323
infostream & endi(infostream &s)
Definition: InfoStream.C:38
Molecule * molecule
Definition: Node.h:176
float Charge
Definition: NamdTypes.h:32
HomePatch * homePatch
double BigReal
Definition: common.h:114
Transform transform
Definition: NamdTypes.h:116
Bool is_atom_gridforced(int atomnum, int gridnum) const
Definition: Molecule.h:1181
int step
Definition: PatchTypes.h:16