Parameters Class Reference

#include <Parameters.h>

List of all members.

Public Member Functions

 Parameters ()
 Parameters (SimParameters *, StringList *f)
 Parameters (Ambertoppar *, BigReal)
void read_parm (Ambertoppar *, BigReal)
 Parameters (const GromacsTopFile *gf, Bool min)
 ~Parameters ()
char * atom_type_name (Index a)
void read_parameter_file (char *)
void read_charmm_parameter_file (char *)
void done_reading_files ()
void done_reading_structure ()
void assign_vdw_index (char *, Atom *)
void assign_bond_index (char *, char *, Bond *)
void assign_angle_index (char *, char *, char *, Angle *, int)
void assign_dihedral_index (char *, char *, char *, char *, Dihedral *, int, int)
void assign_improper_index (char *, char *, char *, char *, Improper *, int)
void assign_crossterm_index (char *, char *, char *, char *, char *, char *, char *, char *, Crossterm *)
void send_Parameters (MOStream *)
void receive_Parameters (MIStream *)
void get_bond_params (Real *k, Real *x0, Index index)
void get_angle_params (Real *k, Real *theta0, Real *k_ub, Real *r_ub, Index index)
int get_improper_multiplicity (Index index)
int get_dihedral_multiplicity (Index index)
void get_improper_params (Real *k, int *n, Real *delta, Index index, int mult)
void get_dihedral_params (Real *k, int *n, Real *delta, Index index, int mult)
void get_vdw_params (Real *sigma, Real *epsilon, Real *sigma14, Real *epsilon14, Index index)
int get_vdw_pair_params (Index ind1, Index ind2, Real *, Real *, Real *, Real *)
int get_num_vdw_params (void)
int get_table_pair_params (Index, Index, int *)
void print_bond_params ()
void print_angle_params ()
void print_dihedral_params ()
void print_improper_params ()
void print_crossterm_params ()
void print_vdw_params ()
void print_vdw_pair_params ()
void print_nbthole_pair_params ()
void print_param_summary ()
void read_ener_table (SimParameters *)
int get_int_table_type (char *)
int read_energy_type (FILE *, const int, BigReal *, const float, const float)
int read_energy_type_cubspline (FILE *, const int, BigReal *, const float, const float)
int read_energy_type_bothcubspline (FILE *, const int, BigReal *, const float, const float)

Public Attributes

BondValuebond_array
AngleValueangle_array
DihedralValuedihedral_array
ImproperValueimproper_array
CrosstermValuecrossterm_array
GromacsPairValuegromacsPair_array
VdwValuevdw_array
NbtholePairValuenbthole_array
int numenerentries
int rowsize
int columnsize
BigRealtable_ener
IndexedVdwPairvdw_pair_tree
IndexedNbtholePairnbthole_pair_tree
IndexedTablePairtab_pair_tree
int tablenumtypes
int NumBondParams
int NumAngleParams
int NumCosAngles
int NumDihedralParams
int NumImproperParams
int NumCrosstermParams
int NumGromacsPairParams
int NumVdwParams
int NumTableParams
int NumVdwParamsAssigned
int NumVdwPairParams
int NumNbtholePairParams
int NumTablePairParams


Detailed Description

Definition at line 214 of file Parameters.h.


Constructor & Destructor Documentation

Parameters::Parameters (  ) 

Definition at line 186 of file Parameters.C.

00186                        {
00187   initialize();
00188 }

Parameters::Parameters ( SimParameters ,
StringList f 
)

Definition at line 249 of file Parameters.C.

References done_reading_files(), f, FALSE, paraCharmm, paraXplor, read_charmm_parameter_file(), read_ener_table(), read_parameter_file(), and simParams.

00250 {
00251   initialize();
00252 
00254   if (simParams->paraTypeXplorOn)
00255   {
00256     paramType = paraXplor;
00257   }
00258   else if (simParams->paraTypeCharmmOn)
00259   {
00260     paramType = paraCharmm;
00261   }
00262   //****** END CHARMM/XPLOR type changes
00263   //Test for cos-based angles
00264   if (simParams->cosAngles) {
00265     cosAngles = true;
00266   } else {
00267     cosAngles = false;
00268   }
00269 
00270   if (simParams->tabulatedEnergies) {
00271           CkPrintf("Working on tables\n");
00272           read_ener_table(simParams);
00273   }
00274 
00275   //****** BEGIN CHARMM/XPLOR type changes
00276   /* Set up AllFilesRead flag to FALSE.  Once all of the files    */
00277   /* have been read in, then this will be set to true and the     */
00278   /* arrays of parameters will be set up        */
00279   AllFilesRead = FALSE;
00280 
00281   if (NULL != f) 
00282   {
00283     do
00284     {
00285       //****** BEGIN CHARMM/XPLOR type changes
00286       if (paramType == paraXplor)
00287       {
00288         read_parameter_file(f->data);
00289       }
00290       else if (paramType == paraCharmm)
00291       {
00292         read_charmm_parameter_file(f->data);
00293       }
00294       //****** END CHARMM/XPLOR type changes
00295       f = f->next;
00296     } while ( f != NULL );
00297 
00298     done_reading_files();
00299   }
00300 
00301 }

Parameters::Parameters ( Ambertoppar ,
BigReal   
)

Definition at line 6340 of file Parameters.C.

06341 {
06342   initialize();
06343 
06344   // Read in parm parameters
06345   read_parm(amber_data,vdw14);
06346 }

Parameters::Parameters ( const GromacsTopFile gf,
Bool  min 
)

Definition at line 6482 of file Parameters.C.

06483 {
06484   initialize();
06485 
06486   // Read in parm parameters
06487   read_parm(gf,min);
06488 }

Parameters::~Parameters (  ) 

Definition at line 313 of file Parameters.C.

References angle_array, bond_array, crossterm_array, dihedral_array, gromacsPair_array, improper_array, nbthole_pair_tree, ResizeArray< Elem >::resize(), ResizeArray< Elem >::size(), tab_pair_tree, vdw_array, and vdw_pair_tree.

00315 {
00316         if (atomTypeNames)
00317           delete [] atomTypeNames;
00318 
00319   if (bondp != NULL)
00320     free_bond_tree(bondp);
00321 
00322   if (anglep != NULL)
00323     free_angle_tree(anglep);
00324 
00325   if (dihedralp != NULL)
00326     free_dihedral_list(dihedralp);
00327 
00328   if (improperp != NULL)
00329     free_improper_list(improperp);
00330 
00331   if (crosstermp != NULL)
00332     free_crossterm_list(crosstermp);
00333 
00334   if (vdwp != NULL)
00335     free_vdw_tree(vdwp);
00336 
00337   if (vdw_pairp != NULL)
00338     free_vdw_pair_list();
00339 
00340   if (nbthole_pairp != NULL)
00341     free_nbthole_pair_list();
00342 
00343   if (bond_array != NULL)
00344     delete [] bond_array;
00345 
00346   if (angle_array != NULL)
00347     delete [] angle_array;
00348 
00349   if (dihedral_array != NULL)
00350     delete [] dihedral_array;
00351 
00352   if (improper_array != NULL)
00353     delete [] improper_array;
00354 
00355   if (crossterm_array != NULL)
00356     delete [] crossterm_array;
00357 
00358   // JLai
00359   if (gromacsPair_array != NULL)
00360     delete [] gromacsPair_array;
00361   // End of JLai
00362 
00363   if (vdw_array != NULL)
00364     delete [] vdw_array;
00365   
00366   if (tab_pair_tree != NULL)
00367     free_table_pair_tree(tab_pair_tree);
00368 
00369   if (vdw_pair_tree != NULL)
00370     free_vdw_pair_tree(vdw_pair_tree);
00371 
00372   if (nbthole_pair_tree != NULL)
00373     free_nbthole_pair_tree(nbthole_pair_tree);
00374 
00375   if (maxDihedralMults != NULL)
00376     delete [] maxDihedralMults;
00377 
00378   if (maxImproperMults != NULL)
00379     delete [] maxImproperMults;
00380 
00381   for( int i = 0; i < error_msgs.size(); ++i ) {
00382     delete [] error_msgs[i];
00383   }
00384   error_msgs.resize(0);
00385 }


Member Function Documentation

void Parameters::assign_angle_index ( char *  ,
char *  ,
char *  ,
Angle ,
int   
)

Definition at line 3849 of file Parameters.C.

References angle::angle_type, angle::atom1, angle_params::atom1name, angle::atom2, angle_params::atom2name, angle::atom3, angle_params::atom3name, angle_params::index, iout, iWARN(), angle_params::left, NAMD_die(), and angle_params::right.

03852 {
03853   struct angle_params *ptr;  //  Current position in tree
03854   int comp_val;      //  value from strcasecmp
03855   int found=0;      //  flag 1->found a match
03856   char tmp_name[15];    //  Temporary atom name
03857 
03858   /*  Check to make sure the files have all been read    */
03859   if (!AllFilesRead)
03860   {
03861     NAMD_die("Tried to assign angle index before all parameter files were read");
03862   }
03863 
03864   /*  We need atom1 < atom3.  If that was not what we were   */
03865   /*  passed, switch them            */
03866   if (strcasecmp(atom1, atom3) > 0)
03867   {
03868     strcpy(tmp_name, atom1);
03869     strcpy(atom1, atom3);
03870     strcpy(atom3, tmp_name);
03871   }
03872 
03873   /*  Start at the top            */
03874   ptr=anglep;
03875 
03876   /*  While we don't have a match and we haven't reached the  */
03877   /*  bottom of the tree, compare values        */
03878   while (!found && (ptr != NULL))
03879   {
03880     comp_val = strcasecmp(atom1, ptr->atom1name);
03881 
03882     if (comp_val == 0)
03883     {
03884       /*  Atom 1 matches, so compare atom 2    */
03885       comp_val = strcasecmp(atom2, ptr->atom2name);
03886       
03887       if (comp_val == 0)
03888       {
03889         /*  Atoms 1&2 match, try atom 3    */
03890         comp_val = strcasecmp(atom3, ptr->atom3name);
03891       }
03892     }
03893 
03894     if (comp_val == 0)
03895     {
03896       /*  Found a match        */
03897       found = 1;
03898       angle_ptr->angle_type = ptr->index;
03899     }
03900     else if (comp_val < 0)
03901     {
03902       /*  Go left          */
03903       ptr=ptr->left;
03904     }
03905     else
03906     {
03907       /*  Go right          */
03908       ptr=ptr->right;
03909     }
03910   }
03911 
03912   /*  Make sure we found a match          */
03913   if (!found)
03914   {
03915     char err_msg[128];
03916 
03917     sprintf(err_msg, "UNABLE TO FIND ANGLE PARAMETERS FOR %s %s %s (ATOMS %i %i %i)",
03918        atom1, atom2, atom3, angle_ptr->atom1+1, angle_ptr->atom2+1, angle_ptr->atom3+1);
03919 
03920     if ( notFoundIndex ) {
03921       angle_ptr->angle_type = notFoundIndex;
03922       iout << iWARN << err_msg << "\n" << endi;
03923       return;
03924     } else NAMD_die(err_msg);
03925   }
03926 
03927   return;
03928 }

void Parameters::assign_bond_index ( char *  ,
char *  ,
Bond  
)

Definition at line 3752 of file Parameters.C.

References bond::atom1, bond_params::atom1name, bond::atom2, bond_params::atom2name, bond::bond_type, bond_params::index, bond_params::left, NAMD_die(), and bond_params::right.

03754 {
03755   struct bond_params *ptr;  //  Current location in tree
03756   int found=0;      //  Flag 1-> found a match
03757   int cmp_code;      //  return code from strcasecmp
03758   char tmp_name[15];    //  Temporary atom name
03759 
03760   /*  Check to make sure the files have all been read    */
03761   if (!AllFilesRead)
03762   {
03763     NAMD_die("Tried to assign bond index before all parameter files were read");
03764   }
03765 
03766   /*  We need atom1 < atom2, so if that's not the way they        */
03767   /*  were passed, flip them          */
03768   if (strcasecmp(atom1, atom2) > 0)
03769   {
03770     strcpy(tmp_name, atom1);
03771     strcpy(atom1, atom2);
03772     strcpy(atom2, tmp_name);
03773   }
03774 
03775   /*  Start at the top            */
03776   ptr=bondp;
03777 
03778   /*  While we haven't found a match and we're not at the end  */
03779   /*  of the tree, compare the bond passed in with the tree  */
03780   while (!found && (ptr!=NULL))
03781   {
03782     cmp_code=strcasecmp(atom1, ptr->atom1name);
03783 
03784     if (cmp_code == 0)
03785     {
03786       cmp_code=strcasecmp(atom2, ptr->atom2name);
03787     }
03788 
03789     if (cmp_code == 0)
03790     {
03791       /*  Found a match        */
03792       found=1;
03793       bond_ptr->bond_type = ptr->index;
03794     }
03795     else if (cmp_code < 0)
03796     {
03797       /*  Go left          */
03798       ptr=ptr->left;
03799     }
03800     else
03801     {
03802       /*  Go right          */
03803       ptr=ptr->right;
03804     }
03805   }
03806 
03807   /*  Check to see if we found anything        */
03808   if (!found)
03809   {
03810     if ((strcmp(atom1, "DRUD")==0 || strcmp(atom2, "DRUD")==0)
03811         && (strcmp(atom1, "X")!=0 && strcmp(atom2, "X")!=0)) {
03812       /* try a wildcard DRUD X match for this Drude bond */
03813       char a1[8] = "DRUD", a2[8] = "X";
03814       return assign_bond_index(a1, a2, bond_ptr);  /* recursive call */
03815     }
03816     else {
03817       char err_msg[128];
03818 
03819       sprintf(err_msg, "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
03820          atom1, atom2, bond_ptr->atom1+1, bond_ptr->atom2+1);
03821       NAMD_die(err_msg);
03822     }
03823   }
03824 
03825   return;
03826 }

void Parameters::assign_crossterm_index ( char *  ,
char *  ,
char *  ,
char *  ,
char *  ,
char *  ,
char *  ,
char *  ,
Crossterm  
)

Definition at line 4160 of file Parameters.C.

References crossterm::atom1, crossterm_params::atom1name, crossterm_params::atom2name, crossterm_params::atom3name, crossterm::atom4, crossterm_params::atom4name, crossterm_params::atom5name, crossterm_params::atom6name, crossterm_params::atom7name, crossterm::atom8, crossterm_params::atom8name, crossterm::crossterm_type, crossterm_params::index, NAMD_die(), and crossterm_params::next.

04163 {
04164   struct crossterm_params *ptr;  //  Current position in list
04165   int found=0;      //  Flag 1->found a match
04166 
04167   /*  Start at the head of the list        */
04168   ptr=crosstermp;
04169 
04170   /*  While we haven't fuond a match and haven't reached the end  */
04171   /*  of the list, keep looking          */
04172   while (!found && (ptr!=NULL))
04173   {
04174     /*  Do a linear search through the linked list of   */
04175     /*  crossterm parameters.  Since the list is arranged    */
04176     /*  with wildcard paramters at the end of the list, we  */
04177     /*  can simply do a linear search and be guaranteed that*/
04178     /*  we will find exact matches before wildcard matches. */
04179     /*  Also, we must check for an exact match, and a match */
04180     /*  in reverse, since they are really the same          */
04181     /*  physically.            */
04182     if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
04183            (strcasecmp(ptr->atom1name, "X")==0) ) &&
04184        ( (strcasecmp(ptr->atom2name, atom2)==0) || 
04185            (strcasecmp(ptr->atom2name, "X")==0) ) &&
04186        ( (strcasecmp(ptr->atom3name, atom3)==0) || 
04187            (strcasecmp(ptr->atom3name, "X")==0) ) &&
04188        ( (strcasecmp(ptr->atom4name, atom4)==0) || 
04189            (strcasecmp(ptr->atom4name, "X")==0) ) )
04190     {
04191       /*  Found an exact match      */
04192       found=1;
04193     }
04194     else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
04195            (strcasecmp(ptr->atom4name, "X")==0) ) &&
04196        ( (strcasecmp(ptr->atom3name, atom2)==0) || 
04197            (strcasecmp(ptr->atom3name, "X")==0) ) &&
04198        ( (strcasecmp(ptr->atom2name, atom3)==0) || 
04199            (strcasecmp(ptr->atom2name, "X")==0) ) &&
04200        ( (strcasecmp(ptr->atom1name, atom4)==0) || 
04201            (strcasecmp(ptr->atom1name, "X")==0) ) )
04202     {
04203       /*  Found a reverse match      */
04204       found=1;
04205     }
04206     if ( ! found ) {
04207       /*  Didn't find a match, go to the next node  */
04208       ptr=ptr->next;
04209       continue;
04210     }
04211     found = 0;
04212     if ( ( (strcasecmp(ptr->atom5name, atom5)==0) || 
04213            (strcasecmp(ptr->atom5name, "X")==0) ) &&
04214        ( (strcasecmp(ptr->atom6name, atom6)==0) || 
04215            (strcasecmp(ptr->atom6name, "X")==0) ) &&
04216        ( (strcasecmp(ptr->atom7name, atom7)==0) || 
04217            (strcasecmp(ptr->atom7name, "X")==0) ) &&
04218        ( (strcasecmp(ptr->atom8name, atom8)==0) || 
04219            (strcasecmp(ptr->atom8name, "X")==0) ) )
04220     {
04221       /*  Found an exact match      */
04222       found=1;
04223     }
04224     else if ( ( (strcasecmp(ptr->atom8name, atom5)==0) || 
04225            (strcasecmp(ptr->atom8name, "X")==0) ) &&
04226        ( (strcasecmp(ptr->atom7name, atom6)==0) || 
04227            (strcasecmp(ptr->atom7name, "X")==0) ) &&
04228        ( (strcasecmp(ptr->atom6name, atom7)==0) || 
04229            (strcasecmp(ptr->atom6name, "X")==0) ) &&
04230        ( (strcasecmp(ptr->atom5name, atom8)==0) || 
04231            (strcasecmp(ptr->atom5name, "X")==0) ) )
04232     {
04233       /*  Found a reverse match      */
04234       found=1;
04235     }
04236     if ( ! found ) {
04237       /*  Didn't find a match, go to the next node  */
04238       ptr=ptr->next;
04239     }
04240   }
04241 
04242   /*  Make sure we found a match          */
04243   if (!found)
04244   {
04245     char err_msg[128];
04246 
04247     sprintf(err_msg, "UNABLE TO FIND CROSSTERM PARAMETERS FOR %s  %s  %s  %s  %s  %s  %s  %s\n"
04248       "(ATOMS %i %i %i %i %i %i %i %i)",
04249       atom1, atom2, atom3, atom4, atom5, atom6, atom7, atom8,
04250       crossterm_ptr->atom1+1, crossterm_ptr->atom1+2, crossterm_ptr->atom1+3, crossterm_ptr->atom4+1,
04251       crossterm_ptr->atom1+5, crossterm_ptr->atom1+6, crossterm_ptr->atom1+7, crossterm_ptr->atom8+1);
04252     
04253     NAMD_die(err_msg);
04254   }
04255 
04256   /*  Assign the constants          */
04257   crossterm_ptr->crossterm_type = ptr->index;
04258 
04259   return;
04260 }

