Parameters Class Reference

#include <Parameters.h>

Public Member Functions
	Parameters ()
	Parameters (SimParameters *, StringList *f)
	Parameters (Ambertoppar *, BigReal)
void	read_parm (Ambertoppar *, BigReal)
	Parameters (AmberParm7Reader::Ambertoppar *, BigReal)
void	read_parm (AmberParm7Reader::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
BondValue *	bond_array
AngleValue *	angle_array
DihedralValue *	dihedral_array
ImproperValue *	improper_array
CrosstermValue *	crossterm_array
GromacsPairValue *	gromacsPair_array
VdwValue *	vdw_array
NbtholePairValue *	nbthole_array
int	numenerentries
int	rowsize
int	columnsize
BigReal *	table_ener
IndexedVdwPair *	vdw_pair_tree
IndexedNbtholePair *	nbthole_pair_tree
IndexedTablePair *	tab_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 218 of file Parameters.h.

Constructor & Destructor Documentation

Parameters::Parameters

(

)

Definition at line 186 of file Parameters.C.

                       {
  initialize();
}

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.

{
  initialize();

  if (simParams->paraTypeXplorOn)
  {
    paramType = paraXplor;
  }
  else if (simParams->paraTypeCharmmOn)
  {
    paramType = paraCharmm;
  }
  //****** END CHARMM/XPLOR type changes
  //Test for cos-based angles
  if (simParams->cosAngles) {
    cosAngles = true;
  } else {
    cosAngles = false;
  }

  if (simParams->tabulatedEnergies) {
          CkPrintf("Working on tables\n");
          read_ener_table(simParams);
  }

  //****** BEGIN CHARMM/XPLOR type changes
  /* Set up AllFilesRead flag to FALSE.  Once all of the files    */
  /* have been read in, then this will be set to true and the     */
  /* arrays of parameters will be set up        */
  AllFilesRead = FALSE;

  if (NULL != f) 
  {
    do
    {
      //****** BEGIN CHARMM/XPLOR type changes
      if (paramType == paraXplor)
      {
        read_parameter_file(f->data);
      }
      else if (paramType == paraCharmm)
      {
        read_charmm_parameter_file(f->data);
      }
      //****** END CHARMM/XPLOR type changes
      f = f->next;
    } while ( f != NULL );

    done_reading_files(simParams->drudeOn && paramType == paraCharmm);
  }

}

Parameters::Parameters	(	Ambertoppar *	amber_data,
		BigReal	vdw14
	)

Definition at line 6344 of file Parameters.C.

{
  initialize();

  // Read in parm parameters
  read_parm(amber_data,vdw14);
}

Parameters::Parameters	(	AmberParm7Reader::Ambertoppar *	amber_data,
		BigReal	vdw14
	)

Definition at line 6479 of file Parameters.C.

{
  initialize();
  
  // Read in parm parameters
  read_parm(amber_data,vdw14);
}

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

Definition at line 6646 of file Parameters.C.

{
  initialize();

  // Read in parm parameters
  read_parm(gf,min);
}

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< T >::resize(), ResizeArray< T >::size(), tab_pair_tree, vdw_array, and vdw_pair_tree.

{
        if (atomTypeNames)
          delete [] atomTypeNames;

  if (bondp != NULL)
    free_bond_tree(bondp);

  if (anglep != NULL)
    free_angle_tree(anglep);

  if (dihedralp != NULL)
    free_dihedral_list(dihedralp);

  if (improperp != NULL)
    free_improper_list(improperp);

  if (crosstermp != NULL)
    free_crossterm_list(crosstermp);

  if (vdwp != NULL)
    free_vdw_tree(vdwp);

  if (vdw_pairp != NULL)
    free_vdw_pair_list();

  if (nbthole_pairp != NULL)
    free_nbthole_pair_list();

  if (bond_array != NULL)
    delete [] bond_array;

  if (angle_array != NULL)
    delete [] angle_array;

  if (dihedral_array != NULL)
    delete [] dihedral_array;

  if (improper_array != NULL)
    delete [] improper_array;

  if (crossterm_array != NULL)
    delete [] crossterm_array;

  // JLai
  if (gromacsPair_array != NULL)
    delete [] gromacsPair_array;
  // End of JLai

  if (vdw_array != NULL)
    delete [] vdw_array;
  
  if (tab_pair_tree != NULL)
    free_table_pair_tree(tab_pair_tree);

  if (vdw_pair_tree != NULL)
    free_vdw_pair_tree(vdw_pair_tree);

  if (nbthole_pair_tree != NULL)
    free_nbthole_pair_tree(nbthole_pair_tree);

  if (maxDihedralMults != NULL)
    delete [] maxDihedralMults;

  if (maxImproperMults != NULL)
    delete [] maxImproperMults;

  for( int i = 0; i < error_msgs.size(); ++i ) {
    delete [] error_msgs[i];
  }
  error_msgs.resize(0);
}

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, endi(), angle_params::index, iout, iWARN(), angle_params::left, NAMD_die(), and angle_params::right.

{
  struct angle_params *ptr;  //  Current position in tree
  int comp_val;      //  value from strcasecmp
  int found=0;      //  flag 1->found a match

  /*  Check to make sure the files have all been read    */
  if (!AllFilesRead)
  {
    NAMD_die("Tried to assign angle index before all parameter files were read");
  }

  /*  We need atom1 < atom3.  If that was not what we were   */
  /*  passed, switch them            */
  if (strcasecmp(atom1, atom3) > 0)
  {
    const char *tmp_name = atom1;
    atom1 = atom3;
    atom3 = tmp_name;
  }

  /*  Start at the top            */
  ptr=anglep;

  /*  While we don't have a match and we haven't reached the  */
  /*  bottom of the tree, compare values        */
  while (!found && (ptr != NULL))
  {
    comp_val = strcasecmp(atom1, ptr->atom1name);

    if (comp_val == 0)
    {
      /*  Atom 1 matches, so compare atom 2    */
      comp_val = strcasecmp(atom2, ptr->atom2name);
      
      if (comp_val == 0)
      {
        /*  Atoms 1&2 match, try atom 3    */
        comp_val = strcasecmp(atom3, ptr->atom3name);
      }
    }

    if (comp_val == 0)
    {
      /*  Found a match        */
      found = 1;
      angle_ptr->angle_type = ptr->index;
    }
    else if (comp_val < 0)
    {
      /*  Go left          */
      ptr=ptr->left;
    }
    else
    {
      /*  Go right          */
      ptr=ptr->right;
    }
  }

  /*  Make sure we found a match          */
  if (!found)
  {
    char err_msg[128];

    sprintf(err_msg, "UNABLE TO FIND ANGLE PARAMETERS FOR %s %s %s (ATOMS %i %i %i)",
       atom1, atom2, atom3, angle_ptr->atom1+1, angle_ptr->atom2+1, angle_ptr->atom3+1);

    if ( notFoundIndex ) {
      angle_ptr->angle_type = notFoundIndex;
      iout << iWARN << err_msg << "\n" << endi;
      return;
    } else NAMD_die(err_msg);
  }

  return;
}

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.

{
  struct bond_params *ptr;  //  Current location in tree
  int found=0;      //  Flag 1-> found a match
  int cmp_code;      //  return code from strcasecmp

  /*  Check to make sure the files have all been read    */
  if (!AllFilesRead)
  {
    NAMD_die("Tried to assign bond index before all parameter files were read");
  }

  /*  We need atom1 < atom2, so if that's not the way they        */
  /*  were passed, flip them          */
  if (strcasecmp(atom1, atom2) > 0)
  {
    const char *tmp_name = atom1;
    atom1 = atom2;
    atom2 = tmp_name;
  }

  /*  Start at the top            */
  ptr=bondp;

  /*  While we haven't found a match and we're not at the end  */
  /*  of the tree, compare the bond passed in with the tree  */
  while (!found && (ptr!=NULL))
  {
    cmp_code=strcasecmp(atom1, ptr->atom1name);

    if (cmp_code == 0)
    {
      cmp_code=strcasecmp(atom2, ptr->atom2name);
    }

    if (cmp_code == 0)
    {
      /*  Found a match        */
      found=1;
      bond_ptr->bond_type = ptr->index;
    }
    else if (cmp_code < 0)
    {
      /*  Go left          */
      ptr=ptr->left;
    }
    else
    {
      /*  Go right          */
      ptr=ptr->right;
    }
  }

  /*  Check to see if we found anything        */
  if (!found)
  {
    if ((strcmp(atom1, "DRUD")==0 || strcmp(atom2, "DRUD")==0)
        && (strcmp(atom1, "X")!=0 && strcmp(atom2, "X")!=0)) {
      /* try a wildcard DRUD X match for this Drude bond */
      char a1[8] = "DRUD", a2[8] = "X";
      return assign_bond_index(a1, a2, bond_ptr);  /* recursive call */
    }
    else {
      char err_msg[128];

      sprintf(err_msg, "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
         atom1, atom2, bond_ptr->atom1+1, bond_ptr->atom2+1);
      NAMD_die(err_msg);
    }
  }

  return;
}

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.

{
  struct crossterm_params *ptr;  //  Current position in list
  int found=0;      //  Flag 1->found a match

  /*  Start at the head of the list        */
  ptr=crosstermp;

  /*  While we haven't fuond a match and haven't reached the end  */
  /*  of the list, keep looking          */
  while (!found && (ptr!=NULL))
  {
    /*  Do a linear search through the linked list of   */
    /*  crossterm parameters.  Since the list is arranged    */
    /*  with wildcard paramters at the end of the list, we  */
    /*  can simply do a linear search and be guaranteed that*/
    /*  we will find exact matches before wildcard matches. */
    /*  Also, we must check for an exact match, and a match */
    /*  in reverse, since they are really the same          */
    /*  physically.            */
    if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
           (strcasecmp(ptr->atom1name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom2name, atom2)==0) || 
           (strcasecmp(ptr->atom2name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom3name, atom3)==0) || 
           (strcasecmp(ptr->atom3name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom4name, atom4)==0) || 
           (strcasecmp(ptr->atom4name, "X")==0) ) )
    {
      /*  Found an exact match      */
      found=1;
    }
    else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
           (strcasecmp(ptr->atom4name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom3name, atom2)==0) || 
           (strcasecmp(ptr->atom3name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom2name, atom3)==0) || 
           (strcasecmp(ptr->atom2name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom1name, atom4)==0) || 
           (strcasecmp(ptr->atom1name, "X")==0) ) )
    {
      /*  Found a reverse match      */
      found=1;
    }
    if ( ! found ) {
      /*  Didn't find a match, go to the next node  */
      ptr=ptr->next;
      continue;
    }
    found = 0;
    if ( ( (strcasecmp(ptr->atom5name, atom5)==0) || 
           (strcasecmp(ptr->atom5name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom6name, atom6)==0) || 
           (strcasecmp(ptr->atom6name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom7name, atom7)==0) || 
           (strcasecmp(ptr->atom7name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom8name, atom8)==0) || 
           (strcasecmp(ptr->atom8name, "X")==0) ) )
    {
      /*  Found an exact match      */
      found=1;
    }
    else if ( ( (strcasecmp(ptr->atom8name, atom5)==0) || 
           (strcasecmp(ptr->atom8name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom7name, atom6)==0) || 
           (strcasecmp(ptr->atom7name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom6name, atom7)==0) || 
           (strcasecmp(ptr->atom6name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom5name, atom8)==0) || 
           (strcasecmp(ptr->atom5name, "X")==0) ) )
    {
      /*  Found a reverse match      */
      found=1;
    }
    if ( ! found ) {
      /*  Didn't find a match, go to the next node  */
      ptr=ptr->next;
    }
  }

  /*  Make sure we found a match          */
  if (!found)
  {
    char err_msg[128];

    sprintf(err_msg, "UNABLE TO FIND CROSSTERM PARAMETERS FOR %s  %s  %s  %s  %s  %s  %s  %s\n"
      "(ATOMS %i %i %i %i %i %i %i %i)",
      atom1, atom2, atom3, atom4, atom5, atom6, atom7, atom8,
      crossterm_ptr->atom1+1, crossterm_ptr->atom1+2, crossterm_ptr->atom1+3, crossterm_ptr->atom4+1,
      crossterm_ptr->atom1+5, crossterm_ptr->atom1+6, crossterm_ptr->atom1+7, crossterm_ptr->atom8+1);
    
    NAMD_die(err_msg);
  }

  /*  Assign the constants          */
  crossterm_ptr->crossterm_type = ptr->index;

  return;
}

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, endi(), dihedral_params::index, iout, iWARN(), dihedral_params::multiplicity, DihedralValue::multiplicity, NAMD_die(), dihedral_params::next, and paraXplor.

