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 6327 of file Parameters.C.

06328 {
06329   initialize();
06330 
06331   // Read in parm parameters
06332   read_parm(amber_data,vdw14);
06333 }

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

Definition at line 6469 of file Parameters.C.

06470 {
06471   initialize();
06472 
06473   // Read in parm parameters
06474   read_parm(gf,min);
06475 }

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 3836 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.

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

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

Definition at line 3739 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.

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

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

Definition at line 4147 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.

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

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

Definition at line 3939 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.

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

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

Definition at line 4048 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.

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

void Parameters::assign_vdw_index ( char *  ,
Atom  
)

Definition at line 3451 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.

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

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   }
02973 
02974   if (NumAngleParams)
02975   {
02976     angle_array = new AngleValue[NumAngleParams];
02977 
02978     if (angle_array == NULL)
02979     {
02980       NAMD_die("memory allocation of angle_array failed!");
02981     }
02982   }
02983 
02984   if (NumDihedralParams)
02985   {
02986     dihedral_array = new DihedralValue[NumDihedralParams];
02987 
02988     if (dihedral_array == NULL)
02989     {
02990       NAMD_die("memory allocation of dihedral_array failed!");
02991     }
02992     memset(dihedral_array, 0, NumDihedralParams*sizeof(DihedralValue));
02993   }
02994 
02995   if (NumImproperParams)
02996   {
02997     improper_array = new ImproperValue[NumImproperParams];
02998 
02999     if (improper_array == NULL)
03000     {
03001       NAMD_die("memory allocation of improper_array failed!");
03002     }
03003     memset(improper_array, 0, NumImproperParams*sizeof(ImproperValue));
03004   }
03005 
03006   if (NumCrosstermParams)
03007   {
03008     crossterm_array = new CrosstermValue[NumCrosstermParams];
03009     memset(crossterm_array, 0, NumCrosstermParams*sizeof(CrosstermValue));
03010   }
03011 
03012   // JLai
03013   if (NumGromacsPairParams)
03014   {
03015     gromacsPair_array = new GromacsPairValue[NumGromacsPairParams];
03016     memset(gromacsPair_array, 0, NumGromacsPairParams*sizeof(GromacsPairValue)); 
03017   }
03018   // End of JLai
03019 
03020   if (NumVdwParams)
03021   {
03022           atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
03023     vdw_array = new VdwValue[NumVdwParams];
03024     
03025     if (vdw_array == NULL)
03026     {
03027       NAMD_die("memory allocation of vdw_array failed!");
03028     }
03029   }
03030   if (NumNbtholePairParams)
03031   {
03032     nbthole_array = new NbtholePairValue[NumNbtholePairParams];
03033 
03034     if(nbthole_array == NULL)
03035     {
03036       NAMD_die("memory allocation of nbthole_array failed!");
03037     }
03038   }
03039   //  Assign indexes to each of the parameters and populate the
03040   //  arrays using the binary trees and linked lists that we have
03041   //  already read in
03042   index_bonds(bondp, 0);
03043   index_angles(anglep, 0);
03044   NumVdwParamsAssigned = index_vdw(vdwp, 0);
03045   index_dihedrals();
03046   index_impropers();
03047   index_crossterms();
03048   convert_nbthole_pairs();
03049   //  Convert the vdw pairs
03050   convert_vdw_pairs();
03051   convert_table_pairs();
03052 }

void Parameters::done_reading_structure (  ) 

Definition at line 4982 of file Parameters.C.

04984 {
04985   if (bondp != NULL)
04986     free_bond_tree(bondp);
04987 
04988   if (anglep != NULL)
04989     free_angle_tree(anglep);
04990 
04991   if (dihedralp != NULL)
04992     free_dihedral_list(dihedralp);
04993 
04994   if (improperp != NULL)
04995     free_improper_list(improperp);
04996 
04997   if (crosstermp != NULL)
04998     free_crossterm_list(crosstermp);
04999 
05000   if (vdwp != NULL)
05001     free_vdw_tree(vdwp);
05002 
05003   //  Free the arrays used to track multiplicity for dihedrals
05004   //  and impropers
05005   if (maxDihedralMults != NULL)
05006     delete [] maxDihedralMults;
05007 
05008   if (maxImproperMults != NULL)
05009     delete [] maxImproperMults;
05010 
05011   bondp=NULL;
05012   anglep=NULL;
05013   dihedralp=NULL;
05014   improperp=NULL;
05015   crosstermp=NULL;
05016   vdwp=NULL;
05017   maxImproperMults=NULL;
05018   maxDihedralMults=NULL;
05019 }

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 7547 of file Parameters.C.

References table_types, and tablenumtypes.