void Parameters::assign_dihedral_index ( char *  ,
char *  ,
char *  ,
char *  ,
Dihedral ,
int  ,
int   
)

Definition at line 3952 of file Parameters.C.

References dihedral::atom1, dihedral_params::atom1name, dihedral_params::atom1wild, dihedral::atom2, dihedral_params::atom2name, dihedral_params::atom2wild, dihedral::atom3, dihedral_params::atom3name, dihedral_params::atom3wild, dihedral::atom4, dihedral_params::atom4name, dihedral_params::atom4wild, dihedral_array, dihedral::dihedral_type, dihedral_params::index, iout, iWARN(), DihedralValue::multiplicity, NAMD_die(), dihedral_params::next, and paraXplor.

03956 {
03957   struct dihedral_params *ptr;  //  Current position in list
03958   int found=0;      //  Flag 1->found a match
03959 
03960   /*  Start at the begining of the list        */
03961   ptr=dihedralp;
03962 
03963   /*  While we haven't found a match and we haven't reached       */
03964   /*  the end of the list, keep looking        */
03965   while (!found && (ptr!=NULL))
03966   {
03967     /*  Do a linear search through the linked list of   */
03968     /*  dihedral parameters.  Since the list is arranged    */
03969     /*  with wildcard paramters at the end of the list, we  */
03970     /*  can simply do a linear search and be guaranteed that*/
03971     /*  we will find exact matches before wildcard matches. */
03972     /*  Also, we must check for an exact match, and a match */
03973     /*  in reverse, since they are really the same          */
03974     /*  physically.            */
03975     if ( ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom1)==0) ) && 
03976          ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom2)==0) ) &&
03977          ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom3)==0) ) &&
03978          ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom4)==0) ) ) 
03979     {
03980       /*  Found an exact match      */
03981       found=1;
03982     }
03983     else if ( ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom1)==0) ) &&
03984               ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom2)==0) ) &&
03985               ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom3)==0) ) &&
03986               ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom4)==0) ) )
03987     {
03988       /*  Found a reverse match      */
03989       found=1;
03990     }
03991     else
03992     {
03993       /*  Didn't find a match, go to the next node  */
03994       ptr=ptr->next;
03995     }
03996   }
03997 
03998   /*  Make sure we found a match          */
03999   if (!found)
04000   {
04001     char err_msg[128];
04002 
04003     sprintf(err_msg, "UNABLE TO FIND DIHEDRAL PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
04004        atom1, atom2, atom3, atom4, dihedral_ptr->atom1+1, dihedral_ptr->atom2+1,
04005        dihedral_ptr->atom3+1, dihedral_ptr->atom4+1);
04006     
04007     if ( notFoundIndex ) {
04008       dihedral_ptr->dihedral_type = notFoundIndex;
04009       iout << iWARN << err_msg << "\n" << endi;
04010       return;
04011     } else NAMD_die(err_msg);
04012   }
04013 
04014  if (paramType == paraXplor) {
04015   //  Check to make sure the number of multiples specified in the psf
04016   //  file doesn't exceed the number of parameters in the parameter
04017   //  files
04018   if (multiplicity > maxDihedralMults[ptr->index])
04019   {
04020     char err_msg[257];
04021 
04022     sprintf(err_msg, "Multiplicity of Paramters for diehedral bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxDihedralMults[ptr->index]);
04023     NAMD_die(err_msg);
04024   }
04025 
04026   //  If the multiplicity from the current bond is larger than that
04027   //  seen in the past, increase the multiplicity for this bond
04028   if (multiplicity > dihedral_array[ptr->index].multiplicity)
04029   {
04030     dihedral_array[ptr->index].multiplicity = multiplicity;
04031   }
04032  }
04033 
04034   dihedral_ptr->dihedral_type = ptr->index;
04035 
04036   return;
04037 }

void Parameters::assign_improper_index ( char *  ,
char *  ,
char *  ,
char *  ,
Improper ,
int   
)

Definition at line 4061 of file Parameters.C.

References improper::atom1, improper_params::atom1name, improper::atom2, improper_params::atom2name, improper::atom3, improper_params::atom3name, improper::atom4, improper_params::atom4name, improper_array, improper::improper_type, improper_params::index, ImproperValue::multiplicity, NAMD_die(), improper_params::next, and paraXplor.

04065 {
04066   struct improper_params *ptr;  //  Current position in list
04067   int found=0;      //  Flag 1->found a match
04068 
04069   /*  Start at the head of the list        */
04070   ptr=improperp;
04071 
04072   /*  While we haven't fuond a match and haven't reached the end  */
04073   /*  of the list, keep looking          */
04074   while (!found && (ptr!=NULL))
04075   {
04076     /*  Do a linear search through the linked list of   */
04077     /*  improper parameters.  Since the list is arranged    */
04078     /*  with wildcard paramters at the end of the list, we  */
04079     /*  can simply do a linear search and be guaranteed that*/
04080     /*  we will find exact matches before wildcard matches. */
04081     /*  Also, we must check for an exact match, and a match */
04082     /*  in reverse, since they are really the same          */
04083     /*  physically.            */
04084     if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
04085            (strcasecmp(ptr->atom1name, "X")==0) ) &&
04086        ( (strcasecmp(ptr->atom2name, atom2)==0) || 
04087            (strcasecmp(ptr->atom2name, "X")==0) ) &&
04088        ( (strcasecmp(ptr->atom3name, atom3)==0) || 
04089            (strcasecmp(ptr->atom3name, "X")==0) ) &&
04090        ( (strcasecmp(ptr->atom4name, atom4)==0) || 
04091            (strcasecmp(ptr->atom4name, "X")==0) ) )
04092     {
04093       /*  Found an exact match      */
04094       found=1;
04095     }
04096     else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
04097            (strcasecmp(ptr->atom4name, "X")==0) ) &&
04098        ( (strcasecmp(ptr->atom3name, atom2)==0) || 
04099            (strcasecmp(ptr->atom3name, "X")==0) ) &&
04100        ( (strcasecmp(ptr->atom2name, atom3)==0) || 
04101            (strcasecmp(ptr->atom2name, "X")==0) ) &&
04102        ( (strcasecmp(ptr->atom1name, atom4)==0) || 
04103            (strcasecmp(ptr->atom1name, "X")==0) ) )
04104     {
04105       /*  Found a reverse match      */
04106       found=1;
04107     }
04108     else
04109     {
04110       /*  Didn't find a match, go to the next node  */
04111       ptr=ptr->next;
04112     }
04113   }
04114 
04115   /*  Make sure we found a match          */
04116   if (!found)
04117   {
04118     char err_msg[128];
04119 
04120     sprintf(err_msg, "UNABLE TO FIND IMPROPER PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
04121        atom1, atom2, atom3, atom4, improper_ptr->atom1+1, improper_ptr->atom2+1,
04122        improper_ptr->atom3+1, improper_ptr->atom4+1);
04123     
04124     NAMD_die(err_msg);
04125   }
04126 
04127  if (paramType == paraXplor) {
04128   //  Check to make sure the number of multiples specified in the psf
04129   //  file doesn't exceed the number of parameters in the parameter
04130   //  files
04131   if (multiplicity > maxImproperMults[ptr->index])
04132   {
04133     char err_msg[257];
04134 
04135     sprintf(err_msg, "Multiplicity of Paramters for improper bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxImproperMults[ptr->index]);
04136     NAMD_die(err_msg);
04137   }
04138 
04139   //  If the multiplicity from the current bond is larger than that
04140   //  seen in the past, increase the multiplicity for this bond
04141   if (multiplicity > improper_array[ptr->index].multiplicity)
04142   {
04143     improper_array[ptr->index].multiplicity = multiplicity;
04144   }
04145  }
04146 
04147   /*  Assign the constants          */
04148   improper_ptr->improper_type = ptr->index;
04149 
04150   return;
04151 }

void Parameters::assign_vdw_index ( char *  ,
Atom  
)

Definition at line 3464 of file Parameters.C.

References vdw_params::atomname, vdw_params::index, iout, iWARN(), vdw_params::left, NAMD_die(), vdw_params::right, and atom_constants::vdw_type.

03466 {
03467   struct vdw_params *ptr;    //  Current position in trees
03468   int found=0;      //  Flag 1->found match
03469   int comp_code;      //  return code from strcasecmp
03470 
03471   /*  Check to make sure the files have all been read    */
03472   if (!AllFilesRead)
03473   {
03474     NAMD_die("Tried to assign vdw index before all parameter files were read");
03475   }
03476 
03477   /*  Start at the top            */
03478   ptr=vdwp;
03479 
03480   /*  While we haven't found a match, and we haven't reached      */
03481   /*  the bottom of the tree, compare the atom passed in with     */
03482   /*  the current value and decide if we have a match, or if not, */
03483   /*  which way to go            */
03484   while (!found && (ptr!=NULL))
03485   {
03486     comp_code = strcasecmp(atomtype, ptr->atomname);
03487 
03488     if (comp_code == 0)
03489     {
03490       /*  Found a match!        */
03491       atom_ptr->vdw_type=ptr->index;
03492       found=1;
03493     }
03494     else if (comp_code < 0)
03495     {
03496       /*  Go to the left        */
03497       ptr=ptr->left;
03498     }
03499     else
03500     {
03501       /*  Go to the right        */
03502       ptr=ptr->right;
03503     }
03504   }
03505 
03506   //****** BEGIN CHARMM/XPLOR type changes
03507   if (!found)
03508   {
03509     // since CHARMM allows wildcards "*" in vdw typenames
03510     // we have to look again if necessary, this way, if
03511     // we already had an exact match, this is never executed
03512     size_t windx;                      //  wildcard index
03513 
03514     /*  Start again at the top                                */
03515     ptr=vdwp;
03516   
03517      while (!found && (ptr!=NULL))
03518      {
03519   
03520        // get index of wildcard wildcard, get index
03521        windx= strcspn(ptr->atomname,"*"); 
03522        if (windx == strlen(ptr->atomname))
03523        {
03524          // there is no wildcard here
03525          comp_code = strcasecmp(atomtype, ptr->atomname);   
03526        }
03527        else
03528        {
03529          comp_code = strncasecmp(atomtype, ptr->atomname, windx); 
03530        }  
03531 
03532        if (comp_code == 0)
03533        {
03534          /*  Found a match!                                */
03535          atom_ptr->vdw_type=ptr->index;
03536          found=1;
03537          char errbuf[100];
03538          sprintf(errbuf,"VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
03539                         atomtype, ptr->atomname);
03540          int i;
03541          for(i=0; i<error_msgs.size(); i++) {
03542            if ( strcmp(errbuf,error_msgs[i]) == 0 ) break;
03543          }
03544          if ( i == error_msgs.size() ) {
03545            char *newbuf = new char[strlen(errbuf)+1];
03546            strcpy(newbuf,errbuf);
03547            error_msgs.add(newbuf);
03548            iout << iWARN << newbuf << "\n" << endi;
03549          }
03550        }
03551        else if (comp_code < 0)
03552        {
03553           /*  Go to the left                                */
03554                 ptr=ptr->left;
03555        }
03556        else
03557        {
03558          /*  Go to the right                                */
03559                 ptr=ptr->right;
03560        }
03561      
03562      }
03563                 
03564   }
03565   //****** END CHARMM/XPLOR type changes
03566 
03567   /*  Make sure we found it          */
03568   if (!found)
03569   {
03570     char err_msg[100];
03571 
03572     sprintf(err_msg, "DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s",
03573        atomtype);
03574     NAMD_die(err_msg);
03575   }
03576 
03577   return;
03578 }

char* Parameters::atom_type_name ( Index  a  )  [inline]

Definition at line 390 of file Parameters.h.

References MAX_ATOMTYPE_CHARS.

00390                                       {
00391           return (atomTypeNames + (a * (MAX_ATOMTYPE_CHARS + 1)));
00392         }

void Parameters::done_reading_files (  ) 

Definition at line 2958 of file Parameters.C.

References angle_array, bond_array, crossterm_array, dihedral_array, gromacsPair_array, improper_array, MAX_ATOMTYPE_CHARS, NAMD_die(), nbthole_array, NumAngleParams, NumBondParams, NumCrosstermParams, NumDihedralParams, NumGromacsPairParams, NumImproperParams, NumNbtholePairParams, NumVdwParams, NumVdwParamsAssigned, TRUE, and vdw_array.

Referenced by Parameters().

02960 {
02961   AllFilesRead = TRUE;
02962 
02963   //  Allocate space for all of the arrays
02964   if (NumBondParams)
02965   {
02966     bond_array = new BondValue[NumBondParams];
02967 
02968     if (bond_array == NULL)
02969     {
02970       NAMD_die("memory allocation of bond_array failed!");
02971     }
02972     memset(bond_array, 0, NumBondParams*sizeof(BondValue));
02973   }
02974 
02975   if (NumAngleParams)
02976   {
02977     angle_array = new AngleValue[NumAngleParams];
02978 
02979     if (angle_array == NULL)
02980     {
02981       NAMD_die("memory allocation of angle_array failed!");
02982     }
02983     memset(angle_array, 0, NumAngleParams*sizeof(AngleValue));
02984     for ( Index i=0; i<NumAngleParams; ++i ) {
02985       angle_array[i].normal = 1;
02986     }
02987   }
02988 
02989   if (NumDihedralParams)
02990   {
02991     dihedral_array = new DihedralValue[NumDihedralParams];
02992 
02993     if (dihedral_array == NULL)
02994     {
02995       NAMD_die("memory allocation of dihedral_array failed!");
02996     }
02997     memset(dihedral_array, 0, NumDihedralParams*sizeof(DihedralValue));
02998   }
02999 
03000   if (NumImproperParams)
03001   {
03002     improper_array = new ImproperValue[NumImproperParams];
03003 
03004     if (improper_array == NULL)
03005     {
03006       NAMD_die("memory allocation of improper_array failed!");
03007     }
03008     memset(improper_array, 0, NumImproperParams*sizeof(ImproperValue));
03009   }
03010 
03011   if (NumCrosstermParams)
03012   {
03013     crossterm_array = new CrosstermValue[NumCrosstermParams];
03014     memset(crossterm_array, 0, NumCrosstermParams*sizeof(CrosstermValue));
03015   }
03016 
03017   // JLai
03018   if (NumGromacsPairParams)
03019   {
03020     gromacsPair_array = new GromacsPairValue[NumGromacsPairParams];
03021     memset(gromacsPair_array, 0, NumGromacsPairParams*sizeof(GromacsPairValue)); 
03022   }
03023   // End of JLai
03024 
03025   if (NumVdwParams)
03026   {
03027           atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
03028     vdw_array = new VdwValue[NumVdwParams];
03029     
03030     if (vdw_array == NULL)
03031     {
03032       NAMD_die("memory allocation of vdw_array failed!");
03033     }
03034   }
03035   if (NumNbtholePairParams)
03036   {
03037     nbthole_array = new NbtholePairValue[NumNbtholePairParams];
03038 
03039     if(nbthole_array == NULL)
03040     {
03041       NAMD_die("memory allocation of nbthole_array failed!");
03042     }
03043   }
03044   //  Assign indexes to each of the parameters and populate the
03045   //  arrays using the binary trees and linked lists that we have
03046   //  already read in
03047 
03048   // Note that if parameters have been overwritten (matching
03049   // atom patterns but different parameter values) the tree
03050   // contains fewer elements than Num...Params would suggest.
03051   // The arrays are initialized above because the end values
03052   // may not be occupied.  Modifying the Num...Params values
03053   // would break backwards compatibility of memopt extraBonds.
03054 
03055   index_bonds(bondp, 0);
03056   index_angles(anglep, 0);
03057   NumVdwParamsAssigned = index_vdw(vdwp, 0);
03058   index_dihedrals();
03059   index_impropers();
03060   index_crossterms();
03061   convert_nbthole_pairs();
03062   //  Convert the vdw pairs
03063   convert_vdw_pairs();
03064   convert_table_pairs();
03065 }