{
  struct dihedral_params *ptr;  //  Current position in list
  int found=0;      //  Flag 1->found a match

  /*  Start at the begining of the list        */
  ptr=dihedralp;

  /*  While we haven't found a match and we haven't reached       */
  /*  the end of the list, keep looking        */
  while (!found && (ptr!=NULL))
  {
    /*  Do a linear search through the linked list of   */
    /*  dihedral parameters.  Since the list is arranged    */
    /*  with wildcard paramters at the end of the list, we  */
    /*  can simply do a linear search and be guaranteed that*/
    /*  we will find exact matches before wildcard matches. */
    /*  Also, we must check for an exact match, and a match */
    /*  in reverse, since they are really the same          */
    /*  physically.            */
    if ( ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom1)==0) ) && 
         ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom2)==0) ) &&
         ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom3)==0) ) &&
         ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom4)==0) ) ) 
    {
      /*  Found an exact match      */
      found=1;
    }
    else if ( ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom1)==0) ) &&
              ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom2)==0) ) &&
              ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom3)==0) ) &&
              ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom4)==0) ) )
    {
      /*  Found a reverse match      */
      found=1;
    }
    else
    {
      /*  Didn't find a match, go to the next node  */
      ptr=ptr->next;
    }
  }

  /*  Make sure we found a match          */
  if (!found)
  {
    char err_msg[128];

    sprintf(err_msg, "UNABLE TO FIND DIHEDRAL PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
       atom1, atom2, atom3, atom4, dihedral_ptr->atom1+1, dihedral_ptr->atom2+1,
       dihedral_ptr->atom3+1, dihedral_ptr->atom4+1);
    
    if ( notFoundIndex ) {
      dihedral_ptr->dihedral_type = notFoundIndex;
      iout << iWARN << err_msg << "\n" << endi;
      return;
    } else NAMD_die(err_msg);
  }

 if (paramType == paraXplor) {
  //  Check to make sure the number of multiples specified in the psf
  //  file doesn't exceed the number of parameters in the parameter
  //  files
  if (multiplicity > maxDihedralMults[ptr->index])
  {
    char err_msg[257];

    sprintf(err_msg, "Multiplicity of Paramters for diehedral bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxDihedralMults[ptr->index]);
    NAMD_die(err_msg);
  }

  //  If the multiplicity from the current bond is larger than that
  //  seen in the past, increase the multiplicity for this bond
  if (multiplicity > dihedral_array[ptr->index].multiplicity)
  {
    dihedral_array[ptr->index].multiplicity = multiplicity;
  }
 }

  dihedral_ptr->dihedral_type = ptr->index;

  return;
}

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, improper_params::multiplicity, ImproperValue::multiplicity, NAMD_die(), improper_params::next, and paraXplor.

{
  struct improper_params *ptr;  //  Current position in list
  int found=0;      //  Flag 1->found a match

  /*  Start at the head of the list        */
  ptr=improperp;

  /*  While we haven't fuond a match and haven't reached the end  */
  /*  of the list, keep looking          */
  while (!found && (ptr!=NULL))
  {
    /*  Do a linear search through the linked list of   */
    /*  improper parameters.  Since the list is arranged    */
    /*  with wildcard paramters at the end of the list, we  */
    /*  can simply do a linear search and be guaranteed that*/
    /*  we will find exact matches before wildcard matches. */
    /*  Also, we must check for an exact match, and a match */
    /*  in reverse, since they are really the same          */
    /*  physically.            */
    if ( ( (strcasecmp(ptr->atom1name, atom1)==0) || 
           (strcasecmp(ptr->atom1name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom2name, atom2)==0) || 
           (strcasecmp(ptr->atom2name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom3name, atom3)==0) || 
           (strcasecmp(ptr->atom3name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom4name, atom4)==0) || 
           (strcasecmp(ptr->atom4name, "X")==0) ) )
    {
      /*  Found an exact match      */
      found=1;
    }
    else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) || 
           (strcasecmp(ptr->atom4name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom3name, atom2)==0) || 
           (strcasecmp(ptr->atom3name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom2name, atom3)==0) || 
           (strcasecmp(ptr->atom2name, "X")==0) ) &&
       ( (strcasecmp(ptr->atom1name, atom4)==0) || 
           (strcasecmp(ptr->atom1name, "X")==0) ) )
    {
      /*  Found a reverse match      */
      found=1;
    }
    else
    {
      /*  Didn't find a match, go to the next node  */
      ptr=ptr->next;
    }
  }

  /*  Make sure we found a match          */
  if (!found)
  {
    char err_msg[128];

    sprintf(err_msg, "UNABLE TO FIND IMPROPER PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
       atom1, atom2, atom3, atom4, improper_ptr->atom1+1, improper_ptr->atom2+1,
       improper_ptr->atom3+1, improper_ptr->atom4+1);
    
    NAMD_die(err_msg);
  }

 if (paramType == paraXplor) {
  //  Check to make sure the number of multiples specified in the psf
  //  file doesn't exceed the number of parameters in the parameter
  //  files
  if (multiplicity > maxImproperMults[ptr->index])
  {
    char err_msg[257];

    sprintf(err_msg, "Multiplicity of Paramters for improper bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxImproperMults[ptr->index]);
    NAMD_die(err_msg);
  }

  //  If the multiplicity from the current bond is larger than that
  //  seen in the past, increase the multiplicity for this bond
  if (multiplicity > improper_array[ptr->index].multiplicity)
  {
    improper_array[ptr->index].multiplicity = multiplicity;
  }
 }

  /*  Assign the constants          */
  improper_ptr->improper_type = ptr->index;

  return;
}

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

Definition at line 3470 of file Parameters.C.

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

Referenced by Molecule::prepare_qm().

{
  struct vdw_params *ptr;    //  Current position in trees
  int found=0;      //  Flag 1->found match
  int comp_code;      //  return code from strcasecmp

  /*  Check to make sure the files have all been read    */
  if (!AllFilesRead)
  {
    NAMD_die("Tried to assign vdw index before all parameter files were read");
  }

  /*  Start at the top            */
  ptr=vdwp;

  /*  While we haven't found a match, and we haven't reached      */
  /*  the bottom of the tree, compare the atom passed in with     */
  /*  the current value and decide if we have a match, or if not, */
  /*  which way to go            */
  while (!found && (ptr!=NULL))
  {
    comp_code = strcasecmp(atomtype, ptr->atomname);

    if (comp_code == 0)
    {
      /*  Found a match!        */
      atom_ptr->vdw_type=ptr->index;
      found=1;
    }
    else if (comp_code < 0)
    {
      /*  Go to the left        */
      ptr=ptr->left;
    }
    else
    {
      /*  Go to the right        */
      ptr=ptr->right;
    }
  }

  //****** BEGIN CHARMM/XPLOR type changes
  if (!found)
  {
    // since CHARMM allows wildcards "*" in vdw typenames
    // we have to look again if necessary, this way, if
    // we already had an exact match, this is never executed
    size_t windx;                      //  wildcard index

    /*  Start again at the top                                */
    ptr=vdwp;
  
     while (!found && (ptr!=NULL))
     {
  
       // get index of wildcard wildcard, get index
       windx= strcspn(ptr->atomname,"*"); 
       if (windx == strlen(ptr->atomname))
       {
         // there is no wildcard here
         comp_code = strcasecmp(atomtype, ptr->atomname);   
       }
       else
       {
         comp_code = strncasecmp(atomtype, ptr->atomname, windx); 
       }  

       if (comp_code == 0)
       {
         /*  Found a match!                                */
         atom_ptr->vdw_type=ptr->index;
         found=1;
         char errbuf[100];
         sprintf(errbuf,"VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
                        atomtype, ptr->atomname);
         int i;
         for(i=0; i<error_msgs.size(); i++) {
           if ( strcmp(errbuf,error_msgs[i]) == 0 ) break;
         }
         if ( i == error_msgs.size() ) {
           char *newbuf = new char[strlen(errbuf)+1];
           strcpy(newbuf,errbuf);
           error_msgs.add(newbuf);
           iout << iWARN << newbuf << "\n" << endi;
         }
       }
       else if (comp_code < 0)
       {
          /*  Go to the left                                */
                ptr=ptr->left;
       }
       else
       {
         /*  Go to the right                                */
                ptr=ptr->right;
       }
     
     }
                
  }
  //****** END CHARMM/XPLOR type changes

  /*  Make sure we found it          */
  if (!found)
  {
    char err_msg[100];

    sprintf(err_msg, "DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s",
       atomtype);
    NAMD_die(err_msg);
  }

  return;
}

char* Parameters::atom_type_name ( Index  a )

inline

Definition at line 396 of file Parameters.h.

References MAX_ATOMTYPE_CHARS.

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

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

{
  AllFilesRead = TRUE;

  if (addDrudeBond) {
    // default definition for Drude bonds if none given
    NumBondParams++;
    add_bond_param("X DRUD 500.0 0.0\n", FALSE);
  }
  //  Allocate space for all of the arrays
  if (NumBondParams)
  {
    bond_array = new BondValue[NumBondParams];

    if (bond_array == NULL)
    {
      NAMD_die("memory allocation of bond_array failed!");
    }
    memset(bond_array, 0, NumBondParams*sizeof(BondValue));
  }

  if (NumAngleParams)
  {
    angle_array = new AngleValue[NumAngleParams];

    if (angle_array == NULL)
    {
      NAMD_die("memory allocation of angle_array failed!");
    }
    memset(angle_array, 0, NumAngleParams*sizeof(AngleValue));
    for ( Index i=0; i<NumAngleParams; ++i ) {
      angle_array[i].normal = 1;
    }
  }

  if (NumDihedralParams)
  {
    dihedral_array = new DihedralValue[NumDihedralParams];

    if (dihedral_array == NULL)
    {
      NAMD_die("memory allocation of dihedral_array failed!");
    }
    memset(dihedral_array, 0, NumDihedralParams*sizeof(DihedralValue));
  }

  if (NumImproperParams)
  {
    improper_array = new ImproperValue[NumImproperParams];

    if (improper_array == NULL)
    {
      NAMD_die("memory allocation of improper_array failed!");
    }
    memset(improper_array, 0, NumImproperParams*sizeof(ImproperValue));
  }

  if (NumCrosstermParams)
  {
    crossterm_array = new CrosstermValue[NumCrosstermParams];
    memset(crossterm_array, 0, NumCrosstermParams*sizeof(CrosstermValue));
  }

  // JLai
  if (NumGromacsPairParams)
  {
    gromacsPair_array = new GromacsPairValue[NumGromacsPairParams];
    memset(gromacsPair_array, 0, NumGromacsPairParams*sizeof(GromacsPairValue)); 
  }
  // End of JLai

  if (NumVdwParams)
  {
          atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
    vdw_array = new VdwValue[NumVdwParams];
    
    if (vdw_array == NULL)
    {
      NAMD_die("memory allocation of vdw_array failed!");
    }
  }
  if (NumNbtholePairParams)
  {
    nbthole_array = new NbtholePairValue[NumNbtholePairParams];

    if(nbthole_array == NULL)
    {
      NAMD_die("memory allocation of nbthole_array failed!");
    }
  }
  //  Assign indexes to each of the parameters and populate the
  //  arrays using the binary trees and linked lists that we have
  //  already read in

  // Note that if parameters have been overwritten (matching
  // atom patterns but different parameter values) the tree
  // contains fewer elements than Num...Params would suggest.
  // The arrays are initialized above because the end values
  // may not be occupied.  Modifying the Num...Params values
  // would break backwards compatibility of memopt extraBonds.

  index_bonds(bondp, 0);
  index_angles(anglep, 0);
  NumVdwParamsAssigned = index_vdw(vdwp, 0);
  index_dihedrals();
  index_impropers();
  index_crossterms();
  convert_nbthole_pairs();
  //  Convert the vdw pairs
  convert_vdw_pairs();
  convert_table_pairs();
}

void Parameters::done_reading_structure

(

)

Definition at line 4999 of file Parameters.C.

{
  if (bondp != NULL)
    free_bond_tree(bondp);

  if (anglep != NULL)
    free_angle_tree(anglep);

  if (dihedralp != NULL)
    free_dihedral_list(dihedralp);

  if (improperp != NULL)
    free_improper_list(improperp);

  if (crosstermp != NULL)
    free_crossterm_list(crosstermp);

  if (vdwp != NULL)
    free_vdw_tree(vdwp);

  //  Free the arrays used to track multiplicity for dihedrals
  //  and impropers
  if (maxDihedralMults != NULL)
    delete [] maxDihedralMults;

  if (maxImproperMults != NULL)
    delete [] maxImproperMults;

  bondp=NULL;
  anglep=NULL;
  dihedralp=NULL;
  improperp=NULL;
  crosstermp=NULL;
  vdwp=NULL;
  maxImproperMults=NULL;
  maxDihedralMults=NULL;
}

void Parameters::get_angle_params ( Real *  k, Real *  theta0, Real *  k_ub, Real *  r_ub, Index  index  )

inline

Definition at line 456 of file Parameters.h.

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

        {
                *k = angle_array[index].k;
                *theta0 = angle_array[index].theta0;
                *k_ub = angle_array[index].k_ub;
                *r_ub = angle_array[index].r_ub;
        }

void Parameters::get_bond_params ( Real *  k, Real *  x0, Index  index  )

inline

Definition at line 450 of file Parameters.h.

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

Referenced by Molecule::print_bonds().

        {
                *k = bond_array[index].k;
                *x0 = bond_array[index].x0;
        }

int Parameters::get_dihedral_multiplicity ( Index  index )

inline

Definition at line 470 of file Parameters.h.

References dihedral_array.

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

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

inline

Definition at line 488 of file Parameters.h.

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

        {
                if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
                {
                        NAMD_die("Bad mult index in Parameters::get_dihedral_params");
                }

                *k = dihedral_array[index].values[mult].k;
                *n = dihedral_array[index].values[mult].n;
                *delta = dihedral_array[index].values[mult].delta;
        }

int Parameters::get_improper_multiplicity ( Index  index )

inline

Definition at line 465 of file Parameters.h.

References improper_array.

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

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

inline

Definition at line 475 of file Parameters.h.

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

        {
                if ( (mult<0) || (mult>MAX_MULTIPLICITY) )
                {
                        NAMD_die("Bad mult index in Parameters::get_improper_params");
                }

                *k = improper_array[index].values[mult].k;
                *n = improper_array[index].values[mult].n;
                *delta = improper_array[index].values[mult].delta;
        }

int Parameters::get_int_table_type

(

char *

tabletype

)

Definition at line 7724 of file Parameters.C.

References table_types, and tablenumtypes.

                                                  {
  for (int i=0; i<tablenumtypes; i++) {
    if (!strncmp(tabletype, table_types[i], 5)) {
      return i;
    }
  }

  return -1;
}

int Parameters::get_num_vdw_params ( void  )

inline

Definition at line 535 of file Parameters.h.