07547                                                   {
07548   for (int i=0; i<tablenumtypes; i++) {
07549     if (!strncmp(tabletype, table_types[i], 5)) {
07550       return i;
07551     }
07552   }
07553 
07554   return -1;
07555 }

int Parameters::get_num_vdw_params ( void   )  [inline]

Definition at line 529 of file Parameters.h.

References NumVdwParamsAssigned.

Referenced by LJTable::LJTable().

00529 { return NumVdwParamsAssigned; }

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

Definition at line 3586 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.

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

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

Definition at line 3661 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().

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

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 4848 of file Parameters.C.

References DebugM, and NumAngleParams.

04849 {
04850   DebugM(3,NumAngleParams << " ANGLE PARAMETERS\n"
04851       << "*****************************************" );
04852   traverse_angle_params(anglep);
04853 }

void Parameters::print_bond_params (  ) 

Definition at line 4830 of file Parameters.C.

References DebugM, and NumBondParams.

04831 {
04832   DebugM(3,NumBondParams << " BOND PARAMETERS\n" \
04833       << "*****************************************"  \
04834       );
04835 
04836   traverse_bond_params(bondp);
04837 }

void Parameters::print_crossterm_params (  ) 

void Parameters::print_dihedral_params (  ) 

Definition at line 4864 of file Parameters.C.

References DebugM, and NumDihedralParams.

04865 {
04866   DebugM(3,NumDihedralParams << " DIHEDRAL PARAMETERS\n" \
04867       << "*****************************************" );
04868 
04869   traverse_dihedral_params(dihedralp);
04870 }

void Parameters::print_improper_params (  ) 

Definition at line 4881 of file Parameters.C.

References DebugM, and NumImproperParams.

04882 {
04883   DebugM(3,NumImproperParams << " IMPROPER PARAMETERS\n" \
04884       << "*****************************************" );
04885 
04886   traverse_improper_params(improperp);
04887 }

void Parameters::print_nbthole_pair_params (  ) 

Definition at line 4932 of file Parameters.C.

References DebugM, and NumNbtholePairParams.

04933 {
04934   DebugM(3,NumNbtholePairParams << " NBTHOLE PAIR PARAMETERS\n" \
04935       << "*****************************************" );
04936 
04937   traverse_nbthole_pair_params(nbthole_pairp);
04938 } 

void Parameters::print_param_summary (  ) 

Definition at line 4949 of file Parameters.C.

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

Referenced by NamdState::loadStructure().

04950 {
04951   iout << iINFO << "SUMMARY OF PARAMETERS:\n" 
04952        << iINFO << NumBondParams << " BONDS\n" 
04953        << iINFO << NumAngleParams << " ANGLES\n" << endi;
04954        if (cosAngles) {
04955          iout << iINFO << "  " << (NumAngleParams - NumCosAngles) << " HARMONIC\n"
04956               << iINFO << "  " << NumCosAngles << " COSINE-BASED\n" << endi;
04957        }
04958   iout << iINFO << NumDihedralParams << " DIHEDRAL\n"
04959        << iINFO << NumImproperParams << " IMPROPER\n"
04960        << iINFO << NumCrosstermParams << " CROSSTERM\n"
04961        << iINFO << NumVdwParams << " VDW\n"
04962        << iINFO << NumVdwPairParams << " VDW_PAIRS\n"
04963        << iINFO << NumNbtholePairParams << " NBTHOLE_PAIRS\n" << endi;
04964 }

void Parameters::print_vdw_pair_params (  ) 

Definition at line 4915 of file Parameters.C.

References DebugM, and NumVdwPairParams.

04916 {
04917   DebugM(3,NumVdwPairParams << " vdW PAIR PARAMETERS\n" \
04918       << "*****************************************" );
04919 
04920   traverse_vdw_pair_params(vdw_pairp);
04921 }

void Parameters::print_vdw_params (  ) 

Definition at line 4898 of file Parameters.C.

References DebugM, and NumVdwParams.

04899 {
04900   DebugM(3,NumVdwParams << " vdW PARAMETERS\n" \
04901       << "*****************************************" );
04902 
04903   traverse_vdw_params(vdwp);
04904 }

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 6711 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().

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

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

Definition at line 7429 of file Parameters.C.

References NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

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

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

Definition at line 6954 of file Parameters.C.

References INDEX, NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

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

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

Definition at line 7190 of file Parameters.C.

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

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

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

void Parameters::receive_Parameters ( MIStream  ) 

Definition at line 5407 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().

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

void Parameters::send_Parameters ( MOStream  ) 

Definition at line 5031 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().

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


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 22 01:17:22 2017 for NAMD by  doxygen 1.4.7