void Parameters::done_reading_structure (  ) 

Definition at line 4995 of file Parameters.C.

04997 {
04998   if (bondp != NULL)
04999     free_bond_tree(bondp);
05000 
05001   if (anglep != NULL)
05002     free_angle_tree(anglep);
05003 
05004   if (dihedralp != NULL)
05005     free_dihedral_list(dihedralp);
05006 
05007   if (improperp != NULL)
05008     free_improper_list(improperp);
05009 
05010   if (crosstermp != NULL)
05011     free_crossterm_list(crosstermp);
05012 
05013   if (vdwp != NULL)
05014     free_vdw_tree(vdwp);
05015 
05016   //  Free the arrays used to track multiplicity for dihedrals
05017   //  and impropers
05018   if (maxDihedralMults != NULL)
05019     delete [] maxDihedralMults;
05020 
05021   if (maxImproperMults != NULL)
05022     delete [] maxImproperMults;
05023 
05024   bondp=NULL;
05025   anglep=NULL;
05026   dihedralp=NULL;
05027   improperp=NULL;
05028   crosstermp=NULL;
05029   vdwp=NULL;
05030   maxImproperMults=NULL;
05031   maxDihedralMults=NULL;
05032 }

void Parameters::get_angle_params ( Real k,
Real theta0,
Real k_ub,
Real r_ub,
Index  index 
) [inline]

Definition at line 450 of file Parameters.h.

References angle_array, AngleValue::k, AngleValue::k_ub, AngleValue::r_ub, and AngleValue::theta0.

00452         {
00453                 *k = angle_array[index].k;
00454                 *theta0 = angle_array[index].theta0;
00455                 *k_ub = angle_array[index].k_ub;
00456                 *r_ub = angle_array[index].r_ub;
00457         }

void Parameters::get_bond_params ( Real k,
Real x0,
Index  index 
) [inline]

Definition at line 444 of file Parameters.h.

References bond_array, BondValue::k, and BondValue::x0.

Referenced by Molecule::print_bonds().

00445         {
00446                 *k = bond_array[index].k;
00447                 *x0 = bond_array[index].x0;
00448         }

int Parameters::get_dihedral_multiplicity ( Index  index  )  [inline]

Definition at line 464 of file Parameters.h.

References dihedral_array.

00465         {
00466                 return(dihedral_array[index].multiplicity);
00467         }

void Parameters::get_dihedral_params ( Real k,
int *  n,
Real delta,
Index  index,
int  mult 
) [inline]

Definition at line 482 of file Parameters.h.

References four_body_consts::delta, dihedral_array, four_body_consts::k, MAX_MULTIPLICITY, four_body_consts::n, NAMD_die(), and DihedralValue::values.

00484         {
00485                 if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
00486                 {
00487                         NAMD_die("Bad mult index in Parameters::get_dihedral_params");
00488                 }
00489 
00490                 *k = dihedral_array[index].values[mult].k;
00491                 *n = dihedral_array[index].values[mult].n;
00492                 *delta = dihedral_array[index].values[mult].delta;
00493         }

int Parameters::get_improper_multiplicity ( Index  index  )  [inline]

Definition at line 459 of file Parameters.h.

References improper_array.

00460         {
00461                 return(improper_array[index].multiplicity);
00462         }

void Parameters::get_improper_params ( Real k,
int *  n,
Real delta,
Index  index,
int  mult 
) [inline]

Definition at line 469 of file Parameters.h.

References four_body_consts::delta, improper_array, four_body_consts::k, MAX_MULTIPLICITY, four_body_consts::n, NAMD_die(), and ImproperValue::values.

00471         {
00472                 if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
00473                 {
00474                         NAMD_die("Bad mult index in Parameters::get_improper_params");
00475                 }
00476 
00477                 *k = improper_array[index].values[mult].k;
00478                 *n = improper_array[index].values[mult].n;
00479                 *delta = improper_array[index].values[mult].delta;
00480         }

int Parameters::get_int_table_type ( char *   ) 

Definition at line 7560 of file Parameters.C.

References table_types, and tablenumtypes.

07560                                                   {
07561   for (int i=0; i<tablenumtypes; i++) {
07562     if (!strncmp(tabletype, table_types[i], 5)) {
07563       return i;
07564     }
07565   }
07566 
07567   return -1;
07568 }

int Parameters::get_num_vdw_params ( void   )  [inline]

Definition at line 529 of file Parameters.h.

References NumVdwParamsAssigned.

Referenced by LJTable::LJTable(), Molecule::receive_Molecule(), and Molecule::send_Molecule().

00529 { return NumVdwParamsAssigned; }

int Parameters::get_table_pair_params ( Index  ,
Index  ,
int *   
)

Definition at line 3599 of file Parameters.C.

References FALSE, indexed_table_pair::ind1, indexed_table_pair::ind2, indexed_table_pair::K, indexed_table_pair::left, indexed_table_pair::right, tab_pair_tree, and TRUE.