References NumVdwParamsAssigned.

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

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

                                                                    {
  IndexedTablePair *ptr;
  Index temp;
  int found=FALSE;

  ptr=tab_pair_tree;
  //
  //  We need the smaller type in ind1, so if it isn't already that 
  //  way, switch them        */
  if (ind1 > ind2)
  {
    temp = ind1;
    ind1 = ind2;
    ind2 = temp;
  }

  /*  While we haven't found a match and we're not at the end  */
  /*  of the tree, compare the bond passed in with the tree  */
  while (!found && (ptr!=NULL))
  {
//    printf("Comparing %i with %i and %i with %i\n", ind1, ptr->ind1, ind2, ptr->ind2);
    if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
    {
       found = TRUE;
    }
    else if ( (ind1 < ptr->ind1) || 
        ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
    {
      /*  Go left          */
      ptr=ptr->left;
    }
    else
    {
      /*  Go right          */
      ptr=ptr->right;
    }
  }

  /*  If we found a match, assign the values      */
  if (found)
  {
    *K = ptr->K;
    return(TRUE);
  }
  else
  {
    return(FALSE);
  }
}

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(), and ComputeNonbondedUtil::select().

{
  IndexedVdwPair *ptr;    //  Current location in tree
  Index temp;      //  Temporary value for swithcing
          // values
  int found=FALSE;    //  Flag 1-> found a match

  ptr=vdw_pair_tree;

  //  We need the smaller type in ind1, so if it isn't already that 
  //  way, switch them        */
  if (ind1 > ind2)
  {
    temp = ind1;
    ind1 = ind2;
    ind2 = temp;
  }

  /*  While we haven't found a match and we're not at the end  */
  /*  of the tree, compare the bond passed in with the tree  */
  while (!found && (ptr!=NULL))
  {
    if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
    {
       found = TRUE;
    }
    else if ( (ind1 < ptr->ind1) || 
        ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
    {
      /*  Go left          */
      ptr=ptr->left;
    }
    else
    {
      /*  Go right          */
      ptr=ptr->right;
    }
  }

  /*  If we found a match, assign the values      */
  if (found)
  {
    *A = ptr->A;
    *B = ptr->B;
    *A14 = ptr->A14;
    *B14 = ptr->B14;

    return(TRUE);
  }
  else
  {
    return(FALSE);
  }
}

void Parameters::get_vdw_params ( Real *  sigma, Real *  epsilon, Real *  sigma14, Real *  epsilon14, Index  index  )

inline

Definition at line 501 of file Parameters.h.

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

Referenced by Molecule::print_atoms(), and ComputeNonbondedUtil::select().

  {
    if ( vdw_array ) {
      *sigma = vdw_array[index].sigma;
      *epsilon = vdw_array[index].epsilon;
      *sigma14 = vdw_array[index].sigma14;
      *epsilon14 = vdw_array[index].epsilon14;
    } else {
      // sigma and epsilon from A and B for each vdw type's interaction with itself
      Real A; Real B; Real A14; Real B14;
      if (get_vdw_pair_params(index, index, &A, &B, &A14, &B14) ) {
        if (A && B) {
          *sigma = sqrt(cbrt(A)) / sqrt(cbrt(B));
          *epsilon = B*B / (4*A);
        }
        else {
          *sigma = 0; *epsilon = 0;
        }
        if (A14 && B14) {
          *sigma14 = sqrt(cbrt(A14)) / sqrt(cbrt(B14));
          *epsilon14 = B14*B14 / (4*A14);
        }
        else {
          *sigma14 = 0; *epsilon14 = 0;
        }
      }
      else {NAMD_die ("Function get_vdw_params failed to derive Lennard-Jones sigma and epsilon from A and B values\n");}
    }
  }

void Parameters::print_angle_params

(

)

Definition at line 4865 of file Parameters.C.

References DebugM, and NumAngleParams.

{
  DebugM(3,NumAngleParams << " ANGLE PARAMETERS\n"
      << "*****************************************" );
  traverse_angle_params(anglep);
}

void Parameters::print_bond_params

(

)

Definition at line 4847 of file Parameters.C.

References DebugM, and NumBondParams.

{
  DebugM(3,NumBondParams << " BOND PARAMETERS\n" \
      << "*****************************************"  \
      );

  traverse_bond_params(bondp);
}

void Parameters::print_crossterm_params

(

)

void Parameters::print_dihedral_params

(

)

Definition at line 4881 of file Parameters.C.

References DebugM, and NumDihedralParams.

{
  DebugM(3,NumDihedralParams << " DIHEDRAL PARAMETERS\n" \
      << "*****************************************" );

  traverse_dihedral_params(dihedralp);
}

void Parameters::print_improper_params

(

)

Definition at line 4898 of file Parameters.C.

References DebugM, and NumImproperParams.

{
  DebugM(3,NumImproperParams << " IMPROPER PARAMETERS\n" \
      << "*****************************************" );

  traverse_improper_params(improperp);
}

void Parameters::print_nbthole_pair_params

(

)

Definition at line 4949 of file Parameters.C.

References DebugM, and NumNbtholePairParams.

{
  DebugM(3,NumNbtholePairParams << " NBTHOLE PAIR PARAMETERS\n" \
      << "*****************************************" );

  traverse_nbthole_pair_params(nbthole_pairp);
} 

void Parameters::print_param_summary

(

)

Definition at line 4966 of file Parameters.C.

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

Referenced by NamdState::loadStructure().

{
  iout << iINFO << "SUMMARY OF PARAMETERS:\n" 
       << iINFO << NumBondParams << " BONDS\n" 
       << iINFO << NumAngleParams << " ANGLES\n" << endi;
       if (cosAngles) {
         iout << iINFO << "  " << (NumAngleParams - NumCosAngles) << " HARMONIC\n"
              << iINFO << "  " << NumCosAngles << " COSINE-BASED\n" << endi;
       }
  iout << iINFO << NumDihedralParams << " DIHEDRAL\n"
       << iINFO << NumImproperParams << " IMPROPER\n"
       << iINFO << NumCrosstermParams << " CROSSTERM\n"
       << iINFO << NumVdwParams << " VDW\n"
       << iINFO << NumVdwPairParams << " VDW_PAIRS\n"
       << iINFO << NumNbtholePairParams << " NBTHOLE_PAIRS\n" << endi;
}

void Parameters::print_vdw_pair_params

(

)

Definition at line 4932 of file Parameters.C.

References DebugM, and NumVdwPairParams.

{
  DebugM(3,NumVdwPairParams << " vdW PAIR PARAMETERS\n" \
      << "*****************************************" );

  traverse_vdw_pair_params(vdw_pairp);
}

void Parameters::print_vdw_params

(

)

Definition at line 4915 of file Parameters.C.

References DebugM, and NumVdwParams.

{
  DebugM(3,NumVdwParams << " vdW PARAMETERS\n" \
      << "*****************************************" );

  traverse_vdw_params(vdwp);
}

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

{
  int  par_type=0;         //  What type of parameter are we currently
                           //  dealing with? (vide infra)
  int  skipline;           //  skip this line?
  int  skipall = 0;        //  skip rest of file;
  char buffer[512];           //  Buffer to store each line of the file
  char first_word[512];           //  First word of the current line
  FILE *pfile;                   //  File descriptor for the parameter file

  /*  Check to make sure that we haven't previously been told     */
  /*  that all the files were read                                */
  if (AllFilesRead)
  {
    NAMD_die("Tried to read another parameter file after being told that all files were read!");
  }

  /*  Try and open the file                                        */
  if ( (pfile = fopen(fname, "r")) == NULL)
  {
    char err_msg[256];

    sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
    NAMD_die(err_msg);
  }

  /*  Keep reading in lines until we hit the EOF                        */
  while (NAMD_read_line(pfile, buffer) != -1)
  {
    /*  Get the first word of the line                        */
    NAMD_find_first_word(buffer, first_word);
    skipline=0;

    /*  First, screen out things that we ignore.                   */   
    /*  blank lines, lines that start with '!' or '*', lines that  */
    /*  start with "END".                                          */
    if (!NAMD_blank_string(buffer) &&
        (strncmp(first_word, "!", 1) != 0) &&
         (strncmp(first_word, "*", 1) != 0) &&
        (strncasecmp(first_word, "END", 3) != 0))
    {
      if ( skipall ) {
        iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
        break;
      }
      /*  Now, determine the apropriate parameter type.   */
      if (strncasecmp(first_word, "bond", 4)==0)
      {
        par_type=1; skipline=1;
      }
      else if (strncasecmp(first_word, "angl", 4)==0)
      {
        par_type=2; skipline=1;
      }
      else if (strncasecmp(first_word, "thet", 4)==0)
      {
        par_type=2; skipline=1;
      }
      else if (strncasecmp(first_word, "dihe", 4)==0)
      {
        par_type=3; skipline=1;
      }
      else if (strncasecmp(first_word, "phi", 3)==0)
      {
        par_type=3; skipline=1;
      }
      else if (strncasecmp(first_word, "impr", 4)==0)
      {
        par_type=4; skipline=1;
      }
      else if (strncasecmp(first_word, "imph", 4)==0)
      {
        par_type=4; skipline=1;
      }
      else if (strncasecmp(first_word, "nbon", 4)==0)
      {
        par_type=5; skipline=1;
      }
      else if (strncasecmp(first_word, "nonb", 4)==0)
      {
        par_type=5; skipline=1;
      }
      else if (strncasecmp(first_word, "nbfi", 4)==0)
      {
        par_type=6; skipline=1;
      }
      else if (strncasecmp(first_word, "hbon", 4)==0)
      {
        par_type=7; skipline=1;
      }
      else if (strncasecmp(first_word, "cmap", 4)==0)
      {
        par_type=8; skipline=1;
      }
      else if (strncasecmp(first_word, "nbta", 4)==0)
      {
        par_type=9; skipline=1;
      }
      else if (strncasecmp(first_word, "thol", 4)==0)
      {
        par_type=10; skipline=1;
      }
      else if (strncasecmp(first_word, "atom", 4)==0)
      {
        par_type=11; skipline=1;
      }
      else if (strncasecmp(first_word, "ioformat", 8)==0)
      {
        skipline=1;
      }
      else if (strncasecmp(first_word, "read", 4)==0)
      {
        skip_stream_read(buffer, pfile);  skipline=1;
      }
      else if (strncasecmp(first_word, "return", 4)==0)
      {
        skipall=1;  skipline=1;
      }
      else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
               (strncasecmp(first_word, "grou", 4) == 0) ||
               (strncasecmp(first_word, "cdie", 4) == 0) ||
               (strncasecmp(first_word, "shif", 4) == 0) ||
               (strncasecmp(first_word, "vgro", 4) == 0) ||
               (strncasecmp(first_word, "vdis", 4) == 0) ||
               (strncasecmp(first_word, "vswi", 4) == 0) ||
               (strncasecmp(first_word, "cutn", 4) == 0) ||
               (strncasecmp(first_word, "ctof", 4) == 0) ||
               (strncasecmp(first_word, "cton", 4) == 0) ||
               (strncasecmp(first_word, "eps", 3) == 0) ||
               (strncasecmp(first_word, "e14f", 4) == 0) ||
               (strncasecmp(first_word, "wmin", 4) == 0) ||
               (strncasecmp(first_word, "aexp", 4) == 0) ||
               (strncasecmp(first_word, "rexp", 4) == 0) ||
               (strncasecmp(first_word, "haex", 4) == 0) ||
               (strncasecmp(first_word, "aaex", 4) == 0) ||
               (strncasecmp(first_word, "noac", 4) == 0) ||
               (strncasecmp(first_word, "hbno", 4) == 0) ||
               (strncasecmp(first_word, "cuth", 4) == 0) ||
               (strncasecmp(first_word, "ctof", 4) == 0) ||
               (strncasecmp(first_word, "cton", 4) == 0) ||
               (strncasecmp(first_word, "cuth", 4) == 0) ||
               (strncasecmp(first_word, "ctof", 4) == 0) ||
               (strncasecmp(first_word, "cton", 4) == 0) ) 
      {
        if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
        {
          char err_msg[512];

          sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
          NAMD_die(err_msg);
        }
        else 
        {
          skipline = 1;
        }
      }        
      else if (par_type == 0)
      {
        /*  This is an unknown paramter.        */
        /*  This is BAD                                */
        char err_msg[512];

        sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
        NAMD_die(err_msg);
      }
    }
    else
    {
      skipline=1;
    }

    if ( (par_type != 0) && (!skipline) )
    {
      /*  Now, call the appropriate function based    */
      /*  on the type of parameter we have                */
      /*  I know, this should really be a switch ...  */
      if (par_type == 1)
      {
        add_bond_param(buffer, TRUE);
        NumBondParams++;
      }
      else if (par_type == 2)
      {
        add_angle_param(buffer);
        NumAngleParams++;
      }
      else if (par_type == 3)
      {
        add_dihedral_param(buffer, pfile);
        NumDihedralParams++;
      }
      else if (par_type == 4)
      {
        add_improper_param(buffer, pfile);
        NumImproperParams++;
      }
      else if (par_type == 5)
      {
        add_vdw_param(buffer);
        NumVdwParams++;
      }
      else if (par_type == 6)
      {
        add_vdw_pair_param(buffer);
        NumVdwPairParams++; 
      }
      else if (par_type == 7)
      {
        add_hb_pair_param(buffer);                  
      }
      else if (par_type == 8)
      {
        add_crossterm_param(buffer, pfile);                  
        NumCrosstermParams++;
      }
      else if (par_type == 9)
      {

        if (table_ener == NULL) {
          continue;
        }

        add_table_pair_param(buffer);                  
        NumTablePairParams++;
      }
      else if (par_type == 10)
      {
        add_nbthole_pair_param(buffer);
        NumNbtholePairParams++;
      }
      else if (par_type == 11)
      {
        if ( strncasecmp(first_word, "mass", 4) != 0 ) {
          char err_msg[512];
          sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
          NAMD_die(err_msg);
        }
      }
      else
      {
        /*  This really should not occour!      */
        /*  This is an internal error.          */
        /*  This is VERY BAD                        */
        char err_msg[512];

        sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
        NAMD_die(err_msg);
      }
    }
  }

  /*  Close the file                                                */
  fclose(pfile);

  return;
}

void Parameters::read_ener_table

(

SimParameters *

simParams

)

Definition at line 6888 of file Parameters.C.

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

Referenced by Parameters().

                                                         {
        char* table_file = simParams->tabulatedEnergiesFile;
  char* interp_type = simParams->tableInterpType;
        FILE* enertable;
        enertable = fopen(table_file, "r");

        if (enertable == NULL) {
                NAMD_die("ERROR: Could not open energy table file!\n");
        }

        char tableline[121];
  char* newtypename;
  int numtypes;
        int atom1;
        int atom2;
        int distbin;
  int readcount;
        Real dist;
        Real energy;
        Real force;
        Real table_spacing;
        Real maxdist;

/* First find the header line and set the variables we need */
        iout << "Beginning table read\n" << endi;
        while(fgets(tableline,120,enertable)) {
                if (strncmp(tableline,"#",1)==0) {
                        continue;
                }
    readcount = sscanf(tableline, "%i %f %f", &numtypes, &table_spacing, &maxdist);
    if (readcount != 3) {
      NAMD_die("ERROR: Couldn't parse table header line\n");
    }
    break;
  }

  if (maxdist < simParams->cutoff) {
    NAMD_die("Tabulated energies must at least extend to the cutoff distance\n");
  }

        if (maxdist > simParams->cutoff) {
                iout << "Info: Reducing max table size to match nonbond cutoff\n";
                maxdist = ceil(simParams->cutoff);
        }

/* Now allocate memory for the table; we know what we should be getting */
        numenerentries = 2 * numtypes * int (namdnearbyint(maxdist/table_spacing) + 1);
        //Set up a full energy lookup table from a file
        //Allocate the table; layout is atom1 x atom2 x distance energy force
        fprintf(stdout, "Table has %i entries\n",numenerentries);
        //iout << "Allocating original energy table\n" << endi;
        if ( table_ener ) delete [] table_ener;
        table_ener = new BigReal[numenerentries];
  if ( table_types ) delete [] table_types;
  table_types = new char*[numtypes];

/* Finally, start reading the table */
  int numtypesread = 0; //Number of types read so far

        while(fgets(tableline,120,enertable)) {
                if (strncmp(tableline,"#",1)==0) {
                        continue;
                }
    if (strncmp(tableline,"TYPE",4)==0) {
      // Read a new type
      newtypename = new char[5];
      int readcount = sscanf(tableline, "%*s %s", newtypename);
      if (readcount != 1) {
        NAMD_die("Couldn't read interaction type from TYPE line\n");
      }
//      printf("Setting table_types[%i] to %s\n", numtypesread, newtypename);
      table_types[numtypesread] = newtypename;

      if (numtypesread == numtypes) {
        NAMD_die("Error: Number of types in table doesn't match table header\n");
      }

      // Read the current energy type with the proper interpolation
      if (!strncasecmp(interp_type, "linear", 6)) {
        if (read_energy_type(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
          char err_msg[512];
          sprintf(err_msg, "Failed to read table energy (linear) type %s\n", newtypename);
          NAMD_die(err_msg);
        }
      } else if (!strncasecmp(interp_type, "cubic", 5)) {
        if (read_energy_type_bothcubspline(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
          char err_msg[512];
          sprintf(err_msg, "Failed to read table energy (cubic) type %s\n", newtypename);
          NAMD_die(err_msg);
        }
      } else {
        NAMD_die("Table interpolation type must be linear or cubic\n");
      }

      numtypesread++;
      continue;
    }
    //char err_msg[512];
    //sprintf(err_msg, "Unrecognized line in energy table file: %s\n", tableline);
    //NAMD_die(err_msg);
  }

  // See if we got what we expected
  if (numtypesread != numtypes) {
    char err_msg[512];
    sprintf(err_msg, "ERROR: Expected %i tabulated energy types but got %i\n", numtypes, numtypesread);
    NAMD_die(err_msg);
  }

  // Move the necessary information into simParams
  simParams->tableNumTypes = numtypes;
  simParams->tableSpacing = table_spacing;
  simParams->tableMaxDist = maxdist;
  tablenumtypes = numtypes;

  /*
xxxxxx
        int numtypes = simParams->tableNumTypes;

        //This parameter controls the table spacing
        BigReal table_spacing = 0.01;
        BigReal maxdist = 20.0;
        if (maxdist > simParams->cutoff) {
                iout << "Info: Reducing max table size to match nonbond cutoff\n";
                maxdist = ceil(simParams->cutoff);
        }

        numenerentries = (numtypes + 1) * numtypes * int (ceil(maxdist/table_spacing));
//      fprintf(stdout, "Table arithmetic: maxdist %f, table_spacing %f, numtypes %i, numentries %i\n", maxdist, table_spacing, numtypes, numenerentries);
        columnsize = 2 * namdnearbyint(maxdist/table_spacing);
        rowsize = numtypes * columnsize;
        //Set up a full energy lookup table from a file
        //Allocate the table; layout is atom1 x atom2 x distance energy force
        fprintf(stdout, "Table has %i entries\n",numenerentries);
        //iout << "Allocating original energy table\n" << endi;
        if ( table_ener ) delete [] table_ener;
        table_ener = new Real[numenerentries];
        //
        //Set sentinel values for uninitialized table energies
        for (int i=0 ; i<numenerentries ; i++) {
                table_ener[i] = 1337.0;
        }
        Real compval = 1337.0;

        //    iout << "Entering table reading\n" << endi;
        //iout << "Finished allocating table\n" << endi;

        //Counter to make sure we read the right number of energy entries
        int numentries = 0;

        //Now, start reading from the file
        char* table_file = simParams->tabulatedEnergiesFile;
        FILE* enertable;
        enertable = fopen(table_file, "r");

        if (enertable == NULL) {
                NAMD_die("ERROR: Could not open energy table file!\n");
        }

        char tableline[121];
        int atom1;
        int atom2;
        int distbin;
        Real dist;
        Real energy;
        Real force;

        iout << "Beginning table read\n" << endi;
        //Read the table entries
        while(fgets(tableline,120,enertable)) {
//              IOut << "Processing line " << tableline << "\n" << endi;
                if (strncmp(tableline,"#",1)==0) {
                        continue;
                }


                sscanf(tableline, "%i %i %f %f %f\n", &atom1, &atom2, &dist, &energy, &force);
                distbin = int(namdnearbyint(dist/table_spacing));
//              fprintf(stdout, "Working on atoms %i and %i at distance %f\n",atom1,atom2,dist);
                if ((atom2 > atom1) || (distbin > int(namdnearbyint(maxdist/table_spacing)))) {
//                      fprintf(stdout,"Rejected\n");
//                      fprintf(stdout, "Error: Throwing out energy line beyond bounds\n");
        //              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)));
                } else {
                        //The magic formula for the number of columns from previous rows
                        //in the triangular matrix is (2ni+i-i**2)/2
                        //Here n is the number of types, and i is atom2
//                      fprintf(stdout, "Input: atom1 %f atom2 %f columnsize %f \n", float(atom1), float(atom2), float(columnsize));
//                      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);
                        int eneraddress = columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2) + columnsize*(atom1-atom2) + 2*distbin - 2;
                        int forceaddress = eneraddress + 1;
//                              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]);
                fflush(stdout);
//                      fprintf(stdout, "Checking for dupes: Looking at: %f %f \n", table_ener[eneraddress], table_ener[forceaddress]);
                        if ((table_ener[eneraddress] == compval && table_ener[forceaddress] == compval)) {
                                numentries++;
                                table_ener[eneraddress] = energy;
                                table_ener[forceaddress] = force;
//                              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]);
                                //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin] = energy;
                                //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin + 1] = force;
//                              fflush(stdout);
                        } else {
//                              fprintf(stdout,"Rejecting duplicate entry\n");
                        }
                }
                //      iout << "Done with line\n"<< endi;
        }

        //    int should = numtypes * numtypes * (maxdist/table_spacing);
        //    cout << "Entries: " << numentries << " should be " << should << "\n" << endi;
//      int exptypes = ceil((numtypes+1) * numtypes * (maxdist/table_spacing));
//fprintf(stdout,"numtypes: %i maxdist: %f table_spacing: %f exptypes: %i\n",numtypes,maxdist,table_spacing);
        if (numentries != int(numenerentries/2)) {
                fprintf(stdout,"ERROR: Expected %i entries but got %i\n",numenerentries/2, numentries);
                NAMD_die("Number of energy table entries did not match expected value\n");
        }
        //      iout << "Done with table\n"<< endi;
        fprintf(stdout, "Numenerentries: %i\n",numenerentries/2);
  */
} /* 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 7606 of file Parameters.C.

References NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

                                                                                                                                          {

  char tableline[120];

  /* Last values read from table */
  BigReal readdist;
  BigReal readener;
  BigReal readforce;
  BigReal lastdist;
  BigReal lastener;
  BigReal lastforce;
  readdist = -1.0;
  readener = 0.0;
  readforce = 0.0;

  /* Keep track of where in the table we are */
  float currdist;
  int distbin;
  currdist = -1.0;
  distbin = -1;

        while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
    printf("At distance %f + %f vs. %f\n", currdist, table_spacing, maxdist);
                if (strncmp(tableline,"#",1)==0) {
                        continue;
                }
    if (strncmp(tableline,"TYPE",4)==0) {
      break;
    }

    // Read an energy line from the table
    lastdist = readdist;
    lastener = readener;
    lastforce = readforce;
    int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
    if (distbin == -1) {
      if (readdist != 0.0) {
        NAMD_die("ERROR: Energy/force tables must start at d=0.0\n");
      } else {
        distbin = 0;
        continue;
      }
    }
   // printf("Read an energy line: %f %f %f\n", readdist, readener, readforce);
    if (readcount != 3) {
      char err_msg[512];
      sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
      NAMD_die(err_msg);
    }

    //Sanity check the current entry
    if (readdist < lastdist) {
      NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
    }

    currdist = lastdist;

    while (currdist <= readdist && distbin <= (int) (namdnearbyint(maxdist / table_spacing))) {
      distbin = (int) (namdnearbyint(currdist / table_spacing));
      int table_loc = 2 * (distbin + (typeindex * (1 + (int) namdnearbyint(maxdist / table_spacing))));
      printf("Doing interpolation for energy between %f %f and %f %f: Dist %f\n", readener, readdist, lastener, lastdist, currdist);
      table_ener[table_loc] = interp_lin(readener, lastener, readdist, lastdist, currdist);
      table_ener[table_loc + 1] = interp_lin(readforce, lastforce, readdist, lastdist, currdist);
      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);
      currdist += table_spacing;
      distbin++;
    }
  }

  // Go back one line, since we should now be into the next TYPE block
  fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR); 

  // Clean up and make sure everything worked ok
  distbin--;
  printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
  if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
  return 0;
}

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

