NAMD
LJTable.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "LJTable.h"
9 #include "Node.h"
10 #include "Parameters.h"
11 // #define DEBUGM
12 #include "Debug.h"
13 #include "Molecule.h"
14 
15 //----------------------------------------------------------------------
17 {
18  Bool soluteScalingOn = Node::Object()->simParameters->soluteScalingOn;
19 
20  if (!soluteScalingOn) {
21  table_dim = Node::Object()->parameters->get_num_vdw_params();
22  } else {
23  int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
24  table_dim = ss_dim + Node::Object()->parameters->get_num_vdw_params();
25  }
26  table_alloc = new char[2*table_dim*table_dim*sizeof(TableEntry) + 31];
27  char *table_align = table_alloc;
28  while ( (long)table_align % 32 ) table_align++;
29  table = (TableEntry *) table_align;
30 
31  for (register int i=0; i < table_dim; i++)
32  for (register int j=i; j < table_dim; j++)
33  {
34  TableEntry *curij = &(table[2*(i*table_dim+j)]);
35  TableEntry *curji = &(table[2*(j*table_dim+i)]);
36  compute_vdw_params(i,j,curij,curij+1);
37 
38  // Copy to transpose entry
39  *curji = *curij;
40  *(curji + 1) = *(curij + 1);
41  }
42 
43 }
44 
45 //----------------------------------------------------------------------
47 {
48  delete [] table_alloc;
49 }
50 
51 //----------------------------------------------------------------------
52 void LJTable::compute_vdw_params(int i, int j,
53  LJTable::TableEntry *cur,
54  LJTable::TableEntry *cur_scaled)
55 {
56  Parameters *params = Node::Object()->parameters;
58  int useGeom = simParams->vdwGeometricSigma;
59  Bool tabulatedEnergies = simParams->tabulatedEnergies;
60  Bool soluteScalingOn = simParams->soluteScalingOn;
61  BigReal soluteScalingFactor = simParams->soluteScalingFactorVdw;
62  unsigned int table_dim_org = params->get_num_vdw_params();
63  int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
64  Real A, B, A14, B14;
65  int K = -1;
66  int *ss_vdw_type = Node::Object()->molecule->ss_vdw_type;
67  // BigReal sigma_max;
68  // We need the A and B parameters for the Van der Waals. These can
69  // be explicitly be specified for this pair or calculated from the
70  // sigma and epsilon values for the two atom types
71 // printf("Looking at interaction of %i with %i\n", i, j);
72  if ( tabulatedEnergies && params->get_table_pair_params(i,j,&K)) {
73 // printf("Making this interaction tabulated. %i %i %i\n", i, j, K);
74 #ifdef NAMD_CUDA
75  NAMD_die("Tabulated energies are not supported in CUDA-enabled NAMD");
76 #endif
77  if ( K < 0 ) NAMD_bug(
78  "LJTable::compute_vdw_params: energy table index is negative");
79 
80  cur->A = -1 - K;
81  cur->B = 0;
82  cur_scaled->A = -1 - K;
83  cur_scaled->B = 0;
84  }
85  else if (params->get_vdw_pair_params(i,j, &A, &B, &A14, &B14))
86  {
87  cur->A = A;
88  cur->B = B;
89  cur_scaled->A = A14;
90  cur_scaled->B = B14;
91 
92  if ( tabulatedEnergies && ( cur->A < 0 || cur_scaled->A < 0 ) )
93  NAMD_die("LJ A is negative with tabulatedEnergies enabled");
94 
95  // BigReal sigma_ij, sigma_ij14;
96 
97  // if ((B == 0) || (A/B < 0)) sigma_ij = 0;
98  // else sigma_ij = pow((BigReal)(A/B),(BigReal)(1./6.));
99 
100  // if ((B14 == 0) || (A14/B14 < 0)) sigma_ij14 = 0;
101  // else sigma_ij14 = pow((BigReal)(A14/B14),(BigReal)(1./6.));
102 
103  // sigma_max = ( sigma_ij > sigma_ij14 ? sigma_ij : sigma_ij14 );
104  }
105  else
106  {
107  // We didn't find explicit parameters for this pair. So instead,
108  // get the parameters for each atom type separately and use them
109  // to calculate the values we need
110  Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
111  Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
112 
113  if (!soluteScalingOn) {
114 
115  params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
116  &epsilon_i14,i);
117  params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14,
118  &epsilon_j14,j);
119  } else {
120  int i_type = (i >= table_dim_org)? ss_vdw_type[i-table_dim_org]:i;
121  int j_type = (j >= table_dim_org)? ss_vdw_type[j-table_dim_org]:j;
122  params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
123  &epsilon_i14,i_type);
124  params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14,
125  &epsilon_j14,j_type);
126  }
127  BigReal sigma_ij =
128  useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
129  BigReal sigma_ij14 =
130  useGeom ? sqrt(sigma_i14*sigma_j14) : 0.5 * (sigma_i14+sigma_j14);
131  BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
132  BigReal epsilon_ij14 = sqrt(epsilon_i14*epsilon_j14);
133 
134  // sigma_max = ( sigma_ij > sigma_ij14 ? sigma_ij : sigma_ij14 );
135 
136  // Calculate sigma^6
137  sigma_ij *= sigma_ij*sigma_ij;
138  sigma_ij *= sigma_ij;
139  sigma_ij14 *= sigma_ij14*sigma_ij14;
140  sigma_ij14 *= sigma_ij14;
141 
142  // Calculate LJ constants A & B
143  cur->B = 4.0 * sigma_ij * epsilon_ij;
144  cur->A = cur->B * sigma_ij;
145  cur_scaled->B = 4.0 * sigma_ij14 * epsilon_ij14;
146  cur_scaled->A = cur_scaled->B * sigma_ij14;
147 
148  if (soluteScalingOn) {
149  if (i >= table_dim_org && i < (table_dim_org+ss_dim) && j < table_dim_org) {
150  cur->A *= sqrt(soluteScalingFactor);
151  cur->B *= sqrt(soluteScalingFactor);
152  cur_scaled->A *= sqrt(soluteScalingFactor);
153  cur_scaled->B *= sqrt(soluteScalingFactor);
154  }
155  if (i < table_dim_org && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
156  cur->A *= sqrt(soluteScalingFactor);
157  cur->B *= sqrt(soluteScalingFactor);
158  cur_scaled->A *= sqrt(soluteScalingFactor);
159  cur_scaled->B *= sqrt(soluteScalingFactor);
160  }
161  if (i >=table_dim_org && i < (table_dim_org+ss_dim) && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
162  cur->A *= soluteScalingFactor;
163  cur->B *= soluteScalingFactor;
164  cur_scaled->A *= soluteScalingFactor;
165  cur_scaled->B *= soluteScalingFactor;
166  }
167  }
168 
169  if ( tabulatedEnergies && ( cur->A < 0 || cur_scaled->A < 0 ) )
170  NAMD_die("LJ A is negative with tabulatedEnergies enabled");
171  }
172  // Calculate exclcut2
173  // cur_scaled->exclcut2 = cur->exclcut2 = 0.64 * sigma_max * sigma_max;
174 
175 }
176 
static Node * Object()
Definition: Node.h:86
Bool vdwGeometricSigma
const BigReal A
SimParameters * simParameters
Definition: Node.h:178
float Real
Definition: common.h:109
void NAMD_bug(const char *err_msg)
Definition: common.C:123
int Bool
Definition: common.h:133
~LJTable(void)
Definition: LJTable.C:46
int get_vdw_pair_params(Index ind1, Index ind2, Real *, Real *, Real *, Real *)
Definition: Parameters.C:3680
void NAMD_die(const char *err_msg)
Definition: common.C:83
BigReal soluteScalingFactorVdw
int get_num_vdw_params(void)
Definition: Parameters.h:531
LJTable(void)
Definition: LJTable.C:16
int get_table_pair_params(Index, Index, int *)
Definition: Parameters.C:3605
Parameters * parameters
Definition: Node.h:177
#define simParams
Definition: Output.C:127
const BigReal B
int * ss_vdw_type
Definition: Molecule.h:457
Bool tabulatedEnergies
void get_vdw_params(Real *sigma, Real *epsilon, Real *sigma14, Real *epsilon14, Index index)
Definition: Parameters.h:497
Molecule * molecule
Definition: Node.h:176
int ss_num_vdw_params
Definition: Molecule.h:456
double BigReal
Definition: common.h:114