03599                                                                     {
03600   IndexedTablePair *ptr;
03601   Index temp;
03602   int found=FALSE;
03603 
03604   ptr=tab_pair_tree;
03605   //
03606   //  We need the smaller type in ind1, so if it isn't already that 
03607   //  way, switch them        */
03608   if (ind1 > ind2)
03609   {
03610     temp = ind1;
03611     ind1 = ind2;
03612     ind2 = temp;
03613   }
03614 
03615   /*  While we haven't found a match and we're not at the end  */
03616   /*  of the tree, compare the bond passed in with the tree  */
03617   while (!found && (ptr!=NULL))
03618   {
03619 //    printf("Comparing %i with %i and %i with %i\n", ind1, ptr->ind1, ind2, ptr->ind2);
03620     if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
03621     {
03622        found = TRUE;
03623     }
03624     else if ( (ind1 < ptr->ind1) || 
03625         ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
03626     {
03627       /*  Go left          */
03628       ptr=ptr->left;
03629     }
03630     else
03631     {
03632       /*  Go right          */
03633       ptr=ptr->right;
03634     }
03635   }
03636 
03637   /*  If we found a match, assign the values      */
03638   if (found)
03639   {
03640     *K = ptr->K;
03641     return(TRUE);
03642   }
03643   else
03644   {
03645     return(FALSE);
03646   }
03647 }

int Parameters::get_vdw_pair_params ( Index  ind1,
Index  ind2,
Real ,
Real ,
Real ,
Real  
)

Definition at line 3674 of file Parameters.C.

References indexed_vdw_pair::A, indexed_vdw_pair::A14, indexed_vdw_pair::B, indexed_vdw_pair::B14, FALSE, indexed_vdw_pair::ind1, indexed_vdw_pair::ind2, indexed_vdw_pair::left, indexed_vdw_pair::right, TRUE, and vdw_pair_tree.

Referenced by get_vdw_params().

03677 {
03678   IndexedVdwPair *ptr;    //  Current location in tree
03679   Index temp;      //  Temporary value for swithcing
03680           // values
03681   int found=FALSE;    //  Flag 1-> found a match
03682 
03683   ptr=vdw_pair_tree;
03684 
03685   //  We need the smaller type in ind1, so if it isn't already that 
03686   //  way, switch them        */
03687   if (ind1 > ind2)
03688   {
03689     temp = ind1;
03690     ind1 = ind2;
03691     ind2 = temp;
03692   }
03693 
03694   /*  While we haven't found a match and we're not at the end  */
03695   /*  of the tree, compare the bond passed in with the tree  */
03696   while (!found && (ptr!=NULL))
03697   {
03698     if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
03699     {
03700        found = TRUE;
03701     }
03702     else if ( (ind1 < ptr->ind1) || 
03703         ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
03704     {
03705       /*  Go left          */
03706       ptr=ptr->left;
03707     }
03708     else
03709     {
03710       /*  Go right          */
03711       ptr=ptr->right;
03712     }
03713   }
03714 
03715   /*  If we found a match, assign the values      */
03716   if (found)
03717   {
03718     *A = ptr->A;
03719     *B = ptr->B;
03720     *A14 = ptr->A14;
03721     *B14 = ptr->B14;
03722 
03723     return(TRUE);
03724   }
03725   else
03726   {
03727     return(FALSE);
03728   }
03729 }

void Parameters::get_vdw_params ( Real sigma,
Real epsilon,
Real sigma14,
Real epsilon14,
Index  index 
) [inline]

Definition at line 495 of file Parameters.h.

References A, B, cbrt, vdw_val::epsilon, vdw_val::epsilon14, get_vdw_pair_params(), NAMD_die(), vdw_val::sigma, vdw_val::sigma14, and vdw_array.

Referenced by Molecule::print_atoms().

00497   {
00498     if ( vdw_array ) {
00499       *sigma = vdw_array[index].sigma;
00500       *epsilon = vdw_array[index].epsilon;
00501       *sigma14 = vdw_array[index].sigma14;
00502       *epsilon14 = vdw_array[index].epsilon14;
00503     } else {
00504       // sigma and epsilon from A and B for each vdw type's interaction with itself
00505       Real A; Real B; Real A14; Real B14;
00506       if (get_vdw_pair_params(index, index, &A, &B, &A14, &B14) ) {
00507         if (A && B) {
00508           *sigma = sqrt(cbrt(A)) / sqrt(cbrt(B));
00509           *epsilon = B*B / (4*A);
00510         }
00511         else {
00512           *sigma = 0; *epsilon = 0;
00513         }
00514         if (A14 && B14) {
00515           *sigma14 = sqrt(cbrt(A14)) / sqrt(cbrt(B14));
00516           *epsilon14 = B14*B14 / (4*A14);
00517         }
00518         else {
00519           *sigma14 = 0; *epsilon14 = 0;
00520         }
00521       }
00522       else {NAMD_die ("Function get_vdw_params failed to derive Lennard-Jones sigma and epsilon from A and B values\n");}
00523     }
00524   }

void Parameters::print_angle_params (  ) 

Definition at line 4861 of file Parameters.C.

References DebugM, and NumAngleParams.

04862 {
04863   DebugM(3,NumAngleParams << " ANGLE PARAMETERS\n"
04864       << "*****************************************" );
04865   traverse_angle_params(anglep);
04866 }

void Parameters::print_bond_params (  ) 

Definition at line 4843 of file Parameters.C.

References DebugM, and NumBondParams.

04844 {
04845   DebugM(3,NumBondParams << " BOND PARAMETERS\n" \
04846       << "*****************************************"  \
04847       );
04848 
04849   traverse_bond_params(bondp);
04850 }

void Parameters::print_crossterm_params (  ) 

void Parameters::print_dihedral_params (  ) 

Definition at line 4877 of file Parameters.C.

References DebugM, and NumDihedralParams.

04878 {
04879   DebugM(3,NumDihedralParams << " DIHEDRAL PARAMETERS\n" \
04880       << "*****************************************" );
04881 
04882   traverse_dihedral_params(dihedralp);
04883 }

void Parameters::print_improper_params (  ) 

Definition at line 4894 of file Parameters.C.

References DebugM, and NumImproperParams.

04895 {
04896   DebugM(3,NumImproperParams << " IMPROPER PARAMETERS\n" \
04897       << "*****************************************" );
04898 
04899   traverse_improper_params(improperp);
04900 }

void Parameters::print_nbthole_pair_params (  ) 

Definition at line 4945 of file Parameters.C.

References DebugM, and NumNbtholePairParams.

04946 {
04947   DebugM(3,NumNbtholePairParams << " NBTHOLE PAIR PARAMETERS\n" \
04948       << "*****************************************" );
04949 
04950   traverse_nbthole_pair_params(nbthole_pairp);
04951 } 

void Parameters::print_param_summary (  ) 

Definition at line 4962 of file Parameters.C.

References iINFO(), iout, NumAngleParams, NumBondParams, NumCosAngles, NumCrosstermParams, NumDihedralParams, NumImproperParams, NumNbtholePairParams, NumVdwPairParams, and NumVdwParams.

Referenced by NamdState::loadStructure().

04963 {
04964   iout << iINFO << "SUMMARY OF PARAMETERS:\n" 
04965        << iINFO << NumBondParams << " BONDS\n" 
04966        << iINFO << NumAngleParams << " ANGLES\n" << endi;
04967        if (cosAngles) {
04968          iout << iINFO << "  " << (NumAngleParams - NumCosAngles) << " HARMONIC\n"
04969               << iINFO << "  " << NumCosAngles << " COSINE-BASED\n" << endi;
04970        }
04971   iout << iINFO << NumDihedralParams << " DIHEDRAL\n"
04972        << iINFO << NumImproperParams << " IMPROPER\n"
04973        << iINFO << NumCrosstermParams << " CROSSTERM\n"
04974        << iINFO << NumVdwParams << " VDW\n"
04975        << iINFO << NumVdwPairParams << " VDW_PAIRS\n"
04976        << iINFO << NumNbtholePairParams << " NBTHOLE_PAIRS\n" << endi;
04977 }

void Parameters::print_vdw_pair_params (  ) 

Definition at line 4928 of file Parameters.C.

References DebugM, and NumVdwPairParams.

04929 {
04930   DebugM(3,NumVdwPairParams << " vdW PAIR PARAMETERS\n" \
04931       << "*****************************************" );
04932 
04933   traverse_vdw_pair_params(vdw_pairp);
04934 }

void Parameters::print_vdw_params (  ) 

Definition at line 4911 of file Parameters.C.

References DebugM, and NumVdwParams.

04912 {
04913   DebugM(3,NumVdwParams << " vdW PARAMETERS\n" \
04914       << "*****************************************" );
04915 
04916   traverse_vdw_params(vdwp);
04917 }

void Parameters::read_charmm_parameter_file ( char *   ) 

Definition at line 539 of file Parameters.C.

References endi(), iout, iWARN(), NAMD_blank_string(), NAMD_die(), NAMD_find_first_word(), NAMD_read_line(), NumAngleParams, NumBondParams, NumCrosstermParams, NumDihedralParams, NumImproperParams, NumNbtholePairParams, NumTablePairParams, NumVdwPairParams, NumVdwParams, and table_ener.

Referenced by Parameters().

00541 {
00542   int  par_type=0;         //  What type of parameter are we currently
00543                            //  dealing with? (vide infra)
00544   int  skipline;           //  skip this line?
00545   int  skipall = 0;        //  skip rest of file;
00546   char buffer[512];           //  Buffer to store each line of the file
00547   char first_word[512];           //  First word of the current line
00548   FILE *pfile;                   //  File descriptor for the parameter file
00549 
00550   /*  Check to make sure that we haven't previously been told     */
00551   /*  that all the files were read                                */
00552   if (AllFilesRead)
00553   {
00554     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00555   }
00556 
00557   /*  Try and open the file                                        */
00558   if ( (pfile = fopen(fname, "r")) == NULL)
00559   {
00560     char err_msg[256];
00561 
00562     sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
00563     NAMD_die(err_msg);
00564   }
00565 
00566   /*  Keep reading in lines until we hit the EOF                        */
00567   while (NAMD_read_line(pfile, buffer) != -1)
00568   {
00569     /*  Get the first word of the line                        */
00570     NAMD_find_first_word(buffer, first_word);
00571     skipline=0;
00572 
00573     /*  First, screen out things that we ignore.                   */   
00574     /*  blank lines, lines that start with '!' or '*', lines that  */
00575     /*  start with "END".                                          */
00576     if (!NAMD_blank_string(buffer) &&
00577         (strncmp(first_word, "!", 1) != 0) &&
00578          (strncmp(first_word, "*", 1) != 0) &&
00579         (strncasecmp(first_word, "END", 3) != 0))
00580     {
00581       if ( skipall ) {
00582         iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
00583         break;
00584       }
00585       /*  Now, determine the apropriate parameter type.   */
00586       if (strncasecmp(first_word, "bond", 4)==0)
00587       {
00588         par_type=1; skipline=1;
00589       }
00590       else if (strncasecmp(first_word, "angl", 4)==0)
00591       {
00592         par_type=2; skipline=1;
00593       }
00594       else if (strncasecmp(first_word, "thet", 4)==0)
00595       {
00596         par_type=2; skipline=1;
00597       }
00598       else if (strncasecmp(first_word, "dihe", 4)==0)
00599       {
00600         par_type=3; skipline=1;
00601       }
00602       else if (strncasecmp(first_word, "phi", 3)==0)
00603       {
00604         par_type=3; skipline=1;
00605       }
00606       else if (strncasecmp(first_word, "impr", 4)==0)
00607       {
00608         par_type=4; skipline=1;
00609       }
00610       else if (strncasecmp(first_word, "imph", 4)==0)
00611       {
00612         par_type=4; skipline=1;
00613       }
00614       else if (strncasecmp(first_word, "nbon", 4)==0)
00615       {
00616         par_type=5; skipline=1;
00617       }
00618       else if (strncasecmp(first_word, "nonb", 4)==0)
00619       {
00620         par_type=5; skipline=1;
00621       }
00622       else if (strncasecmp(first_word, "nbfi", 4)==0)
00623       {
00624         par_type=6; skipline=1;
00625       }
00626       else if (strncasecmp(first_word, "hbon", 4)==0)
00627       {
00628         par_type=7; skipline=1;
00629       }
00630       else if (strncasecmp(first_word, "cmap", 4)==0)
00631       {
00632         par_type=8; skipline=1;
00633       }
00634       else if (strncasecmp(first_word, "nbta", 4)==0)
00635       {
00636         par_type=9; skipline=1;
00637       }
00638       else if (strncasecmp(first_word, "thol", 4)==0)
00639       {
00640         par_type=10; skipline=1;
00641       }
00642       else if (strncasecmp(first_word, "atom", 4)==0)
00643       {
00644         par_type=11; skipline=1;
00645       }
00646       else if (strncasecmp(first_word, "ioformat", 8)==0)
00647       {
00648         skipline=1;
00649       }
00650       else if (strncasecmp(first_word, "read", 4)==0)
00651       {
00652         skip_stream_read(buffer, pfile);  skipline=1;
00653       }
00654       else if (strncasecmp(first_word, "return", 4)==0)
00655       {
00656         skipall=1;  skipline=1;
00657       }
00658       else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
00659                (strncasecmp(first_word, "grou", 4) == 0) ||
00660                (strncasecmp(first_word, "cdie", 4) == 0) ||
00661                (strncasecmp(first_word, "shif", 4) == 0) ||
00662                (strncasecmp(first_word, "vgro", 4) == 0) ||
00663                (strncasecmp(first_word, "vdis", 4) == 0) ||
00664                (strncasecmp(first_word, "vswi", 4) == 0) ||
00665                (strncasecmp(first_word, "cutn", 4) == 0) ||
00666                (strncasecmp(first_word, "ctof", 4) == 0) ||
00667                (strncasecmp(first_word, "cton", 4) == 0) ||
00668                (strncasecmp(first_word, "eps", 3) == 0) ||
00669                (strncasecmp(first_word, "e14f", 4) == 0) ||
00670                (strncasecmp(first_word, "wmin", 4) == 0) ||
00671                (strncasecmp(first_word, "aexp", 4) == 0) ||
00672                (strncasecmp(first_word, "rexp", 4) == 0) ||
00673                (strncasecmp(first_word, "haex", 4) == 0) ||
00674                (strncasecmp(first_word, "aaex", 4) == 0) ||
00675                (strncasecmp(first_word, "noac", 4) == 0) ||
00676                (strncasecmp(first_word, "hbno", 4) == 0) ||
00677                (strncasecmp(first_word, "cuth", 4) == 0) ||
00678                (strncasecmp(first_word, "ctof", 4) == 0) ||
00679                (strncasecmp(first_word, "cton", 4) == 0) ||
00680                (strncasecmp(first_word, "cuth", 4) == 0) ||
00681                (strncasecmp(first_word, "ctof", 4) == 0) ||
00682                (strncasecmp(first_word, "cton", 4) == 0) ) 
00683       {
00684         if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
00685         {
00686           char err_msg[512];
00687 
00688           sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00689           NAMD_die(err_msg);
00690         }
00691         else 
00692         {
00693           skipline = 1;
00694         }
00695       }        
00696       else if (par_type == 0)
00697       {
00698         /*  This is an unknown paramter.        */
00699         /*  This is BAD                                */
00700         char err_msg[512];
00701 
00702         sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00703         NAMD_die(err_msg);
00704       }
00705     }
00706     else
00707     {
00708       skipline=1;
00709     }
00710 
00711     if ( (par_type != 0) && (!skipline) )
00712     {
00713       /*  Now, call the appropriate function based    */
00714       /*  on the type of parameter we have                */
00715       /*  I know, this should really be a switch ...  */
00716       if (par_type == 1)
00717       {
00718         add_bond_param(buffer);
00719         NumBondParams++;
00720       }
00721       else if (par_type == 2)
00722       {
00723         add_angle_param(buffer);
00724         NumAngleParams++;
00725       }
00726       else if (par_type == 3)
00727       {
00728         add_dihedral_param(buffer, pfile);
00729         NumDihedralParams++;
00730       }
00731       else if (par_type == 4)
00732       {
00733         add_improper_param(buffer, pfile);
00734         NumImproperParams++;
00735       }
00736       else if (par_type == 5)
00737       {
00738         add_vdw_param(buffer);
00739         NumVdwParams++;
00740       }
00741       else if (par_type == 6)
00742       {
00743         add_vdw_pair_param(buffer);
00744         NumVdwPairParams++; 
00745       }
00746       else if (par_type == 7)
00747       {
00748         add_hb_pair_param(buffer);                  
00749       }
00750       else if (par_type == 8)
00751       {
00752         add_crossterm_param(buffer, pfile);                  
00753         NumCrosstermParams++;
00754       }
00755       else if (par_type == 9)
00756       {
00757 
00758         if (table_ener == NULL) {
00759           continue;
00760         }
00761 
00762         add_table_pair_param(buffer);                  
00763         NumTablePairParams++;
00764       }
00765       else if (par_type == 10)
00766       {
00767         add_nbthole_pair_param(buffer);
00768         NumNbtholePairParams++;
00769       }
00770       else if (par_type == 11)
00771       {
00772         if ( strncasecmp(first_word, "mass", 4) != 0 ) {
00773           char err_msg[512];
00774           sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00775           NAMD_die(err_msg);
00776         }
00777       }
00778       else
00779       {
00780         /*  This really should not occour!      */
00781         /*  This is an internal error.          */
00782         /*  This is VERY BAD                        */
00783         char err_msg[512];
00784 
00785         sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00786         NAMD_die(err_msg);
00787       }
00788     }
00789   }
00790 
00791   /*  Close the file                                                */
00792   fclose(pfile);
00793 
00794   return;
00795 }

void Parameters::read_ener_table ( SimParameters  ) 

Definition at line 6724 of file Parameters.C.

References iout, NAMD_die(), namdnearbyint, numenerentries, read_energy_type(), read_energy_type_bothcubspline(), simParams, table_ener, table_spacing, table_types, and tablenumtypes.

Referenced by Parameters().

06724                                                          {
06725         char* table_file = simParams->tabulatedEnergiesFile;
06726   char* interp_type = simParams->tableInterpType;
06727         FILE* enertable;
06728         enertable = fopen(table_file, "r");
06729 
06730         if (enertable == NULL) {
06731                 NAMD_die("ERROR: Could not open energy table file!\n");
06732         }
06733 
06734         char tableline[121];
06735   char* newtypename;
06736   int numtypes;
06737         int atom1;
06738         int atom2;
06739         int distbin;
06740   int readcount;
06741         Real dist;
06742         Real energy;
06743         Real force;
06744         Real table_spacing;
06745         Real maxdist;
06746 
06747 /* First find the header line and set the variables we need */
06748         iout << "Beginning table read\n" << endi;
06749         while(fgets(tableline,120,enertable)) {
06750                 if (strncmp(tableline,"#",1)==0) {
06751                         continue;
06752                 }
06753     readcount = sscanf(tableline, "%i %f %f", &numtypes, &table_spacing, &maxdist);
06754     if (readcount != 3) {
06755       NAMD_die("ERROR: Couldn't parse table header line\n");
06756     }
06757     break;
06758   }
06759 
06760   if (maxdist < simParams->cutoff) {
06761     NAMD_die("Tabulated energies must at least extend to the cutoff distance\n");
06762   }
06763 
06764         if (maxdist > simParams->cutoff) {
06765                 iout << "Info: Reducing max table size to match nonbond cutoff\n";
06766                 maxdist = ceil(simParams->cutoff);
06767         }
06768 
06769 /* Now allocate memory for the table; we know what we should be getting */
06770         numenerentries = 2 * numtypes * int (namdnearbyint(maxdist/table_spacing) + 1);
06771         //Set up a full energy lookup table from a file
06772         //Allocate the table; layout is atom1 x atom2 x distance energy force
06773         fprintf(stdout, "Table has %i entries\n",numenerentries);
06774         //iout << "Allocating original energy table\n" << endi;
06775         if ( table_ener ) delete [] table_ener;
06776         table_ener = new BigReal[numenerentries];
06777   if ( table_types ) delete [] table_types;
06778   table_types = new char*[numtypes];
06779 
06780 /* Finally, start reading the table */
06781   int numtypesread = 0; //Number of types read so far
06782 
06783         while(fgets(tableline,120,enertable)) {
06784                 if (strncmp(tableline,"#",1)==0) {
06785                         continue;
06786                 }
06787     if (strncmp(tableline,"TYPE",4)==0) {
06788       // Read a new type
06789       newtypename = new char[5];
06790       int readcount = sscanf(tableline, "%*s %s", newtypename);
06791       if (readcount != 1) {
06792         NAMD_die("Couldn't read interaction type from TYPE line\n");
06793       }
06794 //      printf("Setting table_types[%i] to %s\n", numtypesread, newtypename);
06795       table_types[numtypesread] = newtypename;
06796 
06797       if (numtypesread == numtypes) {
06798         NAMD_die("Error: Number of types in table doesn't match table header\n");
06799       }
06800 
06801       // Read the current energy type with the proper interpolation
06802       if (!strncasecmp(interp_type, "linear", 6)) {
06803         if (read_energy_type(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
06804           char err_msg[512];
06805           sprintf(err_msg, "Failed to read table energy (linear) type %s\n", newtypename);
06806           NAMD_die(err_msg);
06807         }
06808       } else if (!strncasecmp(interp_type, "cubic", 5)) {
06809         if (read_energy_type_bothcubspline(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
06810           char err_msg[512];
06811           sprintf(err_msg, "Failed to read table energy (cubic) type %s\n", newtypename);
06812           NAMD_die(err_msg);
06813         }
06814       } else {
06815         NAMD_die("Table interpolation type must be linear or cubic\n");
06816       }
06817 
06818       numtypesread++;
06819       continue;
06820     }
06821     //char err_msg[512];
06822     //sprintf(err_msg, "Unrecognized line in energy table file: %s\n", tableline);
06823     //NAMD_die(err_msg);
06824   }
06825 
06826   // See if we got what we expected
06827   if (numtypesread != numtypes) {
06828     char err_msg[512];
06829     sprintf(err_msg, "ERROR: Expected %i tabulated energy types but got %i\n", numtypes, numtypesread);
06830     NAMD_die(err_msg);
06831   }
06832 
06833   // Move the necessary information into simParams
06834   simParams->tableNumTypes = numtypes;
06835   simParams->tableSpacing = table_spacing;
06836   simParams->tableMaxDist = maxdist;
06837   tablenumtypes = numtypes;
06838 
06839   /*
06840 xxxxxx
06841         int numtypes = simParams->tableNumTypes;
06842 
06843         //This parameter controls the table spacing
06844         BigReal table_spacing = 0.01;
06845         BigReal maxdist = 20.0;
06846         if (maxdist > simParams->cutoff) {
06847                 iout << "Info: Reducing max table size to match nonbond cutoff\n";
06848                 maxdist = ceil(simParams->cutoff);
06849         }
06850 
06851         numenerentries = (numtypes + 1) * numtypes * int (ceil(maxdist/table_spacing));
06852 //      fprintf(stdout, "Table arithmetic: maxdist %f, table_spacing %f, numtypes %i, numentries %i\n", maxdist, table_spacing, numtypes, numenerentries);
06853         columnsize = 2 * namdnearbyint(maxdist/table_spacing);
06854         rowsize = numtypes * columnsize;
06855         //Set up a full energy lookup table from a file
06856         //Allocate the table; layout is atom1 x atom2 x distance energy force
06857         fprintf(stdout, "Table has %i entries\n",numenerentries);
06858         //iout << "Allocating original energy table\n" << endi;
06859         if ( table_ener ) delete [] table_ener;
06860         table_ener = new Real[numenerentries];
06861         //
06862         //Set sentinel values for uninitialized table energies
06863         for (int i=0 ; i<numenerentries ; i++) {
06864                 table_ener[i] = 1337.0;
06865         }
06866         Real compval = 1337.0;
06867 
06868         //    iout << "Entering table reading\n" << endi;
06869         //iout << "Finished allocating table\n" << endi;
06870 
06871         //Counter to make sure we read the right number of energy entries
06872         int numentries = 0;
06873 
06874         //Now, start reading from the file
06875         char* table_file = simParams->tabulatedEnergiesFile;
06876         FILE* enertable;
06877         enertable = fopen(table_file, "r");
06878 
06879         if (enertable == NULL) {
06880                 NAMD_die("ERROR: Could not open energy table file!\n");
06881         }
06882 
06883         char tableline[121];
06884         int atom1;
06885         int atom2;
06886         int distbin;
06887         Real dist;
06888         Real energy;
06889         Real force;
06890 
06891         iout << "Beginning table read\n" << endi;
06892         //Read the table entries
06893         while(fgets(tableline,120,enertable)) {
06894 //              IOut << "Processing line " << tableline << "\n" << endi;
06895                 if (strncmp(tableline,"#",1)==0) {
06896                         continue;
06897                 }
06898 
06899 
06900                 sscanf(tableline, "%i %i %f %f %f\n", &atom1, &atom2, &dist, &energy, &force);
06901                 distbin = int(namdnearbyint(dist/table_spacing));
06902 //              fprintf(stdout, "Working on atoms %i and %i at distance %f\n",atom1,atom2,dist);
06903                 if ((atom2 > atom1) || (distbin > int(namdnearbyint(maxdist/table_spacing)))) {
06904 //                      fprintf(stdout,"Rejected\n");
06905 //                      fprintf(stdout, "Error: Throwing out energy line beyond bounds\n");
06906         //              fprintf(stdout, "Bad line: Atom 1: %i Atom 2: %i Distance Bin: %i Max Distance Bin: %i \n", atom1, atom2, distbin, int(namdnearbyint(maxdist/table_spacing)));
06907                 } else {
06908                         //The magic formula for the number of columns from previous rows
06909                         //in the triangular matrix is (2ni+i-i**2)/2
06910                         //Here n is the number of types, and i is atom2
06911 //                      fprintf(stdout, "Input: atom1 %f atom2 %f columnsize %f \n", float(atom1), float(atom2), float(columnsize));
06912 //                      fprintf(stdout, "Testing arithmetic: Part1: %i Part2: %i Part3: %i Total: %i\n", columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2), columnsize*(atom1-atom2), 2*distbin, columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2) + columnsize*(atom1-atom2) + 2*distbin - 2);
06913                         int eneraddress = columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2) + columnsize*(atom1-atom2) + 2*distbin - 2;
06914                         int forceaddress = eneraddress + 1;
06915 //                              fprintf(stdout, "Tableloc: %i Atom 1: %i Atom 2: %i Distance Bin: %i Energy: %f Force: %f\n", eneraddress, atom1, atom2, distbin, table_ener[eneraddress], table_ener[forceaddress]);
06916                 fflush(stdout);
06917 //                      fprintf(stdout, "Checking for dupes: Looking at: %f %f \n", table_ener[eneraddress], table_ener[forceaddress]);
06918                         if ((table_ener[eneraddress] == compval && table_ener[forceaddress] == compval)) {
06919                                 numentries++;
06920                                 table_ener[eneraddress] = energy;
06921                                 table_ener[forceaddress] = force;
06922 //                              fprintf(stdout, "Tableloc: %i Atom 1: %i Atom 2: %i Distance Bin: %i Energy: %f Force: %f\n", eneraddress, atom1, atom2, distbin, table_ener[eneraddress], table_ener[forceaddress]);
06923                                 //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin] = energy;
06924                                 //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin + 1] = force;
06925 //                              fflush(stdout);
06926                         } else {
06927 //                              fprintf(stdout,"Rejecting duplicate entry\n");
06928                         }
06929                 }
06930                 //      iout << "Done with line\n"<< endi;
06931         }
06932 
06933         //    int should = numtypes * numtypes * (maxdist/table_spacing);
06934         //    cout << "Entries: " << numentries << " should be " << should << "\n" << endi;
06935 //      int exptypes = ceil((numtypes+1) * numtypes * (maxdist/table_spacing));
06936 //fprintf(stdout,"numtypes: %i maxdist: %f table_spacing: %f exptypes: %i\n",numtypes,maxdist,table_spacing);
06937         if (numentries != int(numenerentries/2)) {
06938                 fprintf(stdout,"ERROR: Expected %i entries but got %i\n",numenerentries/2, numentries);
06939                 NAMD_die("Number of energy table entries did not match expected value\n");
06940         }
06941         //      iout << "Done with table\n"<< endi;
06942         fprintf(stdout, "Numenerentries: %i\n",numenerentries/2);
06943   */
06944 } /* END of function read_ener_table */ 

int Parameters::read_energy_type ( FILE *  ,
const   int,
BigReal ,
const   float,
const   float 
)

Definition at line 7442 of file Parameters.C.

References NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

07442                                                                                                                                           {
07443 
07444   char tableline[120];
07445 
07446   /* Last values read from table */
07447   BigReal readdist;
07448   BigReal readener;
07449   BigReal readforce;
07450   BigReal lastdist;
07451   BigReal lastener;
07452   BigReal lastforce;
07453   readdist = -1.0;
07454   readener = 0.0;
07455   readforce = 0.0;
07456 
07457   /* Keep track of where in the table we are */
07458   float currdist;
07459   int distbin;
07460   currdist = -1.0;
07461   distbin = -1;
07462 
07463         while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
07464     printf("At distance %f + %f vs. %f\n", currdist, table_spacing, maxdist);
07465                 if (strncmp(tableline,"#",1)==0) {
07466                         continue;
07467                 }
07468     if (strncmp(tableline,"TYPE",4)==0) {
07469       break;
07470     }
07471 
07472     // Read an energy line from the table
07473     lastdist = readdist;
07474     lastener = readener;
07475     lastforce = readforce;
07476     int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
07477     if (distbin == -1) {
07478       if (readdist != 0.0) {
07479         NAMD_die("ERROR: Energy/force tables must start at d=0.0\n");
07480       } else {
07481         distbin = 0;
07482         continue;
07483       }
07484     }
07485    // printf("Read an energy line: %f %f %f\n", readdist, readener, readforce);
07486     if (readcount != 3) {
07487       char err_msg[512];
07488       sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
07489       NAMD_die(err_msg);
07490     }
07491 
07492     //Sanity check the current entry
07493     if (readdist < lastdist) {
07494       NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
07495     }
07496 
07497     currdist = lastdist;
07498 
07499     while (currdist <= readdist && distbin <= (int) (namdnearbyint(maxdist / table_spacing))) {
07500       distbin = (int) (namdnearbyint(currdist / table_spacing));
07501       int table_loc = 2 * (distbin + (typeindex * (1 + (int) namdnearbyint(maxdist / table_spacing))));
07502       printf("Doing interpolation for energy between %f %f and %f %f: Dist %f\n", readener, readdist, lastener, lastdist, currdist);
07503       table_ener[table_loc] = interp_lin(readener, lastener, readdist, lastdist, currdist);
07504       table_ener[table_loc + 1] = interp_lin(readforce, lastforce, readdist, lastdist, currdist);
07505       printf("Adding energy/force entry: %f %f in distbin %i (distance %f) to address %i/%i\n", table_ener[table_loc], table_ener[table_loc + 1], distbin, currdist, table_loc, table_loc + 1);
07506       currdist += table_spacing;
07507       distbin++;
07508     }
07509   }
07510 
07511   // Go back one line, since we should now be into the next TYPE block
07512   fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR); 
07513 
07514   // Clean up and make sure everything worked ok
07515   distbin--;
07516   printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
07517   if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
07518   return 0;
07519 }

int Parameters::read_energy_type_bothcubspline ( FILE *  ,
const   int,
BigReal ,
const   float,
const   float 
)

Definition at line 6967 of file Parameters.C.

References INDEX, NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

06967                                                                                                                                                         {
06968 
06969   char tableline[120];
06970   int i,j;
06971 
06972   /* Last values read from table */
06973   BigReal readdist;
06974   BigReal readener;
06975   BigReal readforce;
06976   BigReal spacing;
06977 //  BigReal readforce;
06978   BigReal lastdist;
06979 //  BigReal lastener;
06980 //  BigReal lastforce;
06981 //  readdist = -1.0;
06982 //  readener = 0.0;
06983 //  readforce = 0.0;
06984 
06985   /* Create arrays for holding the input values */
06986   std::vector<BigReal>  dists;
06987   std::vector<BigReal> enervalues;
06988   std::vector<BigReal> forcevalues;
06989   int numentries = 0;
06990 
06991 
06992   /* Keep track of where in the table we are */
06993   BigReal currdist;
06994   int distbin;
06995   currdist = 0.0;
06996   lastdist = -1.0;
06997   distbin = 0;
06998 
06999   // Read all of the values first -- we'll interpolate later
07000         while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
07001                 if (strncmp(tableline,"#",1)==0) {
07002                         continue;
07003                 }
07004     if (strncmp(tableline,"TYPE",4)==0) {
07005       fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR); 
07006       break;
07007     }
07008 
07009     // Read an energy line from the table
07010     int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
07011 
07012     //printf("Read an energy line: %g %g %g\n", readdist, readener, readforce);
07013     if (readcount != 3) {
07014       char err_msg[512];
07015       sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
07016       NAMD_die(err_msg);
07017     }
07018 
07019     //Sanity check the current entry
07020     if (readdist < lastdist) {
07021       NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
07022     }
07023 
07024     lastdist = readdist;
07025     dists.push_back(readdist);
07026     enervalues.push_back(readener);
07027     forcevalues.push_back(readforce);
07028     numentries++;
07029   }
07030 
07031   // Check the spacing and make sure it is uniform
07032   if (dists[0] != 0.0) {
07033     NAMD_die("ERROR: First data point for energy table must be at r=0\n");
07034   }
07035   spacing = dists[1] - dists[0];
07036   for (i=2; i<(numentries - 1); i++) {
07037     BigReal myspacing;
07038     myspacing = dists[i] - dists[i-1];
07039     if (fabs(myspacing - spacing) > 0.00001) {
07040       printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
07041       NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
07042     }
07043   }
07044 
07045 // Perform cubic spline interpolation to get the energies and forces
07046 
07047   /* allocate spline coefficient matrix */
07048   // xe and xf / be and bf for energies and forces, respectively
07049   double* m = new double[numentries*numentries];
07050   double* xe = new double[numentries];
07051   double* xf = new double[numentries];
07052   double* be = new double[numentries];
07053   double* bf = new double[numentries];
07054 
07055   be[0] = 3 * (enervalues[1] - enervalues[0]);
07056   for (i=1; i < (numentries - 1); i++) {
07057 //    printf("Control point %i at %f\n", i, enervalues[i]);
07058     be[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
07059 //    printf("be is %f\n", be[i]);
07060   }
07061   be[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
07062 
07063   bf[0] = 3 * (forcevalues[1] - forcevalues[0]);
07064   for (i=1; i < (numentries - 1); i++) {
07065 //    printf("Control point %i at %f\n", i, forcevalues[i]);
07066     bf[i] = 3 * (forcevalues[i+1] - forcevalues[i-1]);
07067 //    printf("bf is %f\n", bf[i]);
07068   }
07069   bf[numentries - 1] = 3 * (forcevalues[numentries - 1] - forcevalues[numentries - 2]);
07070 
07071   memset(m,0,numentries*numentries*sizeof(double));
07072 
07073   /* initialize spline coefficient matrix */
07074   m[0] = 2;
07075   for (i = 1;  i < numentries;  i++) {
07076     m[INDEX(numentries,i-1,i)] = 1;
07077     m[INDEX(numentries,i,i-1)] = 1;
07078     m[INDEX(numentries,i,i)] = 4;
07079   }
07080   m[INDEX(numentries,numentries-1,numentries-1)] = 2;
07081 
07082   /* Now we need to solve the equation M X = b for X */
07083   // Do this simultaneously for energy and force -- ONLY because we use the same matrix
07084 
07085   //Put m in upper triangular form and apply corresponding operations to b
07086   for (i=0; i<numentries; i++) {
07087     // zero the ith column in all rows below i
07088     const BigReal baseval = m[INDEX(numentries,i,i)];
07089     for (j=i+1; j<numentries; j++) {
07090       const BigReal myval = m[INDEX(numentries,j,i)];
07091       const BigReal factor = myval / baseval;
07092 
07093       for (int k=i; k<numentries; k++) {
07094         const BigReal subval = m[INDEX(numentries,i,k)];
07095         m[INDEX(numentries,j,k)] -= (factor * subval);
07096       }
07097 
07098       be[j] -= (factor * be[i]);
07099       bf[j] -= (factor * bf[i]);
07100 
07101     }
07102   }
07103 
07104   //Now work our way diagonally up from the bottom right to assign values to X
07105   for (i=numentries-1; i>=0; i--) {
07106 
07107     //Subtract the effects of previous columns
07108     for (j=i+1; j<numentries; j++) {
07109       be[i] -= ( m[INDEX(numentries,i,j)] * xe[j] );
07110       bf[i] -= ( m[INDEX(numentries,i,j)] * xf[j] );
07111     }
07112 
07113     xe[i] = be[i] / m[INDEX(numentries,i,i)];
07114     xf[i] = bf[i] / m[INDEX(numentries,i,i)];
07115 
07116   }
07117 
07118   // We now have the coefficient information we need to make the table
07119   // Now interpolate on each interval we want
07120 
07121   distbin = 0;
07122   int entriesperseg = (int) ceil(spacing / table_spacing);
07123   int table_prefix = 2 * typeindex * (int) (namdnearbyint(maxdist / table_spacing) + 1);
07124 
07125   for (i=0; i<numentries-1; i++) {
07126     BigReal Ae,Be,Ce,De;
07127     BigReal Af,Bf,Cf,Df;
07128     currdist = dists[i];
07129 
07130 //    printf("Interpolating on interval %i\n", i);
07131 
07132     // Set up the coefficients for this section
07133     Ae = enervalues[i];
07134     Be = xe[i];
07135     Ce = 3 * (enervalues[i+1] - enervalues[i]) - (2 * xe[i]) - xe[i+1];
07136     De = 2 * (enervalues[i] - enervalues[i+1]) + xe[i] + xe[i+1];
07137 
07138     Af = forcevalues[i];
07139     Bf = xf[i];
07140     Cf = 3 * (forcevalues[i+1] - forcevalues[i]) - (2 * xf[i]) - xf[i+1];
07141     Df = 2 * (forcevalues[i] - forcevalues[i+1]) + xf[i] + xf[i+1];
07142 
07143     // Go over the region of interest and fill in the table
07144     for (j=0; j<entriesperseg; j++) {
07145       const BigReal mydist = currdist + (j * table_spacing);
07146       const BigReal mydistfrac = (float) j / (entriesperseg - 1);
07147       distbin = (int) namdnearbyint(mydist / table_spacing);
07148       if (distbin >= (int) namdnearbyint(maxdist / table_spacing)) break;
07149       BigReal energy;
07150       BigReal force;
07151 
07152       energy = Ae + (Be * mydistfrac) + (Ce * mydistfrac * mydistfrac) + (De * mydistfrac * mydistfrac * mydistfrac);
07153       force = Af + (Bf * mydistfrac) + (Cf * mydistfrac * mydistfrac) + (Df * mydistfrac * mydistfrac * mydistfrac);
07154 
07155 //      printf("Adding energy/force entry %f / %f in bins %i / %i for distance %f (%f)\n", energy, force, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1), mydist, mydistfrac);
07156       table_ener[table_prefix + 2 * distbin] = energy;
07157       table_ener[table_prefix + 2 * distbin + 1] = force;
07158       distbin++;
07159     }
07160     if (currdist >= maxdist) break;
07161   }
07162 
07163   //The procedure above leaves out the last entry -- add it explicitly
07164   distbin = (int) namdnearbyint(maxdist / table_spacing);
07165 //  printf("Adding energy/force entry %f / %f in bins %i / %i\n", enervalues[numentries - 1], 0.0, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1));
07166   table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
07167   table_ener[table_prefix + 2 * distbin + 1] = 0.0;
07168   distbin++;
07169 
07170 
07171   // Clean up and make sure everything worked ok
07172   delete m;
07173   delete xe;
07174   delete xf;
07175   delete be;
07176   delete bf;
07177   distbin--;
07178   printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
07179   if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
07180   return 0;
07181 } /* end read_energy_type_bothcubspline */

int Parameters::read_energy_type_cubspline ( FILE *  ,
const   int,
BigReal ,
const   float,
const   float 
)

Definition at line 7203 of file Parameters.C.

References INDEX, NAMD_die(), namdnearbyint, and x.

07203                                                                                                                                                     {
07204 
07205   char tableline[120];
07206   int i,j;
07207 
07208   /* Last values read from table */
07209   BigReal readdist;
07210   BigReal readener;
07211   BigReal spacing;
07212 //  BigReal readforce;
07213   BigReal lastdist;
07214 //  BigReal lastener;
07215 //  BigReal lastforce;
07216 //  readdist = -1.0;
07217 //  readener = 0.0;
07218 //  readforce = 0.0;
07219 
07220   /* Create arrays for holding the input values */
07221   std::vector<BigReal>  dists;
07222   std::vector<BigReal> enervalues;
07223   int numentries = 0;
07224 
07225 
07226   /* Keep track of where in the table we are */
07227   BigReal currdist;
07228   int distbin;
07229   currdist = 0.0;
07230   lastdist = -1.0;
07231   distbin = 0;
07232 
07233   // Read all of the values first -- we'll interpolate later
07234         while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
07235                 if (strncmp(tableline,"#",1)==0) {
07236                         continue;
07237                 }
07238     if (strncmp(tableline,"TYPE",4)==0) {
07239       fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR); 
07240       break;
07241     }
07242 
07243     // Read an energy line from the table
07244     int readcount = sscanf(tableline, "%lf %lf", &readdist, &readener);
07245 
07246    // printf("Read an energy line: %f %f %f\n", readdist, readener, readforce);
07247     if (readcount != 2) {
07248       char err_msg[512];
07249       sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
07250       NAMD_die(err_msg);
07251     }
07252 
07253     //Sanity check the current entry
07254     if (readdist < lastdist) {
07255       NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
07256     }
07257 
07258     lastdist = readdist;
07259     dists.push_back(readdist);
07260     enervalues.push_back(readener);
07261     numentries++;
07262   }
07263 
07264   // Check the spacing and make sure it is uniform
07265   if (dists[0] != 0.0) {
07266     NAMD_die("ERROR: First data point for energy table must be at r=0\n");
07267   }
07268   spacing = dists[1] - dists[0];
07269   for (i=2; i<(numentries - 1); i++) {
07270     BigReal myspacing;
07271     myspacing = dists[i] - dists[i-1];
07272     if (fabs(myspacing - spacing) > 0.00001) {
07273       printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
07274       NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
07275     }
07276   }
07277 
07278 // Perform cubic spline interpolation to get the energies and forces
07279 
07280   /* allocate spline coefficient matrix */
07281   double* m = new double[numentries*numentries];
07282   double* x = new double[numentries];
07283   double* b = new double[numentries];
07284 
07285   b[0] = 3 * (enervalues[1] - enervalues[0]);
07286   for (i=1; i < (numentries - 1); i++) {
07287     printf("Control point %i at %f\n", i, enervalues[i]);
07288     b[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
07289     printf("b is %f\n", b[i]);
07290   }
07291   b[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
07292 
07293   memset(m,0,numentries*numentries*sizeof(double));
07294 
07295   /* initialize spline coefficient matrix */
07296   m[0] = 2;
07297   for (i = 1;  i < numentries;  i++) {
07298     m[INDEX(numentries,i-1,i)] = 1;
07299     m[INDEX(numentries,i,i-1)] = 1;
07300     m[INDEX(numentries,i,i)] = 4;
07301   }
07302   m[INDEX(numentries,numentries-1,numentries-1)] = 2;
07303 
07304   /* Now we need to solve the equation M X = b for X */
07305 
07306   printf("Solving the matrix equation: \n");
07307 
07308   for (i=0; i<numentries; i++) {
07309     printf(" ( ");
07310     for (j=0; j<numentries; j++) {
07311       printf(" %6.3f,", m[INDEX(numentries, i, j)]);
07312     }
07313     printf(" )  ( D%-3i )  =  ( %6.3f )\n", i, b[i]);
07314   }
07315 
07316   //Put m in upper triangular form and apply corresponding operations to b
07317   for (i=0; i<numentries; i++) {
07318     // zero the ith column in all rows below i
07319     const BigReal baseval = m[INDEX(numentries,i,i)];
07320     for (j=i+1; j<numentries; j++) {
07321       const BigReal myval = m[INDEX(numentries,j,i)];
07322       const BigReal factor = myval / baseval;
07323 
07324       for (int k=i; k<numentries; k++) {
07325         const BigReal subval = m[INDEX(numentries,i,k)];
07326         m[INDEX(numentries,j,k)] -= (factor * subval);
07327       }
07328 
07329       b[j] -= (factor * b[i]);
07330 
07331     }
07332   }
07333 
07334   printf(" In upper diagonal form, equation is:\n");
07335   for (i=0; i<numentries; i++) {
07336     printf(" ( ");
07337     for (j=0; j<numentries; j++) {
07338       printf(" %6.3f,", m[INDEX(numentries, i, j)]);
07339     }
07340     printf(" )  ( D%-3i )  =  ( %6.3f )\n", i, b[i]);
07341   }
07342 
07343   //Now work our way diagonally up from the bottom right to assign values to X
07344   for (i=numentries-1; i>=0; i--) {
07345 
07346     //Subtract the effects of previous columns
07347     for (j=i+1; j<numentries; j++) {
07348       b[i] -= ( m[INDEX(numentries,i,j)] * x[j] );
07349     }
07350 
07351     x[i] = b[i] / m[INDEX(numentries,i,i)];
07352 
07353   }
07354 
07355   printf(" Solution vector is:\n\t(");
07356   for (i=0; i<numentries; i++) {
07357     printf(" %6.3f ", x[i]);
07358   }
07359   printf(" ) \n");
07360 
07361   // We now have the coefficient information we need to make the table
07362   // Now interpolate on each interval we want
07363 
07364   distbin = 0;
07365   int entriesperseg = (int) ceil(spacing / table_spacing);
07366   int table_prefix = 2 * typeindex * (int) (namdnearbyint(maxdist / table_spacing) + 1);
07367 
07368   for (i=0; i<numentries-1; i++) {
07369     BigReal A,B,C,D;
07370     currdist = dists[i];
07371 
07372     printf("Interpolating on interval %i\n", i);
07373 
07374     // Set up the coefficients for this section
07375     A = enervalues[i];
07376     B = x[i];
07377     C = 3 * (enervalues[i+1] - enervalues[i]) - (2 * x[i]) - x[i+1];
07378     D = 2 * (enervalues[i] - enervalues[i+1]) + x[i] + x[i+1];
07379 
07380     printf("Coefficients for this interval: %f %f %f %f\n", A, B, C, D);
07381 
07382     // Go over the region of interest and fill in the table
07383     for (j=0; j<entriesperseg; j++) {
07384       const BigReal mydist = currdist + (j * table_spacing);
07385       const BigReal mydistfrac = (float) j / (entriesperseg - 1);
07386       distbin = (int) namdnearbyint(mydist / table_spacing);
07387       if (distbin >= (int) namdnearbyint(maxdist / table_spacing)) break;
07388       BigReal energy;
07389       BigReal force;
07390 
07391       energy = A + (B * mydistfrac) + (C * mydistfrac * mydistfrac) + (D * mydistfrac * mydistfrac * mydistfrac);
07392       force = B + (2 * C * mydistfrac) + (3 * D * mydistfrac * mydistfrac);
07393       // Multiply force by 1 / (interval length)
07394       force *= (1.0 / spacing);
07395 
07396       printf("Adding energy/force entry %f / %f in bins %i / %i for distance %f (%f)\n", energy, force, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1), mydist, mydistfrac);
07397       table_ener[table_prefix + 2 * distbin] = energy;
07398       table_ener[table_prefix + 2 * distbin + 1] = force;
07399       distbin++;
07400     }
07401     if (currdist >= maxdist) break;
07402   }
07403 
07404   //The procedure above leaves out the last entry -- add it explicitly
07405   distbin = (int) namdnearbyint(maxdist / table_spacing);
07406   printf("Adding energy/force entry %f / %f in bins %i / %i\n", enervalues[numentries - 1], 0.0, (table_prefix + 2 * distbin), (table_prefix + 2 * distbin + 1));
07407   table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
07408   table_ener[table_prefix + 2 * distbin + 1] = 0.0;
07409   distbin++;
07410 
07411 
07412   // Clean up and make sure everything worked ok
07413   delete m;
07414   delete x;
07415   delete b;
07416   distbin--;
07417   printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
07418   if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
07419   return 0;
07420 } /* end read_energy_type_cubspline */