Definition at line 7131 of file Parameters.C.

References INDEX, NAMD_die(), and namdnearbyint.

Referenced by read_ener_table().

                                                                                                                                                        {

  char tableline[120];
  int i,j;

  /* Last values read from table */
  BigReal readdist;
  BigReal readener;
  BigReal readforce;
  BigReal spacing;
//  BigReal readforce;
  BigReal lastdist;
//  BigReal lastener;
//  BigReal lastforce;
//  readdist = -1.0;
//  readener = 0.0;
//  readforce = 0.0;

  /* Create arrays for holding the input values */
  std::vector<BigReal>  dists;
  std::vector<BigReal> enervalues;
  std::vector<BigReal> forcevalues;
  int numentries = 0;


  /* Keep track of where in the table we are */
  BigReal currdist;
  int distbin;
  currdist = 0.0;
  lastdist = -1.0;
  distbin = 0;

  // Read all of the values first -- we'll interpolate later
        while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
                if (strncmp(tableline,"#",1)==0) {
                        continue;
                }
    if (strncmp(tableline,"TYPE",4)==0) {
      fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR); 
      break;
    }

    // Read an energy line from the table
    int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);

    //printf("Read an energy line: %g %g %g\n", readdist, readener, readforce);
    if (readcount != 3) {
      char err_msg[512];
      sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
      NAMD_die(err_msg);
    }

    //Sanity check the current entry
    if (readdist < lastdist) {
      NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
    }

    lastdist = readdist;
    dists.push_back(readdist);
    enervalues.push_back(readener);
    forcevalues.push_back(readforce);
    numentries++;
  }

  // Check the spacing and make sure it is uniform
  if (dists[0] != 0.0) {
    NAMD_die("ERROR: First data point for energy table must be at r=0\n");
  }
  spacing = dists[1] - dists[0];
  for (i=2; i<(numentries - 1); i++) {
    BigReal myspacing;
    myspacing = dists[i] - dists[i-1];
    if (fabs(myspacing - spacing) > 0.00001) {
      printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
      NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
    }
  }

