LJTable.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "LJTable.h"
00009 #include "Node.h"
00010 #include "Parameters.h"
00011 // #define DEBUGM
00012 #include "Debug.h"
00013 #include "Molecule.h"
00014 
00015 //----------------------------------------------------------------------  
00016 LJTable::LJTable()
00017 {
00018   Bool soluteScalingOn = Node::Object()->simParameters->soluteScalingOn;
00019 
00020   if (!soluteScalingOn) {
00021   table_dim = Node::Object()->parameters->get_num_vdw_params();
00022   } else {
00023   int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
00024   table_dim = ss_dim + Node::Object()->parameters->get_num_vdw_params();
00025   }
00026   table_alloc = new char[2*table_dim*table_dim*sizeof(TableEntry) + 31];
00027   char *table_align = table_alloc;
00028   while ( (long)table_align % 32 ) table_align++;
00029   table = (TableEntry *) table_align;
00030 
00031   for (register int i=0; i < table_dim; i++)
00032     for (register int j=i; j < table_dim; j++)
00033     {
00034       TableEntry *curij = &(table[2*(i*table_dim+j)]);
00035       TableEntry *curji = &(table[2*(j*table_dim+i)]);
00036       compute_vdw_params(i,j,curij,curij+1);
00037 
00038       // Copy to transpose entry
00039       *curji = *curij;
00040       *(curji + 1) = *(curij + 1);
00041     }
00042 
00043 }
00044 
00045 //----------------------------------------------------------------------  
00046 LJTable::~LJTable()
00047 {
00048   delete [] table_alloc;
00049 }
00050 
00051 //----------------------------------------------------------------------
00052 void LJTable::compute_vdw_params(int i, int j,
00053                                  LJTable::TableEntry *cur, 
00054                                  LJTable::TableEntry *cur_scaled)
00055 {
00056   Parameters *params = Node::Object()->parameters;
00057   SimParameters *simParams = Node::Object()->simParameters;
00058   int useGeom = simParams->vdwGeometricSigma;
00059   Bool tabulatedEnergies = simParams->tabulatedEnergies;
00060   Bool soluteScalingOn = simParams->soluteScalingOn;
00061   BigReal soluteScalingFactor = simParams->soluteScalingFactorVdw;
00062   unsigned int table_dim_org = params->get_num_vdw_params();
00063   int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
00064   Real A, B, A14, B14;
00065   int K = -1;
00066   int *ss_vdw_type = Node::Object()->molecule->ss_vdw_type;
00067   // BigReal sigma_max;
00068   //  We need the A and B parameters for the Van der Waals.  These can
00069   //  be explicitly be specified for this pair or calculated from the
00070   //  sigma and epsilon values for the two atom types
00071 //  printf("Looking at interaction of  %i with %i\n", i, j);
00072   if ( tabulatedEnergies && params->get_table_pair_params(i,j,&K)) {
00073 //    printf("Making this interaction tabulated. %i %i %i\n", i, j, K);
00074 #ifdef NAMD_CUDA
00075     NAMD_die("Tabulated energies are not supported in CUDA-enabled NAMD");
00076 #endif
00077     if ( K < 0 ) NAMD_bug(
00078         "LJTable::compute_vdw_params: energy table index is negative");
00079 
00080     cur->A = -1 - K;
00081     cur->B = 0;
00082     cur_scaled->A = -1 - K;
00083     cur_scaled->B = 0;
00084   }
00085   else if (params->get_vdw_pair_params(i,j, &A, &B, &A14, &B14))
00086   {
00087     cur->A = A;
00088     cur->B = B;
00089     cur_scaled->A = A14;
00090     cur_scaled->B = B14;
00091 
00092     if ( tabulatedEnergies && ( cur->A < 0 || cur_scaled->A < 0 ) )
00093       NAMD_die("LJ A is negative with tabulatedEnergies enabled");
00094 
00095     // BigReal sigma_ij, sigma_ij14;
00096 
00097     // if ((B == 0) || (A/B < 0)) sigma_ij = 0;
00098     // else sigma_ij = pow((BigReal)(A/B),(BigReal)(1./6.));
00099 
00100     // if ((B14 == 0) || (A14/B14 < 0)) sigma_ij14 = 0;
00101     // else sigma_ij14 = pow((BigReal)(A14/B14),(BigReal)(1./6.));
00102 
00103     // sigma_max = ( sigma_ij > sigma_ij14 ? sigma_ij : sigma_ij14 );
00104   }
00105   else
00106   {
00107     //  We didn't find explicit parameters for this pair. So instead,
00108     //  get the parameters for each atom type separately and use them
00109     //  to calculate the values we need
00110     Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
00111     Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;
00112 
00113    if (!soluteScalingOn) {
00114 
00115     params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
00116                                        &epsilon_i14,i);
00117     params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14, 
00118                                        &epsilon_j14,j);
00119    } else {
00120    int i_type = (i >= table_dim_org)? ss_vdw_type[i-table_dim_org]:i;
00121    int j_type = (j >= table_dim_org)? ss_vdw_type[j-table_dim_org]:j;
00122    params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
00123                                        &epsilon_i14,i_type);
00124    params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14,
00125                                        &epsilon_j14,j_type);
00126    }    
00127     BigReal sigma_ij =
00128        useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
00129     BigReal sigma_ij14 =
00130        useGeom ? sqrt(sigma_i14*sigma_j14) : 0.5 * (sigma_i14+sigma_j14);
00131     BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
00132     BigReal epsilon_ij14 = sqrt(epsilon_i14*epsilon_j14);
00133 
00134     // sigma_max = ( sigma_ij > sigma_ij14 ? sigma_ij : sigma_ij14 );
00135 
00136     //  Calculate sigma^6
00137     sigma_ij *= sigma_ij*sigma_ij;
00138     sigma_ij *= sigma_ij;
00139     sigma_ij14 *= sigma_ij14*sigma_ij14;
00140     sigma_ij14 *= sigma_ij14;
00141     
00142     //  Calculate LJ constants A & B
00143     cur->B = 4.0 * sigma_ij * epsilon_ij;
00144     cur->A = cur->B * sigma_ij;
00145     cur_scaled->B = 4.0 * sigma_ij14 * epsilon_ij14;
00146     cur_scaled->A = cur_scaled->B * sigma_ij14;
00147 
00148     if (soluteScalingOn) {
00149      if (i >= table_dim_org && i < (table_dim_org+ss_dim) && j < table_dim_org) {
00150         cur->A *= sqrt(soluteScalingFactor);
00151         cur->B *= sqrt(soluteScalingFactor);
00152         cur_scaled->A *= sqrt(soluteScalingFactor);
00153         cur_scaled->B *= sqrt(soluteScalingFactor);
00154      }
00155      if (i < table_dim_org && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
00156         cur->A *= sqrt(soluteScalingFactor);
00157         cur->B *= sqrt(soluteScalingFactor);
00158         cur_scaled->A *= sqrt(soluteScalingFactor);
00159         cur_scaled->B *= sqrt(soluteScalingFactor);
00160       }
00161      if (i >=table_dim_org && i < (table_dim_org+ss_dim) && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
00162         cur->A *= soluteScalingFactor;
00163         cur->B *= soluteScalingFactor;
00164         cur_scaled->A *= soluteScalingFactor;
00165         cur_scaled->B *= soluteScalingFactor;
00166      }
00167     }
00168 
00169     if ( tabulatedEnergies && ( cur->A < 0 || cur_scaled->A < 0 ) )
00170       NAMD_die("LJ A is negative with tabulatedEnergies enabled");
00171   }
00172   //  Calculate exclcut2
00173   // cur_scaled->exclcut2 = cur->exclcut2 = 0.64 * sigma_max * sigma_max;
00174 
00175 }
00176 

Generated on Tue May 22 01:17:14 2018 for NAMD by  doxygen 1.4.7