void Parameters::read_parameter_file ( char *   ) 

Definition at line 404 of file Parameters.C.

References Fclose(), Fopen(), NAMD_blank_string(), NAMD_die(), NAMD_find_first_word(), NAMD_read_line(), NumAngleParams, NumBondParams, NumDihedralParams, NumImproperParams, NumTablePairParams, NumVdwPairParams, NumVdwParams, and table_ener.

Referenced by Parameters().

00406 {
00407   char buffer[512];  //  Buffer to store each line of the file
00408   char first_word[512];  //  First word of the current line
00409   FILE *pfile;    //  File descriptor for the parameter file
00410 
00411   /*  Check to make sure that we haven't previously been told     */
00412   /*  that all the files were read        */
00413   if (AllFilesRead)
00414   {
00415     NAMD_die("Tried to read another parameter file after being told that all files were read!");
00416   }
00417 
00418   /*  Try and open the file          */
00419   if ( (pfile = Fopen(fname, "r")) == NULL)
00420   {
00421     char err_msg[256];
00422 
00423     sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
00424     NAMD_die(err_msg);
00425   }
00426 
00427   /*  Keep reading in lines until we hit the EOF      */
00428   while (NAMD_read_line(pfile, buffer) != -1)
00429   {
00430     /*  Get the first word of the line      */
00431     NAMD_find_first_word(buffer, first_word);
00432 
00433     /*  First, screen out things that we ignore, such as    */
00434     /*  blank lines, lines that start with '!', lines that  */
00435     /*  start with "REMARK", lines that start with set",    */
00436     /*  and most of the HBOND parameters which include      */
00437     /*  AEXP, REXP, HAEX, AAEX, but not the HBOND statement */
00438     /*  which is parsed.                                    */
00439     if ((buffer[0] != '!') && 
00440         !NAMD_blank_string(buffer) &&
00441         (strncasecmp(first_word, "REMARK", 6) != 0) &&
00442         (strcasecmp(first_word, "set")!=0) &&
00443         (strncasecmp(first_word, "AEXP", 4) != 0) &&
00444         (strncasecmp(first_word, "REXP", 4) != 0) &&
00445         (strncasecmp(first_word, "HAEX", 4) != 0) &&
00446         (strncasecmp(first_word, "AAEX", 4) != 0) &&
00447         (strncasecmp(first_word, "NBOND", 5) != 0) &&
00448         (strncasecmp(first_word, "CUTNB", 5) != 0) &&
00449         (strncasecmp(first_word, "END", 3) != 0) &&
00450         (strncasecmp(first_word, "CTONN", 5) != 0) &&
00451         (strncasecmp(first_word, "EPS", 3) != 0) &&
00452         (strncasecmp(first_word, "VSWI", 4) != 0) &&
00453         (strncasecmp(first_word, "NBXM", 4) != 0) &&
00454         (strncasecmp(first_word, "INHI", 4) != 0) )
00455     {
00456       /*  Now, call the appropriate function based    */
00457       /*  on the type of parameter we have    */
00458       if (strncasecmp(first_word, "bond", 4)==0)
00459       {
00460         add_bond_param(buffer);
00461         NumBondParams++;
00462       }
00463       else if (strncasecmp(first_word, "angl", 4)==0)
00464       {
00465         add_angle_param(buffer);
00466         NumAngleParams++;
00467       }
00468       else if (strncasecmp(first_word, "dihe", 4)==0)
00469       {
00470         add_dihedral_param(buffer, pfile);
00471         NumDihedralParams++;
00472       }
00473       else if (strncasecmp(first_word, "impr", 4)==0)
00474       {
00475         add_improper_param(buffer, pfile);
00476         NumImproperParams++;
00477       }
00478       else if (strncasecmp(first_word, "nonb", 4)==0)
00479       {
00480         add_vdw_param(buffer);
00481         NumVdwParams++; 
00482       }
00483       else if (strncasecmp(first_word, "nbfi", 4)==0)
00484       {
00485         add_vdw_pair_param(buffer);
00486         NumVdwPairParams++; 
00487       }
00488       else if (strncasecmp(first_word, "nbta", 4)==0)
00489       {
00490 
00491         if (table_ener == NULL) {
00492           continue;
00493         }
00494 
00495         add_table_pair_param(buffer);
00496         NumTablePairParams++; 
00497       }
00498       else if (strncasecmp(first_word, "hbon", 4)==0)
00499       {
00500         add_hb_pair_param(buffer);
00501       }
00502       else
00503       {
00504         /*  This is an unknown paramter.        */
00505         /*  This is BAD        */
00506         char err_msg[512];
00507 
00508         sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
00509            fname, buffer);
00510         NAMD_die(err_msg);
00511       }
00512     }
00513   }
00514 
00515   /*  Close the file            */
00516   Fclose(pfile);
00517 
00518   return;
00519 }