// Perform cubic spline interpolation to get the energies and forces

  /* allocate spline coefficient matrix */
  // xe and xf / be and bf for energies and forces, respectively
  double* m = new double[numentries*numentries];
  double* xe = new double[numentries];
  double* xf = new double[numentries];
  double* be = new double[numentries];
  double* bf = new double[numentries];

  be[0] = 3 * (enervalues[1] - enervalues[0]);
  for (i=1; i < (numentries - 1); i++) {
//    printf("Control point %i at %f\n", i, enervalues[i]);
    be[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
//    printf("be is %f\n", be[i]);
  }
  be[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);

  bf[0] = 3 * (forcevalues[1] - forcevalues[0]);
  for (i=1; i < (numentries - 1); i++) {
//    printf("Control point %i at %f\n", i, forcevalues[i]);
    bf[i] = 3 * (forcevalues[i+1] - forcevalues[i-1]);
//    printf("bf is %f\n", bf[i]);
  }
  bf[numentries - 1] = 3 * (forcevalues[numentries - 1] - forcevalues[numentries - 2]);

  memset(m,0,numentries*numentries*sizeof(double));

  /* initialize spline coefficient matrix */
  m[0] = 2;
  for (i = 1;  i < numentries;  i++) {
    m[INDEX(numentries,i-1,i)] = 1;
    m[INDEX(numentries,i,i-1)] = 1;
    m[INDEX(numentries,i,i)] = 4;
  }
  m[INDEX(numentries,numentries-1,numentries-1)] = 2;

  /* Now we need to solve the equation M X = b for X */
  // Do this simultaneously for energy and force -- ONLY because we use the same matrix

  //Put m in upper triangular form and apply corresponding operations to b
  for (i=0; i<numentries; i++) {
    // zero the ith column in all rows below i
    const BigReal baseval = m[INDEX(numentries,i,i)];
    for (j=i+1; j<numentries; j++) {
      const BigReal myval = m[INDEX(numentries,j,i)];
      const BigReal factor = myval / baseval;

      for (int k=i; k<numentries; k++) {
        const BigReal subval = m[INDEX(numentries,i,k)];
        m[INDEX(numentries,j,k)] -= (factor * subval);
      }

      be[j] -= (factor * be[i]);
      bf[j] -= (factor * bf[i]);

    }
  }

  //Now work our way diagonally up from the bottom right to assign values to X
  for (i=numentries-1; i>=0; i--) {

    //Subtract the effects of previous columns
    for (j=i+1; j<numentries; j++) {
      be[i] -= ( m[INDEX(numentries,i,j)] * xe[j] );
      bf[i] -= ( m[INDEX(numentries,i,j)] * xf[j] );
    }

    xe[i] = be[i] / m[INDEX(numentries,i,i)];
    xf[i] = bf[i] / m[INDEX(numentries,i,i)];

  }

  // We now have the coefficient information we need to make the table
  // Now interpolate on each interval we want

  distbin = 0;
  int entriesperseg = (int) ceil(spacing / table_spacing);
  int table_prefix = 2 * typeindex * (int) (namdnearbyint(maxdist / table_spacing) + 1);

  for (i=0; i<numentries-1; i++) {
    BigReal Ae,Be,Ce,De;
    BigReal Af,Bf,Cf,Df;
    currdist = dists[i];

//    printf("Interpolating on interval %i\n", i);

    // Set up the coefficients for this section
    Ae = enervalues[i];
    Be = xe[i];
    Ce = 3 * (enervalues[i+1] - enervalues[i]) - (2 * xe[i]) - xe[i+1];
    De = 2 * (enervalues[i] - enervalues[i+1]) + xe[i] + xe[i+1];

    Af = forcevalues[i];
    Bf = xf[i];
    Cf = 3 * (forcevalues[i+1] - forcevalues[i]) - (2 * xf[i]) - xf[i+1];
    Df = 2 * (forcevalues[i] - forcevalues[i+1]) + xf[i] + xf[i+1];

    // Go over the region of interest and fill in the table
    for (j=0; j<entriesperseg; j++) {
      const BigReal mydist = currdist + (j * table_spacing);
      const BigReal mydistfrac = (float) j / (entriesperseg - 1);
      distbin = (int) namdnearbyint(mydist / table_spacing);
      if (distbin >= (int) namdnearbyint(maxdist / table_spacing)) break;
      BigReal energy;
      BigReal force;

      energy = Ae + (Be * mydistfrac) + (Ce * mydistfrac * mydistfrac) + (De * mydistfrac * mydistfrac * mydistfrac);
      force = Af + (Bf * mydistfrac) + (Cf * mydistfrac * mydistfrac) + (Df * mydistfrac * mydistfrac * mydistfrac);

//      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);
      table_ener[table_prefix + 2 * distbin] = energy;
      table_ener[table_prefix + 2 * distbin + 1] = force;
      distbin++;
    }
    if (currdist >= maxdist) break;
  }

  //The procedure above leaves out the last entry -- add it explicitly
  distbin = (int) namdnearbyint(maxdist / table_spacing);
//  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));
  table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
  table_ener[table_prefix + 2 * distbin + 1] = 0.0;
  distbin++;


  // Clean up and make sure everything worked ok
  delete m;
  delete xe;
  delete xf;
  delete be;
  delete bf;
  distbin--;
  printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
  if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
  return 0;
} /* 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 7367 of file Parameters.C.

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

                                                                                                                                                    {

  char tableline[120];
  int i,j;

  /* Last values read from table */
  BigReal readdist;
  BigReal readener;
  BigReal spacing;
//  BigReal readforce;
  BigReal lastdist;
//  BigReal lastener;
//  BigReal lastforce;
//  readdist = -1.0;
//  readener = 0.0;
//  readforce = 0.0;

  /* Create arrays for holding the input values */
  std::vector<BigReal>  dists;
  std::vector<BigReal> enervalues;
  int numentries = 0;


  /* Keep track of where in the table we are */
  BigReal currdist;
  int distbin;
  currdist = 0.0;
  lastdist = -1.0;
  distbin = 0;

  // Read all of the values first -- we'll interpolate later
        while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
                if (strncmp(tableline,"#",1)==0) {
                        continue;
                }
    if (strncmp(tableline,"TYPE",4)==0) {
      fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR); 
      break;
    }

    // Read an energy line from the table
    int readcount = sscanf(tableline, "%lf %lf", &readdist, &readener);

   // printf("Read an energy line: %f %f %f\n", readdist, readener, readforce);
    if (readcount != 2) {
      char err_msg[512];
      sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
      NAMD_die(err_msg);
    }

    //Sanity check the current entry
    if (readdist < lastdist) {
      NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
    }

    lastdist = readdist;
    dists.push_back(readdist);
    enervalues.push_back(readener);
    numentries++;
  }

  // Check the spacing and make sure it is uniform
  if (dists[0] != 0.0) {
    NAMD_die("ERROR: First data point for energy table must be at r=0\n");
  }
  spacing = dists[1] - dists[0];
  for (i=2; i<(numentries - 1); i++) {
    BigReal myspacing;
    myspacing = dists[i] - dists[i-1];
    if (fabs(myspacing - spacing) > 0.00001) {
      printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
      NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
    }
  }

// Perform cubic spline interpolation to get the energies and forces

  /* allocate spline coefficient matrix */
  double* m = new double[numentries*numentries];
  double* x = new double[numentries];
  double* b = new double[numentries];

  b[0] = 3 * (enervalues[1] - enervalues[0]);
  for (i=1; i < (numentries - 1); i++) {
    printf("Control point %i at %f\n", i, enervalues[i]);
    b[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
    printf("b is %f\n", b[i]);
  }
  b[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);

  memset(m,0,numentries*numentries*sizeof(double));

  /* initialize spline coefficient matrix */
  m[0] = 2;
  for (i = 1;  i < numentries;  i++) {
    m[INDEX(numentries,i-1,i)] = 1;
    m[INDEX(numentries,i,i-1)] = 1;
    m[INDEX(numentries,i,i)] = 4;
  }
  m[INDEX(numentries,numentries-1,numentries-1)] = 2;

  /* Now we need to solve the equation M X = b for X */

  printf("Solving the matrix equation: \n");

  for (i=0; i<numentries; i++) {
    printf(" ( ");
    for (j=0; j<numentries; j++) {
      printf(" %6.3f,", m[INDEX(numentries, i, j)]);
    }
    printf(" )  ( D%-3i )  =  ( %6.3f )\n", i, b[i]);
  }

  //Put m in upper triangular form and apply corresponding operations to b
  for (i=0; i<numentries; i++) {
    // zero the ith column in all rows below i
    const BigReal baseval = m[INDEX(numentries,i,i)];
    for (j=i+1; j<numentries; j++) {
      const BigReal myval = m[INDEX(numentries,j,i)];
      const BigReal factor = myval / baseval;

      for (int k=i; k<numentries; k++) {
        const BigReal subval = m[INDEX(numentries,i,k)];
        m[INDEX(numentries,j,k)] -= (factor * subval);
      }

      b[j] -= (factor * b[i]);

    }
  }

  printf(" In upper diagonal form, equation is:\n");
  for (i=0; i<numentries; i++) {
    printf(" ( ");
    for (j=0; j<numentries; j++) {
      printf(" %6.3f,", m[INDEX(numentries, i, j)]);
    }
    printf(" )  ( D%-3i )  =  ( %6.3f )\n", i, b[i]);
  }

  //Now work our way diagonally up from the bottom right to assign values to X
  for (i=numentries-1; i>=0; i--) {

    //Subtract the effects of previous columns
    for (j=i+1; j<numentries; j++) {
      b[i] -= ( m[INDEX(numentries,i,j)] * x[j] );
    }

    x[i] = b[i] / m[INDEX(numentries,i,i)];

  }

  printf(" Solution vector is:\n\t(");
  for (i=0; i<numentries; i++) {
    printf(" %6.3f ", x[i]);
  }
  printf(" ) \n");

  // We now have the coefficient information we need to make the table
  // Now interpolate on each interval we want

  distbin = 0;
  int entriesperseg = (int) ceil(spacing / table_spacing);
  int table_prefix = 2 * typeindex * (int) (namdnearbyint(maxdist / table_spacing) + 1);

  for (i=0; i<numentries-1; i++) {
    BigReal A,B,C,D;
    currdist = dists[i];

    printf("Interpolating on interval %i\n", i);

    // Set up the coefficients for this section
    A = enervalues[i];
    B = x[i];
    C = 3 * (enervalues[i+1] - enervalues[i]) - (2 * x[i]) - x[i+1];
    D = 2 * (enervalues[i] - enervalues[i+1]) + x[i] + x[i+1];

    printf("Coefficients for this interval: %f %f %f %f\n", A, B, C, D);

    // Go over the region of interest and fill in the table
    for (j=0; j<entriesperseg; j++) {
      const BigReal mydist = currdist + (j * table_spacing);
      const BigReal mydistfrac = (float) j / (entriesperseg - 1);
      distbin = (int) namdnearbyint(mydist / table_spacing);
      if (distbin >= (int) namdnearbyint(maxdist / table_spacing)) break;
      BigReal energy;
      BigReal force;

      energy = A + (B * mydistfrac) + (C * mydistfrac * mydistfrac) + (D * mydistfrac * mydistfrac * mydistfrac);
      force = B + (2 * C * mydistfrac) + (3 * D * mydistfrac * mydistfrac);
      // Multiply force by 1 / (interval length)
      force *= (1.0 / spacing);

      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);
      table_ener[table_prefix + 2 * distbin] = energy;
      table_ener[table_prefix + 2 * distbin + 1] = force;
      distbin++;
    }
    if (currdist >= maxdist) break;
  }

  //The procedure above leaves out the last entry -- add it explicitly
  distbin = (int) namdnearbyint(maxdist / table_spacing);
  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));
  table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
  table_ener[table_prefix + 2 * distbin + 1] = 0.0;
  distbin++;


  // Clean up and make sure everything worked ok
  delete m;
  delete x;
  delete b;
  distbin--;
  printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
  if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
  return 0;
} /* 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().

{
  char buffer[512];  //  Buffer to store each line of the file
  char first_word[512];  //  First word of the current line
  FILE *pfile;    //  File descriptor for the parameter file

  /*  Check to make sure that we haven't previously been told     */
  /*  that all the files were read        */
  if (AllFilesRead)
  {
    NAMD_die("Tried to read another parameter file after being told that all files were read!");
  }

  /*  Try and open the file          */
  if ( (pfile = Fopen(fname, "r")) == NULL)
  {
    char err_msg[256];

    sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
    NAMD_die(err_msg);
  }

  /*  Keep reading in lines until we hit the EOF      */
  while (NAMD_read_line(pfile, buffer) != -1)
  {
    /*  Get the first word of the line      */
    NAMD_find_first_word(buffer, first_word);

    /*  First, screen out things that we ignore, such as    */
    /*  blank lines, lines that start with '!', lines that  */
    /*  start with "REMARK", lines that start with set",    */
    /*  and most of the HBOND parameters which include      */
    /*  AEXP, REXP, HAEX, AAEX, but not the HBOND statement */
    /*  which is parsed.                                    */
    if ((buffer[0] != '!') && 
        !NAMD_blank_string(buffer) &&
        (strncasecmp(first_word, "REMARK", 6) != 0) &&
        (strcasecmp(first_word, "set")!=0) &&
        (strncasecmp(first_word, "AEXP", 4) != 0) &&
        (strncasecmp(first_word, "REXP", 4) != 0) &&
        (strncasecmp(first_word, "HAEX", 4) != 0) &&
        (strncasecmp(first_word, "AAEX", 4) != 0) &&
        (strncasecmp(first_word, "NBOND", 5) != 0) &&
        (strncasecmp(first_word, "CUTNB", 5) != 0) &&
        (strncasecmp(first_word, "END", 3) != 0) &&
        (strncasecmp(first_word, "CTONN", 5) != 0) &&
        (strncasecmp(first_word, "EPS", 3) != 0) &&
        (strncasecmp(first_word, "VSWI", 4) != 0) &&
        (strncasecmp(first_word, "NBXM", 4) != 0) &&
        (strncasecmp(first_word, "INHI", 4) != 0) )
    {
      /*  Now, call the appropriate function based    */
      /*  on the type of parameter we have    */
      if (strncasecmp(first_word, "bond", 4)==0)
      {
        add_bond_param(buffer, TRUE);
        NumBondParams++;
      }
      else if (strncasecmp(first_word, "angl", 4)==0)
      {
        add_angle_param(buffer);
        NumAngleParams++;
      }
      else if (strncasecmp(first_word, "dihe", 4)==0)
      {
        add_dihedral_param(buffer, pfile);
        NumDihedralParams++;
      }
      else if (strncasecmp(first_word, "impr", 4)==0)
      {
        add_improper_param(buffer, pfile);
        NumImproperParams++;
      }
      else if (strncasecmp(first_word, "nonb", 4)==0)
      {
        add_vdw_param(buffer);
        NumVdwParams++; 
      }
      else if (strncasecmp(first_word, "nbfi", 4)==0)
      {
        add_vdw_pair_param(buffer);
        NumVdwPairParams++; 
      }
      else if (strncasecmp(first_word, "nbta", 4)==0)
      {

        if (table_ener == NULL) {
          continue;
        }

        add_table_pair_param(buffer);
        NumTablePairParams++; 
      }
      else if (strncasecmp(first_word, "hbon", 4)==0)
      {
        add_hb_pair_param(buffer);
      }
      else
      {
        /*  This is an unknown paramter.        */
        /*  This is BAD        */
        char err_msg[512];

        sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
           fname, buffer);
        NAMD_die(err_msg);
      }
    }
  }

  /*  Close the file            */
  Fclose(pfile);

  return;
}

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(), BondValue::k, AngleValue::k, four_body_consts::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, DihedralValue::values, values, vdw_pair_tree, and BondValue::x0.

