31 #define MIN_DEBUG_LEVEL 3
35 #define INDEX(ncols,i,j) ((i)*ncols + (j))
190 void Parameters::initialize() {
217 maxDihedralMults=NULL;
218 maxImproperMults=NULL;
271 CkPrintf(
"Working on tables\n");
279 AllFilesRead =
FALSE;
296 }
while ( f != NULL );
317 delete [] atomTypeNames;
320 free_bond_tree(bondp);
323 free_angle_tree(anglep);
325 if (dihedralp != NULL)
326 free_dihedral_list(dihedralp);
328 if (improperp != NULL)
329 free_improper_list(improperp);
331 if (crosstermp != NULL)
332 free_crossterm_list(crosstermp);
337 if (vdw_pairp != NULL)
338 free_vdw_pair_list();
340 if (nbthole_pairp != NULL)
341 free_nbthole_pair_list();
375 if (maxDihedralMults != NULL)
376 delete [] maxDihedralMults;
378 if (maxImproperMults != NULL)
379 delete [] maxImproperMults;
381 for(
int i = 0; i < error_msgs.
size(); ++i ) {
382 delete [] error_msgs[i];
408 char first_word[512];
415 NAMD_die(
"Tried to read another parameter file after being told that all files were read!");
419 if ( (pfile =
Fopen(fname,
"r")) == NULL)
423 sprintf(err_msg,
"UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
439 if ((buffer[0] !=
'!') &&
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) )
458 if (strncasecmp(first_word,
"bond", 4)==0)
460 add_bond_param(buffer,
TRUE);
463 else if (strncasecmp(first_word,
"angl", 4)==0)
465 add_angle_param(buffer);
468 else if (strncasecmp(first_word,
"dihe", 4)==0)
470 add_dihedral_param(buffer, pfile);
473 else if (strncasecmp(first_word,
"impr", 4)==0)
475 add_improper_param(buffer, pfile);
478 else if (strncasecmp(first_word,
"nonb", 4)==0)
480 add_vdw_param(buffer);
483 else if (strncasecmp(first_word,
"nbfi", 4)==0)
485 add_vdw_pair_param(buffer);
488 else if (strncasecmp(first_word,
"nbta", 4)==0)
495 add_table_pair_param(buffer);
498 else if (strncasecmp(first_word,
"hbon", 4)==0)
500 add_hb_pair_param(buffer);
508 sprintf(err_msg,
"UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
547 char first_word[512];
554 NAMD_die(
"Tried to read another parameter file after being told that all files were read!");
558 if ( (pfile = fopen(fname,
"r")) == NULL)
562 sprintf(err_msg,
"UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
577 (strncmp(first_word,
"!", 1) != 0) &&
578 (strncmp(first_word,
"*", 1) != 0) &&
579 (strncasecmp(first_word,
"END", 3) != 0))
582 iout <<
iWARN <<
"SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" <<
endi;
586 if (strncasecmp(first_word,
"bond", 4)==0)
588 par_type=1; skipline=1;
590 else if (strncasecmp(first_word,
"angl", 4)==0)
592 par_type=2; skipline=1;
594 else if (strncasecmp(first_word,
"thet", 4)==0)
596 par_type=2; skipline=1;
598 else if (strncasecmp(first_word,
"dihe", 4)==0)
600 par_type=3; skipline=1;
602 else if (strncasecmp(first_word,
"phi", 3)==0)
604 par_type=3; skipline=1;
606 else if (strncasecmp(first_word,
"impr", 4)==0)
608 par_type=4; skipline=1;
610 else if (strncasecmp(first_word,
"imph", 4)==0)
612 par_type=4; skipline=1;
614 else if (strncasecmp(first_word,
"nbon", 4)==0)
616 par_type=5; skipline=1;
618 else if (strncasecmp(first_word,
"nonb", 4)==0)
620 par_type=5; skipline=1;
622 else if (strncasecmp(first_word,
"nbfi", 4)==0)
624 par_type=6; skipline=1;
626 else if (strncasecmp(first_word,
"hbon", 4)==0)
628 par_type=7; skipline=1;
630 else if (strncasecmp(first_word,
"cmap", 4)==0)
632 par_type=8; skipline=1;
634 else if (strncasecmp(first_word,
"nbta", 4)==0)
636 par_type=9; skipline=1;
638 else if (strncasecmp(first_word,
"thol", 4)==0)
640 par_type=10; skipline=1;
642 else if (strncasecmp(first_word,
"atom", 4)==0)
644 par_type=11; skipline=1;
646 else if (strncasecmp(first_word,
"ioformat", 8)==0)
650 else if (strncasecmp(first_word,
"read", 4)==0)
652 skip_stream_read(buffer, pfile); skipline=1;
654 else if (strncasecmp(first_word,
"return", 4)==0)
656 skipall=1; skipline=1;
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) )
684 if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
688 sprintf(err_msg,
"ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
696 else if (par_type == 0)
702 sprintf(err_msg,
"UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
711 if ( (par_type != 0) && (!skipline) )
718 add_bond_param(buffer,
TRUE);
721 else if (par_type == 2)
723 add_angle_param(buffer);
726 else if (par_type == 3)
728 add_dihedral_param(buffer, pfile);
731 else if (par_type == 4)
733 add_improper_param(buffer, pfile);
736 else if (par_type == 5)
738 add_vdw_param(buffer);
741 else if (par_type == 6)
743 add_vdw_pair_param(buffer);
746 else if (par_type == 7)
748 add_hb_pair_param(buffer);
750 else if (par_type == 8)
752 add_crossterm_param(buffer, pfile);
755 else if (par_type == 9)
762 add_table_pair_param(buffer);
765 else if (par_type == 10)
767 add_nbthole_pair_param(buffer);
770 else if (par_type == 11)
772 if ( strncasecmp(first_word,
"mass", 4) != 0 ) {
774 sprintf(err_msg,
"UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
785 sprintf(err_msg,
"INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
800 void Parameters::skip_stream_read(
char *buf, FILE *fd) {
803 char first_word[513];
806 int rval = sscanf(buf,
"%s %s", s1, s2);
809 sprintf(err_msg,
"BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
812 if ( ! strncasecmp(s2,
"PARA",4) )
return;
814 iout <<
iINFO <<
"SKIPPING " << s2 <<
" SECTION IN STREAM FILE\n" <<
endi;
821 (strncmp(first_word,
"!", 1) != 0) &&
822 (strncmp(first_word,
"*", 1) != 0) &&
823 (strncasecmp(first_word,
"END", 3) == 0) )
return;
842 void Parameters::add_bond_param(
const char *buf,
Bool overwrite)
857 read_count=sscanf(buf,
"%*s %s %s %f %f\n", atom1name, atom2name,
858 &forceconstant, &distance);
863 read_count=sscanf(buf,
"%s %s %f %f\n", atom1name, atom2name,
864 &forceconstant, &distance);
875 sprintf(err_msg,
"BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
879 sprintf(err_msg,
"BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
887 if (new_node == NULL)
889 NAMD_die(
"memory allocation failed in Parameters::add_bond_param\n");
896 if (strcasecmp(atom1name, atom2name) < 0)
912 new_node->
left = NULL;
913 new_node->
right = NULL;
917 bondp=add_to_bond_tree(new_node, bondp, overwrite);
957 if (compare_code == 0)
963 if (compare_code == 0)
974 iout <<
"\n" <<
iWARN <<
"DUPLICATE BOND ENTRY FOR "
997 if (compare_code < 0)
999 tree->
left = add_to_bond_tree(new_node, tree->
left, overwrite);
1003 tree->
right = add_to_bond_tree(new_node, tree->
right, overwrite);
1023 void Parameters::add_angle_param(
char *buf)
1042 read_count=sscanf(buf,
"%*s %s %s %s %f %f UB %f %f\n",
1043 atom1name, atom2name, atom3name, &forceconstant, &angle,
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,
1056 read_count=sscanf(buf,
"%s %s %s %f %f %f %f\n",
1057 atom1name, atom2name, atom3name, &forceconstant, &angle,
1065 if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
1071 sprintf(err_msg,
"BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1075 sprintf(err_msg,
"BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1083 if (new_node == NULL)
1085 NAMD_die(
"memory allocation failed in Parameters::add_angle_param");
1091 if (strcasecmp(atom1name, atom3name) < 0)
1109 if (strcasecmp(
"cos",norm)==0) {
1121 if (read_count == 7)
1124 if (new_node->
normal == 0) {
1125 NAMD_die(
"ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
1132 new_node->
k_ub = 0.0;
1133 new_node->
r_ub = 0.0;
1136 new_node->
left = NULL;
1137 new_node->
right = NULL;
1140 anglep = add_to_angle_tree(new_node, anglep);
1179 if (compare_code == 0)
1184 if (compare_code == 0)
1187 compare_code = strcasecmp(new_node->
atom3name,
1190 if (compare_code == 0)
1203 iout <<
"\n" <<
iWARN <<
"DUPLICATE ANGLE ENTRY FOR "
1207 <<
"\nPREVIOUS VALUES k="
1209 << tree->
angle <<
" k_ub="
1210 << tree->
k_ub <<
" r_ub="
1212 <<
"\n USING VALUES k="
1214 << new_node->
angle <<
" k_ub="
1215 << new_node->
k_ub <<
" r_ub=" << new_node->
r_ub
1236 if (compare_code < 0)
1238 tree->
left = add_to_angle_tree(new_node, tree->
left);
1242 tree->
right = add_to_angle_tree(new_node, tree->
right);
1262 void Parameters::add_dihedral_param(
char *buf, FILE *fd)
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);
1291 read_count=sscanf(buf,
"%s %s %s %s %f %d %f\n",
1292 atom1name, atom2name, atom3name, atom4name,
1293 &forceconstant, &periodicity, &phase_shift);
1297 if ( (read_count != 4) && (read_count != 8) && (paramType ==
paraXplor) )
1301 sprintf(err_msg,
"BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1304 else if ( (read_count != 7) && (paramType ==
paraCharmm) )
1308 sprintf(err_msg,
"BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1312 if ( (read_count == 4) && (paramType ==
paraXplor) )
1315 read_count=sscanf(buf,
"%*s %*s %*s %*s %*s %f %d %f\n",
1316 &forceconstant, &periodicity, &phase_shift);
1319 if (read_count != 3)
1323 sprintf(err_msg,
"BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1334 sprintf(err_msg,
"Multiple dihedral with multiplicity of %d greater than max of %d",
1342 if (new_node == NULL)
1344 NAMD_die(
"memory allocation failed in Parameters::add_dihedral_param\n");
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");
1360 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
1361 new_node->
values[0].
k = forceconstant;
1362 new_node->
values[0].
n = periodicity;
1365 new_node->
next = NULL;
1368 if (multiplicity > 1)
1388 NAMD_die(
"EOF encoutner in middle of multiple dihedral");
1391 read_count=sscanf(buffer,
"%f %d %f\n",
1392 &forceconstant, &periodicity, &phase_shift);
1394 if (read_count != 3)
1398 sprintf(err_msg,
"BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
1402 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
1403 new_node->
values[i].
k = forceconstant;
1404 new_node->
values[i].
n = periodicity;
1413 add_to_dihedral_list(new_node);
1417 add_to_charmm_dihedral_list(new_node);
1443 void Parameters::add_to_dihedral_list(
1455 if (dihedralp == NULL)
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;}
1496 iout <<
"\n" <<
iWARN <<
"DUPLICATE DIHEDRAL ENTRY FOR "
1501 <<
"\nPREVIOUS VALUES MULTIPLICITY " << ptr->
multiplicity <<
"\n";
1515 <<
" n=" << new_node->
values[i].
n
1553 tail->
next=new_node;
1561 new_node->
next=dihedralp;
1592 void Parameters::add_to_charmm_dihedral_list(
1609 if (dihedralp == NULL)
1623 int same_as_last = 0;
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)))
1666 if (!same_as_last) {
1667 iout <<
"\n" <<
iWARN <<
"DUPLICATE DIHEDRAL ENTRY FOR "
1672 <<
"\nPREVIOUS VALUES MULTIPLICITY: " << ptr->
multiplicity <<
"\n";
1678 if (!same_as_last) {
1685 iout <<
iWARN <<
"IDENTICAL PERIODICITY! REPLACING OLD VALUES BY: \n";
1704 sprintf(err_msg,
"Multiple dihedral with multiplicity of %d greater than max of %d",
1748 tail->
next=new_node;
1757 new_node->
next=dihedralp;
1781 void Parameters::add_improper_param(
char *buf, FILE *fd)
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);
1810 read_count=sscanf(buf,
"%s %s %s %s %f %d %f\n",
1811 atom1name, atom2name, atom3name, atom4name,
1812 &forceconstant, &periodicity, &phase_shift);
1816 if ( (read_count != 4) && (read_count != 8) && (paramType ==
paraXplor) )
1820 sprintf(err_msg,
"BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
1823 else if ( (read_count != 7) && (paramType ==
paraCharmm) )
1827 sprintf(err_msg,
"BAD IMPROPER FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
1831 if ( (read_count == 4) && (paramType ==
paraXplor) )
1834 read_count=sscanf(buf,
"%*s %*s %*s %*s %*s %f %d %f\n",
1835 &forceconstant, &periodicity, &phase_shift);
1838 if (read_count != 3)
1842 sprintf(err_msg,
"BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1853 sprintf(err_msg,
"Multiple improper with multiplicity of %d greater than max of %d",
1861 if (new_node == NULL)
1863 NAMD_die(
"memory allocation failed in Parameters::add_improper_param");
1873 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
1874 new_node->
values[0].
k = forceconstant;
1875 new_node->
values[0].
n = periodicity;
1878 new_node->
next = NULL;
1881 if (multiplicity > 1)
1902 NAMD_die(
"EOF encoutner in middle of multiple improper");
1906 read_count=sscanf(buffer,
"%f %d %f\n",
1907 &forceconstant, &periodicity, &phase_shift);
1909 if (read_count != 3)
1913 sprintf(err_msg,
"BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
1917 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
1918 new_node->
values[i].
k = forceconstant;
1919 new_node->
values[i].
n = periodicity;
1925 add_to_improper_list(new_node);
1949 void Parameters::add_to_improper_list(
struct improper_params *new_node)
1960 if (improperp == NULL)
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;}
2001 iout <<
"\n" <<
iWARN <<
"DUPLICATE IMPROPER DIHEDRAL ENTRY FOR "
2006 <<
"\nPREVIOUS VALUES MULTIPLICITY " << ptr->
multiplicity <<
"\n";
2015 iout <<
"\n" <<
"USING VALUES MULTIPLICITY " << new_node->
multiplicity <<
"\n";
2020 <<
" n=" << new_node->
values[i].
n
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) )
2057 tail->
next=new_node;
2065 new_node->
next=improperp;
2086 void Parameters::add_crossterm_param(
char *buf, FILE *fd)
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,
2109 if ( (read_count != 9) || dimension < 1 || dimension > 1000 )
2113 sprintf(err_msg,
"BAD CMAP FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2131 new_node->
next = NULL;
2137 while ( nread < nterms ) {
2141 if (ret_code == 0) {
2148 if (ret_code == 0) {
2153 if (ret_code != 0) {
2154 NAMD_die(
"EOF encoutner in middle of CMAP");
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);
2168 nread += read_count;
2170 if (read_count == 0 || nread > nterms) {
2173 sprintf(err_msg,
"BAD CMAP FORMAT IN PARAMETER FILE\nLINE=*%s*\n", buffer);
2179 add_to_crossterm_list(new_node);
2214 if (crosstermp == NULL)
2216 crosstermp=new_node;
2245 for (i=0; i<nvals; i++)
2247 if (ptr->
values[i] != new_node->
values[i]) {echoWarn=1;
break;}
2253 iout <<
"\n" <<
iWARN <<
"DUPLICATE CMAP ENTRY FOR "
2261 << ptr->
atom8name <<
", USING NEW VALUES\n";
2269 new_node->
values = tmpvalues;
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) )
2296 tail->
next=new_node;
2304 new_node->
next=crosstermp;
2305 crosstermp=new_node;
2324 void Parameters::add_vdw_param(
char *buf)
2341 read_count=sscanf(buf,
"%*s %s %f %f %f %f\n", atomname,
2342 &epsilon, &sigma, &epsilon14, &sigma14);
2347 read_count=sscanf(buf,
"%s %*f %f %f %*f %f %f\n", atomname,
2348 &epsilon, &sigma, &epsilon14, &sigma14);
2352 if ((read_count != 5) && (paramType ==
paraXplor))
2356 sprintf(err_msg,
"BAD vdW FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
2359 else if ( ((read_count != 5) && (read_count != 3)) && (paramType ==
paraCharmm))
2363 sprintf(err_msg,
"BAD vdW FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
2371 sqrt26=pow(2.,(1./6.));
2372 sigma=2.*sigma/sqrt26;
2374 if (read_count == 3)
2382 sigma14=2.*sigma14/sqrt26;
2387 if ( epsilon < 0. || epsilon14 < 0. ) {
2388 iout <<
iWARN <<
"Ignoring VDW parameter with negative epsilon:\n"
2389 << buf <<
"\n" <<
endi;
2396 if (new_node == NULL)
2398 NAMD_die(
"memory allocation failed in Parameters::add_vdw_param");
2402 strcpy(new_node->
atomname, atomname);
2408 new_node->
left = NULL;
2409 new_node->
right = NULL;
2412 vdwp=add_to_vdw_tree(new_node, vdwp);
2447 if (compare_code==0)
2458 <<
"\nPREVIOUS VALUES sigma=" << tree->
sigma
2459 <<
" epsilon=" << tree->
epsilon
2460 <<
" sigma14=" << tree->
sigma14
2462 <<
"\n USING VALUES sigma=" << new_node->
sigma
2463 <<
" epsilon=" << new_node->
epsilon
2464 <<
" sigma14=" << new_node->
sigma14
2482 if (compare_code < 0)
2484 tree->
left = add_to_vdw_tree(new_node, tree->
left);
2488 tree->
right = add_to_vdw_tree(new_node, tree->
right);
2507 void Parameters::add_table_pair_param(
char *buf)
2518 read_count=sscanf(buf,
"%s %s %s\n", atom1name,
2519 atom2name, tabletype);
2522 if ((read_count != 3))
2526 sprintf(err_msg,
"BAD TABLE PAIR FORMAT IN PARAMETER FILE\nLINE=*%s*", buf);
2533 if (new_node == NULL)
2535 NAMD_die(
"memory allocation failed in Parameters::add_table_pair_param\n");
2546 sprintf(err_msg,
"Couldn't find table parameters for table interaction %s!\n", tabletype);
2553 new_node->
next = NULL;
2557 add_to_table_pair_list(new_node);
2575 void Parameters::add_vdw_pair_param(
char *buf)
2591 read_count=sscanf(buf,
"%*s %s %s %f %f %f %f\n", atom1name,
2592 atom2name, &A, &B, &A14, &B14);
2596 Real well, rmin, well14, rmin14;
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.);
2608 if ((read_count != 6) && (paramType ==
paraXplor))
2612 sprintf(err_msg,
"BAD vdW PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
2615 if ((read_count != 4) && (read_count != 6) && (paramType ==
paraCharmm))
2619 sprintf(err_msg,
"BAD vdW PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2627 if (new_node == NULL)
2629 NAMD_die(
"memory allocation failed in Parameters::add_vdw_pair_param\n");
2641 new_node->
next = NULL;
2644 add_to_vdw_pair_list(new_node);
2662 void Parameters::add_nbthole_pair_param(
char *buf)
2675 read_count=sscanf(buf,
"%s %s %f\n", atom1name,
2676 atom2name, &tholeij);
2680 if ((read_count != 3) && (paramType ==
paraCharmm))
2684 sprintf(err_msg,
"BAD NBTHOLE PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2692 if (new_node == NULL)
2694 NAMD_die(
"memory allocation failed in Parameters::nbthole_pair_param\n");
2703 new_node->
next = NULL;
2706 add_to_nbthole_pair_list(new_node);
2724 void Parameters::add_hb_pair_param(
char *buf)
2736 if (sscanf(buf,
"%*s %s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
2738 sprintf(err_msg,
"BAD HBOND PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
2743 if (sscanf(buf,
"%s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
2745 sprintf(err_msg,
"BAD HBOND PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
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;
2780 if (table_pairp == NULL)
2782 table_pairp = new_node;
2795 if (compare_code == 0)
2800 if (compare_code==0)
2805 iout <<
iWARN <<
"DUPLICATE TABLE PAIR ENTRY FOR "
2823 tail->
next = new_node;
2839 void Parameters::add_to_vdw_pair_list(
struct vdw_pair_params *new_node)
2848 if (vdw_pairp == NULL)
2850 vdw_pairp = new_node;
2863 if (compare_code == 0)
2868 if (compare_code==0)
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))
2878 iout <<
iWARN <<
"DUPLICATE vdW PAIR ENTRY FOR "
2881 <<
"\nPREVIOUS VALUES A=" << ptr->
A
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
2908 tail->
next = new_node;
2933 if (nbthole_pairp == NULL)
2935 nbthole_pairp = new_node;
2940 ptr = nbthole_pairp;
2942 tail->
next = new_node;
2962 AllFilesRead =
TRUE;
2967 add_bond_param(
"X DRUD 500.0 0.0\n",
FALSE);
2976 NAMD_die(
"memory allocation of bond_array failed!");
2987 NAMD_die(
"memory allocation of angle_array failed!");
3001 NAMD_die(
"memory allocation of dihedral_array failed!");
3012 NAMD_die(
"memory allocation of improper_array failed!");
3036 if (vdw_array == NULL)
3038 NAMD_die(
"memory allocation of vdw_array failed!");
3047 NAMD_die(
"memory allocation of nbthole_array failed!");
3061 index_bonds(bondp, 0);
3062 index_angles(anglep, 0);
3067 convert_nbthole_pairs();
3069 convert_vdw_pairs();
3070 convert_table_pairs();
3097 if (tree->
left != NULL)
3099 index=index_bonds(tree->
left, index);
3109 if (tree->
right != NULL)
3111 index=index_bonds(tree->
right, index);
3141 if (tree->
left != NULL)
3143 index=index_angles(tree->
left, index);
3159 if (tree->
right != NULL)
3161 index=index_angles(tree->
right, index);
3178 void Parameters::index_dihedrals()
3203 if (maxDihedralMults == NULL)
3205 NAMD_die(
"memory allocation failed in Parameters::index_dihedrals()");
3262 void Parameters::index_impropers()
3286 if (maxImproperMults == NULL)
3288 NAMD_die(
"memory allocation failed in Parameters::index_impropers()");
3348 void Parameters::index_crossterms()
3365 NAMD_die(
"Sorry, only CMAP dimension of 24 is supported");
3369 for (i=0; i<N; i++) {
3370 for (j=0; j<N; j++) {
3375 for (i=0; i<N; i++) {
3417 if (tree->
left != NULL)
3419 index=index_vdw(tree->
left, index);
3440 if (tree->
right != NULL)
3442 index=index_vdw(tree->
right, index);
3467 #ifdef MEM_OPT_VERSION
3480 NAMD_die(
"Tried to assign vdw index before all parameter files were read");
3490 while (!found && (ptr!=NULL))
3492 comp_code = strcasecmp(atomtype, ptr->
atomname);
3500 else if (comp_code < 0)
3523 while (!found && (ptr!=NULL))
3528 if (windx == strlen(ptr->
atomname))
3531 comp_code = strcasecmp(atomtype, ptr->
atomname);
3535 comp_code = strncasecmp(atomtype, ptr->
atomname, windx);
3544 sprintf(errbuf,
"VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
3547 for(i=0; i<error_msgs.size(); i++) {
3548 if ( strcmp(errbuf,error_msgs[i]) == 0 )
break;
3550 if ( i == error_msgs.size() ) {
3551 char *newbuf =
new char[strlen(errbuf)+1];
3552 strcpy(newbuf,errbuf);
3553 error_msgs.add(newbuf);
3557 else if (comp_code < 0)
3578 sprintf(err_msg,
"DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s",
3623 while (!found && (ptr!=NULL))
3626 if ( (ind1 == ptr->
ind1) && (ind2 == ptr->
ind2) )
3630 else if ( (ind1 < ptr->ind1) ||
3631 ( (ind1==ptr->
ind1) && (ind2 < ptr->ind2) ) )
3702 while (!found && (ptr!=NULL))
3704 if ( (ind1 == ptr->
ind1) && (ind2 == ptr->
ind2) )
3708 else if ( (ind1 < ptr->ind1) ||
3709 ( (ind1==ptr->
ind1) && (ind2 < ptr->ind2) ) )
3768 NAMD_die(
"Tried to assign bond index before all parameter files were read");
3773 if (strcasecmp(atom1, atom2) > 0)
3775 const char *tmp_name = atom1;
3785 while (!found && (ptr!=NULL))
3787 cmp_code=strcasecmp(atom1, ptr->
atom1name);
3791 cmp_code=strcasecmp(atom2, ptr->
atom2name);
3800 else if (cmp_code < 0)
3815 if ((strcmp(atom1,
"DRUD")==0 || strcmp(atom2,
"DRUD")==0)
3816 && (strcmp(atom1,
"X")!=0 && strcmp(atom2,
"X")!=0)) {
3818 char a1[8] =
"DRUD", a2[8] =
"X";
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);
3855 Angle *angle_ptr,
int notFoundIndex)
3865 NAMD_die(
"Tried to assign angle index before all parameter files were read");
3870 if (strcasecmp(atom1, atom3) > 0)
3872 const char *tmp_name = atom1;
3882 while (!found && (ptr != NULL))
3884 comp_val = strcasecmp(atom1, ptr->
atom1name);
3889 comp_val = strcasecmp(atom2, ptr->
atom2name);
3894 comp_val = strcasecmp(atom3, ptr->
atom3name);
3904 else if (comp_val < 0)
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);
3924 if ( notFoundIndex ) {
3957 const char *atom4,
Dihedral *dihedral_ptr,
3958 int multiplicity,
int notFoundIndex)
3969 while (!found && (ptr!=NULL))
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);
4011 if ( notFoundIndex ) {
4022 if (multiplicity > maxDihedralMults[ptr->
index])
4026 sprintf(err_msg,
"Multiplicity of Paramters for diehedral bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxDihedralMults[ptr->
index]);
4066 const char *atom4,
Improper *improper_ptr,
4078 while (!found && (ptr!=NULL))
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) ) )
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) ) )
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);
4135 if (multiplicity > maxImproperMults[ptr->
index])
4139 sprintf(err_msg,
"Multiplicity of Paramters for improper bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxImproperMults[ptr->
index]);
4165 const char *atom4,
const char *atom5,
const char *atom6,
const char *atom7,
4166 const char *atom8,
Crossterm *crossterm_ptr)
4176 while (!found && (ptr!=NULL))
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) ) )
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) ) )
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) ) )
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) ) )
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);
4281 void Parameters::free_bond_tree(
struct bond_params *bond_ptr)
4284 if (bond_ptr->
left != NULL)
4286 free_bond_tree(bond_ptr->
left);
4289 if (bond_ptr->
right != NULL)
4291 free_bond_tree(bond_ptr->
right);
4314 void Parameters::free_angle_tree(
struct angle_params *angle_ptr)
4317 if (angle_ptr->
left != NULL)
4319 free_angle_tree(angle_ptr->
left);
4322 if (angle_ptr->
right != NULL)
4324 free_angle_tree(angle_ptr->
right);
4440 void Parameters::free_vdw_tree(
struct vdw_params *vdw_ptr)
4443 if (vdw_ptr->
left != NULL)
4445 free_vdw_tree(vdw_ptr->
left);
4448 if (vdw_ptr->
right != NULL)
4450 free_vdw_tree(vdw_ptr->
right);
4467 void Parameters::free_vdw_pair_list()
4494 void Parameters::free_nbthole_pair_list()
4509 nbthole_pairp = NULL;
4520 if (table_pair_ptr->
left != NULL)
4522 free_table_pair_tree(table_pair_ptr->
left);
4525 if (table_pair_ptr->
right != NULL)
4527 free_table_pair_tree(table_pair_ptr->
right);
4530 delete table_pair_ptr;
4550 void Parameters::free_vdw_pair_tree(
IndexedVdwPair *vdw_pair_ptr)
4553 if (vdw_pair_ptr->
left != NULL)
4555 free_vdw_pair_tree(vdw_pair_ptr->
left);
4558 if (vdw_pair_ptr->
right != NULL)
4560 free_vdw_pair_tree(vdw_pair_ptr->
right);
4563 delete vdw_pair_ptr;
4586 if (nbthole_pair_ptr->
left != NULL)
4588 free_nbthole_pair_tree(nbthole_pair_ptr->
left);
4591 if (nbthole_pair_ptr->
right != NULL)
4593 free_nbthole_pair_tree(nbthole_pair_ptr->
right);
4596 delete nbthole_pair_ptr;
4615 void Parameters::traverse_bond_params(
struct bond_params *tree)
4621 if (tree->
left != NULL)
4623 traverse_bond_params(tree->
left);
4630 if (tree->
right != NULL)
4632 traverse_bond_params(tree->
right);
4650 void Parameters::traverse_angle_params(
struct angle_params *tree)
4656 if (tree->
left != NULL)
4658 traverse_angle_params(tree->
left);
4665 if (tree->
right != NULL)
4667 traverse_angle_params(tree->
right);
4685 void Parameters::traverse_dihedral_params(
struct dihedral_params *list)
4690 while (list != NULL)
4694 <<
" " << list->
atom4name <<
" index=" \
4701 <<
" n=" << list->
values[i].
n \
4723 void Parameters::traverse_improper_params(
struct improper_params *list)
4728 while (list != NULL)
4732 <<
" " << list->
atom4name <<
" index=" \
4739 <<
" n=" << list->
values[i].
n \
4762 void Parameters::traverse_vdw_params(
struct vdw_params *tree)
4768 if (tree->
left != NULL)
4770 traverse_vdw_params(tree->
left);
4774 <<
" sigma=" << tree->
sigma <<
" epsilon=" << \
4776 <<
" epsilon 1-4=" << tree->
epsilon14 << std::endl);
4778 if (tree->
right != NULL)
4780 traverse_vdw_params(tree->
right);
4797 void Parameters::traverse_vdw_pair_params(
struct vdw_pair_params *list)
4805 <<
" B=" << list->
B <<
" A 1-4=" \
4806 << list->
A14 <<
" B 1-4=" << list->
B14 \
4809 traverse_vdw_pair_params(list->
next);
4834 traverse_nbthole_pair_params(list->
next);
4850 <<
"*****************************************" \
4853 traverse_bond_params(bondp);
4868 <<
"*****************************************" );
4869 traverse_angle_params(anglep);
4884 <<
"*****************************************" );
4886 traverse_dihedral_params(dihedralp);
4901 <<
"*****************************************" );
4903 traverse_improper_params(improperp);
4918 <<
"*****************************************" );
4920 traverse_vdw_params(vdwp);
4935 <<
"*****************************************" );
4937 traverse_vdw_pair_params(vdw_pairp);
4952 <<
"*****************************************" );
4954 traverse_nbthole_pair_params(nbthole_pairp);
4968 iout <<
iINFO <<
"SUMMARY OF PARAMETERS:\n"
5003 free_bond_tree(bondp);
5006 free_angle_tree(anglep);
5008 if (dihedralp != NULL)
5009 free_dihedral_list(dihedralp);
5011 if (improperp != NULL)
5012 free_improper_list(improperp);
5014 if (crosstermp != NULL)
5015 free_crossterm_list(crosstermp);
5018 free_vdw_tree(vdwp);
5022 if (maxDihedralMults != NULL)
5023 delete [] maxDihedralMults;
5025 if (maxImproperMults != NULL)
5026 delete [] maxImproperMults;
5034 maxImproperMults=NULL;
5035 maxDihedralMults=NULL;
5050 Real *a1, *a2, *a3, *a4;
5071 if ( (a1 == NULL) || (a2 == NULL) )
5073 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5082 msg->
put(NumBondParams, a1)->
put(NumBondParams, a2);
5099 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) ||
5100 (a4 == NULL) || (i1 == NULL))
5102 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5114 msg->
put(NumAngleParams, a1)->
put(NumAngleParams, a2);
5115 msg->
put(NumAngleParams, a3)->
put(NumAngleParams, a4);
5116 msg->
put(NumAngleParams, i1);
5135 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5136 (deltavals == NULL) )
5138 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5147 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5149 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5165 msg->
put(NumDihedralParams, i1);
5169 msg->
put(NumDihedralParams, kvals[i]);
5170 msg->
put(NumDihedralParams, nvals[i]);
5171 msg->
put(NumDihedralParams, deltavals[i]);
5175 delete [] deltavals[i];
5181 delete [] deltavals;
5194 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5195 (deltavals == NULL) )
5197 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5206 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5208 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5224 msg->
put(NumImproperParams, i1);
5228 msg->
put(NumImproperParams, kvals[i]);
5229 msg->
put(NumImproperParams, nvals[i]);
5230 msg->
put(NumImproperParams, deltavals[i]);
5234 delete [] deltavals[i];
5240 delete [] deltavals;
5261 if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
5262 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5270 msg->
put(NumGromacsPairParams,pairC6);
5271 msg->
put(NumGromacsPairParams,pairC12);
5305 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) )
5307 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5319 msg->
put(NumVdwParams, a1);
5320 msg->
put(NumVdwParams, a2);
5321 msg->
put(NumVdwParams, a3);
5322 msg->
put(NumVdwParams, a4);
5342 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
5343 (i1 == NULL) || (i2 == NULL) )
5345 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5348 vdw_pair_to_arrays(i1, i2, a1, a2, a3, a4, 0,
vdw_pair_tree);
5367 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (i1 == NULL) || (i2 == NULL) )
5369 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5383 msg->
put(NumNbtholePairParams, i1)->
put(NumNbtholePairParams, i2);
5384 msg->
put(NumNbtholePairParams, a1);
5385 msg->
put(NumNbtholePairParams, a2)->
put(NumNbtholePairParams, a3);
5398 if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
5400 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5428 Real *a1, *a2, *a3, *a4;
5446 if ( (
bond_array == NULL) || (a1 == NULL) || (a2 == NULL) )
5448 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
5476 if ( (
angle_array == NULL) || (a1 == NULL) || (a2 == NULL) ||
5477 (a3 == NULL) || (a4 == NULL) || (i1 == NULL))
5479 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
5516 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5519 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5528 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5530 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5559 delete [] deltavals[i];
5565 delete [] deltavals;
5579 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
5582 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5591 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
5593 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5622 delete [] deltavals[i];
5628 delete [] deltavals;
5654 if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
5655 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
5680 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
5699 if ( (vdw_array==NULL) || (a1==NULL) || (a2==NULL) || (a3==NULL)
5700 || (a4==NULL) || (atomTypeNames==NULL))
5702 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
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];
5737 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
5738 (i1 == NULL) || (i2 == NULL) )
5740 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5754 if (new_node == NULL)
5756 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
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;
5791 if ( (
nbthole_array == NULL) || (a1 == NULL) || (a2 == NULL) || (a3 == NULL)
5792 || (i1 == NULL) || (i2 == NULL) )
5794 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
5830 if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
5832 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
5843 if (new_tab_node == NULL)
5845 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
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;
5866 AllFilesRead =
TRUE;
5883 void Parameters::convert_vdw_pairs()
5886 #ifdef MEM_OPT_VERSION
5891 Index index1, index2;
5903 if (new_node == NULL)
5905 NAMD_die(
"memory allocation failed in Parameters::convert_vdw_pairs");
5915 if (index1 > index2)
5917 new_node->
ind1 = index2;
5918 new_node->
ind2 = index1;
5922 new_node->
ind1 = index1;
5923 new_node->
ind2 = index2;
5926 new_node->
A = ptr->
A;
5927 new_node->
B = ptr->
B;
5928 new_node->
A14 = ptr->
A14;
5929 new_node->
B14 = ptr->
B14;
5931 new_node->
left = NULL;
5932 new_node->
right = NULL;
5961 void Parameters::convert_nbthole_pairs()
5964 #ifdef MEM_OPT_VERSION
5969 Index index1, index2;
5973 ptr = nbthole_pairp;
5981 if (new_node == NULL)
5983 NAMD_die(
"memory allocation failed in Parameters::convert_nbthole_pairs");
5993 if (index1 > index2)
5995 new_node->
ind1 = index2;
5996 new_node->
ind2 = index1;
6000 new_node->
ind1 = index1;
6001 new_node->
ind2 = index2;
6008 new_node->
left = NULL;
6009 new_node->
right = NULL;
6022 nbthole_pairp = NULL;
6038 void Parameters::convert_table_pairs()
6041 #ifdef MEM_OPT_VERSION
6046 Index index1, index2;
6058 if (new_node == NULL)
6060 NAMD_die(
"memory allocation failed in Parameters::convert_table_pairs");
6072 if (index1 > index2)
6074 new_node->
ind1 = index2;
6075 new_node->
ind2 = index1;
6079 new_node->
ind1 = index1;
6080 new_node->
ind2 = index2;
6083 new_node->
K = ptr->
K;
6085 new_node->
left = NULL;
6086 new_node->
right = NULL;
6124 if ( (new_node->
ind1 < tree->
ind1) ||
6127 tree->
left = add_to_indexed_table_pairs(new_node, tree->
left);
6131 tree->
right = add_to_indexed_table_pairs(new_node, tree->
right);
6158 if ( (new_node->
ind1 < tree->
ind1) ||
6161 tree->
left = add_to_indexed_vdw_pairs(new_node, tree->
left);
6165 tree->
right = add_to_indexed_vdw_pairs(new_node, tree->
right);
6192 if ( (new_node->
ind1 < tree->
ind1) ||
6195 tree->
left = add_to_indexed_nbthole_pairs(new_node, tree->
left);
6199 tree->
right = add_to_indexed_nbthole_pairs(new_node, tree->
right);
6227 int Parameters::vdw_pair_to_arrays(
int *ind1_array,
int *ind2_array,
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;
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);
6272 int Parameters::nbthole_pair_to_arrays(
int *ind1_array,
int *ind2_array,
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;
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);
6313 int Parameters::table_pair_to_arrays(
int *ind1_array,
int *ind2_array,
6321 ind1_array[arr_index] = tree->
ind1;
6322 ind2_array[arr_index] = tree->
ind2;
6323 K[arr_index] = tree->
K;
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);
6349 read_parm(amber_data,vdw14);
6368 int i,j,ico,numtype,mul;
6372 NAMD_die(
"No data read from parm file yet!");
6379 NAMD_die(
"memory allocation of bond_array failed!");
6391 NAMD_die(
"memory allocation of angle_array failed!");
6413 NAMD_die(
"memory allocation of dihedral_array failed!");
6420 {
char err_msg[128];
6421 sprintf(err_msg,
"The periodicity of dihedral # %d is zero!", i+1);
6427 if (amber_data->
Pn[i] > 0)
6432 {
char err_msg[181];
6433 sprintf(err_msg,
"Multiple dihedral with multiplicity of %d greater than max of %d",
6439 NAMD_die(
"Negative periodicity in the last dihedral entry!");
6446 for (i=0; i<numtype; ++i)
6447 for (j=i; j<numtype; ++j)
6449 if (new_node == NULL)
6450 NAMD_die(
"memory allocation of vdw_pair_tree failed!");
6453 new_node->
left = new_node->
right = NULL;
6459 ico = amber_data->
Cno[numtype*i+j];
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;
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";
6471 NAMD_die(
"Encounter non-zero 10-12 H-bond term!");
6484 read_parm(amber_data,vdw14);
6489 int i,j,ico,numtype,mul;
6493 NAMD_die(
"No data read from parm file yet!");
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.");
6507 NAMD_die(
"memory allocation of bond_array failed!");
6519 NAMD_die(
"memory allocation of angle_array failed!");
6541 NAMD_die(
"memory allocation of dihedral_array failed!");
6548 {
char err_msg[128];
6549 sprintf(err_msg,
"The periodicity of dihedral # %d is zero!", i+1);
6555 if (amber_data->
Pn[i] > 0)
6560 {
char err_msg[181];
6561 sprintf(err_msg,
"Multiple dihedral with multiplicity of %d greater than max of %d",
6567 NAMD_die(
"Negative periodicity in the last dihedral entry!");
6580 sprintf(err_msg,
"The table of %d crossterm type has a resolution(%d) other than %d!",
6586 for (
int j = 0; j < N; ++j) {
6587 for (
int k = 0; k < N; ++k) {
6592 for (
int j = 0; j < N; ++j) {
6607 for (i=0; i<numtype; ++i)
6608 for (j=i; j<numtype; ++j)
6610 if (new_node == NULL)
6611 NAMD_die(
"memory allocation of vdw_pair_tree failed!");
6614 new_node->
left = new_node->
right = NULL;
6620 ico = amber_data->
Cno[numtype*i+j];
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;
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";
6632 NAMD_die(
"Encounter non-zero 10-12 H-bond term!");
6682 NAMD_die(
"memory allocation of bond_array failed!");
6687 NAMD_die(
"memory allocation of dihedral_array failed!");
6692 NAMD_die(
"memory allocation of angle_array failed!");
6750 NAMD_die(
"I can't do RB dihedrals with MAX_MULTIPLICITY < 5");
6756 if(j%2 == 1) c[j] = -c[j];
6781 test2 = c[0]+1/2.*c[2]+3/8.*c[4];
6786 if(fabs(test1-test2) > 0.0001)
6787 NAMD_die(
"Internal error: failed to handle RB dihedrals");
6804 NAMD_die(
"unknown dihedral type found");
6814 for (i=0; i<numtype; i++) {
6815 for (j=i; j<numtype; j++) {
6819 if (new_node == NULL)
6820 NAMD_die(
"memory allocation of vdw_pair_tree failed!");
6823 new_node->
left = new_node->
right = NULL;
6826 &(new_node->
B14), &(new_node->
A14));
6832 if(min && ( fabs(new_node->
A) < 1.0 )) {
6836 "Adding small LJ repulsion term to some atoms.\n" <<
endi;
6848 Real *pairC6,*pairC12;
6865 const_cast<GromacsTopFile*
>(gf)->getPairLJArrays2(atom1, atom2, pairC6, pairC12);
6892 enertable = fopen(table_file,
"r");
6894 if (enertable == NULL) {
6895 NAMD_die(
"ERROR: Could not open energy table file!\n");
6898 char tableline[121];
6912 iout <<
"Beginning table read\n" <<
endi;
6913 while(fgets(tableline,120,enertable)) {
6914 if (strncmp(tableline,
"#",1)==0) {
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");
6924 if (maxdist < simParams->cutoff) {
6925 NAMD_die(
"Tabulated energies must at least extend to the cutoff distance\n");
6928 if (maxdist > simParams->
cutoff) {
6929 iout <<
"Info: Reducing max table size to match nonbond cutoff\n";
6930 maxdist = ceil(simParams->
cutoff);
6945 int numtypesread = 0;
6947 while(fgets(tableline,120,enertable)) {
6948 if (strncmp(tableline,
"#",1)==0) {
6951 if (strncmp(tableline,
"TYPE",4)==0) {
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");
6961 if (numtypesread == numtypes) {
6962 NAMD_die(
"Error: Number of types in table doesn't match table header\n");
6966 if (!strncasecmp(interp_type,
"linear", 6)) {
6969 sprintf(err_msg,
"Failed to read table energy (linear) type %s\n", newtypename);
6972 }
else if (!strncasecmp(interp_type,
"cubic", 5)) {
6975 sprintf(err_msg,
"Failed to read table energy (cubic) type %s\n", newtypename);
6979 NAMD_die(
"Table interpolation type must be linear or cubic\n");
6991 if (numtypesread != numtypes) {
6993 sprintf(err_msg,
"ERROR: Expected %i tabulated energy types but got %i\n", numtypes, numtypesread);
7133 char tableline[120];
7150 std::vector<BigReal> dists;
7151 std::vector<BigReal> enervalues;
7152 std::vector<BigReal> forcevalues;
7164 while(fgets(tableline,120,enertable) && distbin <= (
int) (
namdnearbyint(maxdist / table_spacing) + 1)) {
7165 if (strncmp(tableline,
"#",1)==0) {
7168 if (strncmp(tableline,
"TYPE",4)==0) {
7169 fseek(enertable, -1 * (
long) strlen(tableline), SEEK_CUR);
7174 int readcount = sscanf(tableline,
"%lf %lf %lf", &readdist, &readener, &readforce);
7177 if (readcount != 3) {
7179 sprintf(err_msg,
"ERROR: Failed to parse table line %s!\n", tableline);
7184 if (readdist < lastdist) {
7185 NAMD_die(
"ERROR: Encountered badly ordered entries in energy table!\n");
7188 lastdist = readdist;
7189 dists.push_back(readdist);
7190 enervalues.push_back(readener);
7191 forcevalues.push_back(readforce);
7196 if (dists[0] != 0.0) {
7197 NAMD_die(
"ERROR: First data point for energy table must be at r=0\n");
7199 spacing = dists[1] - dists[0];
7200 for (i=2; i<(numentries - 1); i++) {
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");
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];
7219 be[0] = 3 * (enervalues[1] - enervalues[0]);
7220 for (i=1; i < (numentries - 1); i++) {
7222 be[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
7225 be[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
7227 bf[0] = 3 * (forcevalues[1] - forcevalues[0]);
7228 for (i=1; i < (numentries - 1); i++) {
7230 bf[i] = 3 * (forcevalues[i+1] - forcevalues[i-1]);
7233 bf[numentries - 1] = 3 * (forcevalues[numentries - 1] - forcevalues[numentries - 2]);
7235 memset(m,0,numentries*numentries*
sizeof(
double));
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;
7244 m[
INDEX(numentries,numentries-1,numentries-1)] = 2;
7250 for (i=0; i<numentries; i++) {
7253 for (j=i+1; j<numentries; j++) {
7255 const BigReal factor = myval / baseval;
7257 for (
int k=i; k<numentries; k++) {
7259 m[
INDEX(numentries,j,k)] -= (factor * subval);
7262 be[j] -= (factor * be[i]);
7263 bf[j] -= (factor * bf[i]);
7269 for (i=numentries-1; i>=0; i--) {
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] );
7277 xe[i] = be[i] / m[
INDEX(numentries,i,i)];
7278 xf[i] = bf[i] / m[
INDEX(numentries,i,i)];
7286 int entriesperseg = (int) ceil(spacing / table_spacing);
7287 int table_prefix = 2 * typeindex * (int) (
namdnearbyint(maxdist / table_spacing) + 1);
7289 for (i=0; i<numentries-1; i++) {
7292 currdist = dists[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];
7302 Af = forcevalues[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];
7308 for (j=0; j<entriesperseg; j++) {
7309 const BigReal mydist = currdist + (j * table_spacing);
7310 const BigReal mydistfrac = (float) j / (entriesperseg - 1);
7312 if (distbin >= (
int)
namdnearbyint(maxdist / table_spacing))
break;
7316 energy = Ae + (Be * mydistfrac) + (Ce * mydistfrac * mydistfrac) + (De * mydistfrac * mydistfrac * mydistfrac);
7317 force = Af + (Bf * mydistfrac) + (Cf * mydistfrac * mydistfrac) + (Df * mydistfrac * mydistfrac * mydistfrac);
7320 table_ener[table_prefix + 2 * distbin] = energy;
7321 table_ener[table_prefix + 2 * distbin + 1] = force;
7324 if (currdist >= maxdist)
break;
7330 table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
7331 table_ener[table_prefix + 2 * distbin + 1] = 0.0;
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;
7369 char tableline[120];
7385 std::vector<BigReal> dists;
7386 std::vector<BigReal> enervalues;
7398 while(fgets(tableline,120,enertable) && distbin <= (
int) (
namdnearbyint(maxdist / table_spacing) + 1)) {
7399 if (strncmp(tableline,
"#",1)==0) {
7402 if (strncmp(tableline,
"TYPE",4)==0) {
7403 fseek(enertable, -1 * (
long) strlen(tableline), SEEK_CUR);
7408 int readcount = sscanf(tableline,
"%lf %lf", &readdist, &readener);
7411 if (readcount != 2) {
7413 sprintf(err_msg,
"ERROR: Failed to parse table line %s!\n", tableline);
7418 if (readdist < lastdist) {
7419 NAMD_die(
"ERROR: Encountered badly ordered entries in energy table!\n");
7422 lastdist = readdist;
7423 dists.push_back(readdist);
7424 enervalues.push_back(readener);
7429 if (dists[0] != 0.0) {
7430 NAMD_die(
"ERROR: First data point for energy table must be at r=0\n");
7432 spacing = dists[1] - dists[0];
7433 for (i=2; i<(numentries - 1); i++) {
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");
7445 double* m =
new double[numentries*numentries];
7446 double*
x =
new double[numentries];
7447 double* b =
new double[numentries];
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]);
7455 b[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
7457 memset(m,0,numentries*numentries*
sizeof(
double));
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;
7466 m[
INDEX(numentries,numentries-1,numentries-1)] = 2;
7470 printf(
"Solving the matrix equation: \n");
7472 for (i=0; i<numentries; i++) {
7474 for (j=0; j<numentries; j++) {
7475 printf(
" %6.3f,", m[
INDEX(numentries, i, j)]);
7477 printf(
" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
7481 for (i=0; i<numentries; i++) {
7484 for (j=i+1; j<numentries; j++) {
7486 const BigReal factor = myval / baseval;
7488 for (
int k=i; k<numentries; k++) {
7490 m[
INDEX(numentries,j,k)] -= (factor * subval);
7493 b[j] -= (factor * b[i]);
7498 printf(
" In upper diagonal form, equation is:\n");
7499 for (i=0; i<numentries; i++) {
7501 for (j=0; j<numentries; j++) {
7502 printf(
" %6.3f,", m[
INDEX(numentries, i, j)]);
7504 printf(
" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
7508 for (i=numentries-1; i>=0; i--) {
7511 for (j=i+1; j<numentries; j++) {
7512 b[i] -= ( m[
INDEX(numentries,i,j)] * x[j] );
7515 x[i] = b[i] / m[
INDEX(numentries,i,i)];
7519 printf(
" Solution vector is:\n\t(");
7520 for (i=0; i<numentries; i++) {
7521 printf(
" %6.3f ", x[i]);
7529 int entriesperseg = (int) ceil(spacing / table_spacing);
7530 int table_prefix = 2 * typeindex * (int) (
namdnearbyint(maxdist / table_spacing) + 1);
7532 for (i=0; i<numentries-1; i++) {
7534 currdist = dists[i];
7536 printf(
"Interpolating on interval %i\n", 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];
7544 printf(
"Coefficients for this interval: %f %f %f %f\n", A, B, C, D);
7547 for (j=0; j<entriesperseg; j++) {
7548 const BigReal mydist = currdist + (j * table_spacing);
7549 const BigReal mydistfrac = (float) j / (entriesperseg - 1);
7551 if (distbin >= (
int)
namdnearbyint(maxdist / table_spacing))
break;
7555 energy = A + (B * mydistfrac) + (C * mydistfrac * mydistfrac) + (D * mydistfrac * mydistfrac * mydistfrac);
7556 force = B + (2 * C * mydistfrac) + (3 * D * mydistfrac * mydistfrac);
7558 force *= (1.0 / spacing);
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;
7565 if (currdist >= maxdist)
break;
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;
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;
7608 char tableline[120];
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) {
7632 if (strncmp(tableline,
"TYPE",4)==0) {
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");
7650 if (readcount != 3) {
7652 sprintf(err_msg,
"ERROR: Failed to parse table line %s!\n", tableline);
7657 if (readdist < lastdist) {
7658 NAMD_die(
"ERROR: Encountered badly ordered entries in energy table!\n");
7661 currdist = lastdist;
7663 while (currdist <= readdist && distbin <= (
int) (
namdnearbyint(maxdist / 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;
7676 fseek(enertable, -1 * (
long) strlen(tableline), SEEK_CUR);
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;
7705 m = (val2 - val1) / (end2 - end1);
7707 val = ((currdist-end1) * m + val1);
struct improper_params * next
GromacsPairValue * gromacsPair_array
void assign_improper_index(const char *, const char *, const char *, const char *, Improper *, int)
std::ostream & iINFO(std::ostream &s)
CrosstermData c[dim][dim]
int NAMD_read_line(FILE *fd, char *buf, int bufsize)
void assign_vdw_index(const char *, Atom *)
void read_charmm_parameter_file(char *)
static char ** table_types
struct indexed_vdw_pair * right
void assign_bond_index(const char *, const char *, Bond *)
void getVDWParams(int typea, int typeb, Real *c6, Real *c12, Real *c6pair, Real *c7) const
struct indexed_table_pair * right
std::ostream & endi(std::ostream &s)
struct nbthole_pair_params * next
std::ostream & iWARN(std::ostream &s)
MIStream * get(char &data)
int read_energy_type_cubspline(FILE *, const int, BigReal *, const float, const float)
DihedralValue * dihedral_array
struct crossterm_params * next
void NAMD_find_first_word(char *source, char *word)
crossterm_params(int dim)
CrosstermValue * crossterm_array
IndexedVdwPair * vdw_pair_tree
int read_energy_type(FILE *, const int, BigReal *, const float, const float)
char tabulatedEnergiesFile[NAMD_FILENAME_BUFFER_SIZE]
#define INDEX(ncols, i, j)
void print_param_summary()
IndexedNbtholePair * nbthole_pair_tree
void read_parameter_file(char *)
FourBodyConsts values[MAX_MULTIPLICITY]
vector< int > CMAPResolution
void getBondParams(int num, Real *b0, Real *kB, int *funct) const
static Units next(Units u)
struct indexed_vdw_pair * left
int get_int_table_type(char *)
struct bond_params * right
void assign_crossterm_index(const char *, const char *, const char *, const char *, const char *, const char *, const char *, const char *, Crossterm *)
int NAMD_blank_string(char *str)
void print_vdw_pair_params()
struct indexed_table_pair * left
void print_dihedral_params()
FILE * Fopen(const char *filename, const char *mode)
struct indexed_vdw_pair IndexedVdwPair
ImproperValue * improper_array
int get_vdw_pair_params(Index ind1, Index ind2, Real *, Real *, Real *, Real *)
void NAMD_die(const char *err_msg)
int getNumDihedralParams() const
NbtholePairValue * nbthole_array
char * atom_type_name(Index a)
FourBodyConsts values[MAX_MULTIPLICITY]
struct indexed_nbthole_pair * right
struct bond_params * left
struct angle_params * left
struct dihedral_params * next
int getNumAtomParams() const
int get_table_pair_params(Index, Index, int *)
void send_Parameters(MOStream *)
void done_reading_files(Bool)
FourBodyConsts values[MAX_MULTIPLICITY]
void print_nbthole_pair_params()
void getDihedralParams(int num, Real *c, int *mult, int *funct) const
vector< vector< _REAL > > CMAPParameter
void crossterm_setup(CrosstermData *table)
void read_ener_table(SimParameters *)
void print_improper_params()
int read_energy_type_bothcubspline(FILE *, const int, BigReal *, const float, const float)
void assign_dihedral_index(const char *, const char *, const char *, const char *, Dihedral *, int, int)
int getNumAngleParams() const
struct table_pair_params * next
void done_reading_structure()
struct vdw_pair_params * next
void assign_angle_index(const char *, const char *, const char *, Angle *, int)
MOStream * put(char data)
int getNumBondParams() const
void NAMD_remove_comment(char *str)
struct indexed_nbthole_pair * left
void print_angle_params()
valT values[SORTTILELISTSKERNEL_ITEMS_PER_THREAD]
IndexedTablePair * tab_pair_tree
struct angle_params * right
FourBodyConsts values[MAX_MULTIPLICITY]
char tableInterpType[128]
struct vdw_params * right
#define MAX_ATOMTYPE_CHARS
void receive_Parameters(MIStream *)
void getAngleParams(int num, Real *th0, Real *kth, int *funct) const