void Parameters::read_parm ( Ambertoppar ,
BigReal   
)

Definition at line 6362 of file Parameters.C.

References indexed_vdw_pair::A, indexed_vdw_pair::A14, angle_array, indexed_vdw_pair::B, indexed_vdw_pair::B14, bond_array, parm::Cn1, parm::Cn2, parm::Cno, parm::data_read, four_body_consts::delta, dihedral_array, parm::HB12, parm::HB6, if(), indexed_vdw_pair::ind1, indexed_vdw_pair::ind2, iout, iWARN(), AngleValue::k_ub, indexed_vdw_pair::left, MAX_MULTIPLICITY, DihedralValue::multiplicity, four_body_consts::n, NAMD_die(), AngleValue::normal, parm::Nptra, parm::Ntypes, parm::Numang, NumAngleParams, parm::Numbnd, NumBondParams, NumDihedralParams, NumVdwPairParams, NumVdwParamsAssigned, parm::Phase, parm::Pk, parm::Pn, AngleValue::r_ub, parm::Req, indexed_vdw_pair::right, parm::Rk, parm::Teq, AngleValue::theta0, parm::Tk, DihedralValue::values, vdw_pair_tree, and BondValue::x0.

06363 { 
06364   int i,j,ico,numtype,mul;
06365   IndexedVdwPair *new_node;
06366 
06367   if (!amber_data->data_read)
06368     NAMD_die("No data read from parm file yet!");
06369 
06370   // Copy bond parameters
06371   NumBondParams = amber_data->Numbnd;
06372   if (NumBondParams)
06373   { bond_array = new BondValue[NumBondParams];
06374     if (bond_array == NULL)
06375       NAMD_die("memory allocation of bond_array failed!");
06376   }
06377   for (i=0; i<NumBondParams; ++i)
06378   { bond_array[i].k = amber_data->Rk[i];
06379     bond_array[i].x0 = amber_data->Req[i];
06380   }
06381 
06382   // Copy angle parameters
06383   NumAngleParams = amber_data->Numang;
06384   if (NumAngleParams)
06385   { angle_array = new AngleValue[NumAngleParams];
06386     if (angle_array == NULL)
06387       NAMD_die("memory allocation of angle_array failed!");
06388   }
06389   for (i=0; i<NumAngleParams; ++i)
06390   { angle_array[i].k = amber_data->Tk[i];
06391     angle_array[i].theta0 = amber_data->Teq[i];
06392     // Amber has no k_ub and r_ub for angle parameters, so they're set to 0
06393     angle_array[i].k_ub = angle_array[i].r_ub = 0;
06394     // All angles are harmonic
06395     angle_array[i].normal = 1;
06396   }
06397 
06398   // Copy dihedral parameters
06399   // Note: If the periodicity is negative, it means the following
06400   //  entry is another term in a multitermed dihedral, and the
06401   //  absolute value is the true periodicity; in this case the
06402   //  following entry in "dihedral_array" should be left empty,
06403   //  NOT be skipped, in order to keep the next dihedral's index
06404   //  unchanged.
06405   NumDihedralParams = amber_data->Nptra;
06406   if (NumDihedralParams)
06407   { dihedral_array = new DihedralValue[amber_data->Nptra];
06408     if (dihedral_array == NULL)
06409       NAMD_die("memory allocation of dihedral_array failed!");
06410   }
06411   mul = 0;
06412   for (i=0; i<NumDihedralParams; ++i)
06413   { dihedral_array[i-mul].values[mul].k = amber_data->Pk[i];
06414     dihedral_array[i-mul].values[mul].n = int(fabs(amber_data->Pn[i])+0.5);
06415     if (dihedral_array[i-mul].values[mul].n == 0)
06416     { char err_msg[128];
06417       sprintf(err_msg, "The periodicity of dihedral # %d is zero!", i+1);
06418       NAMD_die(err_msg);
06419     }
06420     dihedral_array[i-mul].values[mul].delta = amber_data->Phase[i];
06421     // If the periodicity is positive, it means the following
06422     // entry is a new dihedral term.
06423     if (amber_data->Pn[i] > 0)
06424     { dihedral_array[i-mul].multiplicity = mul+1;
06425       mul = 0;
06426     }
06427     else if (++mul >= MAX_MULTIPLICITY)
06428     { char err_msg[181];
06429       sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
06430          mul+1, MAX_MULTIPLICITY);
06431       NAMD_die(err_msg);
06432     }
06433   }
06434   if (mul > 0)
06435     NAMD_die("Negative periodicity in the last dihedral entry!");
06436 
06437   // Copy VDW parameters: AMBER explicitly gives VDW parameters between every
06438   // 2 atom types, so the data are copied to vdw_pair_tree
06439   // In AMBER, all 1-4 VDW interactions are divided by factor vdw14
06440   NumVdwParamsAssigned = numtype = amber_data->Ntypes; // Number of atom types
06441   NumVdwPairParams = numtype * (numtype+1) / 2;
06442   for (i=0; i<numtype; ++i)
06443     for (j=i; j<numtype; ++j)
06444     { new_node = new IndexedVdwPair;
06445       if (new_node == NULL)
06446         NAMD_die("memory allocation of vdw_pair_tree failed!");
06447       new_node->ind1 = i;
06448       new_node->ind2 = j;
06449       new_node->left = new_node->right = NULL;
06450       // ico is the index of interaction between atom type i and j into
06451       // the parameter arrays. If it's positive, the interaction is
06452       // 6-12 VDW; otherwise it's 10-12 H-bond interaction. NAMD doesn't
06453       // have 10-12 term, so if ico is negative, then the 10-12
06454       // coefficients must be 0, otherwise die.
06455       ico = amber_data->Cno[numtype*i+j];
06456       if (ico>0)
06457       { new_node->A = amber_data->Cn1[ico-1];
06458         new_node->A14 = new_node->A / vdw14;
06459         new_node->B = amber_data->Cn2[ico-1];
06460         new_node->B14 = new_node->B / vdw14;
06461       }
06462       else if (amber_data->HB12[abs(ico)-1]==0.0 && amber_data->HB6[abs(ico)-1]==0.0)
06463       { new_node->A = new_node->A14 = new_node->B = new_node->B14 = 0.0;
06464         iout << iWARN << "Encounter 10-12 H-bond term\n";
06465       }
06466       else
06467         NAMD_die("Encounter non-zero 10-12 H-bond term!");
06468       // Add the node to the binary tree
06469       vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
06470     }
06471 }

void Parameters::receive_Parameters ( MIStream  ) 

Definition at line 5420 of file Parameters.C.

References indexed_vdw_pair::A, indexed_vdw_pair::A14, nbthole_pair_value::alphai, nbthole_pair_value::alphaj, angle_array, indexed_vdw_pair::B, indexed_vdw_pair::B14, bond_array, crossterm_array, four_body_consts::delta, dihedral_array, CrosstermValue::dim, MIStream::get(), gromacsPair_array, improper_array, indexed_table_pair::ind1, indexed_vdw_pair::ind1, indexed_table_pair::ind2, nbthole_pair_value::ind2, indexed_vdw_pair::ind2, indexed_table_pair::K, AngleValue::k_ub, indexed_table_pair::left, indexed_vdw_pair::left, MAX_ATOMTYPE_CHARS, MAX_MULTIPLICITY, four_body_consts::n, NAMD_die(), nbthole_array, AngleValue::normal, NumAngleParams, NumBondParams, NumCrosstermParams, NumDihedralParams, numenerentries, NumGromacsPairParams, NumImproperParams, NumNbtholePairParams, NumTablePairParams, NumVdwPairParams, NumVdwParams, NumVdwParamsAssigned, GromacsPairValue::pairC12, AngleValue::r_ub, indexed_table_pair::right, indexed_vdw_pair::right, tab_pair_tree, table_ener, AngleValue::theta0, nbthole_pair_value::tholeij, TRUE, ImproperValue::values, DihedralValue::values, vdw_array, vdw_pair_tree, and BondValue::x0.

Referenced by Node::resendMolecule().