{ 
  int i,j,ico,numtype,mul;
  IndexedVdwPair *new_node;

  if (!amber_data->data_read)
    NAMD_die("No data read from parm file yet!");

  // Copy bond parameters
  NumBondParams = amber_data->Numbnd;
  if (NumBondParams)
  { bond_array = new BondValue[NumBondParams];
    if (bond_array == NULL)
      NAMD_die("memory allocation of bond_array failed!");
  }
  for (i=0; i<NumBondParams; ++i)
  { bond_array[i].k = amber_data->Rk[i];
    bond_array[i].x0 = amber_data->Req[i];
  }

  // Copy angle parameters
  NumAngleParams = amber_data->Numang;
  if (NumAngleParams)
  { angle_array = new AngleValue[NumAngleParams];
    if (angle_array == NULL)
      NAMD_die("memory allocation of angle_array failed!");
  }
  for (i=0; i<NumAngleParams; ++i)
  { angle_array[i].k = amber_data->Tk[i];
    angle_array[i].theta0 = amber_data->Teq[i];
    // Amber has no k_ub and r_ub for angle parameters, so they're set to 0
    angle_array[i].k_ub = angle_array[i].r_ub = 0;
    // All angles are harmonic
    angle_array[i].normal = 1;
  }

  // Copy dihedral parameters
  // Note: If the periodicity is negative, it means the following
  //  entry is another term in a multitermed dihedral, and the
  //  absolute value is the true periodicity; in this case the
  //  following entry in "dihedral_array" should be left empty,
  //  NOT be skipped, in order to keep the next dihedral's index
  //  unchanged.
  NumDihedralParams = amber_data->Nptra;
  if (NumDihedralParams)
  { dihedral_array = new DihedralValue[amber_data->Nptra];
    if (dihedral_array == NULL)
      NAMD_die("memory allocation of dihedral_array failed!");
  }
  mul = 0;
  for (i=0; i<NumDihedralParams; ++i)
  { dihedral_array[i-mul].values[mul].k = amber_data->Pk[i];
    dihedral_array[i-mul].values[mul].n = int(fabs(amber_data->Pn[i])+0.5);
    if (dihedral_array[i-mul].values[mul].n == 0)
    { char err_msg[128];
      sprintf(err_msg, "The periodicity of dihedral # %d is zero!", i+1);
      NAMD_die(err_msg);
    }
    dihedral_array[i-mul].values[mul].delta = amber_data->Phase[i];
    // If the periodicity is positive, it means the following
    // entry is a new dihedral term.
    if (amber_data->Pn[i] > 0)
    { dihedral_array[i-mul].multiplicity = mul+1;
      mul = 0;
    }
    else if (++mul >= MAX_MULTIPLICITY)
    { char err_msg[181];
      sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
         mul+1, MAX_MULTIPLICITY);
      NAMD_die(err_msg);
    }
  }
  if (mul > 0)
    NAMD_die("Negative periodicity in the last dihedral entry!");

  // Copy VDW parameters: AMBER explicitly gives VDW parameters between every
  // 2 atom types, so the data are copied to vdw_pair_tree
  // In AMBER, all 1-4 VDW interactions are divided by factor vdw14
  NumVdwParamsAssigned = numtype = amber_data->Ntypes; // Number of atom types
  NumVdwPairParams = numtype * (numtype+1) / 2;
  for (i=0; i<numtype; ++i)
    for (j=i; j<numtype; ++j)
    { new_node = new IndexedVdwPair;
      if (new_node == NULL)
        NAMD_die("memory allocation of vdw_pair_tree failed!");
      new_node->ind1 = i;
      new_node->ind2 = j;
      new_node->left = new_node->right = NULL;
      // ico is the index of interaction between atom type i and j into
      // the parameter arrays. If it's positive, the interaction is
      // 6-12 VDW; otherwise it's 10-12 H-bond interaction. NAMD doesn't
      // have 10-12 term, so if ico is negative, then the 10-12
      // coefficients must be 0, otherwise die.
      ico = amber_data->Cno[numtype*i+j];
      if (ico>0)
      { new_node->A = amber_data->Cn1[ico-1];
        new_node->A14 = new_node->A / vdw14;
        new_node->B = amber_data->Cn2[ico-1];
        new_node->B14 = new_node->B / vdw14;
      }
      else if (amber_data->HB12[abs(ico)-1]==0.0 && amber_data->HB6[abs(ico)-1]==0.0)
      { new_node->A = new_node->A14 = new_node->B = new_node->B14 = 0.0;
        iout << iWARN << "Encounter 10-12 H-bond term\n";
      }
      else
        NAMD_die("Encounter non-zero 10-12 H-bond term!");
      // Add the node to the binary tree
      vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
    }
}

void Parameters::read_parm	(	AmberParm7Reader::Ambertoppar *	amber_data,
		BigReal	vdw14
	)

Definition at line 6487 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, CrosstermValue::c, AmberParm7Reader::Ambertoppar::CMAPParameter, AmberParm7Reader::Ambertoppar::CMAPResolution, AmberParm7Reader::Ambertoppar::CMAPTypeCount, AmberParm7Reader::Ambertoppar::Cn1, AmberParm7Reader::Ambertoppar::Cn2, AmberParm7Reader::Ambertoppar::Cno, crossterm_array, crossterm_setup(), CrosstermData::d00, four_body_consts::delta, dihedral_array, CrosstermValue::dim, AmberParm7Reader::Ambertoppar::HasCMAP, AmberParm7Reader::Ambertoppar::HasData, AmberParm7Reader::Ambertoppar::HB10, AmberParm7Reader::Ambertoppar::HB12, indexed_vdw_pair::ind1, indexed_vdw_pair::ind2, iout, AmberParm7Reader::Ambertoppar::IsCharmmFF, iWARN(), BondValue::k, AngleValue::k, four_body_consts::k, AngleValue::k_ub, indexed_vdw_pair::left, MAX_MULTIPLICITY, DihedralValue::multiplicity, four_body_consts::n, NAMD_die(), AngleValue::normal, AmberParm7Reader::Ambertoppar::Nptra, AmberParm7Reader::Ambertoppar::Ntypes, AmberParm7Reader::Ambertoppar::Numang, NumAngleParams, AmberParm7Reader::Ambertoppar::Numbnd, NumBondParams, NumCrosstermParams, NumDihedralParams, NumVdwPairParams, NumVdwParamsAssigned, AmberParm7Reader::Ambertoppar::Phase, AmberParm7Reader::Ambertoppar::Pk, AmberParm7Reader::Ambertoppar::Pn, AngleValue::r_ub, AmberParm7Reader::Ambertoppar::Req, indexed_vdw_pair::right, AmberParm7Reader::Ambertoppar::Rk, AmberParm7Reader::Ambertoppar::Teq, AngleValue::theta0, AmberParm7Reader::Ambertoppar::Tk, DihedralValue::values, values, vdw_pair_tree, and BondValue::x0.

{
  int i,j,ico,numtype,mul;
  IndexedVdwPair *new_node;

  if (!amber_data->HasData) {
    NAMD_die("No data read from parm file yet!");
  }

  // Check if we are reading a CHARMBER file
  if (amber_data->IsCharmmFF) {
    NAMD_die("You are using a CHARMM force field in AMBER format generated by CHAMBER or ParmEd.\n"
             "We do not support this format. Please consider using the native format (PSF) for CHARMM force field.");
  }

  // Copy bond parameters
  NumBondParams = amber_data->Numbnd;
  if (NumBondParams)
  { bond_array = new BondValue[NumBondParams];
    if (bond_array == NULL)
      NAMD_die("memory allocation of bond_array failed!");
  }
  for (i=0; i<NumBondParams; ++i)
  { bond_array[i].k = amber_data->Rk[i];
    bond_array[i].x0 = amber_data->Req[i];
  }

  // Copy angle parameters
  NumAngleParams = amber_data->Numang;
  if (NumAngleParams)
  { angle_array = new AngleValue[NumAngleParams];
    if (angle_array == NULL)
      NAMD_die("memory allocation of angle_array failed!");
  }
  for (i=0; i<NumAngleParams; ++i)
  { angle_array[i].k = amber_data->Tk[i];
    angle_array[i].theta0 = amber_data->Teq[i];
    // Amber has no k_ub and r_ub for angle parameters, so they're set to 0
    angle_array[i].k_ub = angle_array[i].r_ub = 0;
    // All angles are harmonic
    angle_array[i].normal = 1;
  }

  // Copy dihedral parameters
  // Note: If the periodicity is negative, it means the following
  //  entry is another term in a multitermed dihedral, and the
  //  absolute value is the true periodicity; in this case the
  //  following entry in "dihedral_array" should be left empty,
  //  NOT be skipped, in order to keep the next dihedral's index
  //  unchanged.
  NumDihedralParams = amber_data->Nptra;
  if (NumDihedralParams)
  { dihedral_array = new DihedralValue[amber_data->Nptra];
    if (dihedral_array == NULL)
      NAMD_die("memory allocation of dihedral_array failed!");
  }
  mul = 0;
  for (i=0; i<NumDihedralParams; ++i)
  { dihedral_array[i-mul].values[mul].k = amber_data->Pk[i];
    dihedral_array[i-mul].values[mul].n = int(fabs(amber_data->Pn[i])+0.5);
    if (dihedral_array[i-mul].values[mul].n == 0)
    { char err_msg[128];
      sprintf(err_msg, "The periodicity of dihedral # %d is zero!", i+1);
      NAMD_die(err_msg);
    }
    dihedral_array[i-mul].values[mul].delta = amber_data->Phase[i];
    // If the periodicity is positive, it means the following
    // entry is a new dihedral term.
    if (amber_data->Pn[i] > 0)
    { dihedral_array[i-mul].multiplicity = mul+1;
      mul = 0;
    }
    else if (++mul >= MAX_MULTIPLICITY)
    { char err_msg[181];
      sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
         mul+1, MAX_MULTIPLICITY);
      NAMD_die(err_msg);
    }
  }
  if (mul > 0)
    NAMD_die("Negative periodicity in the last dihedral entry!");

  // Copy crossterms
  if (amber_data->HasCMAP) {
    NumCrosstermParams = amber_data->CMAPTypeCount;
    if (NumCrosstermParams > 0) {
      crossterm_array = new CrosstermValue[NumCrosstermParams];
      for (i = 0; i < NumCrosstermParams; ++i) {
        // check the resolutions at first
        const int N = amber_data->CMAPResolution[i];
        // NAMD does not support a resolution number other than CrosstermValue::dim - 1
        if (N != CrosstermValue::dim - 1) {
          char err_msg[128];
          sprintf(err_msg, "The table of %d crossterm type has a resolution(%d) other than %d!",
                  i+1, N, CrosstermValue::dim - 1);
          NAMD_die(err_msg);
        }
        // code swipe from Parameters::index_crossterms()
        int l = 0;
        for (int j = 0; j < N; ++j) {
          for (int k = 0; k < N; ++k) {
            crossterm_array[i].c[j][k].d00 = amber_data->CMAPParameter[i][l];
            ++l;
          }
        }
        for (int j = 0; j < N; ++j) {
          crossterm_array[i].c[j][N].d00 = crossterm_array[i].c[j][0].d00;
          crossterm_array[i].c[N][j].d00 = crossterm_array[i].c[0][j].d00;
        }
        crossterm_array[i].c[N][N] = crossterm_array[i].c[0][0];
        crossterm_setup(&crossterm_array[i].c[0][0]);
      }
    }
  }

  // Copy VDW parameters: AMBER explicitly gives VDW parameters between every
  // 2 atom types, so the data are copied to vdw_pair_tree
  // In AMBER, all 1-4 VDW interactions are divided by factor vdw14
  NumVdwParamsAssigned = numtype = amber_data->Ntypes; // Number of atom types
  NumVdwPairParams = numtype * (numtype+1) / 2;
  for (i=0; i<numtype; ++i)
    for (j=i; j<numtype; ++j)
    { new_node = new IndexedVdwPair;
      if (new_node == NULL)
        NAMD_die("memory allocation of vdw_pair_tree failed!");
      new_node->ind1 = i;
      new_node->ind2 = j;
      new_node->left = new_node->right = NULL;
      // ico is the index of interaction between atom type i and j into
      // the parameter arrays. If it's positive, the interaction is
      // 6-12 VDW; otherwise it's 10-12 H-bond interaction. NAMD doesn't
      // have 10-12 term, so if ico is negative, then the 10-12
      // coefficients must be 0, otherwise die.
      ico = amber_data->Cno[numtype*i+j];
      if (ico>0)
      { new_node->A = amber_data->Cn1[ico-1];
        new_node->A14 = new_node->A / vdw14;
        new_node->B = amber_data->Cn2[ico-1];
        new_node->B14 = new_node->B / vdw14;
      }
      else if (amber_data->HB12[abs(ico)-1]==0.0 && amber_data->HB10[abs(ico)-1]==0.0)
      { new_node->A = new_node->A14 = new_node->B = new_node->B14 = 0.0;
        iout << iWARN << "Encounter 10-12 H-bond term\n";
      }
      else
        NAMD_die("Encounter non-zero 10-12 H-bond term!");
      // Add the node to the binary tree
      vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
    }
}

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_vdw_pair::ind1, nbthole_pair_value::ind1, indexed_table_pair::ind1, indexed_vdw_pair::ind2, nbthole_pair_value::ind2, indexed_table_pair::ind2, BondValue::k, AngleValue::k, four_body_consts::k, indexed_table_pair::K, AngleValue::k_ub, indexed_vdw_pair::left, indexed_table_pair::left, MAX_ATOMTYPE_CHARS, MAX_MULTIPLICITY, DihedralValue::multiplicity, ImproperValue::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_vdw_pair::right, indexed_table_pair::right, tab_pair_tree, table_ener, AngleValue::theta0, nbthole_pair_value::tholeij, TRUE, DihedralValue::values, ImproperValue::values, vdw_array, vdw_pair_tree, and BondValue::x0.

