NAMD
ComputeNonbondedCUDAExcl.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
9 #include "Molecule.h"
10 #include "Parameters.h"
11 #include "LJTable.h"
12 #include "Node.h"
13 #include "ReductionMgr.h"
14 #include "Lattice.h"
15 #include "PressureProfile.h"
16 #include "Debug.h"
17 
18 
19 // static initialization
24 
26  (Molecule* mol, int* count, int32*** byatom, Exclusion** structarray)
27 {
28 #ifdef MEM_OPT_VERSION
29  NAMD_die("Should not be called in ExclElem::getMoleculePointers in memory optimized version!");
30 #else
31  *count = mol->numExclusions;
32  *byatom = mol->exclusionsByAtom;
33  *structarray = mol->exclusions;
34 #endif
35 }
36 
37 void ExclElem::getParameterPointers(Parameters *p, const int **v) {
38  *v = 0;
39 }
40 
41 void ExclElem::computeForce(ExclElem *tuples, int ntuple, BigReal *reduction,
42  BigReal *pressureProfileData)
43 {
44  const Lattice & lattice = tuples[0].p[0]->p->lattice;
45  const Flags &flags = tuples[0].p[0]->p->flags;
46  if ( ! flags.doNonbonded ) return;
47  const int doFull = flags.doFullElectrostatics;
48  const int doEnergy = flags.doEnergy;
49 
50  for ( int ituple=0; ituple<ntuple; ++ituple ) {
51  const ExclElem &tup = tuples[ituple];
52  enum { size = 2 };
53  const AtomID (&atomID)[size](tup.atomID);
54  const int (&localIndex)[size](tup.localIndex);
55  TuplePatchElem * const(&p)[size](tup.p);
56  const Real (&scale)(tup.scale);
57  const int (&modified)(tup.modified);
58 
59  const CompAtom &p_i = p[0]->x[localIndex[0]];
60  const CompAtom &p_j = p[1]->x[localIndex[1]];
61 
62  // compute vectors between atoms and their distances
63  const Vector r12 = lattice.delta(p_i.position, p_j.position);
64  BigReal r2 = r12.length2();
65 
66  if ( r2 > cutoff2 ) continue;
67 
68  if ( modified && r2 < 1.0 ) r2 = 1.0; // match CUDA interpolation
69 
70  r2 += r2_delta;
71 
72  union { double f; int64 i; } r2i;
73  r2i.f = r2;
74  const int r2_delta_expc = 64 * (r2_delta_exp - 1023);
75  int table_i = (r2i.i >> (32+14)) + r2_delta_expc; // table_i >= 0
76 
77  const BigReal* const table_four_i = table_noshort + 16*table_i;
78 
79  BigReal diffa = r2 - r2_table[table_i];
80 
81  BigReal fast_a = 0., fast_b = 0., fast_c = 0., fast_d = 0.;
82  BigReal slow_a, slow_b, slow_c, slow_d;
83 
84  if ( modified ) {
85 
87  ljTable->table_row(p_i.vdwType) + 2 * p_j.vdwType;
88 
89  // modified - normal = correction
90  const BigReal A = scaling * ( (lj_pars+1)->A - lj_pars->A );
91  const BigReal B = scaling * ( (lj_pars+1)->B - lj_pars->B );
92 
93  BigReal vdw_d = A * table_four_i[0] - B * table_four_i[4];
94  BigReal vdw_c = A * table_four_i[1] - B * table_four_i[5];
95  BigReal vdw_b = A * table_four_i[2] - B * table_four_i[6];
96  BigReal vdw_a = A * table_four_i[3] - B * table_four_i[7];
97 
98  const BigReal kqq = (1.0 - scale14) *
99  COULOMB * p_i.charge * p_j.charge * scaling * dielectric_1;
100 
101  fast_a = kqq * fast_table[4*table_i+0]; // not used!
102  fast_b = 2. * kqq * fast_table[4*table_i+1];
103  fast_c = 4. * kqq * fast_table[4*table_i+2];
104  fast_d = 6. * kqq * fast_table[4*table_i+3];
105 
106  if ( doFull ) {
107  slow_a = kqq * slow_table[4*table_i+3]; // not used!
108  slow_b = 2. * kqq * slow_table[4*table_i+2];
109  slow_c = 4. * kqq * slow_table[4*table_i+1];
110  slow_d = 6. * kqq * slow_table[4*table_i+0];
111  }
112 
113  if ( doEnergy ) {
114  reduction[vdwEnergyIndex] -=
115  ( ( diffa * (1./6.)*vdw_d + 0.25*vdw_c ) * diffa
116  + 0.5*vdw_b ) * diffa + vdw_a;
117  reduction[electEnergyIndex] -=
118  ( ( diffa * (1./6.)*fast_d + 0.25*fast_c ) * diffa
119  + 0.5*fast_b ) * diffa + fast_a;
120  if ( doFull ) {
121  reduction[fullElectEnergyIndex] -=
122  ( ( diffa * (1./6.)*slow_d + 0.25*slow_c ) * diffa
123  + 0.5*slow_b ) * diffa + slow_a;
124  }
125  }
126 
127  fast_a += vdw_a;
128  fast_b += vdw_b;
129  fast_c += vdw_c;
130  fast_d += vdw_d;
131 
132  } else if ( doFull ) { // full exclusion
133 
134  const BigReal kqq =
135  COULOMB * p_i.charge * p_j.charge * scaling * dielectric_1;
136 
137  slow_d = kqq * ( table_four_i[8] - table_four_i[12] );
138  slow_c = kqq * ( table_four_i[9] - table_four_i[13] );
139  slow_b = kqq * ( table_four_i[10] - table_four_i[14] );
140  slow_a = kqq * ( table_four_i[11] - table_four_i[15] ); // not used!
141 
142  if ( doEnergy ) {
143  reduction[fullElectEnergyIndex] -=
144  ( ( diffa * (1./6.)*slow_d + 0.25*slow_c ) * diffa
145  + 0.5*slow_b ) * diffa + slow_a;
146  }
147  }
148 
149  register BigReal fast_dir =
150  (diffa * fast_d + fast_c) * diffa + fast_b;
151 
152  const Force f12 = fast_dir * r12;
153 
154  // Now add the forces to each force vector
155  p[0]->r->f[Results::nbond][localIndex[0]] += f12;
156  p[1]->r->f[Results::nbond][localIndex[1]] -= f12;
157 
158  // reduction[nonbondedEnergyIndex] += energy;
159  reduction[virialIndex_XX] += f12.x * r12.x;
160  reduction[virialIndex_XY] += f12.x * r12.y;
161  reduction[virialIndex_XZ] += f12.x * r12.z;
162  reduction[virialIndex_YX] += f12.y * r12.x;
163  reduction[virialIndex_YY] += f12.y * r12.y;
164  reduction[virialIndex_YZ] += f12.y * r12.z;
165  reduction[virialIndex_ZX] += f12.z * r12.x;
166  reduction[virialIndex_ZY] += f12.z * r12.y;
167  reduction[virialIndex_ZZ] += f12.z * r12.z;
168 
169  if ( doFull ) {
170  register BigReal slow_dir =
171  (diffa * slow_d + slow_c) * diffa + slow_b;
172 
173  const Force slow_f12 = slow_dir * r12;
174 
175  p[0]->r->f[Results::slow][localIndex[0]] += slow_f12;
176  p[1]->r->f[Results::slow][localIndex[1]] -= slow_f12;
177 
178  // reduction[nonbondedEnergyIndex] += energy;
179  reduction[slowVirialIndex_XX] += slow_f12.x * r12.x;
180  reduction[slowVirialIndex_XY] += slow_f12.x * r12.y;
181  reduction[slowVirialIndex_XZ] += slow_f12.x * r12.z;
182  reduction[slowVirialIndex_YX] += slow_f12.y * r12.x;
183  reduction[slowVirialIndex_YY] += slow_f12.y * r12.y;
184  reduction[slowVirialIndex_YZ] += slow_f12.y * r12.z;
185  reduction[slowVirialIndex_ZX] += slow_f12.z * r12.x;
186  reduction[slowVirialIndex_ZY] += slow_f12.z * r12.y;
187  reduction[slowVirialIndex_ZZ] += slow_f12.z * r12.z;
188  }
189 
190  }
191 }
192 
194 {
195  reduction->item(REDUCTION_ELECT_ENERGY) += data[electEnergyIndex];
196  reduction->item(REDUCTION_LJ_ENERGY) += data[vdwEnergyIndex];
198  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
199  ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,slowVirialIndex);
200 }
201 
static BigReal * fast_table
static int pressureProfileAtomTypes
register BigReal fast_dir
short int32
Definition: dumpdcd.c:24
int AtomID
Definition: NamdTypes.h:29
Lattice & lattice
Definition: Patch.h:126
#define ADD_TENSOR(R, RL, D, DL)
Definition: ReductionMgr.h:32
const BigReal A
Definition: Vector.h:64
static void submitReductionData(BigReal *, SubmitReduction *)
static void computeForce(ExclElem *, int, BigReal *, BigReal *)
float Real
Definition: common.h:109
#define COULOMB
Definition: common.h:46
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define table_four_i
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
if(ComputeNonbondedUtil::goMethod==2)
static void getParameterPointers(Parameters *, const int **)
#define p_j
BigReal vdw_b
#define lj_pars
Flags flags
Definition: Patch.h:127
BigReal vdw_c
Charge charge
Definition: NamdTypes.h:54
static BigReal * table_noshort
int doEnergy
Definition: PatchTypes.h:20
int doFullElectrostatics
Definition: PatchTypes.h:23
Force * f[maxNumForces]
Definition: PatchTypes.h:67
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
BigReal x
Definition: Vector.h:66
int doNonbonded
Definition: PatchTypes.h:22
void NAMD_die(const char *err_msg)
Definition: common.C:85
TuplePatchElem * p[size]
static int pressureProfileSlabs
AtomID atomID[size]
static BigReal * slow_table
BigReal length2(void) const
Definition: Vector.h:173
long long int64
Definition: common.h:34
short vdwType
Definition: NamdTypes.h:55
BigReal vdw_d
BigReal y
Definition: Vector.h:66
static const LJTable * ljTable
const BigReal B
BigReal vdw_a
static BigReal pressureProfileThickness
static BigReal pressureProfileMin
const TableEntry * table_row(unsigned int i) const
Definition: LJTable.h:31
static void getMoleculePointers(Molecule *, int *, int32 ***, Exclusion **)
double BigReal
Definition: common.h:114
int numExclusions
Definition: Molecule.h:571