05422 {
05423   int i, j;      //  Loop counters
05424   Real *a1, *a2, *a3, *a4;  //  Temporary arrays to get data from message in
05425   int *i1, *i2, *i3;      //  Temporary int array to get data from message in
05426   IndexedVdwPair *new_node;  //  New node for vdw pair param tree
05427   IndexedTablePair *new_tab_node;
05428   Real **kvals;      //  Force constant values for dihedrals and impropers
05429   int **nvals;      //  Periodicity values for dihedrals and impropers
05430   Real **deltavals;    //  Phase shift values for dihedrals and impropers
05431   BigReal *pairC6, *pairC12;  // JLai
05432 
05433   //  Get the bonded parameters
05434   msg->get(NumBondParams);
05435 
05436   if (NumBondParams)
05437   {
05438     bond_array = new BondValue[NumBondParams];
05439     a1 = new Real[NumBondParams];
05440     a2 = new Real[NumBondParams];
05441 
05442     if ( (bond_array == NULL) || (a1 == NULL) || (a2 == NULL) )
05443     {
05444       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05445     }
05446 
05447     msg->get(NumBondParams, a1);
05448     msg->get(NumBondParams, a2);
05449 
05450     for (i=0; i<NumBondParams; i++)
05451     {
05452       bond_array[i].k = a1[i];
05453       bond_array[i].x0 = a2[i];
05454     }
05455 
05456     delete [] a1;
05457     delete [] a2;
05458   }
05459 
05460   //  Get the angle parameters
05461   msg->get(NumAngleParams);
05462 
05463   if (NumAngleParams)
05464   {
05465     angle_array = new AngleValue[NumAngleParams];
05466     a1 = new Real[NumAngleParams];
05467     a2 = new Real[NumAngleParams];
05468     a3 = new Real[NumAngleParams];
05469     a4 = new Real[NumAngleParams];
05470     i1 = new int[NumAngleParams];
05471 
05472     if ( (angle_array == NULL) || (a1 == NULL) || (a2 == NULL) ||
05473          (a3 == NULL) || (a4 == NULL) || (i1 == NULL))
05474     {
05475       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05476     }
05477 
05478     msg->get(NumAngleParams, a1);
05479     msg->get(NumAngleParams, a2);
05480     msg->get(NumAngleParams, a3);
05481     msg->get(NumAngleParams, a4);
05482     msg->get(NumAngleParams, i1);
05483 
05484     for (i=0; i<NumAngleParams; i++)
05485     {
05486       angle_array[i].k = a1[i];
05487       angle_array[i].theta0 = a2[i];
05488       angle_array[i].k_ub = a3[i];
05489       angle_array[i].r_ub = a4[i];
05490       angle_array[i].normal = i1[i];
05491     }
05492 
05493     delete [] a1;
05494     delete [] a2;
05495     delete [] a3;
05496     delete [] a4;
05497     delete [] i1;
05498   }
05499 
05500   //  Get the dihedral parameters
05501   msg->get(NumDihedralParams);
05502 
05503   if (NumDihedralParams)
05504   {
05505     dihedral_array = new DihedralValue[NumDihedralParams];
05506 
05507     i1 = new int[NumDihedralParams];
05508     kvals = new Real *[MAX_MULTIPLICITY];
05509     nvals = new int *[MAX_MULTIPLICITY];
05510     deltavals = new Real *[MAX_MULTIPLICITY];
05511 
05512     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
05513          (deltavals == NULL) || (dihedral_array == NULL) )
05514     {
05515       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05516     }
05517 
05518     for (i=0; i<MAX_MULTIPLICITY; i++)
05519     {
05520       kvals[i] = new Real[NumDihedralParams];
05521       nvals[i] = new int[NumDihedralParams];
05522       deltavals[i] = new Real[NumDihedralParams];
05523 
05524       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05525       {
05526         NAMD_die("memory allocation failed in Parameters::send_Parameters");
05527       }
05528     }
05529 
05530     msg->get(NumDihedralParams, i1);
05531 
05532     for (i=0; i<MAX_MULTIPLICITY; i++)
05533     {
05534       msg->get(NumDihedralParams, kvals[i]);
05535       msg->get(NumDihedralParams, nvals[i]);
05536       msg->get(NumDihedralParams, deltavals[i]);
05537     }
05538 
05539     for (i=0; i<NumDihedralParams; i++)
05540     {
05541       dihedral_array[i].multiplicity = i1[i];
05542 
05543       for (j=0; j<MAX_MULTIPLICITY; j++)
05544       {
05545         dihedral_array[i].values[j].k = kvals[j][i];
05546         dihedral_array[i].values[j].n = nvals[j][i];
05547         dihedral_array[i].values[j].delta = deltavals[j][i];
05548       }
05549     }
05550 
05551     for (i=0; i<MAX_MULTIPLICITY; i++)
05552     {
05553       delete [] kvals[i];
05554       delete [] nvals[i];
05555       delete [] deltavals[i];
05556     }
05557 
05558     delete [] i1;
05559     delete [] kvals;
05560     delete [] nvals;
05561     delete [] deltavals;
05562   }
05563 
05564   //  Get the improper parameters
05565   msg->get(NumImproperParams);
05566 
05567   if (NumImproperParams)
05568   {
05569     improper_array = new ImproperValue[NumImproperParams];
05570     i1 = new int[NumImproperParams];
05571     kvals = new Real *[MAX_MULTIPLICITY];
05572     nvals = new int *[MAX_MULTIPLICITY];
05573     deltavals = new Real *[MAX_MULTIPLICITY];
05574 
05575     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
05576          (deltavals == NULL) || (improper_array==NULL) )
05577     {
05578       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05579     }
05580 
05581     for (i=0; i<MAX_MULTIPLICITY; i++)
05582     {
05583       kvals[i] = new Real[NumImproperParams];
05584       nvals[i] = new int[NumImproperParams];
05585       deltavals[i] = new Real[NumImproperParams];
05586 
05587       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05588       {
05589         NAMD_die("memory allocation failed in Parameters::send_Parameters");
05590       }
05591     }
05592 
05593     msg->get(NumImproperParams,i1);
05594 
05595     for (i=0; i<MAX_MULTIPLICITY; i++)
05596     {
05597       msg->get(NumImproperParams,kvals[i]);
05598       msg->get(NumImproperParams,nvals[i]);
05599       msg->get(NumImproperParams,deltavals[i]);
05600     }
05601 
05602     for (i=0; i<NumImproperParams; i++)
05603     {
05604       improper_array[i].multiplicity = i1[i];
05605 
05606       for (j=0; j<MAX_MULTIPLICITY; j++)
05607       {
05608         improper_array[i].values[j].k = kvals[j][i];
05609         improper_array[i].values[j].n = nvals[j][i];
05610         improper_array[i].values[j].delta = deltavals[j][i];
05611       }
05612     }
05613 
05614     for (i=0; i<MAX_MULTIPLICITY; i++)
05615     {
05616       delete [] kvals[i];
05617       delete [] nvals[i];
05618       delete [] deltavals[i];
05619     }
05620 
05621     delete [] i1;
05622     delete [] kvals;
05623     delete [] nvals;
05624     delete [] deltavals;
05625   }
05626 
05627   //  Get the crossterm parameters
05628   msg->get(NumCrosstermParams);
05629 
05630   if (NumCrosstermParams)
05631   {
05632     crossterm_array = new CrosstermValue[NumCrosstermParams];
05633 
05634     for (i=0; i<NumCrosstermParams; ++i) {
05635       int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
05636       msg->get(nvals,&crossterm_array[i].c[0][0].d00);
05637     }
05638   }
05639   
05640   // Get GromacsPairs parameters
05641   // JLai
05642   msg->get(NumGromacsPairParams);
05643   
05644   if (NumGromacsPairParams)
05645   {
05646       gromacsPair_array = new GromacsPairValue[NumGromacsPairParams];
05647       pairC6 = new BigReal[NumGromacsPairParams];
05648       pairC12 = new BigReal[NumGromacsPairParams];
05649       
05650       if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
05651           NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05652       }
05653     
05654       msg->get(NumGromacsPairParams,pairC6);
05655       msg->get(NumGromacsPairParams,pairC12);
05656       
05657       for (i=0; i<NumGromacsPairParams; ++i) {
05658           gromacsPair_array[i].pairC6 = pairC6[i];
05659           gromacsPair_array[i].pairC12 = pairC12[i];
05660       }
05661 
05662       delete [] pairC6;
05663       delete [] pairC12;
05664   }
05665   // JLai
05666 
05667   //Get the energy table
05668   msg->get(numenerentries);
05669   if (numenerentries > 0) {
05670     //fprintf(stdout, "Getting tables\n");
05671     //fprintf(infofile, "Trying to allocate table\n");
05672     table_ener = new BigReal[numenerentries];
05673     //fprintf(infofile, "Finished table allocation\n");
05674     if (table_ener==NULL)
05675     {
05676       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05677     }
05678 
05679     msg->get(numenerentries, table_ener);
05680   }
05681 
05682   //  Get the vdw parameters
05683   msg->get(NumVdwParams);
05684   msg->get(NumVdwParamsAssigned);
05685 
05686   if (NumVdwParams)
05687   {
05688           atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
05689     vdw_array = new VdwValue[NumVdwParams];
05690     a1 = new Real[NumVdwParams];
05691     a2 = new Real[NumVdwParams];
05692     a3 = new Real[NumVdwParams];
05693     a4 = new Real[NumVdwParams];
05694 
05695     if ( (vdw_array==NULL) || (a1==NULL) || (a2==NULL) || (a3==NULL)
05696              || (a4==NULL) || (atomTypeNames==NULL))
05697     {
05698       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05699     }
05700 
05701     msg->get(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
05702     msg->get(NumVdwParams, a1);
05703     msg->get(NumVdwParams, a2);
05704     msg->get(NumVdwParams, a3);
05705     msg->get(NumVdwParams, a4);
05706 
05707     for (i=0; i<NumVdwParams; i++)
05708     {
05709       vdw_array[i].sigma = a1[i];
05710       vdw_array[i].epsilon = a2[i];
05711       vdw_array[i].sigma14 = a3[i];
05712       vdw_array[i].epsilon14 = a4[i];
05713     }
05714 
05715     delete [] a1;
05716     delete [] a2;
05717     delete [] a3;
05718     delete [] a4;
05719   }
05720   
05721   //  Get the vdw_pair_parameters
05722   msg->get(NumVdwPairParams);
05723   
05724   if (NumVdwPairParams)
05725   {
05726     a1 = new Real[NumVdwPairParams];
05727     a2 = new Real[NumVdwPairParams];
05728     a3 = new Real[NumVdwPairParams];
05729     a4 = new Real[NumVdwPairParams];
05730     i1 = new int[NumVdwPairParams];
05731     i2 = new int[NumVdwPairParams];
05732 
05733     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
05734          (i1 == NULL) || (i2 == NULL) )
05735     {
05736       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05737     }
05738     
05739     msg->get(NumVdwPairParams, i1);
05740     msg->get(NumVdwPairParams, i2);
05741     msg->get(NumVdwPairParams, a1);
05742     msg->get(NumVdwPairParams, a2);
05743     msg->get(NumVdwPairParams, a3);
05744     msg->get(NumVdwPairParams, a4);
05745     
05746     for (i=0; i<NumVdwPairParams; i++)
05747     {
05748       new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
05749       
05750       if (new_node == NULL)
05751       {
05752          NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05753       }
05754       
05755       new_node->ind1 = i1[i];
05756       new_node->ind2 = i2[i];
05757       new_node->A = a1[i];
05758       new_node->A14 = a2[i];
05759       new_node->B = a3[i];
05760       new_node->B14 = a4[i];
05761       new_node->left = NULL;
05762       new_node->right = NULL;
05763       
05764       vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
05765     }
05766     
05767     delete [] i1;
05768     delete [] i2;
05769     delete [] a1;
05770     delete [] a2;
05771     delete [] a3;
05772     delete [] a4;
05773   }
05774  
05775   //  Get the nbthole_pair_parameters
05776   msg->get(NumNbtholePairParams); 
05777     
05778   if (NumNbtholePairParams)
05779   {
05780     nbthole_array = new NbtholePairValue[NumNbtholePairParams];
05781     a1 = new Real[NumNbtholePairParams];
05782     a2 = new Real[NumNbtholePairParams];
05783     a3 = new Real[NumNbtholePairParams];
05784     i1 = new int[NumNbtholePairParams];
05785     i2 = new int[NumNbtholePairParams];
05786 
05787     if ( (nbthole_array == NULL) || (a1 == NULL) || (a2 == NULL) || (a3 == NULL)
05788          || (i1 == NULL) || (i2 == NULL) )
05789     {
05790       NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05791     }
05792 
05793     msg->get(NumNbtholePairParams, i1);
05794     msg->get(NumNbtholePairParams, i2);
05795     msg->get(NumNbtholePairParams, a1);
05796     msg->get(NumNbtholePairParams, a2);
05797     msg->get(NumNbtholePairParams, a3);
05798 
05799     for (i=0; i<NumNbtholePairParams; i++)
05800     {
05801 
05802       nbthole_array[i].ind1 = i1[i];
05803       nbthole_array[i].ind2 = i2[i];
05804       nbthole_array[i].alphai = a1[i];
05805       nbthole_array[i].alphaj = a2[i];
05806       nbthole_array[i].tholeij = a3[i];
05807 
05808     }
05809 
05810     delete [] i1;
05811     delete [] i2;
05812     delete [] a1;
05813     delete [] a2;
05814     delete [] a3;
05815   }
05816 
05817   //  Get the table_pair_parameters
05818   msg->get(NumTablePairParams);
05819   
05820   if (NumTablePairParams)
05821   {
05822     i1 = new int[NumTablePairParams];
05823     i2 = new int[NumTablePairParams];
05824     i3 = new int[NumTablePairParams];
05825 
05826     if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05827     {
05828       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05829     }
05830     
05831     msg->get(NumTablePairParams, i1);
05832     msg->get(NumTablePairParams, i2);
05833     msg->get(NumTablePairParams, i3);
05834     
05835     for (i=0; i<NumTablePairParams; i++)
05836     {
05837       new_tab_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
05838       
05839       if (new_tab_node == NULL)
05840       {
05841          NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05842       }
05843       
05844 //      printf("Adding new table pair with ind1 %i ind2 %i k %i\n", i1[i], i2[i],i3[i]);
05845       new_tab_node->ind1 = i1[i];
05846       new_tab_node->ind2 = i2[i];
05847       new_tab_node->K = i3[i];
05848       new_tab_node->left = NULL;
05849       new_tab_node->right = NULL;
05850       
05851       tab_pair_tree = add_to_indexed_table_pairs(new_tab_node, tab_pair_tree);
05852     }
05853     
05854     delete [] i1;
05855     delete [] i2;
05856     delete [] i3;
05857   }
05858   
05859   // receive the hydrogen bond parameters
05860   // hbondParams.receive_message(msg);
05861 
05862   AllFilesRead = TRUE;
05863 
05864   delete msg;
05865 }

void Parameters::send_Parameters ( MOStream  ) 

Definition at line 5044 of file Parameters.C.

References nbthole_pair_value::alphai, nbthole_pair_value::alphaj, angle_array, bond_array, crossterm_array, four_body_consts::delta, dihedral_array, CrosstermValue::dim, MOStream::end(), vdw_val::epsilon, vdw_val::epsilon14, gromacsPair_array, improper_array, nbthole_pair_value::ind2, AngleValue::k_ub, MAX_ATOMTYPE_CHARS, MAX_MULTIPLICITY, four_body_consts::n, NAMD_die(), nbthole_array, nbthole_pair_tree, AngleValue::normal, NumAngleParams, NumBondParams, NumCrosstermParams, NumDihedralParams, numenerentries, NumGromacsPairParams, NumImproperParams, NumNbtholePairParams, NumTablePairParams, NumVdwPairParams, NumVdwParams, NumVdwParamsAssigned, GromacsPairValue::pairC12, MOStream::put(), AngleValue::r_ub, vdw_val::sigma14, tab_pair_tree, table_ener, AngleValue::theta0, nbthole_pair_value::tholeij, ImproperValue::values, DihedralValue::values, vdw_array, vdw_pair_tree, and BondValue::x0.

Referenced by Node::resendMolecule().