Referenced by Node::resendMolecule().

{
  int i, j;      //  Loop counters
  Real *a1, *a2, *a3, *a4;  //  Temporary arrays to get data from message in
  int *i1, *i2, *i3;      //  Temporary int array to get data from message in
  IndexedVdwPair *new_node;  //  New node for vdw pair param tree
  IndexedTablePair *new_tab_node;
  Real **kvals;      //  Force constant values for dihedrals and impropers
  int **nvals;      //  Periodicity values for dihedrals and impropers
  Real **deltavals;    //  Phase shift values for dihedrals and impropers
  BigReal *pairC6, *pairC12;  // JLai

  //  Get the bonded parameters
  msg->get(NumBondParams);

  if (NumBondParams)
  {
    bond_array = new BondValue[NumBondParams];
    a1 = new Real[NumBondParams];
    a2 = new Real[NumBondParams];

    if ( (bond_array == NULL) || (a1 == NULL) || (a2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::receive_Parameters");
    }

    msg->get(NumBondParams, a1);
    msg->get(NumBondParams, a2);

    for (i=0; i<NumBondParams; i++)
    {
      bond_array[i].k = a1[i];
      bond_array[i].x0 = a2[i];
    }

    delete [] a1;
    delete [] a2;
  }

  //  Get the angle parameters
  msg->get(NumAngleParams);

  if (NumAngleParams)
  {
    angle_array = new AngleValue[NumAngleParams];
    a1 = new Real[NumAngleParams];
    a2 = new Real[NumAngleParams];
    a3 = new Real[NumAngleParams];
    a4 = new Real[NumAngleParams];
    i1 = new int[NumAngleParams];

    if ( (angle_array == NULL) || (a1 == NULL) || (a2 == NULL) ||
         (a3 == NULL) || (a4 == NULL) || (i1 == NULL))
    {
      NAMD_die("memory allocation failed in Parameters::receive_Parameters");
    }

    msg->get(NumAngleParams, a1);
    msg->get(NumAngleParams, a2);
    msg->get(NumAngleParams, a3);
    msg->get(NumAngleParams, a4);
    msg->get(NumAngleParams, i1);

    for (i=0; i<NumAngleParams; i++)
    {
      angle_array[i].k = a1[i];
      angle_array[i].theta0 = a2[i];
      angle_array[i].k_ub = a3[i];
      angle_array[i].r_ub = a4[i];
      angle_array[i].normal = i1[i];
    }

    delete [] a1;
    delete [] a2;
    delete [] a3;
    delete [] a4;
    delete [] i1;
  }

  //  Get the dihedral parameters
  msg->get(NumDihedralParams);

  if (NumDihedralParams)
  {
    dihedral_array = new DihedralValue[NumDihedralParams];

    i1 = new int[NumDihedralParams];
    kvals = new Real *[MAX_MULTIPLICITY];
    nvals = new int *[MAX_MULTIPLICITY];
    deltavals = new Real *[MAX_MULTIPLICITY];

    if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
         (deltavals == NULL) || (dihedral_array == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      kvals[i] = new Real[NumDihedralParams];
      nvals[i] = new int[NumDihedralParams];
      deltavals[i] = new Real[NumDihedralParams];

      if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
      {
        NAMD_die("memory allocation failed in Parameters::send_Parameters");
      }
    }

    msg->get(NumDihedralParams, i1);

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      msg->get(NumDihedralParams, kvals[i]);
      msg->get(NumDihedralParams, nvals[i]);
      msg->get(NumDihedralParams, deltavals[i]);
    }

    for (i=0; i<NumDihedralParams; i++)
    {
      dihedral_array[i].multiplicity = i1[i];

      for (j=0; j<MAX_MULTIPLICITY; j++)
      {
        dihedral_array[i].values[j].k = kvals[j][i];
        dihedral_array[i].values[j].n = nvals[j][i];
        dihedral_array[i].values[j].delta = deltavals[j][i];
      }
    }

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      delete [] kvals[i];
      delete [] nvals[i];
      delete [] deltavals[i];
    }

    delete [] i1;
    delete [] kvals;
    delete [] nvals;
    delete [] deltavals;
  }

  //  Get the improper parameters
  msg->get(NumImproperParams);

  if (NumImproperParams)
  {
    improper_array = new ImproperValue[NumImproperParams];
    i1 = new int[NumImproperParams];
    kvals = new Real *[MAX_MULTIPLICITY];
    nvals = new int *[MAX_MULTIPLICITY];
    deltavals = new Real *[MAX_MULTIPLICITY];

    if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
         (deltavals == NULL) || (improper_array==NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      kvals[i] = new Real[NumImproperParams];
      nvals[i] = new int[NumImproperParams];
      deltavals[i] = new Real[NumImproperParams];

      if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
      {
        NAMD_die("memory allocation failed in Parameters::send_Parameters");
      }
    }

    msg->get(NumImproperParams,i1);

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      msg->get(NumImproperParams,kvals[i]);
      msg->get(NumImproperParams,nvals[i]);
      msg->get(NumImproperParams,deltavals[i]);
    }

    for (i=0; i<NumImproperParams; i++)
    {
      improper_array[i].multiplicity = i1[i];

      for (j=0; j<MAX_MULTIPLICITY; j++)
      {
        improper_array[i].values[j].k = kvals[j][i];
        improper_array[i].values[j].n = nvals[j][i];
        improper_array[i].values[j].delta = deltavals[j][i];
      }
    }

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      delete [] kvals[i];
      delete [] nvals[i];
      delete [] deltavals[i];
    }

    delete [] i1;
    delete [] kvals;
    delete [] nvals;
    delete [] deltavals;
  }

  //  Get the crossterm parameters
  msg->get(NumCrosstermParams);

  if (NumCrosstermParams)
  {
    crossterm_array = new CrosstermValue[NumCrosstermParams];

    for (i=0; i<NumCrosstermParams; ++i) {
      int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
      msg->get(nvals,&crossterm_array[i].c[0][0].d00);
    }
  }
  
  // Get GromacsPairs parameters
  // JLai
  msg->get(NumGromacsPairParams);
  
  if (NumGromacsPairParams)
  {
      gromacsPair_array = new GromacsPairValue[NumGromacsPairParams];
      pairC6 = new BigReal[NumGromacsPairParams];
      pairC12 = new BigReal[NumGromacsPairParams];
      
      if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
          NAMD_die("memory allocation failed in Parameters::receive_Parameters");
      }
    
      msg->get(NumGromacsPairParams,pairC6);
      msg->get(NumGromacsPairParams,pairC12);
      
      for (i=0; i<NumGromacsPairParams; ++i) {
          gromacsPair_array[i].pairC6 = pairC6[i];
          gromacsPair_array[i].pairC12 = pairC12[i];
      }

      delete [] pairC6;
      delete [] pairC12;
  }
  // JLai

  //Get the energy table
  msg->get(numenerentries);
  if (numenerentries > 0) {
    //fprintf(stdout, "Getting tables\n");
    //fprintf(infofile, "Trying to allocate table\n");
    table_ener = new BigReal[numenerentries];
    //fprintf(infofile, "Finished table allocation\n");
    if (table_ener==NULL)
    {
      NAMD_die("memory allocation failed in Parameters::receive_Parameters");
    }

    msg->get(numenerentries, table_ener);
  }

  //  Get the vdw parameters
  msg->get(NumVdwParams);
  msg->get(NumVdwParamsAssigned);

  if (NumVdwParams)
  {
          atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
    vdw_array = new VdwValue[NumVdwParams];
    a1 = new Real[NumVdwParams];
    a2 = new Real[NumVdwParams];
    a3 = new Real[NumVdwParams];
    a4 = new Real[NumVdwParams];

    if ( (vdw_array==NULL) || (a1==NULL) || (a2==NULL) || (a3==NULL)
             || (a4==NULL) || (atomTypeNames==NULL))
    {
      NAMD_die("memory allocation failed in Parameters::receive_Parameters");
    }

    msg->get(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
    msg->get(NumVdwParams, a1);
    msg->get(NumVdwParams, a2);
    msg->get(NumVdwParams, a3);
    msg->get(NumVdwParams, a4);

    for (i=0; i<NumVdwParams; i++)
    {
      vdw_array[i].sigma = a1[i];
      vdw_array[i].epsilon = a2[i];
      vdw_array[i].sigma14 = a3[i];
      vdw_array[i].epsilon14 = a4[i];
    }

    delete [] a1;
    delete [] a2;
    delete [] a3;
    delete [] a4;
  }
  
  //  Get the vdw_pair_parameters
  msg->get(NumVdwPairParams);
  
  if (NumVdwPairParams)
  {
    a1 = new Real[NumVdwPairParams];
    a2 = new Real[NumVdwPairParams];
    a3 = new Real[NumVdwPairParams];
    a4 = new Real[NumVdwPairParams];
    i1 = new int[NumVdwPairParams];
    i2 = new int[NumVdwPairParams];

    if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
         (i1 == NULL) || (i2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }
    
    msg->get(NumVdwPairParams, i1);
    msg->get(NumVdwPairParams, i2);
    msg->get(NumVdwPairParams, a1);
    msg->get(NumVdwPairParams, a2);
    msg->get(NumVdwPairParams, a3);
    msg->get(NumVdwPairParams, a4);
    
    for (i=0; i<NumVdwPairParams; i++)
    {
      new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
      
      if (new_node == NULL)
      {
         NAMD_die("memory allocation failed in Parameters::receive_Parameters");
      }
      
      new_node->ind1 = i1[i];
      new_node->ind2 = i2[i];
      new_node->A = a1[i];
      new_node->A14 = a2[i];
      new_node->B = a3[i];
      new_node->B14 = a4[i];
      new_node->left = NULL;
      new_node->right = NULL;
      
      vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
    }
    
    delete [] i1;
    delete [] i2;
    delete [] a1;
    delete [] a2;
    delete [] a3;
    delete [] a4;
  }
 
  //  Get the nbthole_pair_parameters
  msg->get(NumNbtholePairParams); 
    
  if (NumNbtholePairParams)
  {
    nbthole_array = new NbtholePairValue[NumNbtholePairParams];
    a1 = new Real[NumNbtholePairParams];
    a2 = new Real[NumNbtholePairParams];
    a3 = new Real[NumNbtholePairParams];
    i1 = new int[NumNbtholePairParams];
    i2 = new int[NumNbtholePairParams];

    if ( (nbthole_array == NULL) || (a1 == NULL) || (a2 == NULL) || (a3 == NULL)
         || (i1 == NULL) || (i2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::receive_Parameters");
    }

    msg->get(NumNbtholePairParams, i1);
    msg->get(NumNbtholePairParams, i2);
    msg->get(NumNbtholePairParams, a1);
    msg->get(NumNbtholePairParams, a2);
    msg->get(NumNbtholePairParams, a3);

    for (i=0; i<NumNbtholePairParams; i++)
    {

      nbthole_array[i].ind1 = i1[i];
      nbthole_array[i].ind2 = i2[i];
      nbthole_array[i].alphai = a1[i];
      nbthole_array[i].alphaj = a2[i];
      nbthole_array[i].tholeij = a3[i];

    }

    delete [] i1;
    delete [] i2;
    delete [] a1;
    delete [] a2;
    delete [] a3;
  }

  //  Get the table_pair_parameters
  msg->get(NumTablePairParams);
  
  if (NumTablePairParams)
  {
    i1 = new int[NumTablePairParams];
    i2 = new int[NumTablePairParams];
    i3 = new int[NumTablePairParams];

    if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }
    
    msg->get(NumTablePairParams, i1);
    msg->get(NumTablePairParams, i2);
    msg->get(NumTablePairParams, i3);
    
    for (i=0; i<NumTablePairParams; i++)
    {
      new_tab_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
      
      if (new_tab_node == NULL)
      {
         NAMD_die("memory allocation failed in Parameters::receive_Parameters");
      }
      
//      printf("Adding new table pair with ind1 %i ind2 %i k %i\n", i1[i], i2[i],i3[i]);
      new_tab_node->ind1 = i1[i];
      new_tab_node->ind2 = i2[i];
      new_tab_node->K = i3[i];
      new_tab_node->left = NULL;
      new_tab_node->right = NULL;
      
      tab_pair_tree = add_to_indexed_table_pairs(new_tab_node, tab_pair_tree);
    }
    
    delete [] i1;
    delete [] i2;
    delete [] i3;
  }
  
  // receive the hydrogen bond parameters
  // hbondParams.receive_message(msg);

  AllFilesRead = TRUE;

  delete msg;
}

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, BondValue::k, AngleValue::k, four_body_consts::k, AngleValue::k_ub, MAX_ATOMTYPE_CHARS, MAX_MULTIPLICITY, DihedralValue::multiplicity, ImproperValue::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, DihedralValue::values, ImproperValue::values, vdw_array, vdw_pair_tree, and BondValue::x0.

