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 (Bool)
void done_reading_structure ()
void assign_vdw_index (const char *, Atom *)
void assign_bond_index (const char *, const char *, Bond *)
void assign_angle_index (const char *, const char *, const char *, Angle *, int)
void assign_dihedral_index (const char *, const char *, const char *, const char *, Dihedral *, int, int)
void assign_improper_index (const char *, const char *, const char *, const char *, Improper *, int)
void assign_crossterm_index (const char *, const char *, const char *, const char *, const char *, const char *, const char *, const 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 simParams,
StringList f 
)

Definition at line 249 of file Parameters.C.

References SimParameters::cosAngles, StringList::data, done_reading_files(), SimParameters::drudeOn, FALSE, StringList::next, paraCharmm, SimParameters::paraTypeCharmmOn, SimParameters::paraTypeXplorOn, paraXplor, read_charmm_parameter_file(), read_ener_table(), read_parameter_file(), and SimParameters::tabulatedEnergies.

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(simParams->drudeOn && paramType == paraCharmm);
00299   }
00300 
00301 }

Parameters::Parameters ( Ambertoppar amber_data,
BigReal  vdw14 
)

Definition at line 6344 of file Parameters.C.

06345 {
06346   initialize();
06347 
06348   // Read in parm parameters
06349   read_parm(amber_data,vdw14);
06350 }

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

Definition at line 6486 of file Parameters.C.

06487 {
06488   initialize();
06489 
06490   // Read in parm parameters
06491   read_parm(gf,min);
06492 }

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 ( const char *  atom1,
const char *  atom2,
const char *  atom3,
Angle angle_ptr,
int  notFoundIndex 
)

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

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

void Parameters::assign_bond_index ( const char *  atom1,
const char *  atom2,
Bond bond_ptr 
)

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

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

void Parameters::assign_crossterm_index ( const char *  atom1,
const char *  atom2,
const char *  atom3,
const char *  atom4,
const char *  atom5,
const char *  atom6,
const char *  atom7,
const char *  atom8,
Crossterm crossterm_ptr 
)

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

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

void Parameters::assign_dihedral_index ( const char *  atom1,
const char *  atom2,
const char *  atom3,
const char *  atom4,
Dihedral dihedral_ptr,
int  multiplicity,
int  notFoundIndex 
)

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

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

void Parameters::assign_improper_index ( const char *  atom1,
const char *  atom2,
const char *  atom3,
const char *  atom4,
Improper improper_ptr,
int  multiplicity 
)

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

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

void Parameters::assign_vdw_index ( const char *  atomtype,
Atom atom_ptr 
)

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

Referenced by Molecule::prepare_qm().

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

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 ( Bool  addDrudeBond  ) 

Definition at line 2959 of file Parameters.C.

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

Referenced by Parameters().

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

void Parameters::done_reading_structure (  ) 

Definition at line 4999 of file Parameters.C.

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

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 *  tabletype  ) 

Definition at line 7564 of file Parameters.C.

References table_types, and tablenumtypes.

07564                                                   {
07565   for (int i=0; i<tablenumtypes; i++) {
07566     if (!strncmp(tabletype, table_types[i], 5)) {
07567       return i;
07568     }
07569   }
07570 
07571   return -1;
07572 }

int Parameters::get_num_vdw_params ( void   )  [inline]

Definition at line 529 of file Parameters.h.

References NumVdwParamsAssigned.

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

00529 { return NumVdwParamsAssigned; }

int Parameters::get_table_pair_params ( Index  ind1,
Index  ind2,
int *  K 
)

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

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

int Parameters::get_vdw_pair_params ( Index  ind1,
Index  ind2,
Real A,
Real B,
Real A14,
Real B14 
)

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

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

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

References DebugM, and NumAngleParams.

04866 {
04867   DebugM(3,NumAngleParams << " ANGLE PARAMETERS\n"
04868       << "*****************************************" );
04869   traverse_angle_params(anglep);
04870 }

void Parameters::print_bond_params (  ) 

Definition at line 4847 of file Parameters.C.

References DebugM, and NumBondParams.

04848 {
04849   DebugM(3,NumBondParams << " BOND PARAMETERS\n" \
04850       << "*****************************************"  \
04851       );
04852 
04853   traverse_bond_params(bondp);
04854 }

void Parameters::print_crossterm_params (  ) 
void Parameters::print_dihedral_params (  ) 

Definition at line 4881 of file Parameters.C.

References DebugM, and NumDihedralParams.

04882 {
04883   DebugM(3,NumDihedralParams << " DIHEDRAL PARAMETERS\n" \
04884       << "*****************************************" );
04885 
04886   traverse_dihedral_params(dihedralp);
04887 }

void Parameters::print_improper_params (  ) 

Definition at line 4898 of file Parameters.C.

References DebugM, and NumImproperParams.

04899 {
04900   DebugM(3,NumImproperParams << " IMPROPER PARAMETERS\n" \
04901       << "*****************************************" );
04902 
04903   traverse_improper_params(improperp);
04904 }