05045 {
05046   Real *a1, *a2, *a3, *a4;        //  Temporary arrays for sending messages
05047   int *i1, *i2, *i3;      //  Temporary int array
05048   int i, j;      //  Loop counters
05049   Real **kvals;      //  Force constant values for dihedrals and impropers
05050   int **nvals;      //  Periodicity values for  dihedrals and impropers
05051   Real **deltavals;    //  Phase shift values for  dihedrals and impropers
05052   BigReal *pairC6, *pairC12; // JLai
05053   /*MOStream *msg=comm_obj->newOutputStream(ALLBUTME, STATICPARAMSTAG, BUFSIZE);
05054   if ( msg == NULL )
05055   {
05056     NAMD_die("memory allocation failed in Parameters::send_Parameters");
05057   }*/
05058 
05059   //  Send the bond parameters
05060   msg->put(NumBondParams);
05061 
05062   if (NumBondParams)
05063   {
05064     a1 = new Real[NumBondParams];
05065     a2 = new Real[NumBondParams];
05066 
05067     if ( (a1 == NULL) || (a2 == NULL) )
05068     {
05069       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05070     }
05071 
05072     for (i=0; i<NumBondParams; i++)
05073     {
05074       a1[i] = bond_array[i].k;
05075       a2[i] = bond_array[i].x0;
05076     }
05077 
05078     msg->put(NumBondParams, a1)->put(NumBondParams, a2);
05079 
05080     delete [] a1;
05081     delete [] a2;
05082   }
05083 
05084   //  Send the angle parameters
05085   msg->put(NumAngleParams);
05086 
05087   if (NumAngleParams)
05088   {
05089     a1 = new Real[NumAngleParams];
05090     a2 = new Real[NumAngleParams];
05091     a3 = new Real[NumAngleParams];
05092     a4 = new Real[NumAngleParams];
05093     i1 = new int[NumAngleParams];
05094 
05095     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) ||
05096          (a4 == NULL) || (i1 == NULL))
05097     {
05098       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05099     }
05100 
05101     for (i=0; i<NumAngleParams; i++)
05102     {
05103       a1[i] = angle_array[i].k;
05104       a2[i] = angle_array[i].theta0;
05105       a3[i] = angle_array[i].k_ub;
05106       a4[i] = angle_array[i].r_ub;
05107       i1[i] = angle_array[i].normal;
05108     }
05109 
05110     msg->put(NumAngleParams, a1)->put(NumAngleParams, a2);
05111     msg->put(NumAngleParams, a3)->put(NumAngleParams, a4);
05112     msg->put(NumAngleParams, i1);
05113 
05114     delete [] a1;
05115     delete [] a2;
05116     delete [] a3;
05117     delete [] a4;
05118     delete [] i1;
05119   }
05120 
05121   //  Send the dihedral parameters
05122   msg->put(NumDihedralParams);
05123 
05124   if (NumDihedralParams)
05125   {
05126     i1 = new int[NumDihedralParams];
05127     kvals = new Real *[MAX_MULTIPLICITY];
05128     nvals = new int *[MAX_MULTIPLICITY];
05129     deltavals = new Real *[MAX_MULTIPLICITY];
05130 
05131     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
05132          (deltavals == NULL) )
05133     {
05134       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05135     }
05136 
05137     for (i=0; i<MAX_MULTIPLICITY; i++)
05138     {
05139       kvals[i] = new Real[NumDihedralParams];
05140       nvals[i] = new int[NumDihedralParams];
05141       deltavals[i] = new Real[NumDihedralParams];
05142 
05143       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05144       {
05145         NAMD_die("memory allocation failed in Parameters::send_Parameters");
05146       }
05147     }
05148 
05149     for (i=0; i<NumDihedralParams; i++)
05150     {
05151       i1[i] = dihedral_array[i].multiplicity;
05152 
05153       for (j=0; j<MAX_MULTIPLICITY; j++)
05154       {
05155         kvals[j][i] = dihedral_array[i].values[j].k;
05156         nvals[j][i] = dihedral_array[i].values[j].n;
05157         deltavals[j][i] = dihedral_array[i].values[j].delta;
05158       }
05159     }
05160 
05161     msg->put(NumDihedralParams, i1);
05162 
05163     for (i=0; i<MAX_MULTIPLICITY; i++)
05164     {
05165       msg->put(NumDihedralParams, kvals[i]);
05166       msg->put(NumDihedralParams, nvals[i]);
05167       msg->put(NumDihedralParams, deltavals[i]);
05168 
05169       delete [] kvals[i];
05170       delete [] nvals[i];
05171       delete [] deltavals[i];
05172     }
05173 
05174     delete [] i1;
05175     delete [] kvals;
05176     delete [] nvals;
05177     delete [] deltavals;
05178   }
05179 
05180   //  Send the improper parameters
05181   msg->put(NumImproperParams);
05182 
05183   if (NumImproperParams)
05184   {
05185     i1 = new int[NumImproperParams];
05186     kvals = new Real *[MAX_MULTIPLICITY];
05187     nvals = new int *[MAX_MULTIPLICITY];
05188     deltavals = new Real *[MAX_MULTIPLICITY];
05189 
05190     if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
05191          (deltavals == NULL) )
05192     {
05193       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05194     }
05195 
05196     for (i=0; i<MAX_MULTIPLICITY; i++)
05197     {
05198       kvals[i] = new Real[NumImproperParams];
05199       nvals[i] = new int[NumImproperParams];
05200       deltavals[i] = new Real[NumImproperParams];
05201 
05202       if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05203       {
05204         NAMD_die("memory allocation failed in Parameters::send_Parameters");
05205       }
05206     }
05207 
05208     for (i=0; i<NumImproperParams; i++)
05209     {
05210       i1[i] = improper_array[i].multiplicity;
05211 
05212       for (j=0; j<MAX_MULTIPLICITY; j++)
05213       {
05214         kvals[j][i] = improper_array[i].values[j].k;
05215         nvals[j][i] = improper_array[i].values[j].n;
05216         deltavals[j][i] = improper_array[i].values[j].delta;
05217       }
05218     }
05219 
05220     msg->put(NumImproperParams, i1);
05221 
05222     for (i=0; i<MAX_MULTIPLICITY; i++)
05223     {
05224       msg->put(NumImproperParams, kvals[i]);
05225       msg->put(NumImproperParams, nvals[i]);
05226       msg->put(NumImproperParams, deltavals[i]);
05227 
05228       delete [] kvals[i];
05229       delete [] nvals[i];
05230       delete [] deltavals[i];
05231     }
05232 
05233     delete [] i1;
05234     delete [] kvals;
05235     delete [] nvals;
05236     delete [] deltavals;
05237   }
05238 
05239   //  Send the crossterm parameters
05240   msg->put(NumCrosstermParams);
05241 
05242   if (NumCrosstermParams)
05243   {
05244     for (i=0; i<NumCrosstermParams; ++i) {
05245       int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
05246       msg->put(nvals,&crossterm_array[i].c[0][0].d00);
05247     }
05248   }
05249   //  Send the GromacsPairs parameters
05250   // JLai
05251   msg->put(NumGromacsPairParams);
05252  
05253   if(NumGromacsPairParams) 
05254   {
05255       pairC6 = new BigReal[NumGromacsPairParams];
05256       pairC12 = new BigReal[NumGromacsPairParams];
05257       if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
05258           NAMD_die("memory allocation failed in Parameters::send_Parameters");
05259       }
05260 
05261       for (i=0; i<NumGromacsPairParams; i++) {
05262           pairC6[i] = gromacsPair_array[i].pairC6;
05263           pairC12[i] = gromacsPair_array[i].pairC12;
05264       }
05265 
05266       msg->put(NumGromacsPairParams,pairC6);
05267       msg->put(NumGromacsPairParams,pairC12);
05268 
05269       delete [] pairC6;
05270       delete [] pairC12;
05271   }
05272   // End of JLai
05273   
05274   //
05275   //Send the energy table parameters
05276   msg->put(numenerentries);
05277 
05278   if (numenerentries) {
05279           /*
05280     b1 = new Real[numenerentries];
05281     if (b1 == NULL) 
05282     {
05283       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05284     }
05285     */
05286     
05287     msg->put(numenerentries, table_ener);
05288   }
05289 
05290   //  Send the vdw parameters
05291   msg->put(NumVdwParams);
05292   msg->put(NumVdwParamsAssigned);
05293 
05294   if (NumVdwParams)
05295   {
05296     a1 = new Real[NumVdwParams];
05297     a2 = new Real[NumVdwParams];
05298     a3 = new Real[NumVdwParams];
05299     a4 = new Real[NumVdwParams];
05300 
05301     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) )
05302     {
05303       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05304     }
05305 
05306     for (i=0; i<NumVdwParams; i++)
05307     {
05308       a1[i] = vdw_array[i].sigma;
05309       a2[i] = vdw_array[i].epsilon;
05310       a3[i] = vdw_array[i].sigma14;
05311       a4[i] = vdw_array[i].epsilon14;
05312     }
05313 
05314     msg->put(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
05315     msg->put(NumVdwParams, a1);
05316     msg->put(NumVdwParams, a2);
05317     msg->put(NumVdwParams, a3);
05318     msg->put(NumVdwParams, a4);
05319 
05320     delete [] a1;
05321     delete [] a2;
05322     delete [] a3;
05323     delete [] a4;
05324   }
05325   
05326   //  Send the vdw pair parameters
05327   msg->put(NumVdwPairParams);
05328   
05329   if (NumVdwPairParams)
05330   {
05331     a1 = new Real[NumVdwPairParams];
05332     a2 = new Real[NumVdwPairParams];
05333     a3 = new Real[NumVdwPairParams];
05334     a4 = new Real[NumVdwPairParams];
05335     i1 = new int[NumVdwPairParams];
05336     i2 = new int[NumVdwPairParams];    
05337 
05338     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
05339          (i1 == NULL) || (i2 == NULL) )
05340     {
05341       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05342     }
05343     
05344     vdw_pair_to_arrays(i1, i2, a1, a2, a3, a4, 0, vdw_pair_tree);
05345     
05346     msg->put(NumVdwPairParams, i1)->put(NumVdwPairParams, i2);
05347     msg->put(NumVdwPairParams, a1);
05348     msg->put(NumVdwPairParams, a2)->put(NumVdwPairParams, a3);
05349     msg->put(NumVdwPairParams, a4);
05350   }
05351 
05352   //  Send the nbthole pair parameters
05353   msg->put(NumNbtholePairParams);
05354 
05355   if (NumNbtholePairParams)
05356   {
05357     a1 = new Real[NumNbtholePairParams];
05358     a2 = new Real[NumNbtholePairParams];
05359     a3 = new Real[NumNbtholePairParams];
05360     i1 = new int[NumNbtholePairParams];
05361     i2 = new int[NumNbtholePairParams];
05362 
05363     if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05364     {
05365       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05366     }
05367 
05368     nbthole_pair_to_arrays(i1, i2, a1, a2, a3, 0, nbthole_pair_tree);
05369 
05370    for (i=0; i<NumNbtholePairParams; i++)
05371    {
05372     nbthole_array[i].ind1 = i1[i];
05373     nbthole_array[i].ind2 = i2[i];
05374     nbthole_array[i].alphai = a1[i];
05375     nbthole_array[i].alphaj = a2[i];
05376     nbthole_array[i].tholeij = a3[i];
05377    }
05378 
05379     msg->put(NumNbtholePairParams, i1)->put(NumNbtholePairParams, i2);
05380     msg->put(NumNbtholePairParams, a1);
05381     msg->put(NumNbtholePairParams, a2)->put(NumNbtholePairParams, a3);
05382   }
05383   
05384   //  Send the table pair parameters
05385   //printf("Pairs: %i\n", NumTablePairParams);
05386   msg->put(NumTablePairParams);
05387   
05388   if (NumTablePairParams)
05389   {
05390     i1 = new int[NumTablePairParams];
05391     i2 = new int[NumTablePairParams];    
05392     i3 = new int[NumTablePairParams];
05393 
05394     if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05395     {
05396       NAMD_die("memory allocation failed in Parameters::send_Parameters");
05397     }
05398     
05399     table_pair_to_arrays(i1, i2, i3, 0, tab_pair_tree);
05400     
05401     msg->put(NumTablePairParams, i1)->put(NumTablePairParams, i2);
05402     msg->put(NumTablePairParams, i3);
05403   }
05404 
05405   // send the hydrogen bond parameters
05406   // hbondParams.create_message(msg);
05407   msg->end();
05408   delete msg;
05409 }


Member Data Documentation

AngleValue* Parameters::angle_array

Definition at line 240 of file Parameters.h.

Referenced by done_reading_files(), get_angle_params(), read_parm(), receive_Parameters(), send_Parameters(), and ~Parameters().

BondValue* Parameters::bond_array

Definition at line 239 of file Parameters.h.

Referenced by done_reading_files(), get_bond_params(), read_parm(), receive_Parameters(), send_Parameters(), and ~Parameters().

int Parameters::columnsize

Definition at line 251 of file Parameters.h.

Referenced by ComputeNonbondedUtil::select().

CrosstermValue* Parameters::crossterm_array

Definition at line 243 of file Parameters.h.

Referenced by done_reading_files(), receive_Parameters(), send_Parameters(), and ~Parameters().

DihedralValue* Parameters::dihedral_array

Definition at line 241 of file Parameters.h.

Referenced by assign_dihedral_index(), done_reading_files(), get_dihedral_multiplicity(), get_dihedral_params(), outputCompressedFile(), read_parm(), receive_Parameters(), send_Parameters(), and ~Parameters().

GromacsPairValue* Parameters::gromacsPair_array

Definition at line 245 of file Parameters.h.

Referenced by done_reading_files(), receive_Parameters(), send_Parameters(), and ~Parameters().

ImproperValue* Parameters::improper_array

Definition at line 242 of file Parameters.h.

Referenced by assign_improper_index(), done_reading_files(), get_improper_multiplicity(), get_improper_params(), outputCompressedFile(), receive_Parameters(), send_Parameters(), and ~Parameters().

NbtholePairValue* Parameters::nbthole_array

Definition at line 248 of file Parameters.h.

Referenced by done_reading_files(), receive_Parameters(), SELF(), and send_Parameters().

IndexedNbtholePair* Parameters::nbthole_pair_tree

Definition at line 254 of file Parameters.h.

Referenced by send_Parameters(), and ~Parameters().

int Parameters::NumAngleParams

Definition at line 258 of file Parameters.h.

Referenced by done_reading_files(), print_angle_params(), print_param_summary(), read_charmm_parameter_file(), read_parameter_file(), read_parm(), receive_Parameters(), and send_Parameters().

int Parameters::NumBondParams

Definition at line 257 of file Parameters.h.

Referenced by done_reading_files(), print_bond_params(), print_param_summary(), read_charmm_parameter_file(), read_parameter_file(), read_parm(), receive_Parameters(), and send_Parameters().

int Parameters::NumCosAngles

Definition at line 259 of file Parameters.h.

Referenced by print_param_summary().

int Parameters::NumCrosstermParams

Definition at line 262 of file Parameters.h.

Referenced by done_reading_files(), print_param_summary(), read_charmm_parameter_file(), receive_Parameters(), and send_Parameters().

int Parameters::NumDihedralParams

Definition at line 260 of file Parameters.h.

Referenced by done_reading_files(), outputCompressedFile(), print_dihedral_params(), print_param_summary(), read_charmm_parameter_file(), read_parameter_file(), read_parm(), receive_Parameters(), and send_Parameters().

int Parameters::numenerentries

Definition at line 249 of file Parameters.h.

Referenced by read_ener_table(), receive_Parameters(), and send_Parameters().

int Parameters::NumGromacsPairParams

Definition at line 264 of file Parameters.h.

Referenced by done_reading_files(), receive_Parameters(), and send_Parameters().

int Parameters::NumImproperParams

Definition at line 261 of file Parameters.h.

Referenced by done_reading_files(), outputCompressedFile(), print_improper_params(), print_param_summary(), read_charmm_parameter_file(), read_parameter_file(), receive_Parameters(), and send_Parameters().

int Parameters::NumNbtholePairParams

Definition at line 270 of file Parameters.h.

Referenced by done_reading_files(), print_nbthole_pair_params(), print_param_summary(), read_charmm_parameter_file(), receive_Parameters(), SELF(), and send_Parameters().

int Parameters::NumTablePairParams

Definition at line 271 of file Parameters.h.

Referenced by read_charmm_parameter_file(), read_parameter_file(), receive_Parameters(), and send_Parameters().

int Parameters::NumTableParams

Definition at line 267 of file Parameters.h.

int Parameters::NumVdwPairParams

Definition at line 269 of file Parameters.h.

Referenced by print_param_summary(), print_vdw_pair_params(), read_charmm_parameter_file(), read_parameter_file(), read_parm(), receive_Parameters(), and send_Parameters().

int Parameters::NumVdwParams

Definition at line 266 of file Parameters.h.

Referenced by done_reading_files(), print_param_summary(), print_vdw_params(), read_charmm_parameter_file(), read_parameter_file(), receive_Parameters(), and send_Parameters().

int Parameters::NumVdwParamsAssigned

Definition at line 268 of file Parameters.h.

Referenced by done_reading_files(), get_num_vdw_params(), read_parm(), receive_Parameters(), and send_Parameters().

int Parameters::rowsize

Definition at line 250 of file Parameters.h.

Referenced by ComputeNonbondedUtil::select().

IndexedTablePair* Parameters::tab_pair_tree

Definition at line 255 of file Parameters.h.

Referenced by get_table_pair_params(), receive_Parameters(), send_Parameters(), and ~Parameters().

BigReal* Parameters::table_ener

Definition at line 252 of file Parameters.h.

Referenced by read_charmm_parameter_file(), read_ener_table(), read_parameter_file(), receive_Parameters(), ComputeNonbondedUtil::select(), and send_Parameters().

int Parameters::tablenumtypes

Definition at line 256 of file Parameters.h.

Referenced by get_int_table_type(), and read_ener_table().

VdwValue* Parameters::vdw_array

Definition at line 247 of file Parameters.h.

Referenced by done_reading_files(), get_vdw_params(), receive_Parameters(), send_Parameters(), and ~Parameters().

IndexedVdwPair* Parameters::vdw_pair_tree

Definition at line 253 of file Parameters.h.

Referenced by get_vdw_pair_params(), read_parm(), receive_Parameters(), send_Parameters(), and ~Parameters().


The documentation for this class was generated from the following files:
Generated on Wed Nov 21 01:17:20 2018 for NAMD by  doxygen 1.4.7