Referenced by Node::resendMolecule().

{
  Real *a1, *a2, *a3, *a4;        //  Temporary arrays for sending messages
  int *i1, *i2, *i3;      //  Temporary int array
  int i, j;      //  Loop counters
  Real **kvals;      //  Force constant values for dihedrals and impropers
  int **nvals;      //  Periodicity values for  dihedrals and impropers
  Real **deltavals;    //  Phase shift values for  dihedrals and impropers
  BigReal *pairC6, *pairC12; // JLai
  /*MOStream *msg=comm_obj->newOutputStream(ALLBUTME, STATICPARAMSTAG, BUFSIZE);
  if ( msg == NULL )
  {
    NAMD_die("memory allocation failed in Parameters::send_Parameters");
  }*/

  //  Send the bond parameters
  msg->put(NumBondParams);

  if (NumBondParams)
  {
    a1 = new Real[NumBondParams];
    a2 = new Real[NumBondParams];

    if ( (a1 == NULL) || (a2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    for (i=0; i<NumBondParams; i++)
    {
      a1[i] = bond_array[i].k;
      a2[i] = bond_array[i].x0;
    }

    msg->put(NumBondParams, a1)->put(NumBondParams, a2);

    delete [] a1;
    delete [] a2;
  }

  //  Send the angle parameters
  msg->put(NumAngleParams);

  if (NumAngleParams)
  {
    a1 = new Real[NumAngleParams];
    a2 = new Real[NumAngleParams];
    a3 = new Real[NumAngleParams];
    a4 = new Real[NumAngleParams];
    i1 = new int[NumAngleParams];

    if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) ||
         (a4 == NULL) || (i1 == NULL))
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    for (i=0; i<NumAngleParams; i++)
    {
      a1[i] = angle_array[i].k;
      a2[i] = angle_array[i].theta0;
      a3[i] = angle_array[i].k_ub;
      a4[i] = angle_array[i].r_ub;
      i1[i] = angle_array[i].normal;
    }

    msg->put(NumAngleParams, a1)->put(NumAngleParams, a2);
    msg->put(NumAngleParams, a3)->put(NumAngleParams, a4);
    msg->put(NumAngleParams, i1);

    delete [] a1;
    delete [] a2;
    delete [] a3;
    delete [] a4;
    delete [] i1;
  }

  //  Send the dihedral parameters
  msg->put(NumDihedralParams);

  if (NumDihedralParams)
  {
    i1 = new int[NumDihedralParams];
    kvals = new Real *[MAX_MULTIPLICITY];
    nvals = new int *[MAX_MULTIPLICITY];
    deltavals = new Real *[MAX_MULTIPLICITY];

    if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
         (deltavals == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      kvals[i] = new Real[NumDihedralParams];
      nvals[i] = new int[NumDihedralParams];
      deltavals[i] = new Real[NumDihedralParams];

      if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
      {
        NAMD_die("memory allocation failed in Parameters::send_Parameters");
      }
    }

    for (i=0; i<NumDihedralParams; i++)
    {
      i1[i] = dihedral_array[i].multiplicity;

      for (j=0; j<MAX_MULTIPLICITY; j++)
      {
        kvals[j][i] = dihedral_array[i].values[j].k;
        nvals[j][i] = dihedral_array[i].values[j].n;
        deltavals[j][i] = dihedral_array[i].values[j].delta;
      }
    }

    msg->put(NumDihedralParams, i1);

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      msg->put(NumDihedralParams, kvals[i]);
      msg->put(NumDihedralParams, nvals[i]);
      msg->put(NumDihedralParams, deltavals[i]);

      delete [] kvals[i];
      delete [] nvals[i];
      delete [] deltavals[i];
    }

    delete [] i1;
    delete [] kvals;
    delete [] nvals;
    delete [] deltavals;
  }

  //  Send the improper parameters
  msg->put(NumImproperParams);

  if (NumImproperParams)
  {
    i1 = new int[NumImproperParams];
    kvals = new Real *[MAX_MULTIPLICITY];
    nvals = new int *[MAX_MULTIPLICITY];
    deltavals = new Real *[MAX_MULTIPLICITY];

    if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) || 
         (deltavals == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      kvals[i] = new Real[NumImproperParams];
      nvals[i] = new int[NumImproperParams];
      deltavals[i] = new Real[NumImproperParams];

      if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
      {
        NAMD_die("memory allocation failed in Parameters::send_Parameters");
      }
    }

    for (i=0; i<NumImproperParams; i++)
    {
      i1[i] = improper_array[i].multiplicity;

      for (j=0; j<MAX_MULTIPLICITY; j++)
      {
        kvals[j][i] = improper_array[i].values[j].k;
        nvals[j][i] = improper_array[i].values[j].n;
        deltavals[j][i] = improper_array[i].values[j].delta;
      }
    }

    msg->put(NumImproperParams, i1);

    for (i=0; i<MAX_MULTIPLICITY; i++)
    {
      msg->put(NumImproperParams, kvals[i]);
      msg->put(NumImproperParams, nvals[i]);
      msg->put(NumImproperParams, deltavals[i]);

      delete [] kvals[i];
      delete [] nvals[i];
      delete [] deltavals[i];
    }

    delete [] i1;
    delete [] kvals;
    delete [] nvals;
    delete [] deltavals;
  }

  //  Send the crossterm parameters
  msg->put(NumCrosstermParams);

  if (NumCrosstermParams)
  {
    for (i=0; i<NumCrosstermParams; ++i) {
      int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
      msg->put(nvals,&crossterm_array[i].c[0][0].d00);
    }
  }
  //  Send the GromacsPairs parameters
  // JLai
  msg->put(NumGromacsPairParams);
 
  if(NumGromacsPairParams) 
  {
      pairC6 = new BigReal[NumGromacsPairParams];
      pairC12 = new BigReal[NumGromacsPairParams];
      if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
          NAMD_die("memory allocation failed in Parameters::send_Parameters");
      }

      for (i=0; i<NumGromacsPairParams; i++) {
          pairC6[i] = gromacsPair_array[i].pairC6;
          pairC12[i] = gromacsPair_array[i].pairC12;
      }

      msg->put(NumGromacsPairParams,pairC6);
      msg->put(NumGromacsPairParams,pairC12);

      delete [] pairC6;
      delete [] pairC12;
  }
  // End of JLai
  
  //
  //Send the energy table parameters
  msg->put(numenerentries);

  if (numenerentries) {
          /*
    b1 = new Real[numenerentries];
    if (b1 == NULL) 
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }
    */
    
    msg->put(numenerentries, table_ener);
  }

  //  Send the vdw parameters
  msg->put(NumVdwParams);
  msg->put(NumVdwParamsAssigned);

  if (NumVdwParams)
  {
    a1 = new Real[NumVdwParams];
    a2 = new Real[NumVdwParams];
    a3 = new Real[NumVdwParams];
    a4 = new Real[NumVdwParams];

    if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    for (i=0; i<NumVdwParams; i++)
    {
      a1[i] = vdw_array[i].sigma;
      a2[i] = vdw_array[i].epsilon;
      a3[i] = vdw_array[i].sigma14;
      a4[i] = vdw_array[i].epsilon14;
    }

    msg->put(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
    msg->put(NumVdwParams, a1);
    msg->put(NumVdwParams, a2);
    msg->put(NumVdwParams, a3);
    msg->put(NumVdwParams, a4);

    delete [] a1;
    delete [] a2;
    delete [] a3;
    delete [] a4;
  }
  
  //  Send the vdw pair parameters
  msg->put(NumVdwPairParams);
  
  if (NumVdwPairParams)
  {
    a1 = new Real[NumVdwPairParams];
    a2 = new Real[NumVdwPairParams];
    a3 = new Real[NumVdwPairParams];
    a4 = new Real[NumVdwPairParams];
    i1 = new int[NumVdwPairParams];
    i2 = new int[NumVdwPairParams];    

    if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) || 
         (i1 == NULL) || (i2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }
    
    vdw_pair_to_arrays(i1, i2, a1, a2, a3, a4, 0, vdw_pair_tree);
    
    msg->put(NumVdwPairParams, i1)->put(NumVdwPairParams, i2);
    msg->put(NumVdwPairParams, a1);
    msg->put(NumVdwPairParams, a2)->put(NumVdwPairParams, a3);
    msg->put(NumVdwPairParams, a4);
  }

  //  Send the nbthole pair parameters
  msg->put(NumNbtholePairParams);

  if (NumNbtholePairParams)
  {
    a1 = new Real[NumNbtholePairParams];
    a2 = new Real[NumNbtholePairParams];
    a3 = new Real[NumNbtholePairParams];
    i1 = new int[NumNbtholePairParams];
    i2 = new int[NumNbtholePairParams];

    if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (i1 == NULL) || (i2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }

    nbthole_pair_to_arrays(i1, i2, a1, a2, a3, 0, nbthole_pair_tree);

   for (i=0; i<NumNbtholePairParams; i++)
   {
    nbthole_array[i].ind1 = i1[i];
    nbthole_array[i].ind2 = i2[i];
    nbthole_array[i].alphai = a1[i];
    nbthole_array[i].alphaj = a2[i];
    nbthole_array[i].tholeij = a3[i];
   }

    msg->put(NumNbtholePairParams, i1)->put(NumNbtholePairParams, i2);
    msg->put(NumNbtholePairParams, a1);
    msg->put(NumNbtholePairParams, a2)->put(NumNbtholePairParams, a3);
  }
  
  //  Send the table pair parameters
  //printf("Pairs: %i\n", NumTablePairParams);
  msg->put(NumTablePairParams);
  
  if (NumTablePairParams)
  {
    i1 = new int[NumTablePairParams];
    i2 = new int[NumTablePairParams];    
    i3 = new int[NumTablePairParams];

    if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
    {
      NAMD_die("memory allocation failed in Parameters::send_Parameters");
    }
    
    table_pair_to_arrays(i1, i2, i3, 0, tab_pair_tree);
    
    msg->put(NumTablePairParams, i1)->put(NumTablePairParams, i2);
    msg->put(NumTablePairParams, i3);
  }

  // send the hydrogen bond parameters
  // hbondParams.create_message(msg);
  msg->end();
  delete msg;
}

Member Data Documentation

AngleValue* Parameters::angle_array

Definition at line 244 of file Parameters.h.

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

BondValue* Parameters::bond_array

Definition at line 243 of file Parameters.h.

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

int Parameters::columnsize

Definition at line 255 of file Parameters.h.

Referenced by ComputeNonbondedUtil::select().

CrosstermValue* Parameters::crossterm_array

Definition at line 247 of file Parameters.h.

Referenced by done_reading_files(), CrosstermElem::getParameterPointers(), read_parm(), receive_Parameters(), send_Parameters(), and ~Parameters().

DihedralValue* Parameters::dihedral_array

Definition at line 245 of file Parameters.h.

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

GromacsPairValue* Parameters::gromacsPair_array

Definition at line 249 of file Parameters.h.

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

ImproperValue* Parameters::improper_array

Definition at line 246 of file Parameters.h.

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

NbtholePairValue* Parameters::nbthole_array

Definition at line 252 of file Parameters.h.

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

IndexedNbtholePair* Parameters::nbthole_pair_tree

Definition at line 258 of file Parameters.h.

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

int Parameters::NumAngleParams

Definition at line 262 of file Parameters.h.

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

int Parameters::NumBondParams

Definition at line 261 of file Parameters.h.

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

int Parameters::NumCosAngles

Definition at line 263 of file Parameters.h.

Referenced by print_param_summary().

int Parameters::NumCrosstermParams

Definition at line 266 of file Parameters.h.

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

int Parameters::NumDihedralParams

Definition at line 264 of file Parameters.h.

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

int Parameters::numenerentries

Definition at line 253 of file Parameters.h.

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

int Parameters::NumGromacsPairParams

Definition at line 268 of file Parameters.h.

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

int Parameters::NumImproperParams

Definition at line 265 of file Parameters.h.

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

int Parameters::NumNbtholePairParams

Definition at line 274 of file Parameters.h.

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

int Parameters::NumTablePairParams

Definition at line 275 of file Parameters.h.

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

int Parameters::NumTableParams

Definition at line 271 of file Parameters.h.

int Parameters::NumVdwPairParams

Definition at line 273 of file Parameters.h.

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

int Parameters::NumVdwParams

Definition at line 270 of file Parameters.h.

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

int Parameters::NumVdwParamsAssigned

Definition at line 272 of file Parameters.h.

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

int Parameters::rowsize

Definition at line 254 of file Parameters.h.

Referenced by ComputeNonbondedUtil::select().

IndexedTablePair* Parameters::tab_pair_tree

Definition at line 259 of file Parameters.h.

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

BigReal* Parameters::table_ener

Definition at line 256 of file Parameters.h.

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

int Parameters::tablenumtypes

Definition at line 260 of file Parameters.h.

Referenced by get_int_table_type(), and read_ener_table().

VdwValue* Parameters::vdw_array

Definition at line 251 of file Parameters.h.

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

IndexedVdwPair* Parameters::vdw_pair_tree

Definition at line 257 of file Parameters.h.

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

---

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

• Parameters.h
• Parameters.C