NAMD
Parameters.C
Go to the documentation of this file.
1 
7 /*
8  The class Parameters is used to hold all of the parameters read
9  in from the parameter files. The class provides a routine to read in
10  parameter files (as many parameter files as desired can be read in) and
11  a series of routines that allow the parameters that have been read in
12  to be queried.
13 */
14 
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <string.h>
18 #include <vector>
19 #ifndef WIN32
20 #include <strings.h>
21 #endif
22 #include "InfoStream.h"
23 #include <charm++.h>
24 #include "Parameters.h"
25 #include "Communicate.h"
26 #include "ConfigList.h"
27 //****** BEGIN CHARMM/XPLOR type changes
28 #include "SimParameters.h"
29 //****** END CHARMM/XPLOR type changes
30 
31 #define MIN_DEBUG_LEVEL 3
32 //#define DEBUGM
33 #include "Debug.h"
34 
35 #define INDEX(ncols,i,j) ((i)*ncols + (j))
36 
37 #define ENABLETABLES
38 
39 static char** table_types;
40 
41 // struct bond_params is used to form a binary tree of bond parameters.
42 // The two atom names are used to determine the order of the nodes in the
43 // tree. atom1name should ALWAYS be lexically before atom2name
44 
46 {
47  char atom1name[11];
48  char atom2name[11];
52  struct bond_params *left;
53  struct bond_params *right;
54 };
55 
56 // struct angle_params is used to form a binary tree of bond parameters.
57 // The three atom names are used to determine the order of the nodes in
58 // the tree. atom1name should ALWAYS be lexically before atom3name
59 
61 {
62  char atom1name[11];
63  char atom2name[11];
64  char atom3name[11];
66  int normal;
71  struct angle_params *left;
73 };
74 
75 // struct dihedral_params is used to form a linked list of the dihedral
76 // parameters. The linked list is arranged in such a way that any
77 // bonds with wildcards are at the end of the list so that a linear
78 // search can be done but we will still find exact matches before
79 // wildcard matches
80 
82 {
83  char atom1name[11];
84  char atom2name[11];
85  char atom3name[11];
86  char atom4name[11];
87  char atom1wild;
88  char atom2wild;
89  char atom3wild;
90  char atom4wild;
95  dihedral_params() { memset(this, 0, sizeof(dihedral_params)); }
96 };
97 
98 // struct improper_params is used to form a linked list of the improper
99 // parameters. The linked list is arranged in such a way that any
100 // bonds with wildcards are at the end of the list so that a linear
101 // search can be done but we will still find exact matches before
102 // wildcard matches
103 
105 {
106  char atom1name[11];
107  char atom2name[11];
108  char atom3name[11];
109  char atom4name[11];
114 };
115 
117 {
118  crossterm_params(int dim) : dimension(dim) {
119  values = new double[dimension*dimension];
120  }
122  delete [] values;
123  }
124  char atom1name[11];
125  char atom2name[11];
126  char atom3name[11];
127  char atom4name[11];
128  char atom5name[11];
129  char atom6name[11];
130  char atom7name[11];
131  char atom8name[11];
132  int dimension; // usually 24
133  double *values; // dimension * dimension data
136 };
137 
138 // struct vdw_params is used to form a binary serach tree of the
139 // vdw paramters for a single atom.
140 
142 {
143  char atomname[11];
149  struct vdw_params *left;
150  struct vdw_params *right;
151 };
152 
153 // struct vdw_pair_params is used to form a linked list of the
154 // vdw parameters for a pair of atoms
155 
157 {
158  char atom1name[11];
159  char atom2name[11];
165 };
166 
168 {
169  char atom1name[11];
170  char atom2name[11];
171  int K;
173 };
174 
176 {
177  char atom1name[11];
178  char atom2name[11];
184 };
185 
187  initialize();
188 }
189 
190 void Parameters::initialize() {
191 
192  paramType = -1;
193 
194  /* Set all the pointers to NULL */
195  atomTypeNames=NULL;
196  bondp=NULL;
197  anglep=NULL;
198  improperp=NULL;
199  dihedralp=NULL;
200  crosstermp=NULL;
201  vdwp=NULL;
202  vdw_pairp=NULL;
203  nbthole_pairp=NULL;
204  table_pairp=NULL;
205  bond_array=NULL;
206  angle_array=NULL;
207  dihedral_array=NULL;
208  improper_array=NULL;
209  crossterm_array=NULL;
210  // JLai
211  gromacsPair_array=NULL;
212  // End of JLai
213  vdw_array=NULL;
214  vdw_pair_tree=NULL;
215  nbthole_pair_tree=NULL;
216  tab_pair_tree=NULL;
217  maxDihedralMults=NULL;
218  maxImproperMults=NULL;
219  table_ener = NULL;
220 
221  /* Set all the counts to 0 */
222  NumBondParams=0;
223  NumAngleParams=0;
227  // JLai
229  // End of JLai
230  NumVdwParams=0;
234  NumCosAngles=0;
235  numenerentries=0;
236 }
237 
238 /************************************************************************/
239 /* */
240 /* FUNCTION Parameters */
241 /* */
242 /* This is the constructor for the class. It simlpy sets the */
243 /* pointers to the list and trees to NULL and the count of all the */
244 /* parameters to 0. */
245 /* The type (format) of the input parameters (Xplor,Charmm) is set here. */
246 /* */
247 /************************************************************************/
248 
250 {
251  initialize();
252 
254  if (simParams->paraTypeXplorOn)
255  {
256  paramType = paraXplor;
257  }
258  else if (simParams->paraTypeCharmmOn)
259  {
260  paramType = paraCharmm;
261  }
262  //****** END CHARMM/XPLOR type changes
263  //Test for cos-based angles
264  if (simParams->cosAngles) {
265  cosAngles = true;
266  } else {
267  cosAngles = false;
268  }
269 
270  if (simParams->tabulatedEnergies) {
271  CkPrintf("Working on tables\n");
272  read_ener_table(simParams);
273  }
274 
275  //****** BEGIN CHARMM/XPLOR type changes
276  /* Set up AllFilesRead flag to FALSE. Once all of the files */
277  /* have been read in, then this will be set to true and the */
278  /* arrays of parameters will be set up */
279  AllFilesRead = FALSE;
280 
281  if (NULL != f)
282  {
283  do
284  {
285  //****** BEGIN CHARMM/XPLOR type changes
286  if (paramType == paraXplor)
287  {
289  }
290  else if (paramType == paraCharmm)
291  {
293  }
294  //****** END CHARMM/XPLOR type changes
295  f = f->next;
296  } while ( f != NULL );
297 
298  done_reading_files(simParams->drudeOn && paramType == paraCharmm);
299  }
300 
301 }
302 /* END OF FUNCTION Parameters */
303 
304 /************************************************************************/
305 /* */
306 /* FUNCTION ~Parameters */
307 /* */
308 /* This is the destructor for this class. It basically just */
309 /* frees all of the memory allocated for the parameters. */
310 /* */
311 /************************************************************************/
312 
314 
315 {
316  if (atomTypeNames)
317  delete [] atomTypeNames;
318 
319  if (bondp != NULL)
320  free_bond_tree(bondp);
321 
322  if (anglep != NULL)
323  free_angle_tree(anglep);
324 
325  if (dihedralp != NULL)
326  free_dihedral_list(dihedralp);
327 
328  if (improperp != NULL)
329  free_improper_list(improperp);
330 
331  if (crosstermp != NULL)
332  free_crossterm_list(crosstermp);
333 
334  if (vdwp != NULL)
335  free_vdw_tree(vdwp);
336 
337  if (vdw_pairp != NULL)
338  free_vdw_pair_list();
339 
340  if (nbthole_pairp != NULL)
341  free_nbthole_pair_list();
342 
343  if (bond_array != NULL)
344  delete [] bond_array;
345 
346  if (angle_array != NULL)
347  delete [] angle_array;
348 
349  if (dihedral_array != NULL)
350  delete [] dihedral_array;
351 
352  if (improper_array != NULL)
353  delete [] improper_array;
354 
355  if (crossterm_array != NULL)
356  delete [] crossterm_array;
357 
358  // JLai
359  if (gromacsPair_array != NULL)
360  delete [] gromacsPair_array;
361  // End of JLai
362 
363  if (vdw_array != NULL)
364  delete [] vdw_array;
365 
366  if (tab_pair_tree != NULL)
367  free_table_pair_tree(tab_pair_tree);
368 
369  if (vdw_pair_tree != NULL)
370  free_vdw_pair_tree(vdw_pair_tree);
371 
372  if (nbthole_pair_tree != NULL)
373  free_nbthole_pair_tree(nbthole_pair_tree);
374 
375  if (maxDihedralMults != NULL)
376  delete [] maxDihedralMults;
377 
378  if (maxImproperMults != NULL)
379  delete [] maxImproperMults;
380 
381  for( int i = 0; i < error_msgs.size(); ++i ) {
382  delete [] error_msgs[i];
383  }
384  error_msgs.resize(0);
385 }
386 /* END OF FUNCTION ~Parameters */
387 
388 /************************************************************************/
389 /* */
390 /* FUNCTION read_paramter_file */
391 /* */
392 /* INPUTS: */
393 /* fname - name of the parameter file to read */
394 /* */
395 /* This function reads in a parameter file and adds the parameters */
396 /* from this file to the current group of parameters. The basic */
397 /* logic of the routine is to read in a line from the file, looks at */
398 /* the first word of the line to determine what kind of parameter we */
399 /* have, and then call the appropriate routine to add the parameter */
400 /* to the parameters that we have. */
401 /* */
402 /************************************************************************/
403 
405 
406 {
407  char buffer[512]; // Buffer to store each line of the file
408  char first_word[512]; // First word of the current line
409  FILE *pfile; // File descriptor for the parameter file
410 
411  /* Check to make sure that we haven't previously been told */
412  /* that all the files were read */
413  if (AllFilesRead)
414  {
415  NAMD_die("Tried to read another parameter file after being told that all files were read!");
416  }
417 
418  /* Try and open the file */
419  if ( (pfile = Fopen(fname, "r")) == NULL)
420  {
421  char err_msg[256];
422 
423  sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
424  NAMD_die(err_msg);
425  }
426 
427  /* Keep reading in lines until we hit the EOF */
428  while (NAMD_read_line(pfile, buffer) != -1)
429  {
430  /* Get the first word of the line */
431  NAMD_find_first_word(buffer, first_word);
432 
433  /* First, screen out things that we ignore, such as */
434  /* blank lines, lines that start with '!', lines that */
435  /* start with "REMARK", lines that start with set", */
436  /* and most of the HBOND parameters which include */
437  /* AEXP, REXP, HAEX, AAEX, but not the HBOND statement */
438  /* which is parsed. */
439  if ((buffer[0] != '!') &&
440  !NAMD_blank_string(buffer) &&
441  (strncasecmp(first_word, "REMARK", 6) != 0) &&
442  (strcasecmp(first_word, "set")!=0) &&
443  (strncasecmp(first_word, "AEXP", 4) != 0) &&
444  (strncasecmp(first_word, "REXP", 4) != 0) &&
445  (strncasecmp(first_word, "HAEX", 4) != 0) &&
446  (strncasecmp(first_word, "AAEX", 4) != 0) &&
447  (strncasecmp(first_word, "NBOND", 5) != 0) &&
448  (strncasecmp(first_word, "CUTNB", 5) != 0) &&
449  (strncasecmp(first_word, "END", 3) != 0) &&
450  (strncasecmp(first_word, "CTONN", 5) != 0) &&
451  (strncasecmp(first_word, "EPS", 3) != 0) &&
452  (strncasecmp(first_word, "VSWI", 4) != 0) &&
453  (strncasecmp(first_word, "NBXM", 4) != 0) &&
454  (strncasecmp(first_word, "INHI", 4) != 0) )
455  {
456  /* Now, call the appropriate function based */
457  /* on the type of parameter we have */
458  if (strncasecmp(first_word, "bond", 4)==0)
459  {
460  add_bond_param(buffer, TRUE);
461  NumBondParams++;
462  }
463  else if (strncasecmp(first_word, "angl", 4)==0)
464  {
465  add_angle_param(buffer);
466  NumAngleParams++;
467  }
468  else if (strncasecmp(first_word, "dihe", 4)==0)
469  {
470  add_dihedral_param(buffer, pfile);
472  }
473  else if (strncasecmp(first_word, "impr", 4)==0)
474  {
475  add_improper_param(buffer, pfile);
477  }
478  else if (strncasecmp(first_word, "nonb", 4)==0)
479  {
480  add_vdw_param(buffer);
481  NumVdwParams++;
482  }
483  else if (strncasecmp(first_word, "nbfi", 4)==0)
484  {
485  add_vdw_pair_param(buffer);
486  NumVdwPairParams++;
487  }
488  else if (strncasecmp(first_word, "nbta", 4)==0)
489  {
490 
491  if (table_ener == NULL) {
492  continue;
493  }
494 
495  add_table_pair_param(buffer);
497  }
498  else if (strncasecmp(first_word, "hbon", 4)==0)
499  {
500  add_hb_pair_param(buffer);
501  }
502  else
503  {
504  /* This is an unknown paramter. */
505  /* This is BAD */
506  char err_msg[512];
507 
508  sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
509  fname, buffer);
510  NAMD_die(err_msg);
511  }
512  }
513  }
514 
515  /* Close the file */
516  Fclose(pfile);
517 
518  return;
519 }
520 /* END OF FUNCTION read_paramter_file */
521 
522 //****** BEGIN CHARMM/XPLOR type changes
523 /************************************************************************/
524 /* */
525 /* FUNCTION read_charmm_paramter_file */
526 /* */
527 /* INPUTS: */
528 /* fname - name of the parameter file to read */
529 /* */
530 /* This function reads in a CAHRMM parameter file and adds the */
531 /* parameters from this file to the current group of parameters. */
532 /* The basic logic of the routine is to first find out what type of */
533 /* parameter we have in the file. Then look at each line in turn */
534 /* and call the appropriate routine to add the parameters until we hit*/
535 /* a new type of parameter or EOF. */
536 /* */
537 /************************************************************************/
538 
540 
541 {
542  int par_type=0; // What type of parameter are we currently
543  // dealing with? (vide infra)
544  int skipline; // skip this line?
545  int skipall = 0; // skip rest of file;
546  char buffer[512]; // Buffer to store each line of the file
547  char first_word[512]; // First word of the current line
548  FILE *pfile; // File descriptor for the parameter file
549 
550  /* Check to make sure that we haven't previously been told */
551  /* that all the files were read */
552  if (AllFilesRead)
553  {
554  NAMD_die("Tried to read another parameter file after being told that all files were read!");
555  }
556 
557  /* Try and open the file */
558  if ( (pfile = fopen(fname, "r")) == NULL)
559  {
560  char err_msg[256];
561 
562  sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
563  NAMD_die(err_msg);
564  }
565 
566  /* Keep reading in lines until we hit the EOF */
567  while (NAMD_read_line(pfile, buffer) != -1)
568  {
569  /* Get the first word of the line */
570  NAMD_find_first_word(buffer, first_word);
571  skipline=0;
572 
573  /* First, screen out things that we ignore. */
574  /* blank lines, lines that start with '!' or '*', lines that */
575  /* start with "END". */
576  if (!NAMD_blank_string(buffer) &&
577  (strncmp(first_word, "!", 1) != 0) &&
578  (strncmp(first_word, "*", 1) != 0) &&
579  (strncasecmp(first_word, "END", 3) != 0))
580  {
581  if ( skipall ) {
582  iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
583  break;
584  }
585  /* Now, determine the apropriate parameter type. */
586  if (strncasecmp(first_word, "bond", 4)==0)
587  {
588  par_type=1; skipline=1;
589  }
590  else if (strncasecmp(first_word, "angl", 4)==0)
591  {
592  par_type=2; skipline=1;
593  }
594  else if (strncasecmp(first_word, "thet", 4)==0)
595  {
596  par_type=2; skipline=1;
597  }
598  else if (strncasecmp(first_word, "dihe", 4)==0)
599  {
600  par_type=3; skipline=1;
601  }
602  else if (strncasecmp(first_word, "phi", 3)==0)
603  {
604  par_type=3; skipline=1;
605  }
606  else if (strncasecmp(first_word, "impr", 4)==0)
607  {
608  par_type=4; skipline=1;
609  }
610  else if (strncasecmp(first_word, "imph", 4)==0)
611  {
612  par_type=4; skipline=1;
613  }
614  else if (strncasecmp(first_word, "nbon", 4)==0)
615  {
616  par_type=5; skipline=1;
617  }
618  else if (strncasecmp(first_word, "nonb", 4)==0)
619  {
620  par_type=5; skipline=1;
621  }
622  else if (strncasecmp(first_word, "nbfi", 4)==0)
623  {
624  par_type=6; skipline=1;
625  }
626  else if (strncasecmp(first_word, "hbon", 4)==0)
627  {
628  par_type=7; skipline=1;
629  }
630  else if (strncasecmp(first_word, "cmap", 4)==0)
631  {
632  par_type=8; skipline=1;
633  }
634  else if (strncasecmp(first_word, "nbta", 4)==0)
635  {
636  par_type=9; skipline=1;
637  }
638  else if (strncasecmp(first_word, "thol", 4)==0)
639  {
640  par_type=10; skipline=1;
641  }
642  else if (strncasecmp(first_word, "atom", 4)==0)
643  {
644  par_type=11; skipline=1;
645  }
646  else if (strncasecmp(first_word, "ioformat", 8)==0)
647  {
648  skipline=1;
649  }
650  else if (strncasecmp(first_word, "read", 4)==0)
651  {
652  skip_stream_read(buffer, pfile); skipline=1;
653  }
654  else if (strncasecmp(first_word, "return", 4)==0)
655  {
656  skipall=1; skipline=1;
657  }
658  else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
659  (strncasecmp(first_word, "grou", 4) == 0) ||
660  (strncasecmp(first_word, "cdie", 4) == 0) ||
661  (strncasecmp(first_word, "shif", 4) == 0) ||
662  (strncasecmp(first_word, "vgro", 4) == 0) ||
663  (strncasecmp(first_word, "vdis", 4) == 0) ||
664  (strncasecmp(first_word, "vswi", 4) == 0) ||
665  (strncasecmp(first_word, "cutn", 4) == 0) ||
666  (strncasecmp(first_word, "ctof", 4) == 0) ||
667  (strncasecmp(first_word, "cton", 4) == 0) ||
668  (strncasecmp(first_word, "eps", 3) == 0) ||
669  (strncasecmp(first_word, "e14f", 4) == 0) ||
670  (strncasecmp(first_word, "wmin", 4) == 0) ||
671  (strncasecmp(first_word, "aexp", 4) == 0) ||
672  (strncasecmp(first_word, "rexp", 4) == 0) ||
673  (strncasecmp(first_word, "haex", 4) == 0) ||
674  (strncasecmp(first_word, "aaex", 4) == 0) ||
675  (strncasecmp(first_word, "noac", 4) == 0) ||
676  (strncasecmp(first_word, "hbno", 4) == 0) ||
677  (strncasecmp(first_word, "cuth", 4) == 0) ||
678  (strncasecmp(first_word, "ctof", 4) == 0) ||
679  (strncasecmp(first_word, "cton", 4) == 0) ||
680  (strncasecmp(first_word, "cuth", 4) == 0) ||
681  (strncasecmp(first_word, "ctof", 4) == 0) ||
682  (strncasecmp(first_word, "cton", 4) == 0) )
683  {
684  if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
685  {
686  char err_msg[512];
687 
688  sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
689  NAMD_die(err_msg);
690  }
691  else
692  {
693  skipline = 1;
694  }
695  }
696  else if (par_type == 0)
697  {
698  /* This is an unknown paramter. */
699  /* This is BAD */
700  char err_msg[512];
701 
702  sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
703  NAMD_die(err_msg);
704  }
705  }
706  else
707  {
708  skipline=1;
709  }
710 
711  if ( (par_type != 0) && (!skipline) )
712  {
713  /* Now, call the appropriate function based */
714  /* on the type of parameter we have */
715  /* I know, this should really be a switch ... */
716  if (par_type == 1)
717  {
718  add_bond_param(buffer, TRUE);
719  NumBondParams++;
720  }
721  else if (par_type == 2)
722  {
723  add_angle_param(buffer);
724  NumAngleParams++;
725  }
726  else if (par_type == 3)
727  {
728  add_dihedral_param(buffer, pfile);
730  }
731  else if (par_type == 4)
732  {
733  add_improper_param(buffer, pfile);
735  }
736  else if (par_type == 5)
737  {
738  add_vdw_param(buffer);
739  NumVdwParams++;
740  }
741  else if (par_type == 6)
742  {
743  add_vdw_pair_param(buffer);
744  NumVdwPairParams++;
745  }
746  else if (par_type == 7)
747  {
748  add_hb_pair_param(buffer);
749  }
750  else if (par_type == 8)
751  {
752  add_crossterm_param(buffer, pfile);
754  }
755  else if (par_type == 9)
756  {
757 
758  if (table_ener == NULL) {
759  continue;
760  }
761 
762  add_table_pair_param(buffer);
764  }
765  else if (par_type == 10)
766  {
767  add_nbthole_pair_param(buffer);
769  }
770  else if (par_type == 11)
771  {
772  if ( strncasecmp(first_word, "mass", 4) != 0 ) {
773  char err_msg[512];
774  sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
775  NAMD_die(err_msg);
776  }
777  }
778  else
779  {
780  /* This really should not occour! */
781  /* This is an internal error. */
782  /* This is VERY BAD */
783  char err_msg[512];
784 
785  sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
786  NAMD_die(err_msg);
787  }
788  }
789  }
790 
791  /* Close the file */
792  fclose(pfile);
793 
794  return;
795 }
796 /* END OF FUNCTION read_charmm_paramter_file */
797 //****** END CHARMM/XPLOR type changes
798 
799 
800 void Parameters::skip_stream_read(char *buf, FILE *fd) {
801 
802  char buffer[513];
803  char first_word[513];
804  char s1[128];
805  char s2[128];
806  int rval = sscanf(buf, "%s %s", s1, s2);
807  if (rval != 2) {
808  char err_msg[512];
809  sprintf(err_msg, "BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
810  NAMD_die(err_msg);
811  }
812  if ( ! strncasecmp(s2,"PARA",4) ) return; // read parameters
813 
814  iout << iINFO << "SKIPPING " << s2 << " SECTION IN STREAM FILE\n" << endi;
815 
816  while (NAMD_read_line(fd, buffer) != -1)
817  {
818  // read until we find "END"
819  NAMD_find_first_word(buffer, first_word);
820  if (!NAMD_blank_string(buffer) &&
821  (strncmp(first_word, "!", 1) != 0) &&
822  (strncmp(first_word, "*", 1) != 0) &&
823  (strncasecmp(first_word, "END", 3) == 0) ) return;
824  }
825 
826 }
827 
828 
829 /************************************************************************/
830 /* */
831 /* FUNCTION add_bond_param */
832 /* */
833 /* INPUTS: */
834 /* buf - Line from parameter file containing bond parameters */
835 /* */
836 /* This function adds a new bond paramter to the binary tree of */
837 /* angle paramters that we have. If a duplicate is found, a warning */
838 /* message is printed and the new parameters are used. */
839 /* */
840 /************************************************************************/
841 
842 void Parameters::add_bond_param(const char *buf, Bool overwrite)
843 
844 {
845  char atom1name[11]; // Atom type for atom 1
846  char atom2name[11]; // Atom type for atom 2
847  Real forceconstant; // Force constant for bond
848  Real distance; // Rest distance for bond
849  int read_count; // Count from sscanf
850  struct bond_params *new_node; // New node in tree
851 
852  //****** BEGIN CHARMM/XPLOR type changes
853  /* Use sscanf to parse up the input line */
854  if (paramType == paraXplor)
855  {
856  /* read XPLOR format */
857  read_count=sscanf(buf, "%*s %s %s %f %f\n", atom1name, atom2name,
858  &forceconstant, &distance);
859  }
860  else if (paramType == paraCharmm)
861  {
862  /* read CHARMM format */
863  read_count=sscanf(buf, "%s %s %f %f\n", atom1name, atom2name,
864  &forceconstant, &distance);
865  }
866  //****** END CHARMM/XPLOR type changes
867 
868  /* Check to make sure we found everything we expeceted */
869  if (read_count != 4)
870  {
871  char err_msg[512];
872 
873  if (paramType == paraXplor)
874  {
875  sprintf(err_msg, "BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
876  }
877  else if (paramType == paraCharmm)
878  {
879  sprintf(err_msg, "BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
880  }
881  NAMD_die(err_msg);
882  }
883 
884  /* Allocate a new node */
885  new_node = new bond_params;
886 
887  if (new_node == NULL)
888  {
889  NAMD_die("memory allocation failed in Parameters::add_bond_param\n");
890  }
891 
892  /* Order the atoms so that the atom that comes alphabetically */
893  /* first is atom 1. Since the bond is symmetric, it doesn't */
894  /* matter physically which atom is first. And this allows the */
895  /* search of the binary tree to be done in a logical manner */
896  if (strcasecmp(atom1name, atom2name) < 0)
897  {
898  strcpy(new_node->atom1name, atom1name);
899  strcpy(new_node->atom2name, atom2name);
900  }
901  else
902  {
903  strcpy(new_node->atom2name, atom1name);
904  strcpy(new_node->atom1name, atom2name);
905  }
906 
907  /* Assign force constant and distance */
908  new_node->forceconstant = forceconstant;
909  new_node->distance = distance;
910 
911  /* Set pointers to null */
912  new_node->left = NULL;
913  new_node->right = NULL;
914 
915  /* Make call to recursive call to actually add the node to the */
916  /* tree */
917  bondp=add_to_bond_tree(new_node, bondp, overwrite);
918 
919  return;
920 }
921 /* END OF FUNCTION add_bond_param */
922 
923 /************************************************************************/
924 /* */
925 /* FUNCTION add_to_bond_tree */
926 /* */
927 /* INPUTS: */
928 /* new_node - Node to add to the tree */
929 /* tree - tree to add the node to */
930 /* */
931 /* OUTPUTS: */
932 /* ths function returns a pointer to the new tree with the node */
933 /* added to it. Most of the time it will be the same pointer as was */
934 /* passed in, but not if the current tree is empty. */
935 /* */
936 /* this is a receursive function that adds a node to the binary */
937 /* tree used to store bond parameters. */
938 /* */
939 /************************************************************************/
940 
941 struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
942  struct bond_params *tree, Bool overwrite)
943 
944 {
945  int compare_code; // Results from strcasecmp
946 
947  /* If the tree is currently empty, then the new tree consists */
948  /* only of the new node */
949  if (tree == NULL)
950  return(new_node);
951 
952  /* Compare the atom1 name from the new node and the head of */
953  /* the tree */
954  compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
955 
956  /* Check to see if they are the same */
957  if (compare_code == 0)
958  {
959  /* The atom 1 names are the same, compare atom 2 */
960  compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
961 
962  /* If atom 1 AND atom 2 are the same, we have a duplicate */
963  if (compare_code == 0)
964  {
965 
966  /* We have a duplicate. So print out a warning*/
967  /* message. Then assign the new values to the */
968  /* tree and free the new_node */
969  //****** BEGIN CHARMM/XPLOR type changes
970  /* we do not care about identical replacement */
971  if ( ((tree->forceconstant != new_node->forceconstant) ||
972  (tree->distance != new_node->distance)) && overwrite)
973  {
974  iout << "\n" << iWARN << "DUPLICATE BOND ENTRY FOR "
975  << new_node->atom1name << "-"
976  << new_node->atom2name
977  << "\nPREVIOUS VALUES k=" << tree->forceconstant
978  << " x0=" << tree->distance
979  << "\n USING VALUES k=" << new_node->forceconstant
980  << " x0=" << new_node->distance
981  << "\n" << endi;
982 
983  tree->forceconstant=new_node->forceconstant;
984  tree->distance=new_node->distance;
985  }
986  //****** END CHARMM/XPLOR type changes
987 
988  delete new_node;
989 
990  return(tree);
991  }
992  }
993 
994  /* We don't have a duplicate, so if the new value is less */
995  /* than the head of the tree, add it to the left child, */
996  /* otherwise add it to the right child */
997  if (compare_code < 0)
998  {
999  tree->left = add_to_bond_tree(new_node, tree->left, overwrite);
1000  }
1001  else
1002  {
1003  tree->right = add_to_bond_tree(new_node, tree->right, overwrite);
1004  }
1005 
1006  return(tree);
1007 }
1008 /* END OF FUNCTION add_to_bond_tree */
1009 
1010 /************************************************************************/
1011 /* */
1012 /* FUNCTION add_angle_param */
1013 /* */
1014 /* INPUTS: */
1015 /* buf - line from paramter file with angle parameters */
1016 /* */
1017 /* this function adds an angle parameter. It parses up the input */
1018 /* line and then adds it to the binary tree used to store the angle */
1019 /* parameters. */
1020 /* */
1021 /************************************************************************/
1022 
1023 void Parameters::add_angle_param(char *buf)
1024 
1025 {
1026  char atom1name[11]; // Type for atom 1
1027  char atom2name[11]; // Type for atom 2
1028  char atom3name[11]; // Type for atom 3
1029  char norm[4]="xxx";
1030  Real forceconstant; // Force constant
1031  Real angle; // Theta 0
1032  Real k_ub; // Urey-Bradley force constant
1033  Real r_ub; // Urey-Bradley distance
1034  int read_count; // count from sscanf
1035  struct angle_params *new_node; // new node in tree
1036 
1037  //****** BEGIN CHARMM/XPLOR type changes
1038  /* parse up the input line with sscanf */
1039  if (paramType == paraXplor)
1040  {
1041  /* read XPLOR format */
1042  read_count=sscanf(buf, "%*s %s %s %s %f %f UB %f %f\n",
1043  atom1name, atom2name, atom3name, &forceconstant, &angle,
1044  &k_ub, &r_ub);
1045  }
1046  else if ((paramType == paraCharmm) && cosAngles) {
1047  read_count=sscanf(buf, "%s %s %s %f %f %3s %f %f\n",
1048  atom1name, atom2name, atom3name, &forceconstant, &angle, norm,
1049  &k_ub, &r_ub);
1050 // printf("%s\n", buf);
1051 // printf("Data: %s %s %s %f %f %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle, norm, k_ub, r_ub);
1052  }
1053  else if (paramType == paraCharmm)
1054  {
1055  /* read CHARMM format */
1056  read_count=sscanf(buf, "%s %s %s %f %f %f %f\n",
1057  atom1name, atom2name, atom3name, &forceconstant, &angle,
1058  &k_ub, &r_ub);
1059 // printf("%s\n", buf);
1060 // printf("Data: %s %s %s %f %f\n", atom1name, atom2name, atom3name, forceconstant, angle);
1061  }
1062  //****** END CHARMM/XPLOR type changes
1063 
1064  /* Check to make sure we got what we expected */
1065  if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
1066  {
1067  char err_msg[512];
1068 
1069  if (paramType == paraXplor)
1070  {
1071  sprintf(err_msg, "BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1072  }
1073  else if (paramType == paraCharmm)
1074  {
1075  sprintf(err_msg, "BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1076  }
1077  NAMD_die(err_msg);
1078  }
1079 
1080  /* Allocate the new node */
1081  new_node = new angle_params;
1082 
1083  if (new_node == NULL)
1084  {
1085  NAMD_die("memory allocation failed in Parameters::add_angle_param");
1086  }
1087 
1088  /* As with the bond, we want the atom type is comes first */
1089  /* alphbetically first between atom 1 and atom 3 to be in */
1090  /* atom 1 so that we can search the tree reliably. */
1091  if (strcasecmp(atom1name, atom3name) < 0)
1092  {
1093  strcpy(new_node->atom1name, atom1name);
1094  strcpy(new_node->atom2name, atom2name);
1095  strcpy(new_node->atom3name, atom3name);
1096  }
1097  else
1098  {
1099  strcpy(new_node->atom3name, atom1name);
1100  strcpy(new_node->atom2name, atom2name);
1101  strcpy(new_node->atom1name, atom3name);
1102  }
1103 
1104  /* Assign the constants and pointer values */
1105  new_node->forceconstant = forceconstant;
1106  new_node->angle = angle;
1107 
1108  if (cosAngles) {
1109  if (strcasecmp("cos",norm)==0) {
1110 // iout << "Info: Using cos mode for angle " << buf << endl;
1111  NumCosAngles++;
1112  new_node->normal = 0;
1113  } else {
1114 // iout << "Info: Using x^2 mode for angle " << buf << endl;
1115  new_node->normal = 1;
1116  }
1117  } else {
1118  new_node->normal = 1;
1119  }
1120 
1121  if (read_count == 7)
1122  {
1123  // Urey-Bradley constants
1124  if (new_node->normal == 0) {
1125  NAMD_die("ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
1126  }
1127  new_node->k_ub = k_ub;
1128  new_node->r_ub = r_ub;
1129  }
1130  else
1131  {
1132  new_node->k_ub = 0.0;
1133  new_node->r_ub = 0.0;
1134  }
1135 
1136  new_node->left = NULL;
1137  new_node->right = NULL;
1138 
1139  /* Insert it into the tree */
1140  anglep = add_to_angle_tree(new_node, anglep);
1141 
1142  return;
1143 }
1144 /* END OF FUNCTION add_angle_param */
1145 
1146 /************************************************************************/
1147 /* */
1148 /* FUNCTION add_to_angle_tree */
1149 /* */
1150 /* INPUTS: */
1151 /* new_node - new node to add to the angle tree */
1152 /* tree - tree to add the node to */
1153 /* */
1154 /* OUTPUTS: */
1155 /* the function returns a pointer to the new tree with the node */
1156 /* added. Most of the time, this will be the same as the value passed*/
1157 /* in, but not in the case where the tree is empty. */
1158 /* */
1159 /* this is a recursive function that adds an angle parameter */
1160 /* to the binary tree storing the angle parameters. If a duplicate */
1161 /* is found, a warning message is printed, the current values in the */
1162 /* tree are replaced with the new values, and the new node is free'd */
1163 /* */
1164 /************************************************************************/
1165 
1166 struct angle_params *Parameters::add_to_angle_tree(struct angle_params *new_node,
1167  struct angle_params *tree)
1168 
1169 {
1170  int compare_code; // Return code from strcasecmp
1171 
1172  /* If the tree is empty, then the new_node is the tree */
1173  if (tree == NULL)
1174  return(new_node);
1175 
1176  /* Compare atom 1 from the new node and the head of the tree */
1177  compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
1178 
1179  if (compare_code == 0)
1180  {
1181  /* Atom 1 is the same, compare atom 2 */
1182  compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
1183 
1184  if (compare_code == 0)
1185  {
1186  /* Atoms 1 & 2 are the same, compare atom 3 */
1187  compare_code = strcasecmp(new_node->atom3name,
1188  tree->atom3name);
1189 
1190  if (compare_code == 0)
1191  {
1192  /* All three atoms were the same, this */
1193  /* is a duplicate. Print a warning */
1194  /* message, replace the current values,*/
1195  /* and free the new node */
1196  //****** BEGIN CHARMM/XPLOR type changes
1197  /* we do not care about identical replacement */
1198  if ((tree->forceconstant != new_node->forceconstant) ||
1199  (tree->angle != new_node->angle) ||
1200  (tree->k_ub != new_node->k_ub) ||
1201  (tree->r_ub != new_node->r_ub) || (tree->normal != new_node->normal))
1202  {
1203  iout << "\n" << iWARN << "DUPLICATE ANGLE ENTRY FOR "
1204  << new_node->atom1name << "-"
1205  << new_node->atom2name << "-"
1206  << new_node->atom3name
1207  << "\nPREVIOUS VALUES k="
1208  << tree->forceconstant << " theta0="
1209  << tree->angle << " k_ub="
1210  << tree->k_ub << " r_ub="
1211  << tree->r_ub
1212  << "\n USING VALUES k="
1213  << new_node->forceconstant << " theta0="
1214  << new_node->angle << " k_ub="
1215  << new_node->k_ub << " r_ub=" << new_node->r_ub
1216  << "\n" << endi;
1217 
1218  tree->forceconstant=new_node->forceconstant;
1219  tree->angle=new_node->angle;
1220  tree->k_ub=new_node->k_ub;
1221  tree->r_ub=new_node->r_ub;
1222  tree->normal=new_node->normal;
1223  }
1224  //****** END CHARMM/XPLOR type changes
1225 
1226  delete new_node;
1227 
1228  return(tree);
1229  }
1230  }
1231  }
1232 
1233  /* Didn't find a duplicate, so if the new_node is smaller */
1234  /* than the current head, add it to the left child. Otherwise */
1235  /* add it to the right child. */
1236  if (compare_code < 0)
1237  {
1238  tree->left = add_to_angle_tree(new_node, tree->left);
1239  }
1240  else
1241  {
1242  tree->right = add_to_angle_tree(new_node, tree->right);
1243  }
1244 
1245  return(tree);
1246 }
1247 /* END OF FUNCTION add_to_angle_tree */
1248 
1249 /************************************************************************/
1250 /* */
1251 /* FUNCTION add_dihedral_param */
1252 /* */
1253 /* INPUTS: */
1254 /* buf - line from paramter file with dihedral parameters */
1255 /* */
1256 /* this function adds an dihedral parameter. It parses up the */
1257 /* input line and then adds it to the binary tree used to store the */
1258 /* dihedral parameters. */
1259 /* */
1260 /************************************************************************/
1261 
1262 void Parameters::add_dihedral_param(char *buf, FILE *fd)
1263 
1264 {
1265  char atom1name[11]; // Type of atom 1
1266  char atom2name[11]; // Type of atom 2
1267  char atom3name[11]; // Type of atom 3
1268  char atom4name[11]; // Type of atom 4
1269  Real forceconstant; // Force constant
1270  int periodicity; // Periodicity
1271  Real phase_shift; // Phase shift
1272  int read_count; // Count from sscanf
1273  struct dihedral_params *new_node; // New node
1274  int multiplicity; // Multiplicity for bonds
1275  int i; // Loop counter
1276  char buffer[513]; // Buffer for new line
1277  int ret_code; // Return code
1278 
1279  //****** BEGIN CHARMM/XPLOR type changes
1280  /* Parse up the input line using sscanf */
1281  if (paramType == paraXplor)
1282  {
1283  /* read XPLOR format */
1284  read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n",
1285  atom1name, atom2name, atom3name, atom4name, &multiplicity,
1286  &forceconstant, &periodicity, &phase_shift);
1287  }
1288  else if (paramType == paraCharmm)
1289  {
1290  /* read CHARMM format */
1291  read_count=sscanf(buf, "%s %s %s %s %f %d %f\n",
1292  atom1name, atom2name, atom3name, atom4name,
1293  &forceconstant, &periodicity, &phase_shift);
1294  multiplicity=1;
1295  }
1296 
1297  if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
1298  {
1299  char err_msg[512];
1300 
1301  sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1302  NAMD_die(err_msg);
1303  }
1304  else if ( (read_count != 7) && (paramType == paraCharmm) )
1305  {
1306  char err_msg[512];
1307 
1308  sprintf(err_msg, "BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1309  NAMD_die(err_msg);
1310  }
1311 
1312  if ( (read_count == 4) && (paramType == paraXplor) )
1313  //****** END CHARMM/XPLOR type changes
1314  {
1315  read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n",
1316  &forceconstant, &periodicity, &phase_shift);
1317 
1318  /* Check to make sure we got what we expected */
1319  if (read_count != 3)
1320  {
1321  char err_msg[512];
1322 
1323  sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1324  NAMD_die(err_msg);
1325  }
1326 
1327  multiplicity = 1;
1328  }
1329 
1330  if (multiplicity > MAX_MULTIPLICITY)
1331  {
1332  char err_msg[181];
1333 
1334  sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
1335  multiplicity, MAX_MULTIPLICITY);
1336  NAMD_die(err_msg);
1337  }
1338 
1339  /* Allocate new node */
1340  new_node = new dihedral_params;
1341 
1342  if (new_node == NULL)
1343  {
1344  NAMD_die("memory allocation failed in Parameters::add_dihedral_param\n");
1345  }
1346 
1347  /* Assign all of the values for this node. Notice that since */
1348  /* the dihedrals and impropers are implemented with a linked */
1349  /* list rather than a binary tree, we don't really care about */
1350  /* the order of the atoms any more */
1351  strcpy(new_node->atom1name, atom1name);
1352  strcpy(new_node->atom2name, atom2name);
1353  strcpy(new_node->atom3name, atom3name);
1354  strcpy(new_node->atom4name, atom4name);
1355  new_node->atom1wild = ! strcasecmp(atom1name, "X");
1356  new_node->atom2wild = ! strcasecmp(atom2name, "X");
1357  new_node->atom3wild = ! strcasecmp(atom3name, "X");
1358  new_node->atom4wild = ! strcasecmp(atom4name, "X");
1359  new_node->multiplicity = multiplicity;
1360  if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1361  new_node->values[0].k = forceconstant;
1362  new_node->values[0].n = periodicity;
1363  new_node->values[0].delta = phase_shift;
1364 
1365  new_node->next = NULL;
1366 
1367  // If the multiplicity is greater than 1, then read in other parameters
1368  if (multiplicity > 1)
1369  {
1370  for (i=1; i<multiplicity; i++)
1371  {
1372  ret_code = NAMD_read_line(fd, buffer);
1373 
1374  // Get rid of comments at the end of a line
1375  if (ret_code == 0)
1376  {
1377  NAMD_remove_comment(buffer);
1378  }
1379 
1380  // Keep reading lines until we get one that isn't blank
1381  while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
1382  {
1383  ret_code = NAMD_read_line(fd, buffer);
1384  }
1385 
1386  if (ret_code != 0)
1387  {
1388  NAMD_die("EOF encoutner in middle of multiple dihedral");
1389  }
1390 
1391  read_count=sscanf(buffer, "%f %d %f\n",
1392  &forceconstant, &periodicity, &phase_shift);
1393 
1394  if (read_count != 3)
1395  {
1396  char err_msg[512];
1397 
1398  sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
1399  NAMD_die(err_msg);
1400  }
1401 
1402  if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1403  new_node->values[i].k = forceconstant;
1404  new_node->values[i].n = periodicity;
1405  new_node->values[i].delta = phase_shift;
1406  }
1407  }
1408 
1409  //****** BEGIN CHARMM/XPLOR type changes
1410  /* Add this node to the list */
1411  if (paramType == paraXplor)
1412  {
1413  add_to_dihedral_list(new_node); // XPLOR
1414  }
1415  else if (paramType == paraCharmm)
1416  {
1417  add_to_charmm_dihedral_list(new_node); // CHARMM
1418  }
1419  //****** END CHARMM/XPLOR type changes
1420 
1421  return;
1422 }
1423 /* END OF FUNCTION add_dihedral_param */
1424 
1425 /************************************************************************/
1426 /* */
1427 /* FUNCTION add_to_dihedral_list */
1428 /* */
1429 /* INPUTS: */
1430 /* new_node - node that is to be added to dihedral_list */
1431 /* */
1432 /* this function adds a new dihedral parameter to the linked list */
1433 /* of dihedral parameters. First, it checks for duplicates. If a */
1434 /* duplicate is found, a warning message is printed, the old values */
1435 /* are replaced with the new values, and the new node is freed. If */
1436 /* Otherwise, the node is added to the list. This list is arranged */
1437 /* so that bods with wildcards are placed at the tail of the list. */
1438 /* This will guarantee that if we just do a linear search, we will */
1439 /* always find an exact match before a wildcard match. */
1440 /* */
1441 /************************************************************************/
1442 
1443 void Parameters::add_to_dihedral_list(
1444  struct dihedral_params *new_node)
1445 
1446 {
1447  static struct dihedral_params *ptr; // position within list
1448  static struct dihedral_params *tail; // Pointer to the end of
1449  // the list so we can add
1450  // entries to the end of the
1451  // list in constant time
1452  int i; // Loop counter
1453 
1454  /* If the list is currently empty, then the new node is the list*/
1455  if (dihedralp == NULL)
1456  {
1457  dihedralp=new_node;
1458  tail=new_node;
1459 
1460  return;
1461  }
1462 
1463  /* The list isn't empty, so check for a duplicate */
1464  ptr=dihedralp;
1465 
1466  while (ptr != NULL)
1467  {
1468  if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1469  (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1470  (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1471  (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1472  ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1473  (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1474  (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1475  (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
1476  {
1477  /* Found a duplicate */
1478  //****** BEGIN CHARMM/XPLOR type changes
1479  /* we do not care about identical replacement */
1480  int echoWarn=0; // echo warning messages ?
1481 
1482  if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
1483 
1484  if (!echoWarn)
1485  {
1486  for (i=0; i<ptr->multiplicity; i++)
1487  {
1488  if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
1489  if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
1490  if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
1491  }
1492  }
1493 
1494  if (echoWarn)
1495  {
1496  iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
1497  << ptr->atom1name << "-"
1498  << ptr->atom2name << "-"
1499  << ptr->atom3name << "-"
1500  << ptr->atom4name
1501  << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
1502 
1503  for (i=0; i<ptr->multiplicity; i++)
1504  {
1505  iout << " k=" << ptr->values[i].k
1506  << " n=" << ptr->values[i].n
1507  << " delta=" << ptr->values[i].delta;
1508  }
1509 
1510  iout << "\nUSING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
1511 
1512  for (i=0; i<new_node->multiplicity; i++)
1513  {
1514  iout << " k=" << new_node->values[i].k
1515  << " n=" << new_node->values[i].n
1516  << " delta=" << new_node->values[i].delta;
1517  }
1518 
1519  iout << endi;
1520 
1521  ptr->multiplicity = new_node->multiplicity;
1522 
1523  for (i=0; i<new_node->multiplicity; i++)
1524  {
1525  ptr->values[i].k = new_node->values[i].k;
1526  ptr->values[i].n = new_node->values[i].n;
1527  ptr->values[i].delta = new_node->values[i].delta;
1528  }
1529 
1530  }
1531  //****** END CHARMM/XPLOR type changes
1532 
1533  delete new_node;
1534 
1535  return;
1536  }
1537 
1538  ptr=ptr->next;
1539  }
1540 
1541  /* Check to see if we have any wildcards. Since specific */
1542  /* entries are to take precedence, we'll put anything without */
1543  /* wildcards at the begining of the list and anything with */
1544  /* wildcards at the end of the list. Then, we can just do a */
1545  /* linear search for a bond and be guaranteed to have specific */
1546  /* entries take precendence over over wildcards */
1547  if ( new_node->atom1wild ||
1548  new_node->atom2wild ||
1549  new_node->atom3wild ||
1550  new_node->atom4wild )
1551  {
1552  /* add to the end of the list */
1553  tail->next=new_node;
1554  tail=new_node;
1555 
1556  return;
1557  }
1558  else
1559  {
1560  /* add to the head of the list */
1561  new_node->next=dihedralp;
1562  dihedralp=new_node;
1563 
1564  return;
1565  }
1566 
1567 }
1568 /* END OF FUNCTION add_to_dihedral_list */
1569 
1570 //****** BEGIN CHARMM/XPLOR type changes
1571 /************************************************************************/
1572 /* */
1573 /* FUNCTION add_to_charmm_dihedral_list */
1574 /* */
1575 /* INPUTS: */
1576 /* new_node - node that is to be added to dihedral_list */
1577 /* */
1578 /* this function adds a new dihedral parameter to the linked list */
1579 /* of dihedral parameters in CHARMM format. */
1580 /* First, it checks for duplicates. If a duplicate is found, a */
1581 /* warning message is printed. If the periodicity is the same as of */
1582 /* a previous dihedral the old values are replaced with the new */
1583 /* values, otherwise, the dihedral is added and the multiplicity is */
1584 /* increased. */
1585 /* Otherwise, the node is added to the list. This list is arranged */
1586 /* so that bonds with wildcards are placed at the tail of the list. */
1587 /* This will guarantee that if we just do a linear search, we will */
1588 /* always find an exact match before a wildcard match. */
1589 /* */
1590 /************************************************************************/
1591 
1592 void Parameters::add_to_charmm_dihedral_list(
1593  struct dihedral_params *new_node)
1594 
1595 {
1596  static struct dihedral_params *ptr; // position within list
1597  static struct dihedral_params *tail; // Pointer to the end of
1598  // the list so we can add
1599  // entries to the end of the
1600  // list in constant time
1601  int i; // Loop counter
1602  int replace; // replace values?
1603 
1604  // keep track of the last dihedral param read to avoid spurious
1605  // error messages.
1606  static struct dihedral_params last_dihedral;
1607 
1608  /* If the list is currently empty, then the new node is the list*/
1609  if (dihedralp == NULL)
1610  {
1611  dihedralp=new_node;
1612  tail=new_node;
1613  memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1614 
1615  return;
1616  }
1617 
1618  /* The list isn't empty, so check for a duplicate */
1619  ptr=dihedralp;
1620 
1621  while (ptr != NULL)
1622  {
1623  int same_as_last = 0;
1624  if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1625  (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1626  (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1627  (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1628  ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1629  (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1630  (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1631  (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) )
1632  )
1633  {
1634  /* Found a duplicate */
1635 
1636  // check for same_as_last. Note: don't believe the echoWarn crap; it controls
1637  // not just whether we print warning messages, but whether we actually change
1638  // values or not!
1639 
1640  if ( ( !strcmp(ptr->atom1name, last_dihedral.atom1name) &&
1641  !strcmp(ptr->atom2name, last_dihedral.atom2name) &&
1642  !strcmp(ptr->atom3name, last_dihedral.atom3name) &&
1643  !strcmp(ptr->atom4name, last_dihedral.atom4name)))
1644  same_as_last = 1;
1645 
1646  //****** BEGIN CHARMM/XPLOR type changes
1647  /* we do not care about identical replacement */
1648  int echoWarn=1; // echo warning messages ?
1649 
1650  // ptr->multiplicity will always be >= new_node->multiplicity
1651  for (i=0; i<ptr->multiplicity; i++)
1652  {
1653  if ((ptr->values[i].k == new_node->values[0].k) &&
1654  (ptr->values[i].n == new_node->values[0].n) &&
1655  (ptr->values[i].delta == new_node->values[0].delta))
1656  {
1657  // found an identical replacement
1658  echoWarn=0;
1659  break;
1660  }
1661 
1662  }
1663 
1664  if (echoWarn)
1665  {
1666  if (!same_as_last) {
1667  iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
1668  << ptr->atom1name << "-"
1669  << ptr->atom2name << "-"
1670  << ptr->atom3name << "-"
1671  << ptr->atom4name
1672  << "\nPREVIOUS VALUES MULTIPLICITY: " << ptr->multiplicity << "\n";
1673  }
1674  replace=0;
1675 
1676  for (i=0; i<ptr->multiplicity; i++)
1677  {
1678  if (!same_as_last) {
1679  iout << " k=" << ptr->values[i].k
1680  << " n=" << ptr->values[i].n
1681  << " delta=" << ptr->values[i].delta << "\n";
1682  }
1683  if (ptr->values[i].n == new_node->values[0].n)
1684  {
1685  iout << iWARN << "IDENTICAL PERIODICITY! REPLACING OLD VALUES BY: \n";
1686  ptr->values[i].k = new_node->values[0].k;
1687  ptr->values[i].delta = new_node->values[0].delta;
1688  iout << " k=" << ptr->values[i].k
1689  << " n=" << ptr->values[i].n
1690  << " delta=" << ptr->values[i].delta<< "\n";
1691  replace=1;
1692  break;
1693  }
1694  }
1695 
1696  if (!replace)
1697  {
1698  ptr->multiplicity += 1;
1699 
1700  if (ptr->multiplicity > MAX_MULTIPLICITY)
1701  {
1702  char err_msg[181];
1703 
1704  sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
1706  NAMD_die(err_msg);
1707  }
1708  if (!same_as_last)
1709  iout << "INCREASING MULTIPLICITY TO: " << ptr->multiplicity << "\n";
1710 
1711  i= ptr->multiplicity - 1;
1712  ptr->values[i].k = new_node->values[0].k;
1713  ptr->values[i].n = new_node->values[0].n;
1714  ptr->values[i].delta = new_node->values[0].delta;
1715 
1716  if (!same_as_last)
1717  iout << " k=" << ptr->values[i].k
1718  << " n=" << ptr->values[i].n
1719  << " delta=" << ptr->values[i].delta<< "\n";
1720  }
1721 
1722  iout << endi;
1723  }
1724  //****** END CHARMM/XPLOR type changes
1725 
1726  memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1727  delete new_node;
1728 
1729  return;
1730  }
1731 
1732  ptr=ptr->next;
1733  }
1734 
1735  /* CHARMM and XPLOR wildcards for dihedrals are luckily the same */
1736  /* Check to see if we have any wildcards. Since specific */
1737  /* entries are to take precedence, we'll put anything without */
1738  /* wildcards at the begining of the list and anything with */
1739  /* wildcards at the end of the list. Then, we can just do a */
1740  /* linear search for a bond and be guaranteed to have specific */
1741  /* entries take precendence over over wildcards */
1742  if ( new_node->atom1wild ||
1743  new_node->atom2wild ||
1744  new_node->atom3wild ||
1745  new_node->atom4wild )
1746  {
1747  /* add to the end of the list */
1748  tail->next=new_node;
1749  tail=new_node;
1750 
1751  memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1752  return;
1753  }
1754  else
1755  {
1756  /* add to the head of the list */
1757  new_node->next=dihedralp;
1758  dihedralp=new_node;
1759 
1760  memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
1761  return;
1762  }
1763 
1764 }
1765 /* END OF FUNCTION add_to_charmm_dihedral_list */
1766 //****** END CHARMM/XPLOR type changes
1767 
1768 /************************************************************************/
1769 /* */
1770 /* FUNCTION add_improper_param */
1771 /* */
1772 /* INPUTS: */
1773 /* buf - line from paramter file with improper parameters */
1774 /* */
1775 /* this function adds an improper parameter. It parses up the */
1776 /* input line and then adds it to the binary tree used to store the */
1777 /* improper parameters. */
1778 /* */
1779 /************************************************************************/
1780 
1781 void Parameters::add_improper_param(char *buf, FILE *fd)
1782 
1783 {
1784  char atom1name[11]; // Atom 1 type
1785  char atom2name[11]; // Atom 2 type
1786  char atom3name[11]; // Atom 3 type
1787  char atom4name[11]; // Atom 4 type
1788  Real forceconstant; // Force constant
1789  int periodicity; // Periodicity
1790  Real phase_shift; // Phase shift
1791  int read_count; // Count from sscanf
1792  struct improper_params *new_node; // New node
1793  int multiplicity; // Multiplicity for bonds
1794  int i; // Loop counter
1795  char buffer[513]; // Buffer for new line
1796  int ret_code; // Return code
1797 
1798  //****** BEGIN CHARMM/XPLOR type changes
1799  /* Parse up the line with sscanf */
1800  if (paramType == paraXplor)
1801  {
1802  /* read XPLOR format */
1803  read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n",
1804  atom1name, atom2name, atom3name, atom4name, &multiplicity,
1805  &forceconstant, &periodicity, &phase_shift);
1806  }
1807  else if (paramType == paraCharmm)
1808  {
1809  /* read CHARMM format */
1810  read_count=sscanf(buf, "%s %s %s %s %f %d %f\n",
1811  atom1name, atom2name, atom3name, atom4name,
1812  &forceconstant, &periodicity, &phase_shift);
1813  multiplicity=1;
1814  }
1815 
1816  if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
1817  {
1818  char err_msg[512];
1819 
1820  sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
1821  NAMD_die(err_msg);
1822  }
1823  else if ( (read_count != 7) && (paramType == paraCharmm) )
1824  {
1825  char err_msg[512];
1826 
1827  sprintf(err_msg, "BAD IMPROPER FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
1828  NAMD_die(err_msg);
1829  }
1830 
1831  if ( (read_count == 4) && (paramType == paraXplor) )
1832  //****** END CHARMM/XPLOR type changes
1833  {
1834  read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n",
1835  &forceconstant, &periodicity, &phase_shift);
1836 
1837  /* Check to make sure we got what we expected */
1838  if (read_count != 3)
1839  {
1840  char err_msg[512];
1841 
1842  sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1843  NAMD_die(err_msg);
1844  }
1845 
1846  multiplicity = 1;
1847  }
1848 
1849  if (multiplicity > MAX_MULTIPLICITY)
1850  {
1851  char err_msg[181];
1852 
1853  sprintf(err_msg, "Multiple improper with multiplicity of %d greater than max of %d",
1854  multiplicity, MAX_MULTIPLICITY);
1855  NAMD_die(err_msg);
1856  }
1857 
1858  /* Allocate a new node */
1859  new_node = new improper_params;
1860 
1861  if (new_node == NULL)
1862  {
1863  NAMD_die("memory allocation failed in Parameters::add_improper_param");
1864  }
1865 
1866  /* Assign the values for this bond. As with the dihedrals, */
1867  /* the atom order doesn't matter */
1868  strcpy(new_node->atom1name, atom1name);
1869  strcpy(new_node->atom2name, atom2name);
1870  strcpy(new_node->atom3name, atom3name);
1871  strcpy(new_node->atom4name, atom4name);
1872  new_node->multiplicity = multiplicity;
1873  if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1874  new_node->values[0].k = forceconstant;
1875  new_node->values[0].n = periodicity;
1876  new_node->values[0].delta = phase_shift;
1877 
1878  new_node->next = NULL;
1879 
1880  // Check to see if this improper has multiple values
1881  if (multiplicity > 1)
1882  {
1883  // Loop through and read the other values
1884  for (i=1; i<multiplicity; i++)
1885  {
1886  ret_code = NAMD_read_line(fd, buffer);
1887 
1888  // Strip off comments at the end of the line
1889  if (ret_code == 0)
1890  {
1891  NAMD_remove_comment(buffer);
1892  }
1893 
1894  // Skip blank lines
1895  while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
1896  {
1897  ret_code = NAMD_read_line(fd, buffer);
1898  }
1899 
1900  if (ret_code != 0)
1901  {
1902  NAMD_die("EOF encoutner in middle of multiple improper");
1903  }
1904 
1905  // Get the values from the line
1906  read_count=sscanf(buffer, "%f %d %f\n",
1907  &forceconstant, &periodicity, &phase_shift);
1908 
1909  if (read_count != 3)
1910  {
1911  char err_msg[512];
1912 
1913  sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
1914  NAMD_die(err_msg);
1915  }
1916 
1917  if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
1918  new_node->values[i].k = forceconstant;
1919  new_node->values[i].n = periodicity;
1920  new_node->values[i].delta = phase_shift;
1921  }
1922  }
1923 
1924  /* Add the paramter to the list */
1925  add_to_improper_list(new_node); // works for both XPLOR & CHARMM
1926 
1927  return;
1928 }
1929 /* END OF FUNCTION add_improper_param */
1930 
1931 /************************************************************************/
1932 /* */
1933 /* FUNCTION add_to_improper_list */
1934 /* */
1935 /* INPUTS: */
1936 /* new_node - node that is to be added to imporper_list */
1937 /* */
1938 /* this function adds a new dihedral parameter to the linked list */
1939 /* of improper parameters. First, it checks for duplicates. If a */
1940 /* duplicate is found, a warning message is printed, the old values */
1941 /* are replaced with the new values, and the new node is freed. If */
1942 /* Otherwise, the node is added to the list. This list is arranged */
1943 /* so that bods with wildcards are placed at the tail of the list. */
1944 /* This will guarantee that if we just do a linear search, we will */
1945 /* always find an exact match before a wildcard match. */
1946 /* */
1947 /************************************************************************/
1948 
1949 void Parameters::add_to_improper_list(struct improper_params *new_node)
1950 
1951 {
1952  int i; // Loop counter
1953  static struct improper_params *ptr; // position within list
1954  static struct improper_params *tail; // Pointer to the end of
1955  // the list so we can add
1956  // entries to the end of the
1957  // list in constant time
1958 
1959  /* If the list is currently empty, then the new node is the list*/
1960  if (improperp == NULL)
1961  {
1962  improperp=new_node;
1963  tail=new_node;
1964 
1965  return;
1966  }
1967 
1968  /* The list isn't empty, so check for a duplicate */
1969  ptr=improperp;
1970 
1971  while (ptr != NULL)
1972  {
1973  if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1974  (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1975  (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1976  (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1977  ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1978  (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1979  (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1980  (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
1981  {
1982  /* Found a duplicate */
1983  //****** BEGIN CHARMM/XPLOR type changes
1984  /* we do not care about identical replacement */
1985  int echoWarn=0; // echo warning messages ?
1986 
1987  if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
1988 
1989  if (!echoWarn)
1990  {
1991  for (i=0; i<ptr->multiplicity; i++)
1992  {
1993  if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
1994  if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
1995  if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
1996  }
1997  }
1998 
1999  if (echoWarn)
2000  {
2001  iout << "\n" << iWARN << "DUPLICATE IMPROPER DIHEDRAL ENTRY FOR "
2002  << ptr->atom1name << "-"
2003  << ptr->atom2name << "-"
2004  << ptr->atom3name << "-"
2005  << ptr->atom4name
2006  << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
2007 
2008  for (i=0; i<ptr->multiplicity; i++)
2009  {
2010  iout << " k=" << ptr->values[i].k
2011  << " n=" << ptr->values[i].n
2012  << " delta=" << ptr->values[i].delta;
2013  }
2014 
2015  iout << "\n" << "USING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
2016 
2017  for (i=0; i<new_node->multiplicity; i++)
2018  {
2019  iout << " k=" << new_node->values[i].k
2020  << " n=" << new_node->values[i].n
2021  << " delta=" << new_node->values[i].delta;
2022  }
2023 
2024  iout << endi;
2025 
2026  ptr->multiplicity = new_node->multiplicity;
2027 
2028  for (i=0; i<new_node->multiplicity; i++)
2029  {
2030  ptr->values[i].k = new_node->values[i].k;
2031  ptr->values[i].n = new_node->values[i].n;
2032  ptr->values[i].delta = new_node->values[i].delta;
2033  }
2034  }
2035  //****** END CHARMM/XPLOR type changes
2036 
2037  delete new_node;
2038 
2039  return;
2040  }
2041 
2042  ptr=ptr->next;
2043  }
2044 
2045  /* Check to see if we have any wildcards. Since specific */
2046  /* entries are to take precedence, we'll put anything without */
2047  /* wildcards at the begining of the list and anything with */
2048  /* wildcards at the end of the list. Then, we can just do a */
2049  /* linear search for a bond and be guaranteed to have specific */
2050  /* entries take precendence over over wildcards */
2051  if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
2052  (strcasecmp(new_node->atom2name, "X") == 0) ||
2053  (strcasecmp(new_node->atom3name, "X") == 0) ||
2054  (strcasecmp(new_node->atom4name, "X") == 0) )
2055  {
2056  /* add to the end of the list */
2057  tail->next=new_node;
2058  tail=new_node;
2059 
2060  return;
2061  }
2062  else
2063  {
2064  /* add to the head of the list */
2065  new_node->next=improperp;
2066  improperp=new_node;
2067 
2068  return;
2069  }
2070 }
2071 /* END OF FUNCTION add_to_improper_list */
2072 
2073 /************************************************************************/
2074 /* */
2075 /* FUNCTION add_crossterm_param */
2076 /* */
2077 /* INPUTS: */
2078 /* buf - line from paramter file with crossterm parameters */
2079 /* */
2080 /* this function adds an crossterm parameter. It parses up the */
2081 /* input line and then adds it to the binary tree used to store the */
2082 /* crossterm parameters. */
2083 /* */
2084 /************************************************************************/
2085 
2086 void Parameters::add_crossterm_param(char *buf, FILE *fd)
2087 
2088 {
2089  char atom1name[11]; // Atom 1 type
2090  char atom2name[11]; // Atom 2 type
2091  char atom3name[11]; // Atom 3 type
2092  char atom4name[11]; // Atom 4 type
2093  char atom5name[11]; // Atom 1 type
2094  char atom6name[11]; // Atom 2 type
2095  char atom7name[11]; // Atom 3 type
2096  char atom8name[11]; // Atom 4 type
2097  int dimension;
2098  int read_count; // Count from sscanf
2099  struct crossterm_params *new_node; // New node
2100  char buffer[513]; // Buffer for new line
2101  int ret_code; // Return code
2102 
2103  /* read CHARMM format */
2104  read_count=sscanf(buf, "%s %s %s %s %s %s %s %s %d\n",
2105  atom1name, atom2name, atom3name, atom4name,
2106  atom5name, atom6name, atom7name, atom8name,
2107  &dimension);
2108 
2109  if ( (read_count != 9) || dimension < 1 || dimension > 1000 )
2110  {
2111  char err_msg[512];
2112 
2113  sprintf(err_msg, "BAD CMAP FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2114  NAMD_die(err_msg);
2115  }
2116 
2117  /* Allocate a new node */
2118  new_node = new crossterm_params(dimension);
2119 
2120  /* Assign the values for this bond. As with the dihedrals, */
2121  /* the atom order doesn't matter */
2122  strcpy(new_node->atom1name, atom1name);
2123  strcpy(new_node->atom2name, atom2name);
2124  strcpy(new_node->atom3name, atom3name);
2125  strcpy(new_node->atom4name, atom4name);
2126  strcpy(new_node->atom5name, atom5name);
2127  strcpy(new_node->atom6name, atom6name);
2128  strcpy(new_node->atom7name, atom7name);
2129  strcpy(new_node->atom8name, atom8name);
2130 
2131  new_node->next = NULL;
2132 
2133  int nterms = dimension * dimension;
2134  int nread = 0;
2135 
2136  // Loop through and read the other values
2137  while ( nread < nterms ) {
2138  ret_code = NAMD_read_line(fd, buffer);
2139 
2140  // Strip off comments at the end of the line
2141  if (ret_code == 0) {
2142  NAMD_remove_comment(buffer);
2143  }
2144 
2145  // Skip blank lines
2146  while ( (ret_code == 0) && (NAMD_blank_string(buffer)) ) {
2147  ret_code = NAMD_read_line(fd, buffer);
2148  if (ret_code == 0) {
2149  NAMD_remove_comment(buffer);
2150  }
2151  }
2152 
2153  if (ret_code != 0) {
2154  NAMD_die("EOF encoutner in middle of CMAP");
2155  }
2156 
2157  // Get the values from the line
2158  read_count=sscanf(buffer, "%lf %lf %lf %lf %lf %lf %lf %lf\n",
2159  new_node->values + nread,
2160  new_node->values + nread+1,
2161  new_node->values + nread+2,
2162  new_node->values + nread+3,
2163  new_node->values + nread+4,
2164  new_node->values + nread+5,
2165  new_node->values + nread+6,
2166  new_node->values + nread+7);
2167 
2168  nread += read_count;
2169 
2170  if (read_count == 0 || nread > nterms) {
2171  char err_msg[512];
2172 
2173  sprintf(err_msg, "BAD CMAP FORMAT IN PARAMETER FILE\nLINE=*%s*\n", buffer);
2174  NAMD_die(err_msg);
2175  }
2176  }
2177 
2178  /* Add the paramter to the list */
2179  add_to_crossterm_list(new_node);
2180 
2181  return;
2182 }
2183 /* END OF FUNCTION add_crossterm_param */
2184 
2185 /************************************************************************/
2186 /* */
2187 /* FUNCTION add_to_crossterm_list */
2188 /* */
2189 /* INPUTS: */
2190 /* new_node - node that is to be added to crossterm_list */
2191 /* */
2192 /* this function adds a new crossterm parameter to the linked list */
2193 /* of crossterm parameters. First, it checks for duplicates. If a */
2194 /* duplicate is found, a warning message is printed, the old values */
2195 /* are replaced with the new values, and the new node is freed. If */
2196 /* Otherwise, the node is added to the list. This list is arranged */
2197 /* so that bods with wildcards are placed at the tail of the list. */
2198 /* This will guarantee that if we just do a linear search, we will */
2199 /* always find an exact match before a wildcard match. */
2200 /* */
2201 /************************************************************************/
2202 
2203 void Parameters::add_to_crossterm_list(struct crossterm_params *new_node)
2204 
2205 {
2206  int i; // Loop counter
2207  static struct crossterm_params *ptr; // position within list
2208  static struct crossterm_params *tail; // Pointer to the end of
2209  // the list so we can add
2210  // entries to the end of the
2211  // list in constant time
2212 
2213  /* If the list is currently empty, then the new node is the list*/
2214  if (crosstermp == NULL)
2215  {
2216  crosstermp=new_node;
2217  tail=new_node;
2218 
2219  return;
2220  }
2221 
2222  /* The list isn't empty, so check for a duplicate */
2223  ptr=crosstermp;
2224 
2225  while (ptr != NULL)
2226  {
2227  if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
2228  (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
2229  (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
2230  (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) &&
2231  (strcasecmp(new_node->atom5name, ptr->atom5name) == 0) &&
2232  (strcasecmp(new_node->atom6name, ptr->atom6name) == 0) &&
2233  (strcasecmp(new_node->atom7name, ptr->atom7name) == 0) &&
2234  (strcasecmp(new_node->atom8name, ptr->atom8name) == 0) ) )
2235  {
2236  /* Found a duplicate */
2237  /* we do not care about identical replacement */
2238  int echoWarn=0; // echo warning messages ?
2239 
2240  if (ptr->dimension != new_node->dimension) {echoWarn=1;}
2241 
2242  if (!echoWarn)
2243  {
2244  int nvals = ptr->dimension * ptr->dimension;
2245  for (i=0; i<nvals; i++)
2246  {
2247  if (ptr->values[i] != new_node->values[i]) {echoWarn=1; break;}
2248  }
2249  }
2250 
2251  if (echoWarn)
2252  {
2253  iout << "\n" << iWARN << "DUPLICATE CMAP ENTRY FOR "
2254  << ptr->atom1name << "-"
2255  << ptr->atom2name << "-"
2256  << ptr->atom3name << "-"
2257  << ptr->atom4name << " "
2258  << ptr->atom5name << "-"
2259  << ptr->atom6name << "-"
2260  << ptr->atom7name << "-"
2261  << ptr->atom8name << ", USING NEW VALUES\n";
2262 
2263  iout << endi;
2264 
2265  ptr->dimension = new_node->dimension;
2266 
2267  BigReal *tmpvalues = ptr->values;
2268  ptr->values = new_node->values;
2269  new_node->values = tmpvalues;
2270  }
2271 
2272  delete new_node;
2273 
2274  return;
2275  }
2276 
2277  ptr=ptr->next;
2278  }
2279 
2280  /* Check to see if we have any wildcards. Since specific */
2281  /* entries are to take precedence, we'll put anything without */
2282  /* wildcards at the begining of the list and anything with */
2283  /* wildcards at the end of the list. Then, we can just do a */
2284  /* linear search for a bond and be guaranteed to have specific */
2285  /* entries take precendence over over wildcards */
2286  if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
2287  (strcasecmp(new_node->atom2name, "X") == 0) ||
2288  (strcasecmp(new_node->atom3name, "X") == 0) ||
2289  (strcasecmp(new_node->atom4name, "X") == 0) ||
2290  (strcasecmp(new_node->atom5name, "X") == 0) ||
2291  (strcasecmp(new_node->atom6name, "X") == 0) ||
2292  (strcasecmp(new_node->atom7name, "X") == 0) ||
2293  (strcasecmp(new_node->atom8name, "X") == 0) )
2294  {
2295  /* add to the end of the list */
2296  tail->next=new_node;
2297  tail=new_node;
2298 
2299  return;
2300  }
2301  else
2302  {
2303  /* add to the head of the list */
2304  new_node->next=crosstermp;
2305  crosstermp=new_node;
2306 
2307  return;
2308  }
2309 }
2310 /* END OF FUNCTION add_to_crossterm_list */
2311 
2312 /************************************************************************/
2313 /* */
2314 /* FUNCTION add_vdw_param */
2315 /* */
2316 /* INPUTS: */
2317 /* buf - line containing the vdw information */
2318 /* */
2319 /* add_vdw_param adds a vdw parameter for an atom to the current */
2320 /* binary tree of values. */
2321 /* */
2322 /************************************************************************/
2323 
2324 void Parameters::add_vdw_param(char *buf)
2325 
2326 {
2327  char atomname[11]; // atom type of paramter
2328  Real sigma; // sigma value for this atom
2329  Real epsilon; // epsilon value for this atom
2330  Real sigma14; // sigma value for 1-4 interactions
2331  Real epsilon14; // epsilon value for 1-4 interactions
2332  Real sqrt26; // 2^(1/6)
2333  int read_count; // count returned by sscanf
2334  struct vdw_params *new_node; // new node for tree
2335 
2336  //****** BEGIN CHARMM/XPLOR type changes
2337  /* Parse up the line with sscanf */
2338  if (paramType == paraXplor)
2339  {
2340  /* read XPLOR format */
2341  read_count=sscanf(buf, "%*s %s %f %f %f %f\n", atomname,
2342  &epsilon, &sigma, &epsilon14, &sigma14);
2343  }
2344  else if (paramType == paraCharmm)
2345  {
2346  /* read CHARMM format */
2347  read_count=sscanf(buf, "%s %*f %f %f %*f %f %f\n", atomname,
2348  &epsilon, &sigma, &epsilon14, &sigma14);
2349  }
2350 
2351  /* Check to make sure we got what we expected */
2352  if ((read_count != 5) && (paramType == paraXplor))
2353  {
2354  char err_msg[512];
2355 
2356  sprintf(err_msg, "BAD vdW FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
2357  NAMD_die(err_msg);
2358  }
2359  else if ( ((read_count != 5) && (read_count != 3)) && (paramType == paraCharmm))
2360  {
2361  char err_msg[512];
2362 
2363  sprintf(err_msg, "BAD vdW FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
2364  NAMD_die(err_msg);
2365  }
2366 
2367  if (paramType == paraCharmm)
2368  {
2369  // convert CHARMM to XPLOR format
2370  epsilon*=-1.;
2371  sqrt26=pow(2.,(1./6.));
2372  sigma=2.*sigma/sqrt26;
2373 
2374  if (read_count == 3)
2375  {
2376  epsilon14=epsilon;
2377  sigma14=sigma;
2378  }
2379  else
2380  {
2381  epsilon14*=-1.;
2382  sigma14=2.*sigma14/sqrt26;
2383  }
2384  }
2385  //****** END CHARMM/XPLOR type changes
2386 
2387  if ( epsilon < 0. || epsilon14 < 0. ) {
2388  iout << iWARN << "Ignoring VDW parameter with negative epsilon:\n"
2389  << buf << "\n" << endi;
2390  return;
2391  }
2392 
2393  /* Allocate a new node */
2394  new_node = new vdw_params;
2395 
2396  if (new_node == NULL)
2397  {
2398  NAMD_die("memory allocation failed in Parameters::add_vdw_param");
2399  }
2400 
2401  /* Assign the values to the new node */
2402  strcpy(new_node->atomname, atomname);
2403  new_node->sigma = sigma;
2404  new_node->sigma14 = sigma14;
2405  new_node->epsilon = epsilon;
2406  new_node->epsilon14 = epsilon14;
2407 
2408  new_node->left = NULL;
2409  new_node->right = NULL;
2410 
2411  /* Add the new node into the tree */
2412  vdwp=add_to_vdw_tree(new_node, vdwp);
2413 
2414  return;
2415 }
2416 /* END OF FUNCTION add_vdw_param */
2417 
2418 /************************************************************************/
2419 /* */
2420 /* FUNCTION add_to_vdw_tree */
2421 /* */
2422 /* INPUTS: */
2423 /* new_node - node to add to tree */
2424 /* tree - tree to add the node to */
2425 /* */
2426 /* OUTPUTS: */
2427 /* the function returns a pointer to the tree with the node added */
2428 /* */
2429 /* this function adds a vdw to the binary tree containing the */
2430 /* parameters. */
2431 /* */
2432 /************************************************************************/
2433 
2434 struct vdw_params *Parameters::add_to_vdw_tree(struct vdw_params *new_node,
2435  struct vdw_params *tree)
2436 
2437 {
2438  int compare_code; // Return code from strcasecmp
2439 
2440  /* If the tree is currently empty, the new node is the tree */
2441  if (tree == NULL)
2442  return(new_node);
2443 
2444  compare_code = strcasecmp(new_node->atomname, tree->atomname);
2445 
2446  /* Check to see if we have a duplicate */
2447  if (compare_code==0)
2448  {
2449  /* We have a duplicate. So print out a warning */
2450  /* message, copy the new values into the current node */
2451  /* of the tree, and then free the new_node */
2452  if ((tree->sigma != new_node->sigma) ||
2453  (tree->epsilon != new_node->epsilon) ||
2454  (tree->sigma14 != new_node->sigma14) ||
2455  (tree->epsilon14 != new_node->epsilon14))
2456  {
2457  iout << iWARN << "DUPLICATE vdW ENTRY FOR " << tree->atomname
2458  << "\nPREVIOUS VALUES sigma=" << tree->sigma
2459  << " epsilon=" << tree->epsilon
2460  << " sigma14=" << tree->sigma14
2461  << " epsilon14=" << tree->epsilon14
2462  << "\n USING VALUES sigma=" << new_node->sigma
2463  << " epsilon=" << new_node->epsilon
2464  << " sigma14=" << new_node->sigma14
2465  << " epsilon14=" << new_node->epsilon14
2466  << "\n" << endi;
2467 
2468  tree->sigma=new_node->sigma;
2469  tree->epsilon=new_node->epsilon;
2470  tree->sigma14=new_node->sigma14;
2471  tree->epsilon14=new_node->epsilon14;
2472  }
2473 
2474  delete new_node;
2475 
2476  return(tree);
2477  }
2478 
2479  /* Otherwise, if the new node is less than the head of */
2480  /* the tree, add it to the left child, and if it is greater */
2481  /* add it to the right child */
2482  if (compare_code < 0)
2483  {
2484  tree->left = add_to_vdw_tree(new_node, tree->left);
2485  }
2486  else
2487  {
2488  tree->right = add_to_vdw_tree(new_node, tree->right);
2489  }
2490 
2491  return(tree);
2492 }
2493 /* END OF FUNCTION add_to_vdw_tree */
2494 
2495 /************************************************************************/
2496 /* */
2497 /* FUNCTION add_table_pair_param */
2498 /* */
2499 /* INPUTS: */
2500 /* buf - line containing the table_pair information */
2501 /* */
2502 /* this function adds a table_pair parameter to the current */
2503 /* parameters. */
2504 /* */
2505 /************************************************************************/
2506 
2507 void Parameters::add_table_pair_param(char *buf)
2508 
2509 {
2510  char atom1name[11]; // Atom 1 name
2511  char atom2name[11]; // Atom 2 name
2512  char tabletype[11]; // Name of pair interaction
2513  int K; // Table entry type for pair
2514  int read_count; // count from sscanf
2515  struct table_pair_params *new_node; // new node
2516 
2517  /* Parse up the input line using sscanf */
2518  read_count=sscanf(buf, "%s %s %s\n", atom1name,
2519  atom2name, tabletype);
2520 
2521  /* Check to make sure we got what we expected */
2522  if ((read_count != 3))
2523  {
2524  char err_msg[512];
2525 
2526  sprintf(err_msg, "BAD TABLE PAIR FORMAT IN PARAMETER FILE\nLINE=*%s*", buf);
2527  NAMD_die(err_msg);
2528  }
2529 
2530  /* Allocate a new node */
2531  new_node = new table_pair_params;
2532 
2533  if (new_node == NULL)
2534  {
2535  NAMD_die("memory allocation failed in Parameters::add_table_pair_param\n");
2536  }
2537 
2538  strcpy(new_node->atom1name, atom1name);
2539  strcpy(new_node->atom2name, atom2name);
2540 
2541  /* Find the proper table type for this pairing */
2542  K = get_int_table_type(tabletype);
2543 // printf("Looking for type %s; got %i\n", tabletype, K);
2544  if (K < 0) {
2545  char err_msg[512];
2546  sprintf(err_msg, "Couldn't find table parameters for table interaction %s!\n", tabletype);
2547  NAMD_die(err_msg);
2548  }
2549 
2550  /* Assign values to this node */
2551  new_node->K = K;
2552 
2553  new_node->next = NULL;
2554 
2555  /* Add this node to the tree */
2556 // printf("Adding pair parameter with index %i\n", K);
2557  add_to_table_pair_list(new_node);
2558 
2559  return;
2560 }
2561 /* END OF FUNCTION add_vdw_par_param */
2562 
2563 /************************************************************************/
2564 /* */
2565 /* FUNCTION add_vdw_pair_param */
2566 /* */
2567 /* INPUTS: */
2568 /* buf - line containing the vdw_pair information */
2569 /* */
2570 /* this function adds a vdw_pair parameter to the current */
2571 /* parameters. */
2572 /* */
2573 /************************************************************************/
2574 
2575 void Parameters::add_vdw_pair_param(char *buf)
2576 
2577 {
2578  char atom1name[11]; // Atom 1 name
2579  char atom2name[11]; // Atom 2 name
2580  Real A; // A value for pair
2581  Real B; // B value for pair
2582  Real A14; // A value for 1-4 ints
2583  Real B14; // B value for 1-4 ints
2584  int read_count; // count from sscanf
2585  struct vdw_pair_params *new_node; // new node
2586 
2587  /* Parse up the input line using sscanf */
2588  if (paramType == paraXplor)
2589  {
2590  /* read XPLOR format */
2591  read_count=sscanf(buf, "%*s %s %s %f %f %f %f\n", atom1name,
2592  atom2name, &A, &B, &A14, &B14);
2593  }
2594  else if (paramType == paraCharmm)
2595  {
2596  Real well, rmin, well14, rmin14;
2597  /* read CHARMM format */
2598  read_count=sscanf(buf, "%s %s %f %f %f %f\n", atom1name,
2599  atom2name, &well, &rmin, &well14, &rmin14);
2600  if ( read_count == 4 ) { well14 = well; rmin14 = rmin; }
2601  A = -1. * well * pow(rmin, 12.);
2602  B = -2. * well * pow(rmin, 6.);
2603  A14 = -1. * well14 * pow(rmin14, 12.);
2604  B14 = -2. * well14 * pow(rmin14, 6.);
2605  }
2606 
2607  /* Check to make sure we got what we expected */
2608  if ((read_count != 6) && (paramType == paraXplor))
2609  {
2610  char err_msg[512];
2611 
2612  sprintf(err_msg, "BAD vdW PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
2613  NAMD_die(err_msg);
2614  }
2615  if ((read_count != 4) && (read_count != 6) && (paramType == paraCharmm))
2616  {
2617  char err_msg[512];
2618 
2619  sprintf(err_msg, "BAD vdW PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2620  NAMD_die(err_msg);
2621  }
2622 
2623 
2624  /* Allocate a new node */
2625  new_node = new vdw_pair_params;
2626 
2627  if (new_node == NULL)
2628  {
2629  NAMD_die("memory allocation failed in Parameters::add_vdw_pair_param\n");
2630  }
2631 
2632  strcpy(new_node->atom1name, atom1name);
2633  strcpy(new_node->atom2name, atom2name);
2634 
2635  /* Assign values to this node */
2636  new_node->A = A;
2637  new_node->A14 = A14;
2638  new_node->B = B;
2639  new_node->B14 = B14;
2640 
2641  new_node->next = NULL;
2642 
2643  /* Add this node to the tree */
2644  add_to_vdw_pair_list(new_node);
2645 
2646  return;
2647 }
2648 /* END OF FUNCTION add_vdw_par_param */
2649 
2650 /************************************************************************/
2651 /* */
2652 /* FUNCTION add_nbthole_pair_param */
2653 /* */
2654 /* INPUTS: */
2655 /* buf - line containing the nbthole_pair information */
2656 /* */
2657 /* this function adds a nbthole_pair parameter to the current */
2658 /* parameters. */
2659 /* */
2660 /************************************************************************/
2661 
2662 void Parameters::add_nbthole_pair_param(char *buf)
2663 
2664 {
2665  char atom1name[11]; // Atom 1 name
2666  char atom2name[11]; // Atom 2 name
2667  Real tholeij; // nonbonded thole pair thole
2668  int read_count; // count from sscanf
2669  struct nbthole_pair_params *new_node; // new node
2670 
2671  /* Parse up the input line using sscanf */
2672  if (paramType == paraCharmm)
2673  {
2674  /* read CHARMM format */
2675  read_count=sscanf(buf, "%s %s %f\n", atom1name,
2676  atom2name, &tholeij);
2677  }
2678 
2679  /* Check to make sure we got what we expected */
2680  if ((read_count != 3) && (paramType == paraCharmm))
2681  {
2682  char err_msg[512];
2683 
2684  sprintf(err_msg, "BAD NBTHOLE PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2685  NAMD_die(err_msg);
2686  }
2687 
2688 
2689  /* Allocate a new node */
2690  new_node = new nbthole_pair_params;
2691 
2692  if (new_node == NULL)
2693  {
2694  NAMD_die("memory allocation failed in Parameters::nbthole_pair_param\n");
2695  }
2696 
2697  strcpy(new_node->atom1name, atom1name);
2698  strcpy(new_node->atom2name, atom2name);
2699 
2700  /* Assign values to this node */
2701  new_node->tholeij = tholeij;
2702 
2703  new_node->next = NULL;
2704 
2705  /* Add this node to the tree */
2706  add_to_nbthole_pair_list(new_node);
2707 
2708  return;
2709 }
2710 /* END OF FUNCTION add_nbthole_par_param */
2711 
2712 /************************************************************************/
2713 /* */
2714 /* FUNCTION add_hb_pair_param */
2715 /* */
2716 /* INPUTS: */
2717 /* buf - line containing the hydrogen bond information */
2718 /* */
2719 /* this function adds data for a hydrogen bond interaction pair */
2720 /* to the hbondParams object. */
2721 /* */
2722 /************************************************************************/
2723 
2724 void Parameters::add_hb_pair_param(char *buf)
2725 
2726 {
2727 #if 0
2728  char a1n[11]; // Atom 1 name
2729  char a2n[11]; // Atom 2 name
2730  Real A, B; // A, B value for pair
2731 
2732  //****** BEGIN CHARMM/XPLOR type changes
2734  /* Parse up the input line using sscanf */
2735  if (paramType == paraXplor) {
2736  if (sscanf(buf, "%*s %s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
2737  char err_msg[512];
2738  sprintf(err_msg, "BAD HBOND PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
2739  NAMD_die(err_msg);
2740  }
2741  }
2742  else if (paramType == paraCharmm) {
2743  if (sscanf(buf, "%s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
2744  char err_msg[512];
2745  sprintf(err_msg, "BAD HBOND PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2746  NAMD_die(err_msg);
2747  }
2748  }
2749  //****** END CHARMM/XPLOR type changes
2750 
2751  /* add data */
2752  if (hbondParams.add_hbond_pair(a1n, a2n, A, B) == FALSE) {
2753  iout << "\n" << iWARN << "Duplicate HBOND parameters for types " << a1n
2754  << " and " << a2n << " found; using latest values." << "\n" << endi;
2755  }
2756 #endif
2757 }
2758 /* END OF FUNCTION add_hb_par_param */
2759 
2760 /************************************************************************/
2761 /* */
2762 /* FUNCTION add_to_table_pair_list */
2763 /* */
2764 /* INPUTS: */
2765 /* new_node - node to be added to list */
2766 /* */
2767 /* This function adds a link to the end of the table_pair_list list */
2768 /* */
2769 /************************************************************************/
2770 
2771 void Parameters::add_to_table_pair_list(struct table_pair_params *new_node)
2772 
2773 {
2774  static struct table_pair_params *tail=NULL;
2775  struct table_pair_params *ptr;
2776  int compare_code;
2777 
2778 
2779  // If the list was empty, then just make the new node the list
2780  if (table_pairp == NULL)
2781  {
2782  table_pairp = new_node;
2783  tail = new_node;
2784  return;
2785  }
2786 
2787  ptr = table_pairp;
2788 
2789  // Now check the list to see if we have a duplicate entry
2790  while (ptr!=NULL)
2791  {
2792  /* Compare atom 1 */
2793  compare_code = strncasecmp(new_node->atom1name, ptr->atom1name, 5);
2794 
2795  if (compare_code == 0)
2796  {
2797  /* Atom 1 is the same, compare atom 2 */
2798  compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
2799 
2800  if (compare_code==0)
2801  {
2802  /* Found a duplicate. Print out a warning */
2803  /* message, assign the values to the current */
2804  /* node in the tree, and then free the new_node*/
2805  iout << iWARN << "DUPLICATE TABLE PAIR ENTRY FOR "
2806  << new_node->atom1name << "-"
2807  << new_node->atom2name
2808  << "\n" << endi;
2809 
2810  ptr->K=new_node->K;
2811 
2812  delete new_node;
2813 
2814  return;
2815  }
2816  }
2817 
2818  ptr = ptr->next;
2819  }
2820 
2821  // We didn't find a duplicate, so add this node to the end
2822  // of the list
2823  tail->next = new_node;
2824  tail = new_node;
2825 }
2826 /* END OF FUNCTION add_to_vdw_pair_list */
2827 
2828 /************************************************************************/
2829 /* */
2830 /* FUNCTION add_to_vdw_pair_list */
2831 /* */
2832 /* INPUTS: */
2833 /* new_node - node to be added to list */
2834 /* */
2835 /* This function adds a link to the end of the vdw_pair_list list */
2836 /* */
2837 /************************************************************************/
2838 
2839 void Parameters::add_to_vdw_pair_list(struct vdw_pair_params *new_node)
2840 
2841 {
2842  static struct vdw_pair_params *tail=NULL;
2843  struct vdw_pair_params *ptr;
2844  int compare_code;
2845 
2846 
2847  // If the list was empty, then just make the new node the list
2848  if (vdw_pairp == NULL)
2849  {
2850  vdw_pairp = new_node;
2851  tail = new_node;
2852  return;
2853  }
2854 
2855  ptr = vdw_pairp;
2856 
2857  // Now check the list to see if we have a duplicate entry
2858  while (ptr!=NULL)
2859  {
2860  /* Compare atom 1 */
2861  compare_code = strcasecmp(new_node->atom1name, ptr->atom1name);
2862 
2863  if (compare_code == 0)
2864  {
2865  /* Atom 1 is the same, compare atom 2 */
2866  compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
2867 
2868  if (compare_code==0)
2869  {
2870  /* Found a duplicate. Print out a warning */
2871  /* message, assign the values to the current */
2872  /* node in the tree, and then free the new_node*/
2873  if ((ptr->A != new_node->A) ||
2874  (ptr->B != new_node->B) ||
2875  (ptr->A14 != new_node->A14) ||
2876  (ptr->B14 != new_node->B14))
2877  {
2878  iout << iWARN << "DUPLICATE vdW PAIR ENTRY FOR "
2879  << new_node->atom1name << "-"
2880  << new_node->atom2name
2881  << "\nPREVIOUS VALUES A=" << ptr->A
2882  << " B=" << ptr->B
2883  << " A14=" << ptr->A14
2884  << " B14" << ptr->B14
2885  << "\n USING VALUES A=" << new_node->A
2886  << " B=" << new_node->B
2887  << " A14=" << new_node->A14
2888  << " B14" << new_node->B14
2889  << "\n" << endi;
2890 
2891  ptr->A=new_node->A;
2892  ptr->B=new_node->B;
2893  ptr->A14=new_node->A14;
2894  ptr->B14=new_node->B14;
2895  }
2896 
2897  delete new_node;
2898 
2899  return;
2900  }
2901  }
2902 
2903  ptr = ptr->next;
2904  }
2905 
2906  // We didn't find a duplicate, so add this node to the end
2907  // of the list
2908  tail->next = new_node;
2909  tail = new_node;
2910 }
2911 /* END OF FUNCTION add_to_vdw_pair_list */
2912 
2913 /************************************************************************/
2914 /* */
2915 /* FUNCTION add_to_nbthole_pair_list */
2916 /* */
2917 /* INPUTS: */
2918 /* new_node - node to be added to list */
2919 /* */
2920 /* This function adds a link to the end of the nbthole_pair_list list */
2921 /* */
2922 /************************************************************************/
2923 
2924 void Parameters::add_to_nbthole_pair_list(struct nbthole_pair_params *new_node)
2925 
2926 {
2927  static struct nbthole_pair_params *tail=NULL;
2928  struct nbthole_pair_params *ptr;
2929  int compare_code;
2930 
2931 
2932  // If the list was empty, then just make the new node the list
2933  if (nbthole_pairp == NULL)
2934  {
2935  nbthole_pairp = new_node;
2936  tail = new_node;
2937  return;
2938  }
2939 
2940  ptr = nbthole_pairp;
2941 
2942  tail->next = new_node;
2943  tail = new_node;
2944 }
2945 /* END OF FUNCTION add_to_nbthole_pair_list */
2946 
2947 /************************************************************************/
2948 /* */
2949 /* FUNCTION done_reading_files */
2950 /* */
2951 /* This function is used to signal the Parameters object that all */
2952 /* of the parameter files have been read. Once the object knows this, */
2953 /* it can set un indexes for all the parameters and transfer the values*/
2954 /* to linear arrays. This will allow constant time access from this */
2955 /* point on. */
2956 /* */
2957 /************************************************************************/
2958 
2960 
2961 {
2962  AllFilesRead = TRUE;
2963 
2964  if (addDrudeBond) {
2965  // default definition for Drude bonds if none given
2966  NumBondParams++;
2967  add_bond_param("X DRUD 500.0 0.0\n", FALSE);
2968  }
2969  // Allocate space for all of the arrays
2970  if (NumBondParams)
2971  {
2973 
2974  if (bond_array == NULL)
2975  {
2976  NAMD_die("memory allocation of bond_array failed!");
2977  }
2978  memset(bond_array, 0, NumBondParams*sizeof(BondValue));
2979  }
2980 
2981  if (NumAngleParams)
2982  {
2984 
2985  if (angle_array == NULL)
2986  {
2987  NAMD_die("memory allocation of angle_array failed!");
2988  }
2989  memset(angle_array, 0, NumAngleParams*sizeof(AngleValue));
2990  for ( Index i=0; i<NumAngleParams; ++i ) {
2991  angle_array[i].normal = 1;
2992  }
2993  }
2994 
2995  if (NumDihedralParams)
2996  {
2998 
2999  if (dihedral_array == NULL)
3000  {
3001  NAMD_die("memory allocation of dihedral_array failed!");
3002  }
3003  memset(dihedral_array, 0, NumDihedralParams*sizeof(DihedralValue));
3004  }
3005 
3006  if (NumImproperParams)
3007  {
3009 
3010  if (improper_array == NULL)
3011  {
3012  NAMD_die("memory allocation of improper_array failed!");
3013  }
3014  memset(improper_array, 0, NumImproperParams*sizeof(ImproperValue));
3015  }
3016 
3017  if (NumCrosstermParams)
3018  {
3021  }
3022 
3023  // JLai
3025  {
3028  }
3029  // End of JLai
3030 
3031  if (NumVdwParams)
3032  {
3033  atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
3035 
3036  if (vdw_array == NULL)
3037  {
3038  NAMD_die("memory allocation of vdw_array failed!");
3039  }
3040  }
3042  {
3044 
3045  if(nbthole_array == NULL)
3046  {
3047  NAMD_die("memory allocation of nbthole_array failed!");
3048  }
3049  }
3050  // Assign indexes to each of the parameters and populate the
3051  // arrays using the binary trees and linked lists that we have
3052  // already read in
3053 
3054  // Note that if parameters have been overwritten (matching
3055  // atom patterns but different parameter values) the tree
3056  // contains fewer elements than Num...Params would suggest.
3057  // The arrays are initialized above because the end values
3058  // may not be occupied. Modifying the Num...Params values
3059  // would break backwards compatibility of memopt extraBonds.
3060 
3061  index_bonds(bondp, 0);
3062  index_angles(anglep, 0);
3063  NumVdwParamsAssigned = index_vdw(vdwp, 0);
3064  index_dihedrals();
3065  index_impropers();
3066  index_crossterms();
3067  convert_nbthole_pairs();
3068  // Convert the vdw pairs
3069  convert_vdw_pairs();
3070  convert_table_pairs();
3071 }
3072 /* END OF FUNCTION done_reading_files */
3073 
3074 /************************************************************************/
3075 /* */
3076 /* FUNCTION index_bonds */
3077 /* */
3078 /* INPUTS: */
3079 /* tree - The tree that is to be indexed */
3080 /* index - index to start with */
3081 /* */
3082 /* This is a recursive routine that will traverse the binary tree */
3083 /* of bond parameters, assigning an index to each one, and copying */
3084 /* the data from the binary tree to the array that will be used from */
3085 /* here on. */
3086 /* */
3087 /************************************************************************/
3088 
3089 Index Parameters::index_bonds(struct bond_params *tree, Index index)
3090 
3091 {
3092  // Tree is empty, do nothing
3093  if (tree==NULL)
3094  return(index);
3095 
3096  // If I have a left subtree, index it first
3097  if (tree->left != NULL)
3098  {
3099  index=index_bonds(tree->left, index);
3100  }
3101 
3102  // Now assign an index to top node and populate array
3103  tree->index = index;
3104  bond_array[index].k = tree->forceconstant;
3105  bond_array[index].x0 = tree->distance;
3106  index++;
3107 
3108  // If I have a right subtree, index it
3109  if (tree->right != NULL)
3110  {
3111  index=index_bonds(tree->right, index);
3112  }
3113 
3114  return(index);
3115 }
3116 /* END OF FUNCTION index_bonds */
3117 
3118 /************************************************************************/
3119 /* */
3120 /* FUNCTION index_angles */
3121 /* */
3122 /* INPUTS: */
3123 /* tree - The tree that is to be indexed */
3124 /* index - index to start with */
3125 /* */
3126 /* This is a recursive routine that will traverse the binary tree */
3127 /* of angle parameters, assigning an index to each one, and copying */
3128 /* the data from the binary tree to the array that will be used from */
3129 /* here on. */
3130 /* */
3131 /************************************************************************/
3132 
3133 Index Parameters::index_angles(struct angle_params *tree, Index index)
3134 
3135 {
3136  // Tree is empty, do nothing
3137  if (tree==NULL)
3138  return(index);
3139 
3140  // If I have a left subtree, index it first
3141  if (tree->left != NULL)
3142  {
3143  index=index_angles(tree->left, index);
3144  }
3145 
3146  // Now assign an index to top node and populate array
3147  tree->index = index;
3148 
3149  angle_array[index].k = tree->forceconstant;
3150  angle_array[index].k_ub = tree->k_ub;
3151  angle_array[index].r_ub = tree->r_ub;
3152  angle_array[index].normal = tree->normal;
3153 
3154  // Convert the angle to radians before storing it
3155  angle_array[index].theta0 = (tree->angle*PI)/180.0;
3156  index++;
3157 
3158  // If I have a right subtree, index it
3159  if (tree->right != NULL)
3160  {
3161  index=index_angles(tree->right, index);
3162  }
3163 
3164  return(index);
3165 }
3166 /* END OF FUNCTION index_angles */
3167 
3168 /************************************************************************/
3169 /* */
3170 /* FUNCTION index_dihedrals */
3171 /* */
3172 /* This function walks down the linked list of dihedral parameters */
3173 /* and assigns an index to each one. It also copies the data from this*/
3174 /* linked list to the arrays that will be used from here on out */
3175 /* */
3176 /************************************************************************/
3177 
3178 void Parameters::index_dihedrals()
3179 
3180 {
3181  struct dihedral_params *ptr; // Current location in list
3182  Index index=0; // Current index value
3183  int i; // Loop counter
3184 
3185  // Allocate an array to hold the multiplicity present in the
3186  // parameter file for each bond. This will be used to check
3187  // the multiplicities that are detected in the psf file
3188 
3189  // This is kind of ugly, but necessary because of the way that
3190  // X-PLOR psf files deal with Charmm22 parameters. The way
3191  // that multiple periodicities are specified is by having
3192  // the bonds appear multiple times in the psf file. This even
3193  // if a bond type has multiple parameters defined, they
3194  // will be used if the bond appears multiple times in the
3195  // psf file. So we need to store the number of parameters
3196  // we have to make sure the psf file doesn't ask for more
3197  // parameters than we really have, and we also need to track
3198  // how many times the bond appears in the psf file so that
3199  // we can decide how many parameters to actually use.
3200  // This is different for CHARMM parameter files as stated below!
3201  maxDihedralMults = new int[NumDihedralParams];
3202 
3203  if (maxDihedralMults == NULL)
3204  {
3205  NAMD_die("memory allocation failed in Parameters::index_dihedrals()");
3206  }
3207 
3208  // Start at the head
3209  ptr = dihedralp;
3210 
3211  while (ptr != NULL)
3212  {
3213  // Copy data to array and assign index
3214 
3215  // Save the multiplicity in another array
3216  maxDihedralMults[index] = ptr->multiplicity;
3217 
3218 
3219  //****** BEGIN CHARMM/XPLOR type changes
3220  if (paramType == paraXplor)
3221  {
3222  // Assign the multiplicity in the actual structure a bogus value
3223  // that we will update in assign_dihedral_index
3225  }
3226  else if (paramType == paraCharmm)
3227  {
3228  // In a CHARMM psf file each dihedral will be only listed once
3229  // even if it has multiple terms. There is no point in comparing
3230  // to the psf information
3232  }
3233  //****** END CHARMM/XPLOR type changes
3234 
3235  for (i=0; i<ptr->multiplicity; i++)
3236  {
3237  dihedral_array[index].values[i].k = ptr->values[i].k;
3238  dihedral_array[index].values[i].n = ptr->values[i].n;
3239 
3240  // Convert the angle to radians before storing it
3241  dihedral_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
3242  }
3243 
3244  ptr->index = index;
3245 
3246  index++;
3247  ptr=ptr->next;
3248  }
3249 }
3250 /* END OF FUNCTION index_dihedrals */
3251 
3252 /************************************************************************/
3253 /* */
3254 /* FUNCTION index_impropers */
3255 /* */
3256 /* This function walks down the linked list of improper parameters */
3257 /* and assigns an index to each one. It also copies the data from this*/
3258 /* linked list to the arrays that will be used from here on out */
3259 /* */
3260 /************************************************************************/
3261 
3262 void Parameters::index_impropers()
3263 
3264 {
3265  struct improper_params *ptr; // Current place in list
3266  Index index=0; // Current index value
3267  int i; // Loop counter
3268 
3269  // Allocate an array to hold the multiplicity present in the
3270  // parameter file for each bond. This will be used to check
3271  // the multiplicities that are detected in the psf file
3272 
3273  // This is kind of ugly, but necessary because of the way that
3274  // X-PLOR psf files deal with Charmm22 parameters. The way
3275  // that multiple periodicities are specified is by having
3276  // the bonds appear multiple times in the psf file. This even
3277  // if a bond type has multiple parameters defined, they
3278  // will be used if the bond appears multiple times in the
3279  // psf file. So we need to store the number of parameters
3280  // we have to make sure the psf file doesn't ask for more
3281  // parameters than we really have, and we also need to track
3282  // how many times the bond appears in the psf file so that
3283  // we can decide how many parameters to actually use.
3284  maxImproperMults = new int[NumImproperParams];
3285 
3286  if (maxImproperMults == NULL)
3287  {
3288  NAMD_die("memory allocation failed in Parameters::index_impropers()");
3289  }
3290 
3291  // Start at the head
3292  ptr = improperp;
3293 
3294  while (ptr != NULL)
3295  {
3296  // Copy data to array and assign index
3297 
3298  // Save the multiplicity in another array
3299  maxImproperMults[index] = ptr->multiplicity;
3300 
3301 
3302  //****** BEGIN CHARMM/XPLOR type changes
3303  if (paramType == paraXplor)
3304  {
3305  // Assign the multiplicity in the actual structure a bogus value
3306  // that we will update in assign_improper_index
3308  }
3309  else if (paramType == paraCharmm)
3310  {
3311  // In a CHARMM psf file each improper will be only listed once
3312  // even if it has multiple terms. There is no point in comparing
3313  // to the psf information
3315  }
3316  //****** END CHARMM/XPLOR type changes
3317 
3318  for (i=0; i<ptr->multiplicity; i++)
3319  {
3320  improper_array[index].values[i].k = ptr->values[i].k;
3321  improper_array[index].values[i].n = ptr->values[i].n;
3322 
3323  // Convert the angle to radians before storing it
3324  improper_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
3325  }
3326 
3327  ptr->index=index;
3328 
3329  index++;
3330  ptr=ptr->next;
3331  }
3332 }
3333 /* END OF FUNCTION index_impropers */
3334 
3335 
3336 /************************************************************************/
3337 /* */
3338 /* FUNCTION index_crossterms */
3339 /* */
3340 /* This function walks down the linked list of crossterm parameters */
3341 /* and assigns an index to each one. It also copies the data from this*/
3342 /* linked list to the arrays that will be used from here on out */
3343 /* */
3344 /************************************************************************/
3345 
3347 
3348 void Parameters::index_crossterms()
3349 
3350 {
3351  struct crossterm_params *ptr; // Current place in list
3352  Index index=0; // Current index value
3353  int i,j,k; // Loop counter
3354 
3355  // Start at the head
3356  ptr = crosstermp;
3357 
3358  while (ptr != NULL)
3359  {
3360  // Copy data to array and assign index
3361 
3362  int N = CrosstermValue::dim - 1;
3363 
3364  if ( ptr->dimension != N ) {
3365  NAMD_die("Sorry, only CMAP dimension of 24 is supported");
3366  }
3367 
3368  k = 0;
3369  for (i=0; i<N; i++) {
3370  for (j=0; j<N; j++) {
3371  crossterm_array[index].c[i][j].d00 = ptr->values[k];
3372  ++k;
3373  }
3374  }
3375  for (i=0; i<N; i++) {
3376  crossterm_array[index].c[i][N].d00 =
3377  crossterm_array[index].c[i][0].d00;
3378  crossterm_array[index].c[N][i].d00 =
3379  crossterm_array[index].c[0][i].d00;
3380  }
3381  crossterm_array[index].c[N][N].d00 =
3382  crossterm_array[index].c[0][0].d00;
3383 
3384  crossterm_setup(&crossterm_array[index].c[0][0]);
3385 
3386  ptr->index=index;
3387 
3388  index++;
3389  ptr=ptr->next;
3390  }
3391 }
3392 /* END OF FUNCTION index_crossterms */
3393 
3394 /************************************************************************/
3395 /* */
3396 /* FUNCTION index_vdw */
3397 /* */
3398 /* INPUTS: */
3399 /* tree - The tree that is to be indexed */
3400 /* index - index to start with */
3401 /* */
3402 /* This is a recursive routine that will traverse the binary tree */
3403 /* of vdw parameters, assigning an index to each one, and copying */
3404 /* the data from the binary tree to the array that will be used from */
3405 /* here on. */
3406 /* */
3407 /************************************************************************/
3408 
3409 Index Parameters::index_vdw(struct vdw_params *tree, Index index)
3410 
3411 {
3412  // If the tree is empty, do nothing
3413  if (tree==NULL)
3414  return(index);
3415 
3416  // If I have a left subtree, populate it first
3417  if (tree->left != NULL)
3418  {
3419  index=index_vdw(tree->left, index);
3420  }
3421 
3422  // Assign the index and copy the data to the array
3423  tree->index = index;
3424 
3425  vdw_array[index].sigma = tree->sigma;
3426  vdw_array[index].epsilon = tree->epsilon;
3427  vdw_array[index].sigma14 = tree->sigma14;
3429 
3430  char *nameloc = atom_type_name(index);
3431  strncpy(nameloc, tree->atomname, MAX_ATOMTYPE_CHARS);
3432  nameloc[MAX_ATOMTYPE_CHARS] = '\0';
3433 
3434 // iout << iWARN << "Parameters: Stored name for type " << index << ": '";
3435 // iout << iWARN << nameloc << "'" << "\n" << endi;
3436 
3437  index++;
3438 
3439  // If I have a right subtree, index it
3440  if (tree->right != NULL)
3441  {
3442  index=index_vdw(tree->right, index);
3443  }
3444 
3445  return(index);
3446 }
3447 /* END OF FUNCTION index_vdw */
3448 
3449 /************************************************************************/
3450 /* */
3451 /* FUNCTION assign_vdw_index */
3452 /* */
3453 /* INPUTS: */
3454 /* atomtype - atom type to find */
3455 /* atom_ptr - pointer to the atom structure to find vdw paramters */
3456 /* for */
3457 /* */
3458 /* OUTPUTS: */
3459 /* the vdw_index field of the atom structure is populated */
3460 /* */
3461 /* This function searches the binary tree of vdw parameters so */
3462 /* that an index can be assigned to this atom. If the parameter is */
3463 /* is found, then the index is assigned. If the parameter is not */
3464 /* found, then NAMD terminates. */
3465 /* */
3466 /************************************************************************/
3467 #ifdef MEM_OPT_VERSION
3468 void Parameters::assign_vdw_index(const char *atomtype, AtomCstInfo *atom_ptr)
3469 #else
3470 void Parameters::assign_vdw_index(const char *atomtype, Atom *atom_ptr)
3471 #endif
3472 {
3473  struct vdw_params *ptr; // Current position in trees
3474  int found=0; // Flag 1->found match
3475  int comp_code; // return code from strcasecmp
3476 
3477  /* Check to make sure the files have all been read */
3478  if (!AllFilesRead)
3479  {
3480  NAMD_die("Tried to assign vdw index before all parameter files were read");
3481  }
3482 
3483  /* Start at the top */
3484  ptr=vdwp;
3485 
3486  /* While we haven't found a match, and we haven't reached */
3487  /* the bottom of the tree, compare the atom passed in with */
3488  /* the current value and decide if we have a match, or if not, */
3489  /* which way to go */
3490  while (!found && (ptr!=NULL))
3491  {
3492  comp_code = strcasecmp(atomtype, ptr->atomname);
3493 
3494  if (comp_code == 0)
3495  {
3496  /* Found a match! */
3497  atom_ptr->vdw_type=ptr->index;
3498  found=1;
3499  }
3500  else if (comp_code < 0)
3501  {
3502  /* Go to the left */
3503  ptr=ptr->left;
3504  }
3505  else
3506  {
3507  /* Go to the right */
3508  ptr=ptr->right;
3509  }
3510  }
3511 
3512  //****** BEGIN CHARMM/XPLOR type changes
3513  if (!found)
3514  {
3515  // since CHARMM allows wildcards "*" in vdw typenames
3516  // we have to look again if necessary, this way, if
3517  // we already had an exact match, this is never executed
3518  size_t windx; // wildcard index
3519 
3520  /* Start again at the top */
3521  ptr=vdwp;
3522 
3523  while (!found && (ptr!=NULL))
3524  {
3525 
3526  // get index of wildcard wildcard, get index
3527  windx= strcspn(ptr->atomname,"*");
3528  if (windx == strlen(ptr->atomname))
3529  {
3530  // there is no wildcard here
3531  comp_code = strcasecmp(atomtype, ptr->atomname);
3532  }
3533  else
3534  {
3535  comp_code = strncasecmp(atomtype, ptr->atomname, windx);
3536  }
3537 
3538  if (comp_code == 0)
3539  {
3540  /* Found a match! */
3541  atom_ptr->vdw_type=ptr->index;
3542  found=1;
3543  char errbuf[100];
3544  sprintf(errbuf,"VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
3545  atomtype, ptr->atomname);
3546  int i;
3547  for(i=0; i<error_msgs.size(); i++) {
3548  if ( strcmp(errbuf,error_msgs[i]) == 0 ) break;
3549  }
3550  if ( i == error_msgs.size() ) {
3551  char *newbuf = new char[strlen(errbuf)+1];
3552  strcpy(newbuf,errbuf);
3553  error_msgs.add(newbuf);
3554  iout << iWARN << newbuf << "\n" << endi;
3555  }
3556  }
3557  else if (comp_code < 0)
3558  {
3559  /* Go to the left */
3560  ptr=ptr->left;
3561  }
3562  else
3563  {
3564  /* Go to the right */
3565  ptr=ptr->right;
3566  }
3567 
3568  }
3569 
3570  }
3571  //****** END CHARMM/XPLOR type changes
3572 
3573  /* Make sure we found it */
3574  if (!found)
3575  {
3576  char err_msg[100];
3577 
3578  sprintf(err_msg, "DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s",
3579  atomtype);
3580  NAMD_die(err_msg);
3581  }
3582 
3583  return;
3584 }
3585 /* END OF FUNCTION assign_vdw_index */
3586 
3587 /************************************************************************
3588  * FUNCTION get_table_pair_params
3589  *
3590  * Inputs:
3591  * atom1 - atom type for atom 1
3592  * atom2 - atom type for atom 2
3593  * K - an integer value for the table type to populate
3594  *
3595  * Outputs:
3596  * If a match is found, K is populated and 1 is returned. Otherwise,
3597  * 0 is returned.
3598  *
3599  * This function finds the proper type index for tabulated nonbonded
3600  * interactions between two atoms. If no such interactions are found,
3601  * the atoms are assumed to interact through standard VDW potentials.
3602  *
3603  ************************************************************************/
3604 
3606  IndexedTablePair *ptr;
3607  Index temp;
3608  int found=FALSE;
3609 
3610  ptr=tab_pair_tree;
3611  //
3612  // We need the smaller type in ind1, so if it isn't already that
3613  // way, switch them */
3614  if (ind1 > ind2)
3615  {
3616  temp = ind1;
3617  ind1 = ind2;
3618  ind2 = temp;
3619  }
3620 
3621  /* While we haven't found a match and we're not at the end */
3622  /* of the tree, compare the bond passed in with the tree */
3623  while (!found && (ptr!=NULL))
3624  {
3625 // printf("Comparing %i with %i and %i with %i\n", ind1, ptr->ind1, ind2, ptr->ind2);
3626  if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
3627  {
3628  found = TRUE;
3629  }
3630  else if ( (ind1 < ptr->ind1) ||
3631  ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
3632  {
3633  /* Go left */
3634  ptr=ptr->left;
3635  }
3636  else
3637  {
3638  /* Go right */
3639  ptr=ptr->right;
3640  }
3641  }
3642 
3643  /* If we found a match, assign the values */
3644  if (found)
3645  {
3646  *K = ptr->K;
3647  return(TRUE);
3648  }
3649  else
3650  {
3651  return(FALSE);
3652  }
3653 }
3654 /* END OF FUNCTION get_table_pair_params */
3655 
3656 /************************************************************************/
3657 /* */
3658 /* FUNCTION get_vdw_pair_params */
3659 /* */
3660 /* INPUTS: */
3661 /* atom1 - atom type for atom 1 */
3662 /* atom2 - atom type for atom 2 */
3663 /* A - A value to populate */
3664 /* B - B value to populate */
3665 /* A14 - A 1-4 value to populate */
3666 /* B14 - B 1-4 value to populate */
3667 /* */
3668 /* OUTPUTS: */
3669 /* If a match is found, A, B, A14, and B14 are all populated and a */
3670 /* 1 is returned. Otherwise, a 0 is returned. */
3671 /* */
3672 /* This function finds a set of vdw_pair paramters. It is given */
3673 /* the two types of atoms involved. This is the only paramter for */
3674 /* which a match is NOT guaranteed. There will only be a match if */
3675 /* there are specific van der waals parameters for the two atom types */
3676 /* involved. */
3677 /* */
3678 /************************************************************************/
3679 
3681  Real *B, Real *A14, Real *B14)
3682 
3683 {
3684  IndexedVdwPair *ptr; // Current location in tree
3685  Index temp; // Temporary value for swithcing
3686  // values
3687  int found=FALSE; // Flag 1-> found a match
3688 
3689  ptr=vdw_pair_tree;
3690 
3691  // We need the smaller type in ind1, so if it isn't already that
3692  // way, switch them */
3693  if (ind1 > ind2)
3694  {
3695  temp = ind1;
3696  ind1 = ind2;
3697  ind2 = temp;
3698  }
3699 
3700  /* While we haven't found a match and we're not at the end */
3701  /* of the tree, compare the bond passed in with the tree */
3702  while (!found && (ptr!=NULL))
3703  {
3704  if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
3705  {
3706  found = TRUE;
3707  }
3708  else if ( (ind1 < ptr->ind1) ||
3709  ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
3710  {
3711  /* Go left */
3712  ptr=ptr->left;
3713  }
3714  else
3715  {
3716  /* Go right */
3717  ptr=ptr->right;
3718  }
3719  }
3720 
3721  /* If we found a match, assign the values */
3722  if (found)
3723  {
3724  *A = ptr->A;
3725  *B = ptr->B;
3726  *A14 = ptr->A14;
3727  *B14 = ptr->B14;
3728 
3729  return(TRUE);
3730  }
3731  else
3732  {
3733  return(FALSE);
3734  }
3735 }
3736 /* END OF FUNCTION get_vdw_pair_params */
3737 
3738 
3739 /************************************************************************/
3740 /* */
3741 /* FUNCTION assign_bond_index */
3742 /* */
3743 /* INPUTS: */
3744 /* atom1 - atom type for atom 1 */
3745 /* atom2 - atom type for atom 2 */
3746 /* bond_ptr - pointer to bond structure to populate */
3747 /* */
3748 /* OUTPUTS: */
3749 /* the structure pointed to by bond_ptr is populated */
3750 /* */
3751 /* This function finds a bond in the binary tree of bond values */
3752 /* and assigns its index. If the bond is found, than the bond_type */
3753 /* field of the bond structure is populated. If the parameter is */
3754 /* not found, NAMD will terminate. */
3755 /* */
3756 /************************************************************************/
3757 
3758 void Parameters::assign_bond_index(const char *atom1, const char *atom2, Bond *bond_ptr)
3759 
3760 {
3761  struct bond_params *ptr; // Current location in tree
3762  int found=0; // Flag 1-> found a match
3763  int cmp_code; // return code from strcasecmp
3764 
3765  /* Check to make sure the files have all been read */
3766  if (!AllFilesRead)
3767  {
3768  NAMD_die("Tried to assign bond index before all parameter files were read");
3769  }
3770 
3771  /* We need atom1 < atom2, so if that's not the way they */
3772  /* were passed, flip them */
3773  if (strcasecmp(atom1, atom2) > 0)
3774  {
3775  const char *tmp_name = atom1;
3776  atom1 = atom2;
3777  atom2 = tmp_name;
3778  }
3779 
3780  /* Start at the top */
3781  ptr=bondp;
3782 
3783  /* While we haven't found a match and we're not at the end */
3784  /* of the tree, compare the bond passed in with the tree */
3785  while (!found && (ptr!=NULL))
3786  {
3787  cmp_code=strcasecmp(atom1, ptr->atom1name);
3788 
3789  if (cmp_code == 0)
3790  {
3791  cmp_code=strcasecmp(atom2, ptr->atom2name);
3792  }
3793 
3794  if (cmp_code == 0)
3795  {
3796  /* Found a match */
3797  found=1;
3798  bond_ptr->bond_type = ptr->index;
3799  }
3800  else if (cmp_code < 0)
3801  {
3802  /* Go left */
3803  ptr=ptr->left;
3804  }
3805  else
3806  {
3807  /* Go right */
3808  ptr=ptr->right;
3809  }
3810  }
3811 
3812  /* Check to see if we found anything */
3813  if (!found)
3814  {
3815  if ((strcmp(atom1, "DRUD")==0 || strcmp(atom2, "DRUD")==0)
3816  && (strcmp(atom1, "X")!=0 && strcmp(atom2, "X")!=0)) {
3817  /* try a wildcard DRUD X match for this Drude bond */
3818  char a1[8] = "DRUD", a2[8] = "X";
3819  return assign_bond_index(a1, a2, bond_ptr); /* recursive call */
3820  }
3821  else {
3822  char err_msg[128];
3823 
3824  sprintf(err_msg, "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
3825  atom1, atom2, bond_ptr->atom1+1, bond_ptr->atom2+1);
3826  NAMD_die(err_msg);
3827  }
3828  }
3829 
3830  return;
3831 }
3832 /* END OF FUNCTION assign_bond_index */
3833 
3834 /************************************************************************/
3835 /* */
3836 /* FUNCTION assign_angle_index */
3837 /* */
3838 /* INPUTS: */
3839 /* atom1 - atom type for atom 1 */
3840 /* atom2 - atom type for atom 2 */
3841 /* atom3 - atom type for atom 3 */
3842 /* angle_ptr - pointer to angle structure to populate */
3843 /* */
3844 /* OUTPUTS: */
3845 /* the structure pointed to by angle_ptr is populated */
3846 /* */
3847 /* This function assigns an angle index to a specific angle. */
3848 /* It searches the binary tree of angle parameters for the appropriate*/
3849 /* values. If they are found, the index is assigned. If they are */
3850 /* not found, then NAMD will terminate. */
3851 /* */
3852 /************************************************************************/
3853 
3854 void Parameters::assign_angle_index(const char *atom1, const char *atom2, const char*atom3,
3855  Angle *angle_ptr, int notFoundIndex)
3856 
3857 {
3858  struct angle_params *ptr; // Current position in tree
3859  int comp_val; // value from strcasecmp
3860  int found=0; // flag 1->found a match
3861 
3862  /* Check to make sure the files have all been read */
3863  if (!AllFilesRead)
3864  {
3865  NAMD_die("Tried to assign angle index before all parameter files were read");
3866  }
3867 
3868  /* We need atom1 < atom3. If that was not what we were */
3869  /* passed, switch them */
3870  if (strcasecmp(atom1, atom3) > 0)
3871  {
3872  const char *tmp_name = atom1;
3873  atom1 = atom3;
3874  atom3 = tmp_name;
3875  }
3876 
3877  /* Start at the top */
3878  ptr=anglep;
3879 
3880  /* While we don't have a match and we haven't reached the */
3881  /* bottom of the tree, compare values */
3882  while (!found && (ptr != NULL))
3883  {
3884  comp_val = strcasecmp(atom1, ptr->atom1name);
3885 
3886  if (comp_val == 0)
3887  {
3888  /* Atom 1 matches, so compare atom 2 */
3889  comp_val = strcasecmp(atom2, ptr->atom2name);
3890 
3891  if (comp_val == 0)
3892  {
3893  /* Atoms 1&2 match, try atom 3 */
3894  comp_val = strcasecmp(atom3, ptr->atom3name);
3895  }
3896  }
3897 
3898  if (comp_val == 0)
3899  {
3900  /* Found a match */
3901  found = 1;
3902  angle_ptr->angle_type = ptr->index;
3903  }
3904  else if (comp_val < 0)
3905  {
3906  /* Go left */
3907  ptr=ptr->left;
3908  }
3909  else
3910  {
3911  /* Go right */
3912  ptr=ptr->right;
3913  }
3914  }
3915 
3916  /* Make sure we found a match */
3917  if (!found)
3918  {
3919  char err_msg[128];
3920 
3921  sprintf(err_msg, "UNABLE TO FIND ANGLE PARAMETERS FOR %s %s %s (ATOMS %i %i %i)",
3922  atom1, atom2, atom3, angle_ptr->atom1+1, angle_ptr->atom2+1, angle_ptr->atom3+1);
3923 
3924  if ( notFoundIndex ) {
3925  angle_ptr->angle_type = notFoundIndex;
3926  iout << iWARN << err_msg << "\n" << endi;
3927  return;
3928  } else NAMD_die(err_msg);
3929  }
3930 
3931  return;
3932 }
3933 /* END OF FUNCTION assign_angle_index */
3934 
3935 /************************************************************************/
3936 /* */
3937 /* FUNCTION assign_dihedral_index */
3938 /* */
3939 /* INPUTS: */
3940 /* atom1 - atom type for atom 1 */
3941 /* atom2 - atom type for atom 2 */
3942 /* atom3 - atom type for atom 3 */
3943 /* atom4 - atom type for atom 4 */
3944 /* dihedral_ptr - pointer to dihedral structure to populate */
3945 /* multiplicity - Multiplicity to assign to this bond */
3946 /* */
3947 /* OUTPUTS: */
3948 /* the structure pointed to by dihedral_ptr is populated */
3949 /* */
3950 /* This function searchs the linked list of dihedral parameters for*/
3951 /* a given bond. If a match is found, the dihedral type is assigned. */
3952 /* If no match is found, NAMD terminates */
3953 /* */
3954 /************************************************************************/
3955 
3956 void Parameters::assign_dihedral_index(const char *atom1, const char *atom2, const char *atom3,
3957  const char *atom4, Dihedral *dihedral_ptr,
3958  int multiplicity, int notFoundIndex)
3959 
3960 {
3961  struct dihedral_params *ptr; // Current position in list
3962  int found=0; // Flag 1->found a match
3963 
3964  /* Start at the begining of the list */
3965  ptr=dihedralp;
3966 
3967  /* While we haven't found a match and we haven't reached */
3968  /* the end of the list, keep looking */
3969  while (!found && (ptr!=NULL))
3970  {
3971  /* Do a linear search through the linked list of */
3972  /* dihedral parameters. Since the list is arranged */
3973  /* with wildcard paramters at the end of the list, we */
3974  /* can simply do a linear search and be guaranteed that*/
3975  /* we will find exact matches before wildcard matches. */
3976  /* Also, we must check for an exact match, and a match */
3977  /* in reverse, since they are really the same */
3978  /* physically. */
3979  if ( ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom1)==0) ) &&
3980  ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom2)==0) ) &&
3981  ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom3)==0) ) &&
3982  ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom4)==0) ) )
3983  {
3984  /* Found an exact match */
3985  found=1;
3986  }
3987  else if ( ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom1)==0) ) &&
3988  ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom2)==0) ) &&
3989  ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom3)==0) ) &&
3990  ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom4)==0) ) )
3991  {
3992  /* Found a reverse match */
3993  found=1;
3994  }
3995  else
3996  {
3997  /* Didn't find a match, go to the next node */
3998  ptr=ptr->next;
3999  }
4000  }
4001 
4002  /* Make sure we found a match */
4003  if (!found)
4004  {
4005  char err_msg[128];
4006 
4007  sprintf(err_msg, "UNABLE TO FIND DIHEDRAL PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
4008  atom1, atom2, atom3, atom4, dihedral_ptr->atom1+1, dihedral_ptr->atom2+1,
4009  dihedral_ptr->atom3+1, dihedral_ptr->atom4+1);
4010 
4011  if ( notFoundIndex ) {
4012  dihedral_ptr->dihedral_type = notFoundIndex;
4013  iout << iWARN << err_msg << "\n" << endi;
4014  return;
4015  } else NAMD_die(err_msg);
4016  }
4017 
4018  if (paramType == paraXplor) {
4019  // Check to make sure the number of multiples specified in the psf
4020  // file doesn't exceed the number of parameters in the parameter
4021  // files
4022  if (multiplicity > maxDihedralMults[ptr->index])
4023  {
4024  char err_msg[257];
4025 
4026  sprintf(err_msg, "Multiplicity of Paramters for diehedral bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxDihedralMults[ptr->index]);
4027  NAMD_die(err_msg);
4028  }
4029 
4030  // If the multiplicity from the current bond is larger than that
4031  // seen in the past, increase the multiplicity for this bond
4032  if (multiplicity > dihedral_array[ptr->index].multiplicity)
4033  {
4035  }
4036  }
4037 
4038  dihedral_ptr->dihedral_type = ptr->index;
4039 
4040  return;
4041 }
4042 /* END OF FUNCTION assign_dihedral_index */
4043 
4044 /************************************************************************/
4045 /* */
4046 /* FUNCTION assign_improper_index */
4047 /* */
4048 /* INPUTS: */
4049 /* atom1 - atom type for atom 1 */
4050 /* atom2 - atom type for atom 2 */
4051 /* atom3 - atom type for atom 3 */
4052 /* atom4 - atom type for atom 4 */
4053 /* improper_ptr - pointer to improper structure to populate */
4054 /* multiplicity - Multiplicity to assign to this bond */
4055 /* */
4056 /* OUTPUTS: */
4057 /* the structure pointed to by improper_ptr is populated */
4058 /* */
4059 /* This function searchs the linked list of improper parameters for*/
4060 /* a given bond. If a match is found, the improper_type is assigned. */
4061 /* If no match is found, NAMD will terminate. */
4062 /* */
4063 /************************************************************************/
4064 
4065 void Parameters::assign_improper_index(const char *atom1, const char *atom2, const char *atom3,
4066  const char *atom4, Improper *improper_ptr,
4067  int multiplicity)
4068 
4069 {
4070  struct improper_params *ptr; // Current position in list
4071  int found=0; // Flag 1->found a match
4072 
4073  /* Start at the head of the list */
4074  ptr=improperp;
4075 
4076  /* While we haven't fuond a match and haven't reached the end */
4077  /* of the list, keep looking */
4078  while (!found && (ptr!=NULL))
4079  {
4080  /* Do a linear search through the linked list of */
4081  /* improper parameters. Since the list is arranged */
4082  /* with wildcard paramters at the end of the list, we */
4083  /* can simply do a linear search and be guaranteed that*/
4084  /* we will find exact matches before wildcard matches. */
4085  /* Also, we must check for an exact match, and a match */
4086  /* in reverse, since they are really the same */
4087  /* physically. */
4088  if ( ( (strcasecmp(ptr->atom1name, atom1)==0) ||
4089  (strcasecmp(ptr->atom1name, "X")==0) ) &&
4090  ( (strcasecmp(ptr->atom2name, atom2)==0) ||
4091  (strcasecmp(ptr->atom2name, "X")==0) ) &&
4092  ( (strcasecmp(ptr->atom3name, atom3)==0) ||
4093  (strcasecmp(ptr->atom3name, "X")==0) ) &&
4094  ( (strcasecmp(ptr->atom4name, atom4)==0) ||
4095  (strcasecmp(ptr->atom4name, "X")==0) ) )
4096  {
4097  /* Found an exact match */
4098  found=1;
4099  }
4100  else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) ||
4101  (strcasecmp(ptr->atom4name, "X")==0) ) &&
4102  ( (strcasecmp(ptr->atom3name, atom2)==0) ||
4103  (strcasecmp(ptr->atom3name, "X")==0) ) &&
4104  ( (strcasecmp(ptr->atom2name, atom3)==0) ||
4105  (strcasecmp(ptr->atom2name, "X")==0) ) &&
4106  ( (strcasecmp(ptr->atom1name, atom4)==0) ||
4107  (strcasecmp(ptr->atom1name, "X")==0) ) )
4108  {
4109  /* Found a reverse match */
4110  found=1;
4111  }
4112  else
4113  {
4114  /* Didn't find a match, go to the next node */
4115  ptr=ptr->next;
4116  }
4117  }
4118 
4119  /* Make sure we found a match */
4120  if (!found)
4121  {
4122  char err_msg[128];
4123 
4124  sprintf(err_msg, "UNABLE TO FIND IMPROPER PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
4125  atom1, atom2, atom3, atom4, improper_ptr->atom1+1, improper_ptr->atom2+1,
4126  improper_ptr->atom3+1, improper_ptr->atom4+1);
4127 
4128  NAMD_die(err_msg);
4129  }
4130 
4131  if (paramType == paraXplor) {
4132  // Check to make sure the number of multiples specified in the psf
4133  // file doesn't exceed the number of parameters in the parameter
4134  // files
4135  if (multiplicity > maxImproperMults[ptr->index])
4136  {
4137  char err_msg[257];
4138 
4139  sprintf(err_msg, "Multiplicity of Paramters for improper bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxImproperMults[ptr->index]);
4140  NAMD_die(err_msg);
4141  }
4142 
4143  // If the multiplicity from the current bond is larger than that
4144  // seen in the past, increase the multiplicity for this bond
4145  if (multiplicity > improper_array[ptr->index].multiplicity)
4146  {
4148  }
4149  }
4150 
4151  /* Assign the constants */
4152  improper_ptr->improper_type = ptr->index;
4153 
4154  return;
4155 }
4156 /* END OF FUNCTION assign_improper_index */
4157 
4158 /************************************************************************/
4159 /* */
4160 /* FUNCTION assign_crossterm_index */
4161 /* */
4162 /************************************************************************/
4163 
4164 void Parameters::assign_crossterm_index(const char *atom1, const char *atom2, const char *atom3,
4165  const char *atom4, const char *atom5, const char *atom6, const char *atom7,
4166  const char *atom8, Crossterm *crossterm_ptr)
4167 {
4168  struct crossterm_params *ptr; // Current position in list
4169  int found=0; // Flag 1->found a match
4170 
4171  /* Start at the head of the list */
4172  ptr=crosstermp;
4173 
4174  /* While we haven't fuond a match and haven't reached the end */
4175  /* of the list, keep looking */
4176  while (!found && (ptr!=NULL))
4177  {
4178  /* Do a linear search through the linked list of */
4179  /* crossterm parameters. Since the list is arranged */
4180  /* with wildcard paramters at the end of the list, we */
4181  /* can simply do a linear search and be guaranteed that*/
4182  /* we will find exact matches before wildcard matches. */
4183  /* Also, we must check for an exact match, and a match */
4184  /* in reverse, since they are really the same */
4185  /* physically. */
4186  if ( ( (strcasecmp(ptr->atom1name, atom1)==0) ||
4187  (strcasecmp(ptr->atom1name, "X")==0) ) &&
4188  ( (strcasecmp(ptr->atom2name, atom2)==0) ||
4189  (strcasecmp(ptr->atom2name, "X")==0) ) &&
4190  ( (strcasecmp(ptr->atom3name, atom3)==0) ||
4191  (strcasecmp(ptr->atom3name, "X")==0) ) &&
4192  ( (strcasecmp(ptr->atom4name, atom4)==0) ||
4193  (strcasecmp(ptr->atom4name, "X")==0) ) )
4194  {
4195  /* Found an exact match */
4196  found=1;
4197  }
4198  else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) ||
4199  (strcasecmp(ptr->atom4name, "X")==0) ) &&
4200  ( (strcasecmp(ptr->atom3name, atom2)==0) ||
4201  (strcasecmp(ptr->atom3name, "X")==0) ) &&
4202  ( (strcasecmp(ptr->atom2name, atom3)==0) ||
4203  (strcasecmp(ptr->atom2name, "X")==0) ) &&
4204  ( (strcasecmp(ptr->atom1name, atom4)==0) ||
4205  (strcasecmp(ptr->atom1name, "X")==0) ) )
4206  {
4207  /* Found a reverse match */
4208  found=1;
4209  }
4210  if ( ! found ) {
4211  /* Didn't find a match, go to the next node */
4212  ptr=ptr->next;
4213  continue;
4214  }
4215  found = 0;
4216  if ( ( (strcasecmp(ptr->atom5name, atom5)==0) ||
4217  (strcasecmp(ptr->atom5name, "X")==0) ) &&
4218  ( (strcasecmp(ptr->atom6name, atom6)==0) ||
4219  (strcasecmp(ptr->atom6name, "X")==0) ) &&
4220  ( (strcasecmp(ptr->atom7name, atom7)==0) ||
4221  (strcasecmp(ptr->atom7name, "X")==0) ) &&
4222  ( (strcasecmp(ptr->atom8name, atom8)==0) ||
4223  (strcasecmp(ptr->atom8name, "X")==0) ) )
4224  {
4225  /* Found an exact match */
4226  found=1;
4227  }
4228  else if ( ( (strcasecmp(ptr->atom8name, atom5)==0) ||
4229  (strcasecmp(ptr->atom8name, "X")==0) ) &&
4230  ( (strcasecmp(ptr->atom7name, atom6)==0) ||
4231  (strcasecmp(ptr->atom7name, "X")==0) ) &&
4232  ( (strcasecmp(ptr->atom6name, atom7)==0) ||
4233  (strcasecmp(ptr->atom6name, "X")==0) ) &&
4234  ( (strcasecmp(ptr->atom5name, atom8)==0) ||
4235  (strcasecmp(ptr->atom5name, "X")==0) ) )
4236  {
4237  /* Found a reverse match */
4238  found=1;
4239  }
4240  if ( ! found ) {
4241  /* Didn't find a match, go to the next node */
4242  ptr=ptr->next;
4243  }
4244  }
4245 
4246  /* Make sure we found a match */
4247  if (!found)
4248  {
4249  char err_msg[128];
4250 
4251  sprintf(err_msg, "UNABLE TO FIND CROSSTERM PARAMETERS FOR %s %s %s %s %s %s %s %s\n"
4252  "(ATOMS %i %i %i %i %i %i %i %i)",
4253  atom1, atom2, atom3, atom4, atom5, atom6, atom7, atom8,
4254  crossterm_ptr->atom1+1, crossterm_ptr->atom1+2, crossterm_ptr->atom1+3, crossterm_ptr->atom4+1,
4255  crossterm_ptr->atom1+5, crossterm_ptr->atom1+6, crossterm_ptr->atom1+7, crossterm_ptr->atom8+1);
4256 
4257  NAMD_die(err_msg);
4258  }
4259 
4260  /* Assign the constants */
4261  crossterm_ptr->crossterm_type = ptr->index;
4262 
4263  return;
4264 }
4265 /* END OF FUNCTION assign_improper_index */
4266 
4267 /************************************************************************/
4268 /* */
4269 /* FUNCTION free_bond_tree */
4270 /* */
4271 /* INPUTS: */
4272 /* bond_ptr - pointer to bond tree to free */
4273 /* */
4274 /* this is a recursive function that is used to free the memory */
4275 /* allocated for a bond paramter tree. It makes recursive calls to */
4276 /* free the left an right subtress, and then frees the head. It is */
4277 /* only called by the destructor */
4278 /* */
4279 /************************************************************************/
4280 
4281 void Parameters::free_bond_tree(struct bond_params *bond_ptr)
4282 
4283 {
4284  if (bond_ptr->left != NULL)
4285  {
4286  free_bond_tree(bond_ptr->left);
4287  }
4288 
4289  if (bond_ptr->right != NULL)
4290  {
4291  free_bond_tree(bond_ptr->right);
4292  }
4293 
4294  delete bond_ptr;
4295 
4296  return;
4297 }
4298 /* END OF FUNCTION free_bond_tree */
4299 
4300 /************************************************************************/
4301 /* */
4302 /* FUNCTION free_angle_tree */
4303 /* */
4304 /* INPUTS: */
4305 /* angle_ptr - pointer to angle tree to free */
4306 /* */
4307 /* this is a recursive function that is used to free the memory */
4308 /* allocated for a angle paramter tree. It makes recursive calls to */
4309 /* free the left an right subtress, and then frees the head. It is */
4310 /* only called by the destructor */
4311 /* */
4312 /************************************************************************/
4313 
4314 void Parameters::free_angle_tree(struct angle_params *angle_ptr)
4315 
4316 {
4317  if (angle_ptr->left != NULL)
4318  {
4319  free_angle_tree(angle_ptr->left);
4320  }
4321 
4322  if (angle_ptr->right != NULL)
4323  {
4324  free_angle_tree(angle_ptr->right);
4325  }
4326 
4327  delete angle_ptr;
4328 
4329  return;
4330 }
4331 /* END OF FUNCTION free_angle_tree */
4332 
4333 /************************************************************************/
4334 /* */
4335 /* FUNCTION free_dihedral_list */
4336 /* */
4337 /* INPUTS: */
4338 /* dih_ptr - pointer to the list to free */
4339 /* */
4340 /* this function frees a linked list of dihedral parameters. It */
4341 /* is only called by the destructor. */
4342 /* */
4343 /************************************************************************/
4344 
4345 void Parameters::free_dihedral_list(struct dihedral_params *dih_ptr)
4346 
4347 {
4348  struct dihedral_params *ptr; // Current position in list
4349  struct dihedral_params *next; // Next position in list
4350 
4351  ptr=dih_ptr;
4352 
4353  while (ptr != NULL)
4354  {
4355  next=ptr->next;
4356  delete ptr;
4357  ptr=next;
4358  }
4359 
4360  return;
4361 }
4362 /* END OF FUNCTION free_dihedral_list */
4363 
4364 /************************************************************************/
4365 /* */
4366 /* FUNCTION free_improper_list */
4367 /* */
4368 /* INPUTS: */
4369 /* imp_ptr - pointer to the list to free */
4370 /* */
4371 /* this function frees a linked list of improper parameters. It */
4372 /* is only called by the destructor. */
4373 /* */
4374 /************************************************************************/
4375 
4376 void Parameters::free_improper_list(struct improper_params *imp_ptr)
4377 
4378 {
4379  struct improper_params *ptr; // Current position in list
4380  struct improper_params *next; // Next position in list
4381 
4382  ptr=imp_ptr;
4383 
4384  while (ptr != NULL)
4385  {
4386  next=ptr->next;
4387  delete ptr;
4388  ptr=next;
4389  }
4390 
4391  return;
4392 }
4393 /* END OF FUNCTION free_improper_list */
4394 
4395 /************************************************************************/
4396 /* */
4397 /* FUNCTION free_crossterm_list */
4398 /* */
4399 /* INPUTS: */
4400 /* imp_ptr - pointer to the list to free */
4401 /* */
4402 /* this function frees a linked list of crossterm parameters. It */
4403 /* is only called by the destructor. */
4404 /* */
4405 /************************************************************************/
4406 
4407 void Parameters::free_crossterm_list(struct crossterm_params *imp_ptr)
4408 
4409 {
4410  struct crossterm_params *ptr; // Current position in list
4411  struct crossterm_params *next; // Next position in list
4412 
4413  ptr=imp_ptr;
4414 
4415  while (ptr != NULL)
4416  {
4417  next=ptr->next;
4418  delete ptr;
4419  ptr=next;
4420  }
4421 
4422  return;
4423 }
4424 /* END OF FUNCTION free_crossterm_list */
4425 
4426 /************************************************************************/
4427 /* */
4428 /* FUNCTION free_vdw_tree */
4429 /* */
4430 /* INPUTS: */
4431 /* vdw_ptr - pointer to vdw tree to free */
4432 /* */
4433 /* this is a recursive function that is used to free the memory */
4434 /* allocated for a vdw paramter tree. It makes recursive calls to */
4435 /* free the left an right subtress, and then frees the head. It is */
4436 /* only called by the destructor */
4437 /* */
4438 /************************************************************************/
4439 
4440 void Parameters::free_vdw_tree(struct vdw_params *vdw_ptr)
4441 
4442 {
4443  if (vdw_ptr->left != NULL)
4444  {
4445  free_vdw_tree(vdw_ptr->left);
4446  }
4447 
4448  if (vdw_ptr->right != NULL)
4449  {
4450  free_vdw_tree(vdw_ptr->right);
4451  }
4452 
4453  delete vdw_ptr;
4454 
4455  return;
4456 }
4457 /* END OF FUNCTION free_vdw_tree */
4458 
4459 /************************************************************************/
4460 /* */
4461 /* FUNCTION free_vdw_pair_list */
4462 /* */
4463 /* This function frees the vdw_pair_list */
4464 /* */
4465 /************************************************************************/
4466 
4467 void Parameters::free_vdw_pair_list()
4468 {
4469  struct vdw_pair_params *ptr, *next;
4470 
4471  ptr=vdw_pairp;
4472 
4473  while (ptr != NULL)
4474  {
4475  next = ptr->next;
4476 
4477  delete ptr;
4478 
4479  ptr = next;
4480  }
4481 
4482  vdw_pairp = NULL;
4483 }
4484 /* END OF FUNCTION free_vdw_pair_list */
4485 
4486 /************************************************************************/
4487 /* */
4488 /* FUNCTION free_nbthole_pair_list */
4489 /* */
4490 /* This function frees the nbthole_pair_list */
4491 /* */
4492 /************************************************************************/
4493 
4494 void Parameters::free_nbthole_pair_list()
4495 {
4496  struct nbthole_pair_params *ptr, *next;
4497 
4498  ptr=nbthole_pairp;
4499 
4500  while (ptr != NULL)
4501  {
4502  next = ptr->next;
4503 
4504  delete ptr;
4505 
4506  ptr = next;
4507  }
4508 
4509  nbthole_pairp = NULL;
4510 }
4511 /* END OF FUNCTION free_vdw_pair_list */
4512 
4513 /************************************************************************
4514  * FUNCTION free_table_pair_tree
4515  *
4516  * Free a table_pair_tree given a pointer to its head
4517  * **********************************************************************/
4518 
4519 void Parameters::free_table_pair_tree(IndexedTablePair *table_pair_ptr) {
4520  if (table_pair_ptr->left != NULL)
4521  {
4522  free_table_pair_tree(table_pair_ptr->left);
4523  }
4524 
4525  if (table_pair_ptr->right != NULL)
4526  {
4527  free_table_pair_tree(table_pair_ptr->right);
4528  }
4529 
4530  delete table_pair_ptr;
4531 
4532  return;
4533 }
4534 
4535 
4536 /************************************************************************/
4537 /* */
4538 /* FUNCTION free_vdw_pair_tree */
4539 /* */
4540 /* INPUTS: */
4541 /* vdw_pair_ptr - pointer to vdw_pair tree to free */
4542 /* */
4543 /* this is a recursive function that is used to free the memory */
4544 /* allocated for a vdw_pair paramter tree. It makes recursive calls */
4545 /* to free the left an right subtress, and then frees the head. It is*/
4546 /* only called by the destructor */
4547 /* */
4548 /************************************************************************/
4549 
4550 void Parameters::free_vdw_pair_tree(IndexedVdwPair *vdw_pair_ptr)
4551 
4552 {
4553  if (vdw_pair_ptr->left != NULL)
4554  {
4555  free_vdw_pair_tree(vdw_pair_ptr->left);
4556  }
4557 
4558  if (vdw_pair_ptr->right != NULL)
4559  {
4560  free_vdw_pair_tree(vdw_pair_ptr->right);
4561  }
4562 
4563  delete vdw_pair_ptr;
4564 
4565  return;
4566 }
4567 /* END OF FUNCTION free_vdw_pair_tree */
4568 
4569 /************************************************************************/
4570 /* */
4571 /* FUNCTION free_nbthole_pair_tree */
4572 /* */
4573 /* INPUTS: */
4574 /* nbthole_pair_ptr - pointer to nbthole_pair tree to free */
4575 /* */
4576 /* this is a recursive function that is used to free the memory */
4577 /* allocated for a nbthole_pair paramter tree. It makes recursive calls */
4578 /* to free the left an right subtress, and then frees the head. It is*/
4579 /* only called by the destructor */
4580 /* */
4581 /************************************************************************/
4582 
4583 void Parameters::free_nbthole_pair_tree(IndexedNbtholePair *nbthole_pair_ptr)
4584 
4585 {
4586  if (nbthole_pair_ptr->left != NULL)
4587  {
4588  free_nbthole_pair_tree(nbthole_pair_ptr->left);
4589  }
4590 
4591  if (nbthole_pair_ptr->right != NULL)
4592  {
4593  free_nbthole_pair_tree(nbthole_pair_ptr->right);
4594  }
4595 
4596  delete nbthole_pair_ptr;
4597 
4598  return;
4599 }
4600 /* END OF FUNCTION free_nbthole_pair_tree */
4601 
4602 /************************************************************************/
4603 /* */
4604 /* FUNCTION traverse_bond_params */
4605 /* */
4606 /* INPUTS: */
4607 /* tree - the bond binary tree to traverse */
4608 /* */
4609 /* This is a recursive call used for debugging purposes that */
4610 /* prints out all the bond paramters in the bond parameter binary */
4611 /* search tree. It is only called by print_bond_params */
4612 /* */
4613 /************************************************************************/
4614 
4615 void Parameters::traverse_bond_params(struct bond_params *tree)
4616 
4617 {
4618  if (tree==NULL)
4619  return;
4620 
4621  if (tree->left != NULL)
4622  {
4623  traverse_bond_params(tree->left);
4624  }
4625 
4626  DebugM(3,"BOND " << tree->atom1name << " " << tree->atom2name \
4627  << " index=" << tree->index << " k=" << tree->forceconstant \
4628  << " x0=" << tree->distance);
4629 
4630  if (tree->right != NULL)
4631  {
4632  traverse_bond_params(tree->right);
4633  }
4634 }
4635 /* END OF FUNCTION traverse_bond_params */
4636 
4637 /************************************************************************/
4638 /* */
4639 /* FUNCTION traverse_angle_params */
4640 /* */
4641 /* INPUTS: */
4642 /* tree - the angle binary tree to traverse */
4643 /* */
4644 /* This is a recursive call used for debugging purposes that */
4645 /* prints out all the angle paramters in the angle parameter binary */
4646 /* search tree. It is only called by print_angle_params */
4647 /* */
4648 /************************************************************************/
4649 
4650 void Parameters::traverse_angle_params(struct angle_params *tree)
4651 
4652 {
4653  if (tree==NULL)
4654  return;
4655 
4656  if (tree->left != NULL)
4657  {
4658  traverse_angle_params(tree->left);
4659  }
4660  DebugM(3,"ANGLE " << tree->atom1name << " " << tree->atom2name \
4661  << " " << tree->atom3name << " index=" << tree->index \
4662  << " k=" << tree->forceconstant << " theta0=" << tree->angle \
4663  );
4664 
4665  if (tree->right != NULL)
4666  {
4667  traverse_angle_params(tree->right);
4668  }
4669 }
4670 /* END OF FUNCTION traverse_angle_params */
4671 
4672 /************************************************************************/
4673 /* */
4674 /* FUNCTION traverse_dihedral_params */
4675 /* */
4676 /* INPUTS: */
4677 /* list - the dihedral linked list to traverse */
4678 /* */
4679 /* This is a call used for debugging purposes that prints out all */
4680 /* the bond paramters in the dihedral parameter linked list. It is */
4681 /* only called by print_dihedral_params. */
4682 /* */
4683 /************************************************************************/
4684 
4685 void Parameters::traverse_dihedral_params(struct dihedral_params *list)
4686 
4687 {
4688  int i;
4689 
4690  while (list != NULL)
4691  {
4692  DebugM(3,"DIHEDRAL " << list->atom1name << " " \
4693  << list->atom2name << " " << list->atom3name \
4694  << " " << list->atom4name << " index=" \
4695  << list->index \
4696  << " multiplicity=" << list->multiplicity << "\n");
4697 
4698  for (i=0; i<list->multiplicity; i++)
4699  {
4700  DebugM(3,"k=" << list->values[i].k \
4701  << " n=" << list->values[i].n \
4702  << " delta=" << list->values[i].delta);
4703  }
4704 
4705  list=list->next;
4706  }
4707 }
4708 /* END OF FUNCTION traverse_dihedral_params */
4709 
4710 /************************************************************************/
4711 /* */
4712 /* FUNCTION traverse_improper_params */
4713 /* */
4714 /* INPUTS: */
4715 /* list - the improper linked list to traverse */
4716 /* */
4717 /* This is a call used for debugging purposes that prints out all */
4718 /* the improper paramters in the improper parameter linked list. It is*/
4719 /* only called by print_improper_params. */
4720 /* */
4721 /************************************************************************/
4722 
4723 void Parameters::traverse_improper_params(struct improper_params *list)
4724 
4725 {
4726  int i;
4727 
4728  while (list != NULL)
4729  {
4730  DebugM(3,"Improper " << list->atom1name << " " \
4731  << list->atom2name << " " << list->atom3name \
4732  << " " << list->atom4name << " index=" \
4733  << list->index \
4734  << " multiplicity=" << list->multiplicity << "\n");
4735 
4736  for (i=0; i<list->multiplicity; i++)
4737  {
4738  DebugM(3,"k=" << list->values[i].k \
4739  << " n=" << list->values[i].n \
4740  << " delta=" << list->values[i].delta);
4741  }
4742 
4743  list=list->next;
4744  }
4745 }
4746 /* END OF FUNCTION traverse_improper_params */
4747 
4748 
4749 /************************************************************************/
4750 /* */
4751 /* FUNCTION traverse_vdw_params */
4752 /* */
4753 /* INPUTS: */
4754 /* tree - the vw binary tree to traverse */
4755 /* */
4756 /* This is a recursive call used for debugging purposes that */
4757 /* prints out all the vdw paramters in the vdw parameter binary */
4758 /* search tree. It is only called by print_vdw_params */
4759 /* */
4760 /************************************************************************/
4761 
4762 void Parameters::traverse_vdw_params(struct vdw_params *tree)
4763 
4764 {
4765  if (tree==NULL)
4766  return;
4767 
4768  if (tree->left != NULL)
4769  {
4770  traverse_vdw_params(tree->left);
4771  }
4772 
4773  DebugM(3,"vdW " << tree->atomname << " index=" << tree->index \
4774  << " sigma=" << tree->sigma << " epsilon=" << \
4775  tree->epsilon << " sigma 1-4=" << tree->sigma14 \
4776  << " epsilon 1-4=" << tree->epsilon14 << std::endl);
4777 
4778  if (tree->right != NULL)
4779  {
4780  traverse_vdw_params(tree->right);
4781  }
4782 }
4783 /* END OF FUNCTION traverse_vdw_params */
4784 
4785 
4786 /************************************************************************/
4787 /* */
4788 /* FUNCTION traverse_vdw_pair_params */
4789 /* */
4790 /* INPUTS: */
4791 /* list - the vdw_pair list to traverse */
4792 /* */
4793 /* This call simply prints out the vdw_pair list */
4794 /* */
4795 /************************************************************************/
4796 
4797 void Parameters::traverse_vdw_pair_params(struct vdw_pair_params *list)
4798 
4799 {
4800  if (list==NULL)
4801  return;
4802 
4803  DebugM(3,"vdW PAIR " << list->atom1name << " " \
4804  << list->atom2name << " A=" << list->A \
4805  << " B=" << list->B << " A 1-4=" \
4806  << list->A14 << " B 1-4=" << list->B14 \
4807  );
4808 
4809  traverse_vdw_pair_params(list->next);
4810 }
4811 /* END OF FUNCTION traverse_vdw_pair_params */
4812 
4813 /************************************************************************/
4814 /* */
4815 /* FUNCTION traverse_nbthole_pair_params */
4816 /* */
4817 /* INPUTS: */
4818 /* list - the nbthole_pair list to traverse */
4819 /* */
4820 /* This call simply prints out the nbthole_pair list */
4821 /* */
4822 /************************************************************************/
4823 
4824 void Parameters::traverse_nbthole_pair_params(struct nbthole_pair_params *list)
4825 
4826 {
4827  if (list==NULL)
4828  return;
4829 
4830  DebugM(3,"NBTHOLE PAIR " << list->atom1name << " " \
4831  << list->atom2name << " tholeij =" << list->tholeij \
4832  );
4833 
4834  traverse_nbthole_pair_params(list->next);
4835 }
4836 /* END OF FUNCTION traverse_nbthole_pair_params */
4837 
4838 /************************************************************************/
4839 /* */
4840 /* FUNCTION print_bond_params */
4841 /* */
4842 /* This is a debugging routine used to print out all the bond */
4843 /* parameters */
4844 /* */
4845 /************************************************************************/
4846 
4848 {
4849  DebugM(3,NumBondParams << " BOND PARAMETERS\n" \
4850  << "*****************************************" \
4851  );
4852 
4853  traverse_bond_params(bondp);
4854 }
4855 
4856 /************************************************************************/
4857 /* */
4858 /* FUNCTION print_angle_params */
4859 /* */
4860 /* This is a debugging routine used to print out all the angle */
4861 /* parameters */
4862 /* */
4863 /************************************************************************/
4864 
4866 {
4867  DebugM(3,NumAngleParams << " ANGLE PARAMETERS\n"
4868  << "*****************************************" );
4869  traverse_angle_params(anglep);
4870 }
4871 
4872 /************************************************************************/
4873 /* */
4874 /* FUNCTION print_dihedral_params */
4875 /* */
4876 /* This is a debugging routine used to print out all the dihedral */
4877 /* parameters */
4878 /* */
4879 /************************************************************************/
4880 
4882 {
4883  DebugM(3,NumDihedralParams << " DIHEDRAL PARAMETERS\n" \
4884  << "*****************************************" );
4885 
4886  traverse_dihedral_params(dihedralp);
4887 }
4888 
4889 /************************************************************************/
4890 /* */
4891 /* FUNCTION print_improper_params */
4892 /* */
4893 /* This is a debugging routine used to print out all the improper */
4894 /* parameters */
4895 /* */
4896 /************************************************************************/
4897 
4899 {
4900  DebugM(3,NumImproperParams << " IMPROPER PARAMETERS\n" \
4901  << "*****************************************" );
4902 
4903  traverse_improper_params(improperp);
4904 }
4905 
4906 /************************************************************************/
4907 /* */
4908 /* FUNCTION print_vdw_params */
4909 /* */
4910 /* This is a debugging routine used to print out all the vdw */
4911 /* parameters */
4912 /* */
4913 /************************************************************************/
4914 
4916 {
4917  DebugM(3,NumVdwParams << " vdW PARAMETERS\n" \
4918  << "*****************************************" );
4919 
4920  traverse_vdw_params(vdwp);
4921 }
4922 
4923 /************************************************************************/
4924 /* */
4925 /* FUNCTION print_vdw_pair_params */
4926 /* */
4927 /* This is a debugging routine used to print out all the vdw_pair */
4928 /* parameters */
4929 /* */
4930 /************************************************************************/
4931 
4933 {
4934  DebugM(3,NumVdwPairParams << " vdW PAIR PARAMETERS\n" \
4935  << "*****************************************" );
4936 
4937  traverse_vdw_pair_params(vdw_pairp);
4938 }
4939 
4940 /************************************************************************/
4941 /* */
4942 /* FUNCTION print_nbthole_pair_params */
4943 /* */
4944 /* This is a debugging routine used to print out all the nbthole_pair */
4945 /* parameters */
4946 /* */
4947 /************************************************************************/
4948 
4950 {
4951  DebugM(3,NumNbtholePairParams << " NBTHOLE PAIR PARAMETERS\n" \
4952  << "*****************************************" );
4953 
4954  traverse_nbthole_pair_params(nbthole_pairp);
4955 }
4956 
4957 /************************************************************************/
4958 /* */
4959 /* FUNCTION print_param_summary */
4960 /* */
4961 /* This function just prints out a brief summary of the paramters */
4962 /* that have been read in. It is intended for debugging purposes */
4963 /* */
4964 /************************************************************************/
4965 
4967 {
4968  iout << iINFO << "SUMMARY OF PARAMETERS:\n"
4969  << iINFO << NumBondParams << " BONDS\n"
4970  << iINFO << NumAngleParams << " ANGLES\n" << endi;
4971  if (cosAngles) {
4972  iout << iINFO << " " << (NumAngleParams - NumCosAngles) << " HARMONIC\n"
4973  << iINFO << " " << NumCosAngles << " COSINE-BASED\n" << endi;
4974  }
4975  iout << iINFO << NumDihedralParams << " DIHEDRAL\n"
4976  << iINFO << NumImproperParams << " IMPROPER\n"
4977  << iINFO << NumCrosstermParams << " CROSSTERM\n"
4978  << iINFO << NumVdwParams << " VDW\n"
4979  << iINFO << NumVdwPairParams << " VDW_PAIRS\n"
4980  << iINFO << NumNbtholePairParams << " NBTHOLE_PAIRS\n" << endi;
4981 }
4982 
4983 
4984 /************************************************************************/
4985 /* */
4986 /* FUNCTION done_reading_structure */
4987 /* */
4988 /* This function is used to tell the Parameters object that the */
4989 /* structure has been read in. This is so that the Parameters object */
4990 /* can now release the binary trees and linked lists that it was using */
4991 /* to search for parameters based on the atom type. From this point */
4992 /* on, only the arrays of parameter data will be used. If this object */
4993 /* resides on any node BUT the master node, it will never even have */
4994 /* these trees and lists. For the master node, this just frees up */
4995 /* some memory for better uses. */
4996 /* */
4997 /************************************************************************/
4998 
5000 
5001 {
5002  if (bondp != NULL)
5003  free_bond_tree(bondp);
5004 
5005  if (anglep != NULL)
5006  free_angle_tree(anglep);
5007 
5008  if (dihedralp != NULL)
5009  free_dihedral_list(dihedralp);
5010 
5011  if (improperp != NULL)
5012  free_improper_list(improperp);
5013 
5014  if (crosstermp != NULL)
5015  free_crossterm_list(crosstermp);
5016 
5017  if (vdwp != NULL)
5018  free_vdw_tree(vdwp);
5019 
5020  // Free the arrays used to track multiplicity for dihedrals
5021  // and impropers
5022  if (maxDihedralMults != NULL)
5023  delete [] maxDihedralMults;
5024 
5025  if (maxImproperMults != NULL)
5026  delete [] maxImproperMults;
5027 
5028  bondp=NULL;
5029  anglep=NULL;
5030  dihedralp=NULL;
5031  improperp=NULL;
5032  crosstermp=NULL;
5033  vdwp=NULL;
5034  maxImproperMults=NULL;
5035  maxDihedralMults=NULL;
5036 }
5037 /* END OF FUNCTION done_reading_structure */
5038 
5039 /************************************************************************/
5040 /* */
5041 /* FUNCTION send_Parameters */
5042 /* */
5043 /* This function is used by the master node to broadcast the */
5044 /* structure Parameters to all the other nodes. */
5045 /* */
5046 /************************************************************************/
5047 
5049 {
5050  Real *a1, *a2, *a3, *a4; // Temporary arrays for sending messages
5051  int *i1, *i2, *i3; // Temporary int array
5052  int i, j; // Loop counters
5053  Real **kvals; // Force constant values for dihedrals and impropers
5054  int **nvals; // Periodicity values for dihedrals and impropers
5055  Real **deltavals; // Phase shift values for dihedrals and impropers
5056  BigReal *pairC6, *pairC12; // JLai
5057  /*MOStream *msg=comm_obj->newOutputStream(ALLBUTME, STATICPARAMSTAG, BUFSIZE);
5058  if ( msg == NULL )
5059  {
5060  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5061  }*/
5062 
5063  // Send the bond parameters
5064  msg->put(NumBondParams);
5065 
5066  if (NumBondParams)
5067  {
5068  a1 = new Real[NumBondParams];
5069  a2 = new Real[NumBondParams];
5070 
5071  if ( (a1 == NULL) || (a2 == NULL) )
5072  {
5073  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5074  }
5075 
5076  for (i=0; i<NumBondParams; i++)
5077  {
5078  a1[i] = bond_array[i].k;
5079  a2[i] = bond_array[i].x0;
5080  }
5081 
5082  msg->put(NumBondParams, a1)->put(NumBondParams, a2);
5083 
5084  delete [] a1;
5085  delete [] a2;
5086  }
5087 
5088  // Send the angle parameters
5089  msg->put(NumAngleParams);
5090 
5091  if (NumAngleParams)
5092  {
5093  a1 = new Real[NumAngleParams];
5094  a2 = new Real[NumAngleParams];
5095  a3 = new Real[NumAngleParams];
5096  a4 = new Real[NumAngleParams];
5097  i1 = new int[NumAngleParams];
5098 
5099  if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) ||
5100  (a4 == NULL) || (i1 == NULL))
5101  {
5102  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5103  }
5104 
5105  for (i=0; i<NumAngleParams; i++)
5106  {
5107  a1[i] = angle_array[i].k;
5108  a2[i] = angle_array[i].theta0;
5109  a3[i] = angle_array[i].k_ub;
5110  a4[i] = angle_array[i].r_ub;
5111  i1[i] = angle_array[i].normal;
5112  }
5113 
5114  msg->put(NumAngleParams, a1)->put(NumAngleParams, a2);
5115  msg->put(NumAngleParams, a3)->put(NumAngleParams, a4);
5116  msg->put(NumAngleParams, i1);
5117 
5118  delete [] a1;
5119  delete [] a2;
5120  delete [] a3;
5121  delete [] a4;
5122  delete [] i1;
5123  }
5124 
5125  // Send the dihedral parameters
5126  msg->put(NumDihedralParams);
5127 
5128  if (NumDihedralParams)
5129  {
5130  i1 = new int[NumDihedralParams];
5131  kvals = new Real *[MAX_MULTIPLICITY];
5132  nvals = new int *[MAX_MULTIPLICITY];
5133  deltavals = new Real *[MAX_MULTIPLICITY];
5134 
5135  if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5136  (deltavals == NULL) )
5137  {
5138  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5139  }
5140 
5141  for (i=0; i<MAX_MULTIPLICITY; i++)
5142  {
5143  kvals[i] = new Real[NumDihedralParams];
5144  nvals[i] = new int[NumDihedralParams];
5145  deltavals[i] = new Real[NumDihedralParams];
5146 
5147  if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5148  {
5149  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5150  }
5151  }
5152 
5153  for (i=0; i<NumDihedralParams; i++)
5154  {
5155  i1[i] = dihedral_array[i].multiplicity;
5156 
5157  for (j=0; j<MAX_MULTIPLICITY; j++)
5158  {
5159  kvals[j][i] = dihedral_array[i].values[j].k;
5160  nvals[j][i] = dihedral_array[i].values[j].n;
5161  deltavals[j][i] = dihedral_array[i].values[j].delta;
5162  }
5163  }
5164 
5165  msg->put(NumDihedralParams, i1);
5166 
5167  for (i=0; i<MAX_MULTIPLICITY; i++)
5168  {
5169  msg->put(NumDihedralParams, kvals[i]);
5170  msg->put(NumDihedralParams, nvals[i]);
5171  msg->put(NumDihedralParams, deltavals[i]);
5172 
5173  delete [] kvals[i];
5174  delete [] nvals[i];
5175  delete [] deltavals[i];
5176  }
5177 
5178  delete [] i1;
5179  delete [] kvals;
5180  delete [] nvals;
5181  delete [] deltavals;
5182  }
5183 
5184  // Send the improper parameters
5185  msg->put(NumImproperParams);
5186 
5187  if (NumImproperParams)
5188  {
5189  i1 = new int[NumImproperParams];
5190  kvals = new Real *[MAX_MULTIPLICITY];
5191  nvals = new int *[MAX_MULTIPLICITY];
5192  deltavals = new Real *[MAX_MULTIPLICITY];
5193 
5194  if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5195  (deltavals == NULL) )
5196  {
5197  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5198  }
5199 
5200  for (i=0; i<MAX_MULTIPLICITY; i++)
5201  {
5202  kvals[i] = new Real[NumImproperParams];
5203  nvals[i] = new int[NumImproperParams];
5204  deltavals[i] = new Real[NumImproperParams];
5205 
5206  if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5207  {
5208  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5209  }
5210  }
5211 
5212  for (i=0; i<NumImproperParams; i++)
5213  {
5214  i1[i] = improper_array[i].multiplicity;
5215 
5216  for (j=0; j<MAX_MULTIPLICITY; j++)
5217  {
5218  kvals[j][i] = improper_array[i].values[j].k;
5219  nvals[j][i] = improper_array[i].values[j].n;
5220  deltavals[j][i] = improper_array[i].values[j].delta;
5221  }
5222  }
5223 
5224  msg->put(NumImproperParams, i1);
5225 
5226  for (i=0; i<MAX_MULTIPLICITY; i++)
5227  {
5228  msg->put(NumImproperParams, kvals[i]);
5229  msg->put(NumImproperParams, nvals[i]);
5230  msg->put(NumImproperParams, deltavals[i]);
5231 
5232  delete [] kvals[i];
5233  delete [] nvals[i];
5234  delete [] deltavals[i];
5235  }
5236 
5237  delete [] i1;
5238  delete [] kvals;
5239  delete [] nvals;
5240  delete [] deltavals;
5241  }
5242 
5243  // Send the crossterm parameters
5244  msg->put(NumCrosstermParams);
5245 
5246  if (NumCrosstermParams)
5247  {
5248  for (i=0; i<NumCrosstermParams; ++i) {
5249  int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
5250  msg->put(nvals,&crossterm_array[i].c[0][0].d00);
5251  }
5252  }
5253  // Send the GromacsPairs parameters
5254  // JLai
5255  msg->put(NumGromacsPairParams);
5256 
5258  {
5259  pairC6 = new BigReal[NumGromacsPairParams];
5260  pairC12 = new BigReal[NumGromacsPairParams];
5261  if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
5262  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5263  }
5264 
5265  for (i=0; i<NumGromacsPairParams; i++) {
5266  pairC6[i] = gromacsPair_array[i].pairC6;
5267  pairC12[i] = gromacsPair_array[i].pairC12;
5268  }
5269 
5270  msg->put(NumGromacsPairParams,pairC6);
5271  msg->put(NumGromacsPairParams,pairC12);
5272 
5273  delete [] pairC6;
5274  delete [] pairC12;
5275  }
5276  // End of JLai
5277 
5278  //
5279  //Send the energy table parameters
5280  msg->put(numenerentries);
5281 
5282  if (numenerentries) {
5283  /*
5284  b1 = new Real[numenerentries];
5285  if (b1 == NULL)
5286  {
5287  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5288  }
5289  */
5290 
5291  msg->put(numenerentries, table_ener);
5292  }
5293 
5294  // Send the vdw parameters
5295  msg->put(NumVdwParams);
5296  msg->put(NumVdwParamsAssigned);
5297 
5298  if (NumVdwParams)
5299  {
5300  a1 = new Real[NumVdwParams];
5301  a2 = new Real[NumVdwParams];
5302  a3 = new Real[NumVdwParams];
5303  a4 = new Real[NumVdwParams];
5304 
5305  if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) )
5306  {
5307  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5308  }
5309 
5310  for (i=0; i<NumVdwParams; i++)
5311  {
5312  a1[i] = vdw_array[i].sigma;
5313  a2[i] = vdw_array[i].epsilon;
5314  a3[i] = vdw_array[i].sigma14;
5315  a4[i] = vdw_array[i].epsilon14;
5316  }
5317 
5318  msg->put(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
5319  msg->put(NumVdwParams, a1);
5320  msg->put(NumVdwParams, a2);
5321  msg->put(NumVdwParams, a3);
5322  msg->put(NumVdwParams, a4);
5323 
5324  delete [] a1;
5325  delete [] a2;
5326  delete [] a3;
5327  delete [] a4;
5328  }
5329 
5330  // Send the vdw pair parameters
5331  msg->put(NumVdwPairParams);
5332 
5333  if (NumVdwPairParams)
5334  {
5335  a1 = new Real[NumVdwPairParams];
5336  a2 = new Real[NumVdwPairParams];
5337  a3 = new Real[NumVdwPairParams];
5338  a4 = new Real[NumVdwPairParams];
5339  i1 = new int[NumVdwPairParams];
5340  i2 = new int[NumVdwPairParams];
5341 
5342  if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
5343  (i1 == NULL) || (i2 == NULL) )
5344  {
5345  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5346  }
5347 
5348  vdw_pair_to_arrays(i1, i2, a1, a2, a3, a4, 0, vdw_pair_tree);
5349 
5350  msg->put(NumVdwPairParams, i1)->put(NumVdwPairParams, i2);
5351  msg->put(NumVdwPairParams, a1);
5352  msg->put(NumVdwPairParams, a2)->put(NumVdwPairParams, a3);
5353  msg->put(NumVdwPairParams, a4);
5354  }
5355 
5356  // Send the nbthole pair parameters
5357  msg->put(NumNbtholePairParams);
5358 
5360  {
5361  a1 = new Real[NumNbtholePairParams];
5362  a2 = new Real[NumNbtholePairParams];
5363  a3 = new Real[NumNbtholePairParams];
5364  i1 = new int[NumNbtholePairParams];
5365  i2 = new int[NumNbtholePairParams];
5366 
5367  if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (i1 == NULL) || (i2 == NULL) )
5368  {
5369  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5370  }
5371 
5372  nbthole_pair_to_arrays(i1, i2, a1, a2, a3, 0, nbthole_pair_tree);
5373 
5374  for (i=0; i<NumNbtholePairParams; i++)
5375  {
5376  nbthole_array[i].ind1 = i1[i];
5377  nbthole_array[i].ind2 = i2[i];
5378  nbthole_array[i].alphai = a1[i];
5379  nbthole_array[i].alphaj = a2[i];
5380  nbthole_array[i].tholeij = a3[i];
5381  }
5382 
5383  msg->put(NumNbtholePairParams, i1)->put(NumNbtholePairParams, i2);
5384  msg->put(NumNbtholePairParams, a1);
5385  msg->put(NumNbtholePairParams, a2)->put(NumNbtholePairParams, a3);
5386  }
5387 
5388  // Send the table pair parameters
5389  //printf("Pairs: %i\n", NumTablePairParams);
5390  msg->put(NumTablePairParams);
5391 
5392  if (NumTablePairParams)
5393  {
5394  i1 = new int[NumTablePairParams];
5395  i2 = new int[NumTablePairParams];
5396  i3 = new int[NumTablePairParams];
5397 
5398  if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
5399  {
5400  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5401  }
5402 
5403  table_pair_to_arrays(i1, i2, i3, 0, tab_pair_tree);
5404 
5405  msg->put(NumTablePairParams, i1)->put(NumTablePairParams, i2);
5406  msg->put(NumTablePairParams, i3);
5407  }
5408 
5409  // send the hydrogen bond parameters
5410  // hbondParams.create_message(msg);
5411  msg->end();
5412  delete msg;
5413 }
5414 
5415 /************************************************************************/
5416 /* */
5417 /* FUNCTION receive_Parameters */
5418 /* */
5419 /* This function is used by all the client processes to receive */
5420 /* the structure parameters from the master node. */
5421 /* */
5422 /************************************************************************/
5423 
5425 
5426 {
5427  int i, j; // Loop counters
5428  Real *a1, *a2, *a3, *a4; // Temporary arrays to get data from message in
5429  int *i1, *i2, *i3; // Temporary int array to get data from message in
5430  IndexedVdwPair *new_node; // New node for vdw pair param tree
5431  IndexedTablePair *new_tab_node;
5432  Real **kvals; // Force constant values for dihedrals and impropers
5433  int **nvals; // Periodicity values for dihedrals and impropers
5434  Real **deltavals; // Phase shift values for dihedrals and impropers
5435  BigReal *pairC6, *pairC12; // JLai
5436 
5437  // Get the bonded parameters
5438  msg->get(NumBondParams);
5439 
5440  if (NumBondParams)
5441  {
5443  a1 = new Real[NumBondParams];
5444  a2 = new Real[NumBondParams];
5445 
5446  if ( (bond_array == NULL) || (a1 == NULL) || (a2 == NULL) )
5447  {
5448  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5449  }
5450 
5451  msg->get(NumBondParams, a1);
5452  msg->get(NumBondParams, a2);
5453 
5454  for (i=0; i<NumBondParams; i++)
5455  {
5456  bond_array[i].k = a1[i];
5457  bond_array[i].x0 = a2[i];
5458  }
5459 
5460  delete [] a1;
5461  delete [] a2;
5462  }
5463 
5464  // Get the angle parameters
5465  msg->get(NumAngleParams);
5466 
5467  if (NumAngleParams)
5468  {
5470  a1 = new Real[NumAngleParams];
5471  a2 = new Real[NumAngleParams];
5472  a3 = new Real[NumAngleParams];
5473  a4 = new Real[NumAngleParams];
5474  i1 = new int[NumAngleParams];
5475 
5476  if ( (angle_array == NULL) || (a1 == NULL) || (a2 == NULL) ||
5477  (a3 == NULL) || (a4 == NULL) || (i1 == NULL))
5478  {
5479  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5480  }
5481 
5482  msg->get(NumAngleParams, a1);
5483  msg->get(NumAngleParams, a2);
5484  msg->get(NumAngleParams, a3);
5485  msg->get(NumAngleParams, a4);
5486  msg->get(NumAngleParams, i1);
5487 
5488  for (i=0; i<NumAngleParams; i++)
5489  {
5490  angle_array[i].k = a1[i];
5491  angle_array[i].theta0 = a2[i];
5492  angle_array[i].k_ub = a3[i];
5493  angle_array[i].r_ub = a4[i];
5494  angle_array[i].normal = i1[i];
5495  }
5496 
5497  delete [] a1;
5498  delete [] a2;
5499  delete [] a3;
5500  delete [] a4;
5501  delete [] i1;
5502  }
5503 
5504  // Get the dihedral parameters
5505  msg->get(NumDihedralParams);
5506 
5507  if (NumDihedralParams)
5508  {
5510 
5511  i1 = new int[NumDihedralParams];
5512  kvals = new Real *[MAX_MULTIPLICITY];
5513  nvals = new int *[MAX_MULTIPLICITY];
5514  deltavals = new Real *[MAX_MULTIPLICITY];
5515 
5516  if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5517  (deltavals == NULL) || (dihedral_array == NULL) )
5518  {
5519  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5520  }
5521 
5522  for (i=0; i<MAX_MULTIPLICITY; i++)
5523  {
5524  kvals[i] = new Real[NumDihedralParams];
5525  nvals[i] = new int[NumDihedralParams];
5526  deltavals[i] = new Real[NumDihedralParams];
5527 
5528  if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5529  {
5530  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5531  }
5532  }
5533 
5534  msg->get(NumDihedralParams, i1);
5535 
5536  for (i=0; i<MAX_MULTIPLICITY; i++)
5537  {
5538  msg->get(NumDihedralParams, kvals[i]);
5539  msg->get(NumDihedralParams, nvals[i]);
5540  msg->get(NumDihedralParams, deltavals[i]);
5541  }
5542 
5543  for (i=0; i<NumDihedralParams; i++)
5544  {
5545  dihedral_array[i].multiplicity = i1[i];
5546 
5547  for (j=0; j<MAX_MULTIPLICITY; j++)
5548  {
5549  dihedral_array[i].values[j].k = kvals[j][i];
5550  dihedral_array[i].values[j].n = nvals[j][i];
5551  dihedral_array[i].values[j].delta = deltavals[j][i];
5552  }
5553  }
5554 
5555  for (i=0; i<MAX_MULTIPLICITY; i++)
5556  {
5557  delete [] kvals[i];
5558  delete [] nvals[i];
5559  delete [] deltavals[i];
5560  }
5561 
5562  delete [] i1;
5563  delete [] kvals;
5564  delete [] nvals;
5565  delete [] deltavals;
5566  }
5567 
5568  // Get the improper parameters
5569  msg->get(NumImproperParams);
5570 
5571  if (NumImproperParams)
5572  {
5574  i1 = new int[NumImproperParams];
5575  kvals = new Real *[MAX_MULTIPLICITY];
5576  nvals = new int *[MAX_MULTIPLICITY];
5577  deltavals = new Real *[MAX_MULTIPLICITY];
5578 
5579  if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5580  (deltavals == NULL) || (improper_array==NULL) )
5581  {
5582  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5583  }
5584 
5585  for (i=0; i<MAX_MULTIPLICITY; i++)
5586  {
5587  kvals[i] = new Real[NumImproperParams];
5588  nvals[i] = new int[NumImproperParams];
5589  deltavals[i] = new Real[NumImproperParams];
5590 
5591  if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5592  {
5593  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5594  }
5595  }
5596 
5597  msg->get(NumImproperParams,i1);
5598 
5599  for (i=0; i<MAX_MULTIPLICITY; i++)
5600  {
5601  msg->get(NumImproperParams,kvals[i]);
5602  msg->get(NumImproperParams,nvals[i]);
5603  msg->get(NumImproperParams,deltavals[i]);
5604  }
5605 
5606  for (i=0; i<NumImproperParams; i++)
5607  {
5608  improper_array[i].multiplicity = i1[i];
5609 
5610  for (j=0; j<MAX_MULTIPLICITY; j++)
5611  {
5612  improper_array[i].values[j].k = kvals[j][i];
5613  improper_array[i].values[j].n = nvals[j][i];
5614  improper_array[i].values[j].delta = deltavals[j][i];
5615  }
5616  }
5617 
5618  for (i=0; i<MAX_MULTIPLICITY; i++)
5619  {
5620  delete [] kvals[i];
5621  delete [] nvals[i];
5622  delete [] deltavals[i];
5623  }
5624 
5625  delete [] i1;
5626  delete [] kvals;
5627  delete [] nvals;
5628  delete [] deltavals;
5629  }
5630 
5631  // Get the crossterm parameters
5632  msg->get(NumCrosstermParams);
5633 
5634  if (NumCrosstermParams)
5635  {
5637 
5638  for (i=0; i<NumCrosstermParams; ++i) {
5639  int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
5640  msg->get(nvals,&crossterm_array[i].c[0][0].d00);
5641  }
5642  }
5643 
5644  // Get GromacsPairs parameters
5645  // JLai
5646  msg->get(NumGromacsPairParams);
5647 
5649  {
5651  pairC6 = new BigReal[NumGromacsPairParams];
5652  pairC12 = new BigReal[NumGromacsPairParams];
5653 
5654  if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
5655  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5656  }
5657 
5658  msg->get(NumGromacsPairParams,pairC6);
5659  msg->get(NumGromacsPairParams,pairC12);
5660 
5661  for (i=0; i<NumGromacsPairParams; ++i) {
5662  gromacsPair_array[i].pairC6 = pairC6[i];
5663  gromacsPair_array[i].pairC12 = pairC12[i];
5664  }
5665 
5666  delete [] pairC6;
5667  delete [] pairC12;
5668  }
5669  // JLai
5670 
5671  //Get the energy table
5672  msg->get(numenerentries);
5673  if (numenerentries > 0) {
5674  //fprintf(stdout, "Getting tables\n");
5675  //fprintf(infofile, "Trying to allocate table\n");
5677  //fprintf(infofile, "Finished table allocation\n");
5678  if (table_ener==NULL)
5679  {
5680  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5681  }
5682 
5683  msg->get(numenerentries, table_ener);
5684  }
5685 
5686  // Get the vdw parameters
5687  msg->get(NumVdwParams);
5688  msg->get(NumVdwParamsAssigned);
5689 
5690  if (NumVdwParams)
5691  {
5692  atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
5694  a1 = new Real[NumVdwParams];
5695  a2 = new Real[NumVdwParams];
5696  a3 = new Real[NumVdwParams];
5697  a4 = new Real[NumVdwParams];
5698 
5699  if ( (vdw_array==NULL) || (a1==NULL) || (a2==NULL) || (a3==NULL)
5700  || (a4==NULL) || (atomTypeNames==NULL))
5701  {
5702  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5703  }
5704 
5705  msg->get(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
5706  msg->get(NumVdwParams, a1);
5707  msg->get(NumVdwParams, a2);
5708  msg->get(NumVdwParams, a3);
5709  msg->get(NumVdwParams, a4);
5710 
5711  for (i=0; i<NumVdwParams; i++)
5712  {
5713  vdw_array[i].sigma = a1[i];
5714  vdw_array[i].epsilon = a2[i];
5715  vdw_array[i].sigma14 = a3[i];
5716  vdw_array[i].epsilon14 = a4[i];
5717  }
5718 
5719  delete [] a1;
5720  delete [] a2;
5721  delete [] a3;
5722  delete [] a4;
5723  }
5724 
5725  // Get the vdw_pair_parameters
5726  msg->get(NumVdwPairParams);
5727 
5728  if (NumVdwPairParams)
5729  {
5730  a1 = new Real[NumVdwPairParams];
5731  a2 = new Real[NumVdwPairParams];
5732  a3 = new Real[NumVdwPairParams];
5733  a4 = new Real[NumVdwPairParams];
5734  i1 = new int[NumVdwPairParams];
5735  i2 = new int[NumVdwPairParams];
5736 
5737  if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
5738  (i1 == NULL) || (i2 == NULL) )
5739  {
5740  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5741  }
5742 
5743  msg->get(NumVdwPairParams, i1);
5744  msg->get(NumVdwPairParams, i2);
5745  msg->get(NumVdwPairParams, a1);
5746  msg->get(NumVdwPairParams, a2);
5747  msg->get(NumVdwPairParams, a3);
5748  msg->get(NumVdwPairParams, a4);
5749 
5750  for (i=0; i<NumVdwPairParams; i++)
5751  {
5752  new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
5753 
5754  if (new_node == NULL)
5755  {
5756  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5757  }
5758 
5759  new_node->ind1 = i1[i];
5760  new_node->ind2 = i2[i];
5761  new_node->A = a1[i];
5762  new_node->A14 = a2[i];
5763  new_node->B = a3[i];
5764  new_node->B14 = a4[i];
5765  new_node->left = NULL;
5766  new_node->right = NULL;
5767 
5768  vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
5769  }
5770 
5771  delete [] i1;
5772  delete [] i2;
5773  delete [] a1;
5774  delete [] a2;
5775  delete [] a3;
5776  delete [] a4;
5777  }
5778 
5779  // Get the nbthole_pair_parameters
5780  msg->get(NumNbtholePairParams);
5781 
5783  {
5785  a1 = new Real[NumNbtholePairParams];
5786  a2 = new Real[NumNbtholePairParams];
5787  a3 = new Real[NumNbtholePairParams];
5788  i1 = new int[NumNbtholePairParams];
5789  i2 = new int[NumNbtholePairParams];
5790 
5791  if ( (nbthole_array == NULL) || (a1 == NULL) || (a2 == NULL) || (a3 == NULL)
5792  || (i1 == NULL) || (i2 == NULL) )
5793  {
5794  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5795  }
5796 
5797  msg->get(NumNbtholePairParams, i1);
5798  msg->get(NumNbtholePairParams, i2);
5799  msg->get(NumNbtholePairParams, a1);
5800  msg->get(NumNbtholePairParams, a2);
5801  msg->get(NumNbtholePairParams, a3);
5802 
5803  for (i=0; i<NumNbtholePairParams; i++)
5804  {
5805 
5806  nbthole_array[i].ind1 = i1[i];
5807  nbthole_array[i].ind2 = i2[i];
5808  nbthole_array[i].alphai = a1[i];
5809  nbthole_array[i].alphaj = a2[i];
5810  nbthole_array[i].tholeij = a3[i];
5811 
5812  }
5813 
5814  delete [] i1;
5815  delete [] i2;
5816  delete [] a1;
5817  delete [] a2;
5818  delete [] a3;
5819  }
5820 
5821  // Get the table_pair_parameters
5822  msg->get(NumTablePairParams);
5823 
5824  if (NumTablePairParams)
5825  {
5826  i1 = new int[NumTablePairParams];
5827  i2 = new int[NumTablePairParams];
5828  i3 = new int[NumTablePairParams];
5829 
5830  if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
5831  {
5832  NAMD_die("memory allocation failed in Parameters::send_Parameters");
5833  }
5834 
5835  msg->get(NumTablePairParams, i1);
5836  msg->get(NumTablePairParams, i2);
5837  msg->get(NumTablePairParams, i3);
5838 
5839  for (i=0; i<NumTablePairParams; i++)
5840  {
5841  new_tab_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
5842 
5843  if (new_tab_node == NULL)
5844  {
5845  NAMD_die("memory allocation failed in Parameters::receive_Parameters");
5846  }
5847 
5848 // printf("Adding new table pair with ind1 %i ind2 %i k %i\n", i1[i], i2[i],i3[i]);
5849  new_tab_node->ind1 = i1[i];
5850  new_tab_node->ind2 = i2[i];
5851  new_tab_node->K = i3[i];
5852  new_tab_node->left = NULL;
5853  new_tab_node->right = NULL;
5854 
5855  tab_pair_tree = add_to_indexed_table_pairs(new_tab_node, tab_pair_tree);
5856  }
5857 
5858  delete [] i1;
5859  delete [] i2;
5860  delete [] i3;
5861  }
5862 
5863  // receive the hydrogen bond parameters
5864  // hbondParams.receive_message(msg);
5865 
5866  AllFilesRead = TRUE;
5867 
5868  delete msg;
5869 }
5870 /* END OF FUNCTION receive_Parameters */
5871 
5872 /************************************************************************/
5873 /* */
5874 /* FUNCTION convert_vdw_pairs */
5875 /* */
5876 /* This function converts the linked list of vdw_pairs indexed by */
5877 /* atom name into a binary search tree of parameters stored by vdw */
5878 /* type index. This tree is what will be used for real when searching */
5879 /* for parameters during the simulation. */
5880 /* */
5881 /************************************************************************/
5882 
5883 void Parameters::convert_vdw_pairs()
5884 
5885 {
5886  #ifdef MEM_OPT_VERSION
5887  AtomCstInfo atom_struct;
5888  #else
5889  Atom atom_struct; // Dummy structure for getting indexes
5890  #endif
5891  Index index1, index2; // Indexes for the two atoms
5892  IndexedVdwPair *new_node; // New node for tree
5893  struct vdw_pair_params *ptr, *next; // Pointers for traversing list
5894 
5895  ptr = vdw_pairp;
5896 
5897  // Go down then entire list and insert each node into the
5898  // binary search tree
5899  while (ptr != NULL)
5900  {
5901  new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
5902 
5903  if (new_node == NULL)
5904  {
5905  NAMD_die("memory allocation failed in Parameters::convert_vdw_pairs");
5906  }
5907 
5908  // Get the vdw indexes for the two atoms. This is kind of a hack
5909  // using the goofy Atom structure, but hey, it works
5910  assign_vdw_index(ptr->atom1name, &atom_struct);
5911  index1 = atom_struct.vdw_type;
5912  assign_vdw_index(ptr->atom2name, &atom_struct);
5913  index2 = atom_struct.vdw_type;
5914 
5915  if (index1 > index2)
5916  {
5917  new_node->ind1 = index2;
5918  new_node->ind2 = index1;
5919  }
5920  else
5921  {
5922  new_node->ind1 = index1;
5923  new_node->ind2 = index2;
5924  }
5925 
5926  new_node->A = ptr->A;
5927  new_node->B = ptr->B;
5928  new_node->A14 = ptr->A14;
5929  new_node->B14 = ptr->B14;
5930 
5931  new_node->left = NULL;
5932  new_node->right = NULL;
5933 
5934  // Add it to the tree
5935  vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
5936 
5937  // Free the current node and move to the next
5938  next = ptr->next;
5939 
5940  delete ptr;
5941 
5942  ptr = next;
5943  }
5944 
5945  vdw_pairp = NULL;
5946 
5947 }
5948 /* END OF FUNCTION convert_vdw_pairs */
5949 
5950 /************************************************************************/
5951 /* */
5952 /* FUNCTION convert_nbthole_pairs */
5953 /* */
5954 /* This function converts the linked list of nbthole_pairs indexed by */
5955 /* atom name into a binary search tree of parameters stored by vdw */
5956 /* type index. This tree is what will be used for real when searching */
5957 /* for parameters during the simulation. */
5958 /* */
5959 /************************************************************************/
5960 
5961 void Parameters::convert_nbthole_pairs()
5962 
5963 {
5964  #ifdef MEM_OPT_VERSION
5965  AtomCstInfo atom_struct;
5966  #else
5967  Atom atom_struct; // Dummy structure for getting indexes
5968  #endif
5969  Index index1, index2; // Indexes for the two atoms
5970  IndexedNbtholePair *new_node; // New node for tree
5971  struct nbthole_pair_params *ptr, *next; // Pointers for traversing list
5972 
5973  ptr = nbthole_pairp;
5974 
5975  // Go down then entire list and insert each node into the
5976  // binary search tree
5977  while (ptr != NULL)
5978  {
5979  new_node = (IndexedNbtholePair *) malloc(sizeof(IndexedNbtholePair));
5980 
5981  if (new_node == NULL)
5982  {
5983  NAMD_die("memory allocation failed in Parameters::convert_nbthole_pairs");
5984  }
5985 
5986  // Get the vdw indexes for the two atoms. This is kind of a hack
5987  // using the goofy Atom structure, but hey, it works
5988  assign_vdw_index(ptr->atom1name, &atom_struct);
5989  index1 = atom_struct.vdw_type;
5990  assign_vdw_index(ptr->atom2name, &atom_struct);
5991  index2 = atom_struct.vdw_type;
5992 
5993  if (index1 > index2)
5994  {
5995  new_node->ind1 = index2;
5996  new_node->ind2 = index1;
5997  }
5998  else
5999  {
6000  new_node->ind1 = index1;
6001  new_node->ind2 = index2;
6002  }
6003 
6004  new_node->alphai = ptr->alphai;
6005  new_node->alphaj = ptr->alphaj;
6006  new_node->tholeij = ptr->tholeij;
6007 
6008  new_node->left = NULL;
6009  new_node->right = NULL;
6010 
6011  // Add it to the tree
6012  nbthole_pair_tree = add_to_indexed_nbthole_pairs(new_node, nbthole_pair_tree);
6013 
6014  // Free the current node and move to the next
6015  next = ptr->next;
6016 
6017  delete ptr;
6018 
6019  ptr = next;
6020  }
6021 
6022  nbthole_pairp = NULL;
6023 
6024 }
6025 /* END OF FUNCTION convert_nbthole_pairs */
6026 
6027 /************************************************************************/
6028 /* */
6029 /* FUNCTION convert_table_pairs */
6030 /* */
6031 /* This function converts the linked list of table_pairs indexed by */
6032 /* atom name into a binary search tree of parameters stored by table */
6033 /* type index. This tree is what will be used for real when searching */
6034 /* for parameters during the simulation. */
6035 /* */
6036 /************************************************************************/
6037 
6038 void Parameters::convert_table_pairs()
6039 
6040 {
6041  #ifdef MEM_OPT_VERSION
6042  AtomCstInfo atom_struct;
6043  #else
6044  Atom atom_struct; // Dummy structure for getting indexes
6045  #endif
6046  Index index1, index2; // Indexes for the two atoms
6047  IndexedTablePair *new_node; // New node for tree
6048  struct table_pair_params *ptr, *next; // Pointers for traversing list
6049 
6050  ptr = table_pairp;
6051 
6052  // Go down then entire list and insert each node into the
6053  // binary search tree
6054  while (ptr != NULL)
6055  {
6056  new_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
6057 
6058  if (new_node == NULL)
6059  {
6060  NAMD_die("memory allocation failed in Parameters::convert_table_pairs");
6061  }
6062 
6063  // Get the vdw indexes for the two atoms. This is kind of a hack
6064  // using the goofy Atom structure, but hey, it works
6065  assign_vdw_index(ptr->atom1name, &atom_struct);
6066  index1 = atom_struct.vdw_type;
6067  assign_vdw_index(ptr->atom2name, &atom_struct);
6068  index2 = atom_struct.vdw_type;
6069 
6070 // printf("Adding new table pair with index1 %i, index2 %i, k %i\n", index1, index2, ptr->K);
6071 
6072  if (index1 > index2)
6073  {
6074  new_node->ind1 = index2;
6075  new_node->ind2 = index1;
6076  }
6077  else
6078  {
6079  new_node->ind1 = index1;
6080  new_node->ind2 = index2;
6081  }
6082 
6083  new_node->K = ptr->K;
6084 
6085  new_node->left = NULL;
6086  new_node->right = NULL;
6087 
6088  // Add it to the tree
6089  tab_pair_tree = add_to_indexed_table_pairs(new_node, tab_pair_tree);
6090 
6091  // Free the current node and move to the next
6092  next = ptr->next;
6093 
6094  delete ptr;
6095 
6096  ptr = next;
6097  }
6098 
6099  table_pairp = NULL;
6100 
6101 }
6102 /* END OF FUNCTION convert_table_pairs */
6103 
6104 /************************************************************************/
6105 /* */
6106 /* FUNCTION add_to_indexed_table_pairs */
6107 /* */
6108 /* INPUTS: */
6109 /* new_node - new node to be added to the tree */
6110 /* tree - tree to add the node to */
6111 /* */
6112 /* This is a recursive function that adds a node to the */
6113 /* binary search tree of table_pair parameters */
6114 /* */
6115 /************************************************************************/
6116 
6117 IndexedTablePair *Parameters::add_to_indexed_table_pairs(IndexedTablePair *new_node,
6118  IndexedTablePair *tree)
6119 
6120 {
6121  if (tree == NULL)
6122  return(new_node);
6123 
6124  if ( (new_node->ind1 < tree->ind1) ||
6125  ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
6126  {
6127  tree->left = add_to_indexed_table_pairs(new_node, tree->left);
6128  }
6129  else
6130  {
6131  tree->right = add_to_indexed_table_pairs(new_node, tree->right);
6132  }
6133 
6134  return(tree);
6135 }
6136 /* END OF FUNCTION add_to_indexed_table_pairs */
6137 
6138 /************************************************************************/
6139 /* */
6140 /* FUNCTION add_to_indexed_vdw_pairs */
6141 /* */
6142 /* INPUTS: */
6143 /* new_node - new node to be added to the tree */
6144 /* tree - tree to add the node to */
6145 /* */
6146 /* This is a recursive function that adds a node to the */
6147 /* binary search tree of vdw_pair parameters */
6148 /* */
6149 /************************************************************************/
6150 
6151 IndexedVdwPair *Parameters::add_to_indexed_vdw_pairs(IndexedVdwPair *new_node,
6152  IndexedVdwPair *tree)
6153 
6154 {
6155  if (tree == NULL)
6156  return(new_node);
6157 
6158  if ( (new_node->ind1 < tree->ind1) ||
6159  ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
6160  {
6161  tree->left = add_to_indexed_vdw_pairs(new_node, tree->left);
6162  }
6163  else
6164  {
6165  tree->right = add_to_indexed_vdw_pairs(new_node, tree->right);
6166  }
6167 
6168  return(tree);
6169 }
6170 /* END OF FUNCTION add_to_indexed_vdw_pairs */
6171 
6172 /************************************************************************/
6173 /* */
6174 /* FUNCTION add_to_indexed_nbthole_pairs */
6175 /* */
6176 /* INPUTS: */
6177 /* new_node - new node to be added to the tree */
6178 /* tree - tree to add the node to */
6179 /* */
6180 /* This is a recursive function that adds a node to the */
6181 /* binary search tree of nbthole_pair parameters */
6182 /* */
6183 /************************************************************************/
6184 
6185 IndexedNbtholePair *Parameters::add_to_indexed_nbthole_pairs(IndexedNbtholePair *new_node,
6186  IndexedNbtholePair *tree)
6187 
6188 {
6189  if (tree == NULL)
6190  return(new_node);
6191 
6192  if ( (new_node->ind1 < tree->ind1) ||
6193  ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
6194  {
6195  tree->left = add_to_indexed_nbthole_pairs(new_node, tree->left);
6196  }
6197  else
6198  {
6199  tree->right = add_to_indexed_nbthole_pairs(new_node, tree->right);
6200  }
6201 
6202  return(tree);
6203 }
6204 /* END OF FUNCTION add_to_indexed_nbthole_pairs */
6205 
6206 /************************************************************************/
6207 /* */
6208 /* FUNCTION vdw_pair_to_arrays */
6209 /* */
6210 /* INPUTS: */
6211 /* ind1_array - Array of index 1 values */
6212 /* ind2_array - Array of index 2 values */
6213 /* A - Array of A values */
6214 /* A14 - Array of A 1-4 values */
6215 /* B - Array of B values */
6216 /* B14 - Array of B 1-4 values */
6217 /* arr_index - current position in arrays */
6218 /* tree - tree to traverse */
6219 /* */
6220 /* This is a recursive function that places all the entries of */
6221 /* the tree passed in into arrays of values. This is done so that */
6222 /* the parameters can be sent from the master node to the other */
6223 /* nodes. */
6224 /* */
6225 /************************************************************************/
6226 
6227 int Parameters::vdw_pair_to_arrays(int *ind1_array, int *ind2_array,
6228  Real *A, Real *A14,
6229  Real *B, Real *B14,
6230  int arr_index, IndexedVdwPair *tree)
6231 
6232 {
6233  if (tree == NULL)
6234  return(arr_index);
6235 
6236  ind1_array[arr_index] = tree->ind1;
6237  ind2_array[arr_index] = tree->ind2;
6238  A[arr_index] = tree->A;
6239  A14[arr_index] = tree->A14;
6240  B[arr_index] = tree->B;
6241  B14[arr_index] = tree->B14;
6242 
6243  arr_index++;
6244 
6245  arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
6246  arr_index, tree->left);
6247  arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
6248  arr_index, tree->right);
6249 
6250  return(arr_index);
6251 }
6252 /* END OF FUNCTION vdw_pair_to_arrays */
6253 
6254 /************************************************************************/
6255 /* */
6256 /* FUNCTION nbthole_pair_to_arrays */
6257 /* */
6258 /* INPUTS: */
6259 /* ind1_array - Array of index 1 values */
6260 /* ind2_array - Array of index 2 values */
6261 /* tholeij - Array of tholeij values */
6262 /* arr_index - current position in arrays */
6263 /* tree - tree to traverse */
6264 /* */
6265 /* This is a recursive function that places all the entries of */
6266 /* the tree passed in into arrays of values. This is done so that */
6267 /* the parameters can be sent from the master node to the other */
6268 /* nodes. */
6269 /* */
6270 /************************************************************************/
6271 
6272 int Parameters::nbthole_pair_to_arrays(int *ind1_array, int *ind2_array,
6273  Real *alphai, Real *alphaj, Real *tholeij,
6274  int arr_index, IndexedNbtholePair *tree)
6275 
6276 {
6277  if (tree == NULL)
6278  return(arr_index);
6279 
6280  ind1_array[arr_index] = tree->ind1;
6281  ind2_array[arr_index] = tree->ind2;
6282  alphai[arr_index] = tree->alphai;
6283  alphaj[arr_index] = tree->alphaj;
6284  tholeij[arr_index] = tree->tholeij;
6285 
6286  arr_index++;
6287 
6288  arr_index = nbthole_pair_to_arrays(ind1_array, ind2_array, alphai,
6289  alphaj, tholeij, arr_index, tree->left);
6290  arr_index = nbthole_pair_to_arrays(ind1_array, ind2_array, alphai,
6291  alphaj, tholeij, arr_index, tree->right);
6292 
6293  return(arr_index);
6294 }
6295 /* END OF FUNCTION nbthole_pair_to_arrays */
6296 
6297 /************************************************************************/
6298 /* */
6299 /* FUNCTION table_pair_to_arrays */
6300 /* */
6301 /* INPUTS: */
6302 /* ind1_array - Array of index 1 values */
6303 /* ind2_array - Array of index 2 values */
6304 /* K - Array of K values */
6305 /* */
6306 /* This is a recursive function that places all the entries of */
6307 /* the tree passed in into arrays of values. This is done so that */
6308 /* the parameters can be sent from the master node to the other */
6309 /* nodes. */
6310 /* */
6311 /************************************************************************/
6312 
6313 int Parameters::table_pair_to_arrays(int *ind1_array, int *ind2_array,
6314  int *K,
6315  int arr_index, IndexedTablePair *tree)
6316 
6317 {
6318  if (tree == NULL)
6319  return(arr_index);
6320 
6321  ind1_array[arr_index] = tree->ind1;
6322  ind2_array[arr_index] = tree->ind2;
6323  K[arr_index] = tree->K;
6324 
6325  arr_index++;
6326 
6327  arr_index = table_pair_to_arrays(ind1_array, ind2_array, K,
6328  arr_index, tree->left);
6329  arr_index = table_pair_to_arrays(ind1_array, ind2_array, K,
6330  arr_index, tree->right);
6331 
6332  return(arr_index);
6333 }
6334 /* END OF FUNCTION table_pair_to_arrays */
6335 
6336 /************************************************************************/
6337 /* */
6338 /* FUNCTION Parameters */
6339 /* */
6340 /* This is the constructor for reading AMBER parameter data */
6341 /* */
6342 /************************************************************************/
6343 
6345 {
6346  initialize();
6347 
6348  // Read in parm parameters
6349  read_parm(amber_data,vdw14);
6350 }
6351 /* END OF FUNCTION Parameters */
6352 
6353 
6354 /************************************************************************/
6355 /* */
6356 /* FUNCTION read_parm */
6357 /* */
6358 /* INPUTS: */
6359 /* amber_data - AMBER data structure */
6360 /* */
6361 /* This function copys AMBER parameter data to the corresponding data */
6362 /* structures */
6363 /* */
6364 /************************************************************************/
6365 
6366 void Parameters::read_parm(Ambertoppar *amber_data, BigReal vdw14)
6367 {
6368  int i,j,ico,numtype,mul;
6369  IndexedVdwPair *new_node;
6370 
6371  if (!amber_data->data_read)
6372  NAMD_die("No data read from parm file yet!");
6373 
6374  // Copy bond parameters
6375  NumBondParams = amber_data->Numbnd;
6376  if (NumBondParams)
6378  if (bond_array == NULL)
6379  NAMD_die("memory allocation of bond_array failed!");
6380  }
6381  for (i=0; i<NumBondParams; ++i)
6382  { bond_array[i].k = amber_data->Rk[i];
6383  bond_array[i].x0 = amber_data->Req[i];
6384  }
6385 
6386  // Copy angle parameters
6387  NumAngleParams = amber_data->Numang;
6388  if (NumAngleParams)
6390  if (angle_array == NULL)
6391  NAMD_die("memory allocation of angle_array failed!");
6392  }
6393  for (i=0; i<NumAngleParams; ++i)
6394  { angle_array[i].k = amber_data->Tk[i];
6395  angle_array[i].theta0 = amber_data->Teq[i];
6396  // Amber has no k_ub and r_ub for angle parameters, so they're set to 0
6397  angle_array[i].k_ub = angle_array[i].r_ub = 0;
6398  // All angles are harmonic
6399  angle_array[i].normal = 1;
6400  }
6401 
6402  // Copy dihedral parameters
6403  // Note: If the periodicity is negative, it means the following
6404  // entry is another term in a multitermed dihedral, and the
6405  // absolute value is the true periodicity; in this case the
6406  // following entry in "dihedral_array" should be left empty,
6407  // NOT be skipped, in order to keep the next dihedral's index
6408  // unchanged.
6409  NumDihedralParams = amber_data->Nptra;
6410  if (NumDihedralParams)
6411  { dihedral_array = new DihedralValue[amber_data->Nptra];
6412  if (dihedral_array == NULL)
6413  NAMD_die("memory allocation of dihedral_array failed!");
6414  }
6415  mul = 0;
6416  for (i=0; i<NumDihedralParams; ++i)
6417  { dihedral_array[i-mul].values[mul].k = amber_data->Pk[i];
6418  dihedral_array[i-mul].values[mul].n = int(fabs(amber_data->Pn[i])+0.5);
6419  if (dihedral_array[i-mul].values[mul].n == 0)
6420  { char err_msg[128];
6421  sprintf(err_msg, "The periodicity of dihedral # %d is zero!", i+1);
6422  NAMD_die(err_msg);
6423  }
6424  dihedral_array[i-mul].values[mul].delta = amber_data->Phase[i];
6425  // If the periodicity is positive, it means the following
6426  // entry is a new dihedral term.
6427  if (amber_data->Pn[i] > 0)
6428  { dihedral_array[i-mul].multiplicity = mul+1;
6429  mul = 0;
6430  }
6431  else if (++mul >= MAX_MULTIPLICITY)
6432  { char err_msg[181];
6433  sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
6434  mul+1, MAX_MULTIPLICITY);
6435  NAMD_die(err_msg);
6436  }
6437  }
6438  if (mul > 0)
6439  NAMD_die("Negative periodicity in the last dihedral entry!");
6440 
6441  // Copy VDW parameters: AMBER explicitly gives VDW parameters between every
6442  // 2 atom types, so the data are copied to vdw_pair_tree
6443  // In AMBER, all 1-4 VDW interactions are divided by factor vdw14
6444  NumVdwParamsAssigned = numtype = amber_data->Ntypes; // Number of atom types
6445  NumVdwPairParams = numtype * (numtype+1) / 2;
6446  for (i=0; i<numtype; ++i)
6447  for (j=i; j<numtype; ++j)
6448  { new_node = new IndexedVdwPair;
6449  if (new_node == NULL)
6450  NAMD_die("memory allocation of vdw_pair_tree failed!");
6451  new_node->ind1 = i;
6452  new_node->ind2 = j;
6453  new_node->left = new_node->right = NULL;
6454  // ico is the index of interaction between atom type i and j into
6455  // the parameter arrays. If it's positive, the interaction is
6456  // 6-12 VDW; otherwise it's 10-12 H-bond interaction. NAMD doesn't
6457  // have 10-12 term, so if ico is negative, then the 10-12
6458  // coefficients must be 0, otherwise die.
6459  ico = amber_data->Cno[numtype*i+j];
6460  if (ico>0)
6461  { new_node->A = amber_data->Cn1[ico-1];
6462  new_node->A14 = new_node->A / vdw14;
6463  new_node->B = amber_data->Cn2[ico-1];
6464  new_node->B14 = new_node->B / vdw14;
6465  }
6466  else if (amber_data->HB12[abs(ico)-1]==0.0 && amber_data->HB6[abs(ico)-1]==0.0)
6467  { new_node->A = new_node->A14 = new_node->B = new_node->B14 = 0.0;
6468  iout << iWARN << "Encounter 10-12 H-bond term\n";
6469  }
6470  else
6471  NAMD_die("Encounter non-zero 10-12 H-bond term!");
6472  // Add the node to the binary tree
6473  vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
6474  }
6475 }
6476 /* END OF FUNCTION read_parm */
6477 
6478 // new AMBER parm7 reader
6480 {
6481  initialize();
6482 
6483  // Read in parm parameters
6484  read_parm(amber_data,vdw14);
6485 }
6486 
6487 void Parameters::read_parm(AmberParm7Reader::Ambertoppar *amber_data, BigReal vdw14)
6488 {
6489  int i,j,ico,numtype,mul;
6490  IndexedVdwPair *new_node;
6491 
6492  if (!amber_data->HasData) {
6493  NAMD_die("No data read from parm file yet!");
6494  }
6495 
6496  // Check if we are reading a CHARMBER file
6497  if (amber_data->IsCharmmFF) {
6498  NAMD_die("You are using a CHARMM force field in AMBER format generated by CHAMBER or ParmEd.\n"
6499  "We do not support this format. Please consider using the native format (PSF) for CHARMM force field.");
6500  }
6501 
6502  // Copy bond parameters
6503  NumBondParams = amber_data->Numbnd;
6504  if (NumBondParams)
6506  if (bond_array == NULL)
6507  NAMD_die("memory allocation of bond_array failed!");
6508  }
6509  for (i=0; i<NumBondParams; ++i)
6510  { bond_array[i].k = amber_data->Rk[i];
6511  bond_array[i].x0 = amber_data->Req[i];
6512  }
6513 
6514  // Copy angle parameters
6515  NumAngleParams = amber_data->Numang;
6516  if (NumAngleParams)
6518  if (angle_array == NULL)
6519  NAMD_die("memory allocation of angle_array failed!");
6520  }
6521  for (i=0; i<NumAngleParams; ++i)
6522  { angle_array[i].k = amber_data->Tk[i];
6523  angle_array[i].theta0 = amber_data->Teq[i];
6524  // Amber has no k_ub and r_ub for angle parameters, so they're set to 0
6525  angle_array[i].k_ub = angle_array[i].r_ub = 0;
6526  // All angles are harmonic
6527  angle_array[i].normal = 1;
6528  }
6529 
6530  // Copy dihedral parameters
6531  // Note: If the periodicity is negative, it means the following
6532  // entry is another term in a multitermed dihedral, and the
6533  // absolute value is the true periodicity; in this case the
6534  // following entry in "dihedral_array" should be left empty,
6535  // NOT be skipped, in order to keep the next dihedral's index
6536  // unchanged.
6537  NumDihedralParams = amber_data->Nptra;
6538  if (NumDihedralParams)
6539  { dihedral_array = new DihedralValue[amber_data->Nptra];
6540  if (dihedral_array == NULL)
6541  NAMD_die("memory allocation of dihedral_array failed!");
6542  }
6543  mul = 0;
6544  for (i=0; i<NumDihedralParams; ++i)
6545  { dihedral_array[i-mul].values[mul].k = amber_data->Pk[i];
6546  dihedral_array[i-mul].values[mul].n = int(fabs(amber_data->Pn[i])+0.5);
6547  if (dihedral_array[i-mul].values[mul].n == 0)
6548  { char err_msg[128];
6549  sprintf(err_msg, "The periodicity of dihedral # %d is zero!", i+1);
6550  NAMD_die(err_msg);
6551  }
6552  dihedral_array[i-mul].values[mul].delta = amber_data->Phase[i];
6553  // If the periodicity is positive, it means the following
6554  // entry is a new dihedral term.
6555  if (amber_data->Pn[i] > 0)
6556  { dihedral_array[i-mul].multiplicity = mul+1;
6557  mul = 0;
6558  }
6559  else if (++mul >= MAX_MULTIPLICITY)
6560  { char err_msg[181];
6561  sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
6562  mul+1, MAX_MULTIPLICITY);
6563  NAMD_die(err_msg);
6564  }
6565  }
6566  if (mul > 0)
6567  NAMD_die("Negative periodicity in the last dihedral entry!");
6568 
6569  // Copy crossterms
6570  if (amber_data->HasCMAP) {
6571  NumCrosstermParams = amber_data->CMAPTypeCount;
6572  if (NumCrosstermParams > 0) {
6574  for (i = 0; i < NumCrosstermParams; ++i) {
6575  // check the resolutions at first
6576  const int N = amber_data->CMAPResolution[i];
6577  // NAMD does not support a resolution number other than CrosstermValue::dim - 1
6578  if (N != CrosstermValue::dim - 1) {
6579  char err_msg[128];
6580  sprintf(err_msg, "The table of %d crossterm type has a resolution(%d) other than %d!",
6581  i+1, N, CrosstermValue::dim - 1);
6582  NAMD_die(err_msg);
6583  }
6584  // code swipe from Parameters::index_crossterms()
6585  int l = 0;
6586  for (int j = 0; j < N; ++j) {
6587  for (int k = 0; k < N; ++k) {
6588  crossterm_array[i].c[j][k].d00 = amber_data->CMAPParameter[i][l];
6589  ++l;
6590  }
6591  }
6592  for (int j = 0; j < N; ++j) {
6593  crossterm_array[i].c[j][N].d00 = crossterm_array[i].c[j][0].d00;
6594  crossterm_array[i].c[N][j].d00 = crossterm_array[i].c[0][j].d00;
6595  }
6596  crossterm_array[i].c[N][N] = crossterm_array[i].c[0][0];
6597  crossterm_setup(&crossterm_array[i].c[0][0]);
6598  }
6599  }
6600  }
6601 
6602  // Copy VDW parameters: AMBER explicitly gives VDW parameters between every
6603  // 2 atom types, so the data are copied to vdw_pair_tree
6604  // In AMBER, all 1-4 VDW interactions are divided by factor vdw14
6605  NumVdwParamsAssigned = numtype = amber_data->Ntypes; // Number of atom types
6606  NumVdwPairParams = numtype * (numtype+1) / 2;
6607  for (i=0; i<numtype; ++i)
6608  for (j=i; j<numtype; ++j)
6609  { new_node = new IndexedVdwPair;
6610  if (new_node == NULL)
6611  NAMD_die("memory allocation of vdw_pair_tree failed!");
6612  new_node->ind1 = i;
6613  new_node->ind2 = j;
6614  new_node->left = new_node->right = NULL;
6615  // ico is the index of interaction between atom type i and j into
6616  // the parameter arrays. If it's positive, the interaction is
6617  // 6-12 VDW; otherwise it's 10-12 H-bond interaction. NAMD doesn't
6618  // have 10-12 term, so if ico is negative, then the 10-12
6619  // coefficients must be 0, otherwise die.
6620  ico = amber_data->Cno[numtype*i+j];
6621  if (ico>0)
6622  { new_node->A = amber_data->Cn1[ico-1];
6623  new_node->A14 = new_node->A / vdw14;
6624  new_node->B = amber_data->Cn2[ico-1];
6625  new_node->B14 = new_node->B / vdw14;
6626  }
6627  else if (amber_data->HB12[abs(ico)-1]==0.0 && amber_data->HB10[abs(ico)-1]==0.0)
6628  { new_node->A = new_node->A14 = new_node->B = new_node->B14 = 0.0;
6629  iout << iWARN << "Encounter 10-12 H-bond term\n";
6630  }
6631  else
6632  NAMD_die("Encounter non-zero 10-12 H-bond term!");
6633  // Add the node to the binary tree
6634  vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
6635  }
6636 }
6637 
6638 /************************************************************************/
6639 /* */
6640 /* FUNCTION Parameters */
6641 /* */
6642 /* This is the constructor for reading GROMACS parameter data */
6643 /* */
6644 /************************************************************************/
6645 
6647 {
6648  initialize();
6649 
6650  // Read in parm parameters
6651  read_parm(gf,min);
6652 }
6653 /* END OF FUNCTION Parameters */
6654 
6655 
6656 /************************************************************************/
6657 /* */
6658 /* FUNCTION read_parm */
6659 /* */
6660 /* INPUTS: */
6661 /* gf - GROMACS topology file */
6662 /* */
6663 /* This function copys GROMACS parameter data to the corresponding data */
6664 /* structures */
6665 /* */
6666 /************************************************************************/
6667 
6668 void Parameters::read_parm(const GromacsTopFile *gf, Bool min)
6669 {
6670  int numtype;
6671  IndexedVdwPair *new_node;
6672  int i,j,funct;
6673  Real test1,test2;
6674 
6675  // Allocate space for all of the arrays first
6679  if (NumBondParams) {
6681  if (bond_array == NULL)
6682  NAMD_die("memory allocation of bond_array failed!");
6683  }
6684  if (NumDihedralParams) {
6686  if (dihedral_array == NULL)
6687  NAMD_die("memory allocation of dihedral_array failed!");
6688  }
6689  if (NumAngleParams) {
6691  if (angle_array == NULL)
6692  NAMD_die("memory allocation of angle_array failed!");
6693  }
6694 
6695  // Copy bond parameters
6696  // XXX Warning: we are discarding the GROMACS function type - since
6697  // NAMD does not let us choose between different spring models.
6698  for (i=0;i<NumBondParams;i++) {
6699  Real x0,k;
6700  gf->getBondParams(i,
6701  &x0, // the bond length
6702  &k, // the spring constant
6703  &funct); // discarded
6704  bond_array[i].x0 = x0;
6705  bond_array[i].k = k;
6706  }
6707 
6708  // Copy angle parameters
6709  // XXX Warning: we are discarding the GROMACS function type here
6710  // too.
6711  for (i=0;i<NumAngleParams;i++) {
6712  Real theta0,k;
6713  gf->getAngleParams(i,
6714  &theta0, // the angle size
6715  &k, // the spring constant
6716  &funct); // discarded
6717  angle_array[i].theta0 = theta0*PI/180;
6718  angle_array[i].k = k;
6719  // Gromacs has no Urey-Bradley angle parameters, so they're set to 0
6720  angle_array[i].k_ub = angle_array[i].r_ub = 0;
6721  angle_array[i].normal = 1;
6722  }
6723 
6724  // Copy dihedral parameters
6725  // Here we use the function type (carefully)
6726  for (i=0; i<NumDihedralParams; ++i) { // iterate over all dihedral angles
6727  Real c[6];
6728  int num_periods; // number of periods in one full rotation
6729  int funct;
6730 
6731  gf->getDihedralParams(i,c,&num_periods,&funct); // get the parameters
6732 
6733  switch(funct) {
6734  case 1:
6735  dihedral_array[i].values[0].delta = c[0]*PI/180; // the phase shift
6736  dihedral_array[i].values[0].k = c[1]; // the spring constant
6737  dihedral_array[i].values[0].n = num_periods; // the periodicity
6739  break;
6740  case 2:
6741  dihedral_array[i].values[0].delta = c[0]*PI/180; // the phase shift
6742  dihedral_array[i].values[0].k = c[1]; // the spring constant
6743  dihedral_array[i].values[0].n = 0; // 'periodicity'=0 for impropers
6745  break;
6746  case 3:
6747 
6748  // first a quick check to make sure this is legal
6749  if(MAX_MULTIPLICITY < 5)
6750  NAMD_die("I can't do RB dihedrals with MAX_MULTIPLICITY < 5");
6752 
6753  // Next we negate every other term, since GROMACS does this
6754  // silly thing with psi = 180 - phi:
6755  for(j=0;j<6;j++) {
6756  if(j%2 == 1) c[j] = -c[j];
6757  }
6758 
6759  // Now fill up all the terms. Each is k(1 + cos(n*x - delta))
6760  // so we first let delta = 0 and let n range from 1-6:
6761  for(j=0;j<5;j++) dihedral_array[i].values[j].delta = 0;
6762  for(j=0;j<5;j++) dihedral_array[i].values[j].n = j+1;
6763 
6764  // and now we have a sum of kn(1+cos(nx))
6765  // Gromacs RB gives you a sum of cn(cos(x))^n, so we have to
6766  // use trigonometry to compute the kn:
6767  dihedral_array[i].values[0].k = 1*c[1] + 3/4.*c[3] + 10/16.*c[5];
6768  dihedral_array[i].values[1].k = 1/2.*c[2] + 4/8.*c[4] ;
6769  dihedral_array[i].values[2].k = 1/4.*c[3] + 5/16.*c[5];
6770  dihedral_array[i].values[3].k = 1/8.*c[4] ;
6771  dihedral_array[i].values[4].k = 1/16.*c[5];
6772 
6773  // That was a lot of math, so we had better check it:
6774  // The constant term (which is missing) is c0 + 1/2 c2 + 3/8 c4
6775  test1 = 0;
6776  for(j=5;j>=0;j--) { // sum (..(c5 cos x + c4) cos x + c3)..) + c1
6777  test1 *= cos(0.5);
6778  test1 += c[j];
6779  }
6780 
6781  test2 = c[0]+1/2.*c[2]+3/8.*c[4];
6782  for(j=0;j<5;j++) { // sum k1 cos x + k2 cos 2x + ...
6783  test2 += dihedral_array[i].values[j].k * cos((j+1)*0.5);
6784  }
6785 
6786  if(fabs(test1-test2) > 0.0001)
6787  NAMD_die("Internal error: failed to handle RB dihedrals");
6788 
6789  // Turn this on to have a look at the values if you *still*
6790  // don't believe that they are right!
6791 
6792  /* iout << iINFO << "k: ";
6793  for(j=0;j<5;j++)
6794  iout << dihedral_array[i].values[j].k << " ";
6795  iout << "\n" << endi;
6796 
6797  iout << iINFO << "c: ";
6798  for(j=0;j<6;j++)
6799  iout << c[j] << " ";
6800  iout << "\n" << endi;*/
6801 
6802  break;
6803  default:
6804  NAMD_die("unknown dihedral type found");
6805  }
6806  }
6807 
6808  // Copy VDW parameters.
6809 
6810  Bool warned=false; // warned the user about extra LJ term yet?
6811 
6812  NumVdwParamsAssigned = numtype = gf->getNumAtomParams(); // # of atom types
6813  NumVdwPairParams = numtype * (numtype+1) / 2;
6814  for (i=0; i<numtype; i++) {
6815  for (j=i; j<numtype; j++) {
6816 
6817  // set up a node to put one set of VDW parameters in
6818  new_node = new IndexedVdwPair;
6819  if (new_node == NULL)
6820  NAMD_die("memory allocation of vdw_pair_tree failed!");
6821  new_node->ind1 = i;
6822  new_node->ind2 = j;
6823  new_node->left = new_node->right = NULL;
6824 
6825  gf->getVDWParams(i,j, &(new_node->B), &(new_node->A),
6826  &(new_node->B14), &(new_node->A14));
6827 
6828  /* If we have any VDW radii equal to zero, atoms can just sit on
6829  each other during minimization. So, we'll give a minimum of
6830  1.0 kcal*A^12 to the LJ-repulsion when we are minimizing.
6831  But a warning should be displayed to the user... */
6832  if(min && ( fabs(new_node->A) < 1.0 )) {
6833  new_node->A = 1.0;
6834  if(!warned) {
6835  iout << iWARN <<
6836  "Adding small LJ repulsion term to some atoms.\n" << endi;
6837  warned=true;
6838  }
6839  }
6840 
6841  vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
6842  }
6843  }
6844 
6845  // JLai
6846  // Allocate space for all of the GromacsPair arrays first
6847  int numPair, numLJPair, numGaussPair;
6848  Real *pairC6,*pairC12; // constants to define LJ potential
6849  int *atom1,*atom2; // atom indices for LJ code
6850  atom1 = 0;
6851  atom2 = 0;
6852  pairC6 = 0;
6853  pairC12 = 0;
6854  numPair = gf->getNumPair();
6855  NumGromacsPairParams = numLJPair = gf->getNumLJPair();
6856  if (numLJPair) {
6857  atom1 = new int[numLJPair];
6858  atom2 = new int[numLJPair];
6859  pairC6 = new Real[numLJPair];
6860  pairC12 = new Real[numLJPair];
6862  }
6863 
6864  // Copy GromacsPair data into gromacsPair array structures
6865  const_cast<GromacsTopFile*>(gf)->getPairLJArrays2(atom1, atom2, pairC6, pairC12);
6866  for (i=0;i<numLJPair;i++) {
6867  gromacsPair_array[i].pairC6 = pairC6[i];
6868  gromacsPair_array[i].pairC12 = pairC12[i];
6869  }
6870  delete [] atom1;
6871  delete [] atom2;
6872  delete [] pairC6;
6873  delete [] pairC12;
6874 }
6875 /* END OF FUNCTION read_parm */
6876 
6877 /************************************************************************/
6878 /* */
6879 /* FUNCTION read_ener_table */
6880 /* */
6881 /* INPUTS: */
6882 /* simParams -- Simulation Parameters */
6883 /* */
6884 /* This function reads energy tables from a file and places them into */
6885 /* memory. */
6886 /************************************************************************/
6887 
6889  char* table_file = simParams->tabulatedEnergiesFile;
6890  char* interp_type = simParams->tableInterpType;
6891  FILE* enertable;
6892  enertable = fopen(table_file, "r");
6893 
6894  if (enertable == NULL) {
6895  NAMD_die("ERROR: Could not open energy table file!\n");
6896  }
6897 
6898  char tableline[121];
6899  char* newtypename;
6900  int numtypes;
6901  int atom1;
6902  int atom2;
6903  int distbin;
6904  int readcount;
6905  Real dist;
6906  Real energy;
6907  Real force;
6908  Real table_spacing;
6909  Real maxdist;
6910 
6911 /* First find the header line and set the variables we need */
6912  iout << "Beginning table read\n" << endi;
6913  while(fgets(tableline,120,enertable)) {
6914  if (strncmp(tableline,"#",1)==0) {
6915  continue;
6916  }
6917  readcount = sscanf(tableline, "%i %f %f", &numtypes, &table_spacing, &maxdist);
6918  if (readcount != 3) {
6919  NAMD_die("ERROR: Couldn't parse table header line\n");
6920  }
6921  break;
6922  }
6923 
6924  if (maxdist < simParams->cutoff) {
6925  NAMD_die("Tabulated energies must at least extend to the cutoff distance\n");
6926  }
6927 
6928  if (maxdist > simParams->cutoff) {
6929  iout << "Info: Reducing max table size to match nonbond cutoff\n";
6930  maxdist = ceil(simParams->cutoff);
6931  }
6932 
6933 /* Now allocate memory for the table; we know what we should be getting */
6934  numenerentries = 2 * numtypes * int (namdnearbyint(maxdist/table_spacing) + 1);
6935  //Set up a full energy lookup table from a file
6936  //Allocate the table; layout is atom1 x atom2 x distance energy force
6937  fprintf(stdout, "Table has %i entries\n",numenerentries);
6938  //iout << "Allocating original energy table\n" << endi;
6939  if ( table_ener ) delete [] table_ener;
6941  if ( table_types ) delete [] table_types;
6942  table_types = new char*[numtypes];
6943 
6944 /* Finally, start reading the table */
6945  int numtypesread = 0; //Number of types read so far
6946 
6947  while(fgets(tableline,120,enertable)) {
6948  if (strncmp(tableline,"#",1)==0) {
6949  continue;
6950  }
6951  if (strncmp(tableline,"TYPE",4)==0) {
6952  // Read a new type
6953  newtypename = new char[5];
6954  int readcount = sscanf(tableline, "%*s %s", newtypename);
6955  if (readcount != 1) {
6956  NAMD_die("Couldn't read interaction type from TYPE line\n");
6957  }
6958 // printf("Setting table_types[%i] to %s\n", numtypesread, newtypename);
6959  table_types[numtypesread] = newtypename;
6960 
6961  if (numtypesread == numtypes) {
6962  NAMD_die("Error: Number of types in table doesn't match table header\n");
6963  }
6964 
6965  // Read the current energy type with the proper interpolation
6966  if (!strncasecmp(interp_type, "linear", 6)) {
6967  if (read_energy_type(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
6968  char err_msg[512];
6969  sprintf(err_msg, "Failed to read table energy (linear) type %s\n", newtypename);
6970  NAMD_die(err_msg);
6971  }
6972  } else if (!strncasecmp(interp_type, "cubic", 5)) {
6973  if (read_energy_type_bothcubspline(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
6974  char err_msg[512];
6975  sprintf(err_msg, "Failed to read table energy (cubic) type %s\n", newtypename);
6976  NAMD_die(err_msg);
6977  }
6978  } else {
6979  NAMD_die("Table interpolation type must be linear or cubic\n");
6980  }
6981 
6982  numtypesread++;
6983  continue;
6984  }
6985  //char err_msg[512];
6986  //sprintf(err_msg, "Unrecognized line in energy table file: %s\n", tableline);
6987  //NAMD_die(err_msg);
6988  }
6989 
6990  // See if we got what we expected
6991  if (numtypesread != numtypes) {
6992  char err_msg[512];
6993  sprintf(err_msg, "ERROR: Expected %i tabulated energy types but got %i\n", numtypes, numtypesread);
6994  NAMD_die(err_msg);
6995  }
6996 
6997  // Move the necessary information into simParams
6998  simParams->tableNumTypes = numtypes;
6999  simParams->tableSpacing = table_spacing;
7000  simParams->tableMaxDist = maxdist;
7001  tablenumtypes = numtypes;
7002 
7003  /*
7004 xxxxxx
7005  int numtypes = simParams->tableNumTypes;
7006 
7007  //This parameter controls the table spacing
7008  BigReal table_spacing = 0.01;
7009  BigReal maxdist = 20.0;
7010  if (maxdist > simParams->cutoff) {
7011  iout << "Info: Reducing max table size to match nonbond cutoff\n";
7012  maxdist = ceil(simParams->cutoff);
7013  }
7014 
7015  numenerentries = (numtypes + 1) * numtypes * int (ceil(maxdist/table_spacing));
7016 // fprintf(stdout, "Table arithmetic: maxdist %f, table_spacing %f, numtypes %i, numentries %i\n", maxdist, table_spacing, numtypes, numenerentries);
7017  columnsize = 2 * namdnearbyint(maxdist/table_spacing);
7018  rowsize = numtypes * columnsize;
7019  //Set up a full energy lookup table from a file
7020  //Allocate the table; layout is atom1 x atom2 x distance energy force
7021  fprintf(stdout, "Table has %i entries\n",numenerentries);
7022  //iout << "Allocating original energy table\n" << endi;
7023  if ( table_ener ) delete [] table_ener;
7024  table_ener = new Real[numenerentries];
7025  //
7026  //Set sentinel values for uninitialized table energies
7027  for (int i=0 ; i<numenerentries ; i++) {
7028  table_ener[i] = 1337.0;
7029  }
7030  Real compval = 1337.0;
7031 
7032  // iout << "Entering table reading\n" << endi;
7033  //iout << "Finished allocating table\n" << endi;
7034 
7035  //Counter to make sure we read the right number of energy entries
7036  int numentries = 0;
7037 
7038  //Now, start reading from the file
7039  char* table_file = simParams->tabulatedEnergiesFile;
7040  FILE* enertable;
7041  enertable = fopen(table_file, "r");
7042 
7043  if (enertable == NULL) {
7044  NAMD_die("ERROR: Could not open energy table file!\n");
7045  }
7046 
7047  char tableline[121];
7048  int atom1;
7049  int atom2;
7050  int distbin;
7051  Real dist;
7052  Real energy;
7053  Real force;
7054 
7055  iout << "Beginning table read\n" << endi;
7056  //Read the table entries
7057  while(fgets(tableline,120,enertable)) {
7058 // IOut << "Processing line " << tableline << "\n" << endi;
7059  if (strncmp(tableline,"#",1)==0) {
7060  continue;
7061  }
7062 
7063 
7064  sscanf(tableline, "%i %i %f %f %f\n", &atom1, &atom2, &dist, &energy, &force);
7065  distbin = int(namdnearbyint(dist/table_spacing));
7066 // fprintf(stdout, "Working on atoms %i and %i at distance %f\n",atom1,atom2,dist);
7067  if ((atom2 > atom1) || (distbin > int(namdnearbyint(maxdist/table_spacing)))) {
7068 // fprintf(stdout,"Rejected\n");
7069 // fprintf(stdout, "Error: Throwing out energy line beyond bounds\n");
7070  // 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)));
7071  } else {
7072  //The magic formula for the number of columns from previous rows
7073  //in the triangular matrix is (2ni+i-i**2)/2
7074  //Here n is the number of types, and i is atom2
7075 // fprintf(stdout, "Input: atom1 %f atom2 %f columnsize %f \n", float(atom1), float(atom2), float(columnsize));
7076 // 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);
7077  int eneraddress = columnsize*((2*numtypes*atom2 + atom2 - atom2*atom2)/2) + columnsize*(atom1-atom2) + 2*distbin - 2;
7078  int forceaddress = eneraddress + 1;
7079 // 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]);
7080  fflush(stdout);
7081 // fprintf(stdout, "Checking for dupes: Looking at: %f %f \n", table_ener[eneraddress], table_ener[forceaddress]);
7082  if ((table_ener[eneraddress] == compval && table_ener[forceaddress] == compval)) {
7083  numentries++;
7084  table_ener[eneraddress] = energy;
7085  table_ener[forceaddress] = force;
7086 // 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]);
7087  //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin] = energy;
7088  //table_ener[rowsize*atom2 + columnsize*atom1 + 2*distbin + 1] = force;
7089 // fflush(stdout);
7090  } else {
7091 // fprintf(stdout,"Rejecting duplicate entry\n");
7092  }
7093  }
7094  // iout << "Done with line\n"<< endi;
7095  }
7096 
7097  // int should = numtypes * numtypes * (maxdist/table_spacing);
7098  // cout << "Entries: " << numentries << " should be " << should << "\n" << endi;
7099 // int exptypes = ceil((numtypes+1) * numtypes * (maxdist/table_spacing));
7100 //fprintf(stdout,"numtypes: %i maxdist: %f table_spacing: %f exptypes: %i\n",numtypes,maxdist,table_spacing);
7101  if (numentries != int(numenerentries/2)) {
7102  fprintf(stdout,"ERROR: Expected %i entries but got %i\n",numenerentries/2, numentries);
7103  NAMD_die("Number of energy table entries did not match expected value\n");
7104  }
7105  // iout << "Done with table\n"<< endi;
7106  fprintf(stdout, "Numenerentries: %i\n",numenerentries/2);
7107  */
7108 } /* END of function read_ener_table */
7109 
7110 /**************************************************************************
7111  * FUNCTION read_energy_type_bothcubspline
7112  *
7113  * Read a single type block from an energy table file, using cubic spline interpolation
7114  * Unlike _cubspline, the forces are interpolated separately
7115  *
7116  * Inputs:
7117  * enertable - File stream positioned at the start of the type block
7118  * typeindex - integer index of current type
7119  * table_ener - pointer to array to be filled with table entries
7120  * table_spacing - Spacing between table points (A)
7121  * maxdist - Longest distance needed in table
7122  *
7123  * Return values:
7124  * 0 on normal exit
7125  * 1 if not enough entries were present to fill out the table
7126  *
7127  * Note: enertable should be left positioned immediately BEFORE the next
7128  * TYPE block starts
7129  * **********************************************************************/
7130 
7131 int Parameters::read_energy_type_bothcubspline(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
7132 
7133  char tableline[120];
7134  int i,j;
7135 
7136  /* Last values read from table */
7137  BigReal readdist;
7138  BigReal readener;
7139  BigReal readforce;
7140  BigReal spacing;
7141 // BigReal readforce;
7142  BigReal lastdist;
7143 // BigReal lastener;
7144 // BigReal lastforce;
7145 // readdist = -1.0;
7146 // readener = 0.0;
7147 // readforce = 0.0;
7148 
7149  /* Create arrays for holding the input values */
7150  std::vector<BigReal> dists;
7151  std::vector<BigReal> enervalues;
7152  std::vector<BigReal> forcevalues;
7153  int numentries = 0;
7154 
7155 
7156  /* Keep track of where in the table we are */
7157  BigReal currdist;
7158  int distbin;
7159  currdist = 0.0;
7160  lastdist = -1.0;
7161  distbin = 0;
7162 
7163  // Read all of the values first -- we'll interpolate later
7164  while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
7165  if (strncmp(tableline,"#",1)==0) {
7166  continue;
7167  }
7168  if (strncmp(tableline,"TYPE",4)==0) {
7169  fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR);
7170  break;
7171  }
7172 
7173  // Read an energy line from the table
7174  int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
7175 
7176  //printf("Read an energy line: %g %g %g\n", readdist, readener, readforce);
7177  if (readcount != 3) {
7178  char err_msg[512];
7179  sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
7180  NAMD_die(err_msg);
7181  }
7182 
7183  //Sanity check the current entry
7184  if (readdist < lastdist) {
7185  NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
7186  }
7187 
7188  lastdist = readdist;
7189  dists.push_back(readdist);
7190  enervalues.push_back(readener);
7191  forcevalues.push_back(readforce);
7192  numentries++;
7193  }
7194 
7195  // Check the spacing and make sure it is uniform
7196  if (dists[0] != 0.0) {
7197  NAMD_die("ERROR: First data point for energy table must be at r=0\n");
7198  }
7199  spacing = dists[1] - dists[0];
7200  for (i=2; i<(numentries - 1); i++) {
7201  BigReal myspacing;
7202  myspacing = dists[i] - dists[i-1];
7203  if (fabs(myspacing - spacing) > 0.00001) {
7204  printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
7205  NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
7206  }
7207  }
7208 
7209 // Perform cubic spline interpolation to get the energies and forces
7210 
7211  /* allocate spline coefficient matrix */
7212  // xe and xf / be and bf for energies and forces, respectively
7213  double* m = new double[numentries*numentries];
7214  double* xe = new double[numentries];
7215  double* xf = new double[numentries];
7216  double* be = new double[numentries];
7217  double* bf = new double[numentries];
7218 
7219  be[0] = 3 * (enervalues[1] - enervalues[0]);
7220  for (i=1; i < (numentries - 1); i++) {
7221 // printf("Control point %i at %f\n", i, enervalues[i]);
7222  be[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
7223 // printf("be is %f\n", be[i]);
7224  }
7225  be[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
7226 
7227  bf[0] = 3 * (forcevalues[1] - forcevalues[0]);
7228  for (i=1; i < (numentries - 1); i++) {
7229 // printf("Control point %i at %f\n", i, forcevalues[i]);
7230  bf[i] = 3 * (forcevalues[i+1] - forcevalues[i-1]);
7231 // printf("bf is %f\n", bf[i]);
7232  }
7233  bf[numentries - 1] = 3 * (forcevalues[numentries - 1] - forcevalues[numentries - 2]);
7234 
7235  memset(m,0,numentries*numentries*sizeof(double));
7236 
7237  /* initialize spline coefficient matrix */
7238  m[0] = 2;
7239  for (i = 1; i < numentries; i++) {
7240  m[INDEX(numentries,i-1,i)] = 1;
7241  m[INDEX(numentries,i,i-1)] = 1;
7242  m[INDEX(numentries,i,i)] = 4;
7243  }
7244  m[INDEX(numentries,numentries-1,numentries-1)] = 2;
7245 
7246  /* Now we need to solve the equation M X = b for X */
7247  // Do this simultaneously for energy and force -- ONLY because we use the same matrix
7248 
7249  //Put m in upper triangular form and apply corresponding operations to b
7250  for (i=0; i<numentries; i++) {
7251  // zero the ith column in all rows below i
7252  const BigReal baseval = m[INDEX(numentries,i,i)];
7253  for (j=i+1; j<numentries; j++) {
7254  const BigReal myval = m[INDEX(numentries,j,i)];
7255  const BigReal factor = myval / baseval;
7256 
7257  for (int k=i; k<numentries; k++) {
7258  const BigReal subval = m[INDEX(numentries,i,k)];
7259  m[INDEX(numentries,j,k)] -= (factor * subval);
7260  }
7261 
7262  be[j] -= (factor * be[i]);
7263  bf[j] -= (factor * bf[i]);
7264 
7265  }
7266  }
7267 
7268  //Now work our way diagonally up from the bottom right to assign values to X
7269  for (i=numentries-1; i>=0; i--) {
7270 
7271  //Subtract the effects of previous columns
7272  for (j=i+1; j<numentries; j++) {
7273  be[i] -= ( m[INDEX(numentries,i,j)] * xe[j] );
7274  bf[i] -= ( m[INDEX(numentries,i,j)] * xf[j] );
7275  }
7276 
7277  xe[i] = be[i] / m[INDEX(numentries,i,i)];
7278  xf[i] = bf[i] / m[INDEX(numentries,i,i)];
7279 
7280  }
7281 
7282  // We now have the coefficient information we need to make the table
7283  // Now interpolate on each interval we want
7284 
7285  distbin = 0;
7286  int entriesperseg = (int) ceil(spacing / table_spacing);
7287  int table_prefix = 2 * typeindex * (int) (namdnearbyint(maxdist / table_spacing) + 1);
7288 
7289  for (i=0; i<numentries-1; i++) {
7290  BigReal Ae,Be,Ce,De;
7291  BigReal Af,Bf,Cf,Df;
7292  currdist = dists[i];
7293 
7294 // printf("Interpolating on interval %i\n", i);
7295 
7296  // Set up the coefficients for this section
7297  Ae = enervalues[i];
7298  Be = xe[i];
7299  Ce = 3 * (enervalues[i+1] - enervalues[i]) - (2 * xe[i]) - xe[i+1];
7300  De = 2 * (enervalues[i] - enervalues[i+1]) + xe[i] + xe[i+1];
7301 
7302  Af = forcevalues[i];
7303  Bf = xf[i];
7304  Cf = 3 * (forcevalues[i+1] - forcevalues[i]) - (2 * xf[i]) - xf[i+1];
7305  Df = 2 * (forcevalues[i] - forcevalues[i+1]) + xf[i] + xf[i+1];
7306 
7307  // Go over the region of interest and fill in the table
7308  for (j=0; j<entriesperseg; j++) {
7309  const BigReal mydist = currdist + (j * table_spacing);
7310  const BigReal mydistfrac = (float) j / (entriesperseg - 1);
7311  distbin = (int) namdnearbyint(mydist / table_spacing);
7312  if (distbin >= (int) namdnearbyint(maxdist / table_spacing)) break;
7313  BigReal energy;
7314  BigReal force;
7315 
7316  energy = Ae + (Be * mydistfrac) + (Ce * mydistfrac * mydistfrac) + (De * mydistfrac * mydistfrac * mydistfrac);
7317  force = Af + (Bf * mydistfrac) + (Cf * mydistfrac * mydistfrac) + (Df * mydistfrac * mydistfrac * mydistfrac);
7318 
7319 // 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);
7320  table_ener[table_prefix + 2 * distbin] = energy;
7321  table_ener[table_prefix + 2 * distbin + 1] = force;
7322  distbin++;
7323  }
7324  if (currdist >= maxdist) break;
7325  }
7326 
7327  //The procedure above leaves out the last entry -- add it explicitly
7328  distbin = (int) namdnearbyint(maxdist / table_spacing);
7329 // 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));
7330  table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
7331  table_ener[table_prefix + 2 * distbin + 1] = 0.0;
7332  distbin++;
7333 
7334 
7335  // Clean up and make sure everything worked ok
7336  delete m;
7337  delete xe;
7338  delete xf;
7339  delete be;
7340  delete bf;
7341  distbin--;
7342  printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
7343  if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
7344  return 0;
7345 } /* end read_energy_type_bothcubspline */
7346 
7347 /**************************************************************************
7348  * FUNCTION read_energy_type_cubspline
7349  *
7350  * Read a single type block from an energy table file, using cubic spline interpolation
7351  *
7352  * Inputs:
7353  * enertable - File stream positioned at the start of the type block
7354  * typeindex - integer index of current type
7355  * table_ener - pointer to array to be filled with table entries
7356  * table_spacing - Spacing between table points (A)
7357  * maxdist - Longest distance needed in table
7358  *
7359  * Return values:
7360  * 0 on normal exit
7361  * 1 if not enough entries were present to fill out the table
7362  *
7363  * Note: enertable should be left positioned immediately BEFORE the next
7364  * TYPE block starts
7365  * **********************************************************************/
7366 
7367 int Parameters::read_energy_type_cubspline(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
7368 
7369  char tableline[120];
7370  int i,j;
7371 
7372  /* Last values read from table */
7373  BigReal readdist;
7374  BigReal readener;
7375  BigReal spacing;
7376 // BigReal readforce;
7377  BigReal lastdist;
7378 // BigReal lastener;
7379 // BigReal lastforce;
7380 // readdist = -1.0;
7381 // readener = 0.0;
7382 // readforce = 0.0;
7383 
7384  /* Create arrays for holding the input values */
7385  std::vector<BigReal> dists;
7386  std::vector<BigReal> enervalues;
7387  int numentries = 0;
7388 
7389 
7390  /* Keep track of where in the table we are */
7391  BigReal currdist;
7392  int distbin;
7393  currdist = 0.0;
7394  lastdist = -1.0;
7395  distbin = 0;
7396 
7397  // Read all of the values first -- we'll interpolate later
7398  while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
7399  if (strncmp(tableline,"#",1)==0) {
7400  continue;
7401  }
7402  if (strncmp(tableline,"TYPE",4)==0) {
7403  fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR);
7404  break;
7405  }
7406 
7407  // Read an energy line from the table
7408  int readcount = sscanf(tableline, "%lf %lf", &readdist, &readener);
7409 
7410  // printf("Read an energy line: %f %f %f\n", readdist, readener, readforce);
7411  if (readcount != 2) {
7412  char err_msg[512];
7413  sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
7414  NAMD_die(err_msg);
7415  }
7416 
7417  //Sanity check the current entry
7418  if (readdist < lastdist) {
7419  NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
7420  }
7421 
7422  lastdist = readdist;
7423  dists.push_back(readdist);
7424  enervalues.push_back(readener);
7425  numentries++;
7426  }
7427 
7428  // Check the spacing and make sure it is uniform
7429  if (dists[0] != 0.0) {
7430  NAMD_die("ERROR: First data point for energy table must be at r=0\n");
7431  }
7432  spacing = dists[1] - dists[0];
7433  for (i=2; i<(numentries - 1); i++) {
7434  BigReal myspacing;
7435  myspacing = dists[i] - dists[i-1];
7436  if (fabs(myspacing - spacing) > 0.00001) {
7437  printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
7438  NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
7439  }
7440  }
7441 
7442 // Perform cubic spline interpolation to get the energies and forces
7443 
7444  /* allocate spline coefficient matrix */
7445  double* m = new double[numentries*numentries];
7446  double* x = new double[numentries];
7447  double* b = new double[numentries];
7448 
7449  b[0] = 3 * (enervalues[1] - enervalues[0]);
7450  for (i=1; i < (numentries - 1); i++) {
7451  printf("Control point %i at %f\n", i, enervalues[i]);
7452  b[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
7453  printf("b is %f\n", b[i]);
7454  }
7455  b[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
7456 
7457  memset(m,0,numentries*numentries*sizeof(double));
7458 
7459  /* initialize spline coefficient matrix */
7460  m[0] = 2;
7461  for (i = 1; i < numentries; i++) {
7462  m[INDEX(numentries,i-1,i)] = 1;
7463  m[INDEX(numentries,i,i-1)] = 1;
7464  m[INDEX(numentries,i,i)] = 4;
7465  }
7466  m[INDEX(numentries,numentries-1,numentries-1)] = 2;
7467 
7468  /* Now we need to solve the equation M X = b for X */
7469 
7470  printf("Solving the matrix equation: \n");
7471 
7472  for (i=0; i<numentries; i++) {
7473  printf(" ( ");
7474  for (j=0; j<numentries; j++) {
7475  printf(" %6.3f,", m[INDEX(numentries, i, j)]);
7476  }
7477  printf(" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
7478  }
7479 
7480  //Put m in upper triangular form and apply corresponding operations to b
7481  for (i=0; i<numentries; i++) {
7482  // zero the ith column in all rows below i
7483  const BigReal baseval = m[INDEX(numentries,i,i)];
7484  for (j=i+1; j<numentries; j++) {
7485  const BigReal myval = m[INDEX(numentries,j,i)];
7486  const BigReal factor = myval / baseval;
7487 
7488  for (int k=i; k<numentries; k++) {
7489  const BigReal subval = m[INDEX(numentries,i,k)];
7490  m[INDEX(numentries,j,k)] -= (factor * subval);
7491  }
7492 
7493  b[j] -= (factor * b[i]);
7494 
7495  }
7496  }
7497 
7498  printf(" In upper diagonal form, equation is:\n");
7499  for (i=0; i<numentries; i++) {
7500  printf(" ( ");
7501  for (j=0; j<numentries; j++) {
7502  printf(" %6.3f,", m[INDEX(numentries, i, j)]);
7503  }
7504  printf(" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
7505  }
7506 
7507  //Now work our way diagonally up from the bottom right to assign values to X
7508  for (i=numentries-1; i>=0; i--) {
7509 
7510  //Subtract the effects of previous columns
7511  for (j=i+1; j<numentries; j++) {
7512  b[i] -= ( m[INDEX(numentries,i,j)] * x[j] );
7513  }
7514 
7515  x[i] = b[i] / m[INDEX(numentries,i,i)];
7516 
7517  }
7518 
7519  printf(" Solution vector is:\n\t(");
7520  for (i=0; i<numentries; i++) {
7521  printf(" %6.3f ", x[i]);
7522  }
7523  printf(" ) \n");
7524 
7525  // We now have the coefficient information we need to make the table
7526  // Now interpolate on each interval we want
7527 
7528  distbin = 0;
7529  int entriesperseg = (int) ceil(spacing / table_spacing);
7530  int table_prefix = 2 * typeindex * (int) (namdnearbyint(maxdist / table_spacing) + 1);
7531 
7532  for (i=0; i<numentries-1; i++) {
7533  BigReal A,B,C,D;
7534  currdist = dists[i];
7535 
7536  printf("Interpolating on interval %i\n", i);
7537 
7538  // Set up the coefficients for this section
7539  A = enervalues[i];
7540  B = x[i];
7541  C = 3 * (enervalues[i+1] - enervalues[i]) - (2 * x[i]) - x[i+1];
7542  D = 2 * (enervalues[i] - enervalues[i+1]) + x[i] + x[i+1];
7543 
7544  printf("Coefficients for this interval: %f %f %f %f\n", A, B, C, D);
7545 
7546  // Go over the region of interest and fill in the table
7547  for (j=0; j<entriesperseg; j++) {
7548  const BigReal mydist = currdist + (j * table_spacing);
7549  const BigReal mydistfrac = (float) j / (entriesperseg - 1);
7550  distbin = (int) namdnearbyint(mydist / table_spacing);
7551  if (distbin >= (int) namdnearbyint(maxdist / table_spacing)) break;
7552  BigReal energy;
7553  BigReal force;
7554 
7555  energy = A + (B * mydistfrac) + (C * mydistfrac * mydistfrac) + (D * mydistfrac * mydistfrac * mydistfrac);
7556  force = B + (2 * C * mydistfrac) + (3 * D * mydistfrac * mydistfrac);
7557  // Multiply force by 1 / (interval length)
7558  force *= (1.0 / spacing);
7559 
7560  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);
7561  table_ener[table_prefix + 2 * distbin] = energy;
7562  table_ener[table_prefix + 2 * distbin + 1] = force;
7563  distbin++;
7564  }
7565  if (currdist >= maxdist) break;
7566  }
7567 
7568  //The procedure above leaves out the last entry -- add it explicitly
7569  distbin = (int) namdnearbyint(maxdist / table_spacing);
7570  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));
7571  table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
7572  table_ener[table_prefix + 2 * distbin + 1] = 0.0;
7573  distbin++;
7574 
7575 
7576  // Clean up and make sure everything worked ok
7577  delete m;
7578  delete x;
7579  delete b;
7580  distbin--;
7581  printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
7582  if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
7583  return 0;
7584 } /* end read_energy_type_cubspline */
7585 
7586 /**************************************************************************
7587  * FUNCTION read_energy_type
7588  *
7589  * Read a single type block from an energy table file
7590  *
7591  * Inputs:
7592  * enertable - File stream positioned at the start of the type block
7593  * typeindex - integer index of current type
7594  * table_ener - pointer to array to be filled with table entries
7595  * table_spacing - Spacing between table points (A)
7596  * maxdist - Longest distance needed in table
7597  *
7598  * Return values:
7599  * 0 on normal exit
7600  * 1 if not enough entries were present to fill out the table
7601  *
7602  * Note: enertable should be left positioned immediately BEFORE the next
7603  * TYPE block starts
7604  * **********************************************************************/
7605 
7606 int Parameters::read_energy_type(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
7607 
7608  char tableline[120];
7609 
7610  /* Last values read from table */
7611  BigReal readdist;
7612  BigReal readener;
7613  BigReal readforce;
7614  BigReal lastdist;
7615  BigReal lastener;
7616  BigReal lastforce;
7617  readdist = -1.0;
7618  readener = 0.0;
7619  readforce = 0.0;
7620 
7621  /* Keep track of where in the table we are */
7622  float currdist;
7623  int distbin;
7624  currdist = -1.0;
7625  distbin = -1;
7626 
7627  while(fgets(tableline,120,enertable) && distbin <= (int) (namdnearbyint(maxdist / table_spacing) + 1)) {
7628  printf("At distance %f + %f vs. %f\n", currdist, table_spacing, maxdist);
7629  if (strncmp(tableline,"#",1)==0) {
7630  continue;
7631  }
7632  if (strncmp(tableline,"TYPE",4)==0) {
7633  break;
7634  }
7635 
7636  // Read an energy line from the table
7637  lastdist = readdist;
7638  lastener = readener;
7639  lastforce = readforce;
7640  int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
7641  if (distbin == -1) {
7642  if (readdist != 0.0) {
7643  NAMD_die("ERROR: Energy/force tables must start at d=0.0\n");
7644  } else {
7645  distbin = 0;
7646  continue;
7647  }
7648  }
7649  // printf("Read an energy line: %f %f %f\n", readdist, readener, readforce);
7650  if (readcount != 3) {
7651  char err_msg[512];
7652  sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
7653  NAMD_die(err_msg);
7654  }
7655 
7656  //Sanity check the current entry
7657  if (readdist < lastdist) {
7658  NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
7659  }
7660 
7661  currdist = lastdist;
7662 
7663  while (currdist <= readdist && distbin <= (int) (namdnearbyint(maxdist / table_spacing))) {
7664  distbin = (int) (namdnearbyint(currdist / table_spacing));
7665  int table_loc = 2 * (distbin + (typeindex * (1 + (int) namdnearbyint(maxdist / table_spacing))));
7666  printf("Doing interpolation for energy between %f %f and %f %f: Dist %f\n", readener, readdist, lastener, lastdist, currdist);
7667  table_ener[table_loc] = interp_lin(readener, lastener, readdist, lastdist, currdist);
7668  table_ener[table_loc + 1] = interp_lin(readforce, lastforce, readdist, lastdist, currdist);
7669  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);
7670  currdist += table_spacing;
7671  distbin++;
7672  }
7673  }
7674 
7675  // Go back one line, since we should now be into the next TYPE block
7676  fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR);
7677 
7678  // Clean up and make sure everything worked ok
7679  distbin--;
7680  printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
7681  if (distbin != (int) (namdnearbyint(maxdist / table_spacing))) return 1;
7682  return 0;
7683 }
7684 
7685 /*********************************************************************
7686  * FUNCTION interp_lin
7687  *
7688  * Perform a linear interpolation to fill in energy/force tables
7689  * This should be replaced in the near future with a better interpolation
7690  *
7691  * Input:
7692  * val1,val2 -- Y Values at the endpoints of the segments we interpolate on
7693  * end1,end2 -- X coordinates of the corresponding endpoints
7694  * currdist -- Distance we want the value at
7695  * ** It is assumed that end2 > end1 **
7696  *
7697  * Output: Returns a floating point value at the requested point
7698  * ********************************************************************/
7699 
7700 BigReal Parameters::interp_lin(BigReal val1, BigReal val2, BigReal end1, BigReal end2, BigReal currdist) {
7701 
7702  BigReal m; //slope of line
7703  BigReal val; // Value at desired point
7704 
7705  m = (val2 - val1) / (end2 - end1);
7706 
7707  val = ((currdist-end1) * m + val1);
7708  return val;
7709 }
7710 
7711 /*************************************************************************
7712  * FUNCTION get_int_table_type
7713  *
7714  * Find and return the integer index of a table type given its name
7715  *
7716  * Input:
7717  * tabletype -- char array containing the name of the type to be looked up
7718  *
7719  * Output:
7720  * Returns an integer index < the total number of types, or -1 if the type could
7721  * not be found
7722  * ************************************************************************/
7723 
7724 int Parameters::get_int_table_type(char* tabletype) {
7725  for (int i=0; i<tablenumtypes; i++) {
7726  if (!strncmp(tabletype, table_types[i], 5)) {
7727  return i;
7728  }
7729  }
7730 
7731  return -1;
7732 }
7733 
7734 
_REAL * Pk
Definition: parm.h:24
struct improper_params * next
Definition: Parameters.C:113
int numenerentries
Definition: Parameters.h:253
Real forceconstant
Definition: Parameters.C:49
int NumTablePairParams
Definition: Parameters.h:275
#define paraCharmm
Definition: Parameters.h:79
char atom8name[11]
Definition: Parameters.C:131
#define namdnearbyint(x)
Definition: common.h:76
GromacsPairValue * gromacsPair_array
Definition: Parameters.h:249
void assign_improper_index(const char *, const char *, const char *, const char *, Improper *, int)
Definition: Parameters.C:4065
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:81
Index index
Definition: Parameters.C:51
Index crossterm_type
Definition: structures.h:89
char atom4name[11]
Definition: Parameters.C:86
int32 atom4
Definition: structures.h:75
_REAL * Teq
Definition: parm.h:24
Index improper_type
Definition: structures.h:76
int NumBondParams
Definition: Parameters.h:261
CrosstermData c[dim][dim]
Definition: Parameters.h:124
void end(void)
Definition: MStream.C:176
Real epsilon
Definition: Parameters.h:153
char atom2name[11]
Definition: Parameters.C:107
char atom5name[11]
Definition: Parameters.C:128
int numGaussPair
int NAMD_read_line(FILE *fd, char *buf, int bufsize)
Definition: strlib.C:38
void assign_vdw_index(const char *, Atom *)
Definition: Parameters.C:3470
void print_bond_params()
Definition: Parameters.C:4847
void read_charmm_parameter_file(char *)
Definition: Parameters.C:539
static char ** table_types
Definition: Parameters.C:39
_REAL * Phase
Definition: parm.h:24
Real distance
Definition: Parameters.C:50
char atom1name[11]
Definition: Parameters.C:124
struct indexed_vdw_pair * right
Definition: Parameters.h:170
int32 atom3
Definition: structures.h:65
void assign_bond_index(const char *, const char *, Bond *)
Definition: Parameters.C:3758
void getVDWParams(int typea, int typeb, Real *c6, Real *c12, Real *c6pair, Real *c7) const
struct indexed_table_pair * right
Definition: Parameters.h:203
int NumVdwParams
Definition: Parameters.h:270
const BigReal A
Real r_ub
Definition: Parameters.h:96
float Real
Definition: common.h:109
#define DebugM(x, y)
Definition: Debug.h:59
char atom2name[11]
Definition: Parameters.C:159
int NumVdwPairParams
Definition: Parameters.h:273
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
#define FALSE
Definition: common.h:118
struct nbthole_pair_params * next
Definition: Parameters.C:183
Real sigma14
Definition: Parameters.h:154
int32 atom8
Definition: structures.h:88
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:82
MIStream * get(char &data)
Definition: MStream.h:29
char atom3name[11]
Definition: Parameters.C:108
int read_energy_type_cubspline(FILE *, const int, BigReal *, const float, const float)
Definition: Parameters.C:7367
Real sigma14
Definition: Parameters.C:146
DihedralValue * dihedral_array
Definition: Parameters.h:245
struct crossterm_params * next
Definition: Parameters.C:135
int32 atom1
Definition: structures.h:81
char atom3name[11]
Definition: Parameters.C:85
_REAL * Tk
Definition: parm.h:24
#define iout
Definition: InfoStream.h:51
void NAMD_find_first_word(char *source, char *word)
Definition: strlib.C:258
crossterm_params(int dim)
Definition: Parameters.C:118
char atomname[11]
Definition: Parameters.C:143
int NumAngleParams
Definition: Parameters.h:262
CrosstermValue * crossterm_array
Definition: Parameters.h:247
char atom1name[11]
Definition: Parameters.C:169
int32 atom4
Definition: structures.h:84
char atom2name[11]
Definition: Parameters.C:170
char atom3name[11]
Definition: Parameters.C:64
IndexedVdwPair * vdw_pair_tree
Definition: Parameters.h:257
int read_energy_type(FILE *, const int, BigReal *, const float, const float)
Definition: Parameters.C:7606
Real epsilon14
Definition: Parameters.C:147
int NumDihedralParams
Definition: Parameters.h:264
char tabulatedEnergiesFile[NAMD_FILENAME_BUFFER_SIZE]
#define INDEX(ncols, i, j)
Definition: Parameters.C:35
int Numbnd
Definition: parm.h:17
void print_param_summary()
Definition: Parameters.C:4966
IndexedNbtholePair * nbthole_pair_tree
Definition: Parameters.h:258
void read_parameter_file(char *)
Definition: Parameters.C:404
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.h:110
int32 atom4
Definition: structures.h:66
AngleValue * angle_array
Definition: Parameters.h:244
struct vdw_params * left
Definition: Parameters.C:149
int32 atom1
Definition: structures.h:48
_REAL * HB6
Definition: parm.h:24
int Index
Definition: structures.h:26
void getBondParams(int num, Real *b0, Real *kB, int *funct) const
int getNumLJPair() const
#define PI
Definition: common.h:83
_REAL * HB12
Definition: parm.h:24
static Units next(Units u)
Definition: ParseOptions.C:48
struct indexed_vdw_pair * left
Definition: Parameters.h:171
int get_int_table_type(char *)
Definition: Parameters.C:7724
int NumNbtholePairParams
Definition: Parameters.h:274
#define MAX_MULTIPLICITY
Definition: Parameters.h:70
struct bond_params * right
Definition: Parameters.C:53
int Ntypes
Definition: parm.h:17
void assign_crossterm_index(const char *, const char *, const char *, const char *, const char *, const char *, const char *, const char *, Crossterm *)
Definition: Parameters.C:4164
char atom7name[11]
Definition: Parameters.C:130
int normal
Definition: Parameters.h:97
Real sigma
Definition: Parameters.h:152
int NAMD_blank_string(char *str)
Definition: strlib.C:222
Index dihedral_type
Definition: structures.h:67
void print_vdw_pair_params()
Definition: Parameters.C:4932
int32 atom2
Definition: structures.h:56
int NumImproperParams
Definition: Parameters.h:265
struct indexed_table_pair * left
Definition: Parameters.h:204
int32 atom1
Definition: structures.h:55
int Bool
Definition: common.h:133
void print_dihedral_params()
Definition: Parameters.C:4881
_REAL * Cn2
Definition: parm.h:24
FILE * Fopen(const char *filename, const char *mode)
Definition: common.C:273
double * values
Definition: Parameters.C:133
int32 atom3
Definition: structures.h:57
char atom1name[11]
Definition: Parameters.C:47
BigReal d00
Definition: Parameters.h:119
struct indexed_vdw_pair IndexedVdwPair
char atom1name[11]
Definition: Parameters.C:62
char atom3name[11]
Definition: Parameters.C:126
int NumCosAngles
Definition: Parameters.h:263
int NumCrosstermParams
Definition: Parameters.h:266
Definition: parm.h:15
ImproperValue * improper_array
Definition: Parameters.h:246
int get_vdw_pair_params(Index ind1, Index ind2, Real *, Real *, Real *, Real *)
Definition: Parameters.C:3680
void NAMD_die(const char *err_msg)
Definition: common.C:85
BigReal tableMaxDist
int getNumDihedralParams() const
int NumVdwParamsAssigned
Definition: Parameters.h:272
char atom1name[11]
Definition: Parameters.C:83
NbtholePairValue * nbthole_array
Definition: Parameters.h:252
Real forceconstant
Definition: Parameters.C:65
char atom4name[11]
Definition: Parameters.C:109
char * atom_type_name(Index a)
Definition: Parameters.h:396
Index angle_type
Definition: structures.h:58
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.h:116
struct indexed_nbthole_pair * right
Definition: Parameters.h:181
struct bond_params * left
Definition: Parameters.C:52
struct angle_params * left
Definition: Parameters.C:71
struct dihedral_params * next
Definition: Parameters.C:94
int getNumAtomParams() const
_REAL * Rk
Definition: parm.h:24
int32 atom2
Definition: structures.h:64
int get_table_pair_params(Index, Index, int *)
Definition: Parameters.C:3605
void send_Parameters(MOStream *)
Definition: Parameters.C:5048
void done_reading_files(Bool)
Definition: Parameters.C:2959
int32 atom1
Definition: structures.h:72
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.C:92
int numLJPair
int data_read
Definition: parm.h:34
void print_nbthole_pair_params()
Definition: Parameters.C:4949
int Fclose(FILE *fout)
Definition: common.C:367
char atom4name[11]
Definition: Parameters.C:127
char atom2name[11]
Definition: Parameters.C:125
VdwValue * vdw_array
Definition: Parameters.h:251
#define simParams
Definition: Output.C:127
Index vdw_type
Definition: structures.h:40
StringList * next
Definition: ConfigList.h:49
void resize(int i)
Definition: ResizeArray.h:84
void print_vdw_params()
Definition: Parameters.C:4915
Index index
Definition: Parameters.C:148
void getDihedralParams(int num, Real *c, int *mult, int *funct) const
char * data
Definition: ConfigList.h:48
#define paraXplor
Definition: Parameters.h:78
vector< vector< _REAL > > CMAPParameter
BigReal * table_ener
Definition: Parameters.h:256
_REAL * Cn1
Definition: parm.h:24
Real x0
Definition: Parameters.h:87
void crossterm_setup(CrosstermData *table)
char atom1name[11]
Definition: Parameters.C:158
BondValue * bond_array
Definition: Parameters.h:243
Index index
Definition: Parameters.C:70
void read_ener_table(SimParameters *)
Definition: Parameters.C:6888
void print_improper_params()
Definition: Parameters.C:4898
const BigReal B
int32 atom2
Definition: structures.h:73
char atom2name[11]
Definition: Parameters.C:63
int read_energy_type_bothcubspline(FILE *, const int, BigReal *, const float, const float)
Definition: Parameters.C:7131
void assign_dihedral_index(const char *, const char *, const char *, const char *, Dihedral *, int, int)
Definition: Parameters.C:3956
int Numang
Definition: parm.h:17
int getNumAngleParams() const
_REAL * Req
Definition: parm.h:24
int NumGromacsPairParams
Definition: Parameters.h:268
char atom2name[11]
Definition: Parameters.C:84
struct table_pair_params * next
Definition: Parameters.C:172
void done_reading_structure()
Definition: Parameters.C:4999
Real epsilon14
Definition: Parameters.h:155
char atom2name[11]
Definition: Parameters.C:48
struct vdw_pair_params * next
Definition: Parameters.C:164
void assign_angle_index(const char *, const char *, const char *, Angle *, int)
Definition: Parameters.C:3854
MOStream * put(char data)
Definition: MStream.h:112
int getNumBondParams() const
Real theta0
Definition: Parameters.h:94
void NAMD_remove_comment(char *str)
Definition: strlib.C:119
struct indexed_nbthole_pair * left
Definition: Parameters.h:182
int tablenumtypes
Definition: Parameters.h:260
int size(void) const
Definition: ResizeArray.h:127
void print_angle_params()
Definition: Parameters.C:4865
Real epsilon
Definition: Parameters.C:145
Real sigma
Definition: Parameters.C:144
valT values[SORTTILELISTSKERNEL_ITEMS_PER_THREAD]
Bool tabulatedEnergies
Index bond_type
Definition: structures.h:50
gridSize x
char atom1name[11]
Definition: Parameters.C:106
IndexedTablePair * tab_pair_tree
Definition: Parameters.h:259
int32 atom3
Definition: structures.h:74
#define TRUE
Definition: common.h:119
char atom6name[11]
Definition: Parameters.C:129
struct angle_params * right
Definition: Parameters.C:72
int getNumPair() const
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.C:111
Real k_ub
Definition: Parameters.h:95
int Nptra
Definition: parm.h:17
Real k
Definition: Parameters.h:86
int32 atom1
Definition: structures.h:63
int32 atom2
Definition: structures.h:49
char tableInterpType[128]
struct vdw_params * right
Definition: Parameters.C:150
#define MAX_ATOMTYPE_CHARS
Definition: Parameters.h:73
int * Cno
Definition: parm.h:27
_REAL * Pn
Definition: parm.h:24
void receive_Parameters(MIStream *)
Definition: Parameters.C:5424
double BigReal
Definition: common.h:114
void getAngleParams(int num, Real *th0, Real *kth, int *funct) const