void Parameters::print_nbthole_pair_params (  ) 

Definition at line 4949 of file Parameters.C.

References DebugM, and NumNbtholePairParams.

04950 {
04951   DebugM(3,NumNbtholePairParams << " NBTHOLE PAIR PARAMETERS\n" \
04952       << "*****************************************" );
04953 
04954   traverse_nbthole_pair_params(nbthole_pairp);
04955 } 

void Parameters::print_param_summary (  ) 

Definition at line 4966 of file Parameters.C.

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

Referenced by NamdState::loadStructure().

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

void Parameters::print_vdw_pair_params (  ) 

Definition at line 4932 of file Parameters.C.

References DebugM, and NumVdwPairParams.

04933 {
04934   DebugM(3,NumVdwPairParams << " vdW PAIR PARAMETERS\n" \
04935       << "*****************************************" );
04936 
04937   traverse_vdw_pair_params(vdw_pairp);
04938 }

void Parameters::print_vdw_params (  ) 

Definition at line 4915 of file Parameters.C.

References DebugM, and NumVdwParams.

04916 {
04917   DebugM(3,NumVdwParams << " vdW PARAMETERS\n" \
04918       << "*****************************************" );
04919 
04920   traverse_vdw_params(vdwp);
04921 }

void Parameters::read_charmm_parameter_file ( char *  fname  ) 

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, table_ener, and TRUE.

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, TRUE);
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 simParams  ) 

Definition at line 6728 of file Parameters.C.

References SimParameters::cutoff, iout, NAMD_die(), namdnearbyint, numenerentries, read_energy_type(), read_energy_type_bothcubspline(), table_ener, table_spacing, table_types, SimParameters::tableInterpType, SimParameters::tableMaxDist, tablenumtypes, SimParameters::tableNumTypes, SimParameters::tableSpacing, and SimParameters::tabulatedEnergiesFile.

Referenced by Parameters().

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

int Parameters::read_energy_type ( FILE *  enertable,
const int  typeindex,
BigReal table_ener,
const float  table_spacing,
const float  maxdist 
)

Definition at line 7446 of file Parameters.C.

References NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

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

int Parameters::read_energy_type_bothcubspline ( FILE *  enertable,
const int  typeindex,
BigReal table_ener,
const float  table_spacing,
const float  maxdist 
)

Definition at line 6971 of file Parameters.C.

References INDEX, NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

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

int Parameters::read_energy_type_cubspline ( FILE *  enertable,
const int  typeindex,
BigReal table_ener,
const float  table_spacing,
const float  maxdist 
)

Definition at line 7207 of file Parameters.C.

References INDEX, NAMD_die(), and namdnearbyint.

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

void Parameters::read_parameter_file ( char *  fname  ) 

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, table_ener, and TRUE.

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, TRUE);
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 amber_data,
BigReal  vdw14 
)

Definition at line 6366 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, indexed_vdw_pair::ind1, indexed_vdw_pair::ind2, iout, iWARN(), four_body_consts::k, AngleValue::k, BondValue::k, 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, values, DihedralValue::values, vdw_pair_tree, and BondValue::x0.

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

void Parameters::receive_Parameters ( MIStream msg  ) 

Definition at line 5424 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, nbthole_pair_value::ind1, indexed_vdw_pair::ind1, indexed_table_pair::ind2, nbthole_pair_value::ind2, indexed_vdw_pair::ind2, indexed_table_pair::K, four_body_consts::k, AngleValue::k, BondValue::k, AngleValue::k_ub, indexed_table_pair::left, indexed_vdw_pair::left, MAX_ATOMTYPE_CHARS, MAX_MULTIPLICITY, ImproperValue::multiplicity, DihedralValue::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, GromacsPairValue::pairC6, 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().

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

void Parameters::send_Parameters ( MOStream msg  ) 

Definition at line 5048 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::ind1, nbthole_pair_value::ind2, four_body_consts::k, AngleValue::k, BondValue::k, AngleValue::k_ub, MAX_ATOMTYPE_CHARS, MAX_MULTIPLICITY, ImproperValue::multiplicity, DihedralValue::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, GromacsPairValue::pairC6, MOStream::put(), AngleValue::r_ub, vdw_val::sigma, 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().

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


Member Data Documentation

Definition at line 251 of file Parameters.h.

Referenced by ComputeNonbondedUtil::select().

Definition at line 248 of file Parameters.h.

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

Definition at line 254 of file Parameters.h.

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

Definition at line 259 of file Parameters.h.

Referenced by print_param_summary().

Definition at line 249 of file Parameters.h.

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

Definition at line 264 of file Parameters.h.

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

Definition at line 267 of file Parameters.h.

Definition at line 250 of file Parameters.h.

Referenced by ComputeNonbondedUtil::select().

Definition at line 256 of file Parameters.h.

Referenced by get_int_table_type(), and read_ener_table().


The documentation for this class was generated from the following files:

Generated on 6 Dec 2019 for NAMD by  doxygen 1.6.1