31 #define MIN_DEBUG_LEVEL 3 35 #define INDEX(ncols,i,j) ((i)*ncols + (j)) 51 #ifndef ACCELERATED_INPUT 72 #ifndef ACCELERATED_INPUT_ANGLES 86 #ifndef ACCELERATED_INPUT_DIHEDRALS 197 void Parameters::initialize() {
203 #ifdef ACCELERATED_INPUT 232 maxDihedralMults=NULL;
233 maxImproperMults=NULL;
286 CkPrintf(
"Working on tables\n");
289 #ifdef ACCELERATED_INPUT 305 AllFilesRead =
FALSE;
322 }
while ( f != NULL );
343 delete [] atomTypeNames;
346 free_bond_tree(bondp);
349 free_angle_tree(anglep);
351 if (dihedralp != NULL)
352 free_dihedral_list(dihedralp);
354 if (improperp != NULL)
355 free_improper_list(improperp);
357 if (crosstermp != NULL)
358 free_crossterm_list(crosstermp);
363 if (vdw_pairp != NULL)
364 free_vdw_pair_list();
366 if (nbthole_pairp != NULL)
367 free_nbthole_pair_list();
401 if (maxDihedralMults != NULL)
402 delete [] maxDihedralMults;
404 if (maxImproperMults != NULL)
405 delete [] maxImproperMults;
407 for(
int i = 0; i < error_msgs.
size(); ++i ) {
408 delete [] error_msgs[i];
434 char first_word[512];
441 NAMD_die(
"Tried to read another parameter file after being told that all files were read!");
445 if ( (pfile =
Fopen(fname,
"r")) == NULL)
449 snprintf(err_msg,
sizeof(err_msg),
450 "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
466 if ((buffer[0] !=
'!') &&
468 (strncasecmp(first_word,
"REMARK", 6) != 0) &&
469 (strcasecmp(first_word,
"set")!=0) &&
470 (strncasecmp(first_word,
"AEXP", 4) != 0) &&
471 (strncasecmp(first_word,
"REXP", 4) != 0) &&
472 (strncasecmp(first_word,
"HAEX", 4) != 0) &&
473 (strncasecmp(first_word,
"AAEX", 4) != 0) &&
474 (strncasecmp(first_word,
"NBOND", 5) != 0) &&
475 (strncasecmp(first_word,
"CUTNB", 5) != 0) &&
476 (strncasecmp(first_word,
"END", 3) != 0) &&
477 (strncasecmp(first_word,
"CTONN", 5) != 0) &&
478 (strncasecmp(first_word,
"EPS", 3) != 0) &&
479 (strncasecmp(first_word,
"VSWI", 4) != 0) &&
480 (strncasecmp(first_word,
"NBXM", 4) != 0) &&
481 (strncasecmp(first_word,
"INHI", 4) != 0) )
485 if (strncasecmp(first_word,
"bond", 4)==0)
487 add_bond_param(buffer,
TRUE);
490 else if (strncasecmp(first_word,
"angl", 4)==0)
492 add_angle_param(buffer);
495 else if (strncasecmp(first_word,
"dihe", 4)==0)
497 add_dihedral_param(buffer, pfile);
500 else if (strncasecmp(first_word,
"impr", 4)==0)
502 add_improper_param(buffer, pfile);
505 else if (strncasecmp(first_word,
"nonb", 4)==0)
507 add_vdw_param(buffer);
510 else if (strncasecmp(first_word,
"nbfi", 4)==0)
512 add_vdw_pair_param(buffer);
515 else if (strncasecmp(first_word,
"nbta", 4)==0)
522 add_table_pair_param(buffer);
525 else if (strncasecmp(first_word,
"hbon", 4)==0)
527 add_hb_pair_param(buffer);
535 snprintf(err_msg,
sizeof(err_msg),
536 "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
575 char first_word[512];
582 NAMD_die(
"Tried to read another parameter file after being told that all files were read!");
586 if ( (pfile = fopen(fname,
"r")) == NULL)
590 snprintf(err_msg,
sizeof(err_msg),
591 "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
606 (strncmp(first_word,
"!", 1) != 0) &&
607 (strncmp(first_word,
"*", 1) != 0) &&
608 (strncasecmp(first_word,
"END", 3) != 0))
611 iout <<
iWARN <<
"SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" <<
endi;
615 if (strncasecmp(first_word,
"bond", 4)==0)
617 par_type=1; skipline=1;
619 else if (strncasecmp(first_word,
"angl", 4)==0)
621 par_type=2; skipline=1;
623 else if (strncasecmp(first_word,
"thet", 4)==0)
625 par_type=2; skipline=1;
627 else if (strncasecmp(first_word,
"dihe", 4)==0)
629 par_type=3; skipline=1;
631 else if (strncasecmp(first_word,
"phi", 3)==0)
633 par_type=3; skipline=1;
635 else if (strncasecmp(first_word,
"impr", 4)==0)
637 par_type=4; skipline=1;
639 else if (strncasecmp(first_word,
"imph", 4)==0)
641 par_type=4; skipline=1;
643 else if (strncasecmp(first_word,
"nbon", 4)==0)
645 par_type=5; skipline=1;
647 else if (strncasecmp(first_word,
"nonb", 4)==0)
649 par_type=5; skipline=1;
651 else if (strncasecmp(first_word,
"nbfi", 4)==0)
653 par_type=6; skipline=1;
655 else if (strncasecmp(first_word,
"hbon", 4)==0)
657 par_type=7; skipline=1;
659 else if (strncasecmp(first_word,
"cmap", 4)==0)
661 par_type=8; skipline=1;
663 else if (strncasecmp(first_word,
"nbta", 4)==0)
665 par_type=9; skipline=1;
667 else if (strncasecmp(first_word,
"thol", 4)==0)
669 par_type=10; skipline=1;
671 else if (strncasecmp(first_word,
"atom", 4)==0)
673 par_type=11; skipline=1;
675 else if (strncasecmp(first_word,
"ioformat", 8)==0)
679 else if (strncasecmp(first_word,
"read", 4)==0)
681 skip_stream_read(buffer, pfile); skipline=1;
683 else if (strncasecmp(first_word,
"return", 4)==0)
685 skipall=1; skipline=1;
687 else if ((strncasecmp(first_word,
"nbxm", 4) == 0) ||
688 (strncasecmp(first_word,
"grou", 4) == 0) ||
689 (strncasecmp(first_word,
"cdie", 4) == 0) ||
690 (strncasecmp(first_word,
"shif", 4) == 0) ||
691 (strncasecmp(first_word,
"vgro", 4) == 0) ||
692 (strncasecmp(first_word,
"vdis", 4) == 0) ||
693 (strncasecmp(first_word,
"vswi", 4) == 0) ||
694 (strncasecmp(first_word,
"cutn", 4) == 0) ||
695 (strncasecmp(first_word,
"ctof", 4) == 0) ||
696 (strncasecmp(first_word,
"cton", 4) == 0) ||
697 (strncasecmp(first_word,
"eps", 3) == 0) ||
698 (strncasecmp(first_word,
"e14f", 4) == 0) ||
699 (strncasecmp(first_word,
"wmin", 4) == 0) ||
700 (strncasecmp(first_word,
"aexp", 4) == 0) ||
701 (strncasecmp(first_word,
"rexp", 4) == 0) ||
702 (strncasecmp(first_word,
"haex", 4) == 0) ||
703 (strncasecmp(first_word,
"aaex", 4) == 0) ||
704 (strncasecmp(first_word,
"noac", 4) == 0) ||
705 (strncasecmp(first_word,
"hbno", 4) == 0) ||
706 (strncasecmp(first_word,
"cuth", 4) == 0) ||
707 (strncasecmp(first_word,
"ctof", 4) == 0) ||
708 (strncasecmp(first_word,
"cton", 4) == 0) ||
709 (strncasecmp(first_word,
"cuth", 4) == 0) ||
710 (strncasecmp(first_word,
"ctof", 4) == 0) ||
711 (strncasecmp(first_word,
"cton", 4) == 0) )
713 if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
717 snprintf(err_msg,
sizeof(err_msg),
718 "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*", fname, buffer);
726 else if (par_type == 0)
732 snprintf(err_msg,
sizeof(err_msg),
733 "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",
743 if ( (par_type != 0) && (!skipline) )
750 add_bond_param(buffer,
TRUE);
753 else if (par_type == 2)
755 add_angle_param(buffer);
758 else if (par_type == 3)
760 add_dihedral_param(buffer, pfile);
763 else if (par_type == 4)
765 add_improper_param(buffer, pfile);
768 else if (par_type == 5)
770 add_vdw_param(buffer);
773 else if (par_type == 6)
775 add_vdw_pair_param(buffer);
778 else if (par_type == 7)
780 add_hb_pair_param(buffer);
782 else if (par_type == 8)
784 add_crossterm_param(buffer, pfile);
787 else if (par_type == 9)
794 add_table_pair_param(buffer);
797 else if (par_type == 10)
799 add_nbthole_pair_param(buffer);
802 else if (par_type == 11)
804 if ( strncasecmp(first_word,
"mass", 4) != 0 ) {
806 snprintf(err_msg,
sizeof(err_msg),
807 "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",
819 snprintf(err_msg,
sizeof(err_msg),
820 "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",
835 void Parameters::skip_stream_read(
char *buf, FILE *fd) {
838 char first_word[513];
841 int rval = sscanf(buf,
"%s %s", s1, s2);
844 snprintf(err_msg,
sizeof(err_msg),
845 "BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
848 if ( ! strncasecmp(s2,
"PARA",4) )
return;
850 iout <<
iINFO <<
"SKIPPING " << s2 <<
" SECTION IN STREAM FILE\n" <<
endi;
857 (strncmp(first_word,
"!", 1) != 0) &&
858 (strncmp(first_word,
"*", 1) != 0) &&
859 (strncasecmp(first_word,
"END", 3) == 0) )
return;
878 void Parameters::add_bond_param(
const char *buf,
Bool overwrite)
886 #ifndef ACCELERATED_INPUT 915 snprintf(err_msg,
sizeof(err_msg),
916 "BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
920 snprintf(err_msg,
sizeof(err_msg),
921 "BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
925 #ifndef ACCELERATED_INPUT 929 if (new_node == NULL)
931 NAMD_die(
"memory allocation failed in Parameters::add_bond_param\n");
940 #ifndef ACCELERATED_INPUT 949 #ifndef ACCELERATED_INPUT 956 #ifndef ACCELERATED_INPUT 962 new_node->left = NULL;
963 new_node->right = NULL;
966 bondp=add_to_bond_tree(new_node, bondp, overwrite);
974 if ( !retval.first && ((retval.second->k != new_node.
k) ||
975 (retval.second->x0 != new_node.
x0)) && overwrite)
977 iout <<
"\n" <<
iWARN <<
"DUPLICATE BOND ENTRY FOR " 980 <<
"\nPREVIOUS VALUES k=" << retval.second->k
981 <<
" x0=" << retval.second->x0
982 <<
"\n USING VALUES k=" << new_node.
k 983 <<
" x0=" << new_node.
x0 986 retval.second->k = new_node.
k;
987 retval.second->x0 = new_node.
x0;
988 retval.second->x1 = new_node.
x1;
1018 #ifndef ACCELERATED_INPUT 1031 if (compare_code == 0)
1037 if (compare_code == 0)
1048 iout <<
"\n" <<
iWARN <<
"DUPLICATE BOND ENTRY FOR " 1071 if (compare_code < 0)
1073 tree->left = add_to_bond_tree(new_node, tree->left, overwrite);
1077 tree->right = add_to_bond_tree(new_node, tree->right, overwrite);
1100 void Parameters::add_angle_param(
char *buf)
1112 #ifndef ACCELERATED_INPUT_ANGLES 1123 read_count=sscanf(buf,
"%*s %s %s %s %f %f UB %f %f\n",
1127 else if ((paramType ==
paraCharmm) && cosAngles) {
1128 read_count=sscanf(buf,
"%s %s %s %f %f %3s %f %f\n",
1137 read_count=sscanf(buf,
"%s %s %s %f %f %f %f\n",
1146 if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
1152 snprintf(err_msg,
sizeof(err_msg),
1153 "BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1157 snprintf(err_msg,
sizeof(err_msg),
1158 "BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1162 #ifndef ACCELERATED_INPUT_ANGLES 1166 if (new_node == NULL)
1168 NAMD_die(
"memory allocation failed in Parameters::add_angle_param");
1176 #ifndef ACCELERATED_INPUT_ANGLES 1186 #ifndef ACCELERATED_INPUT_ANGLES 1194 #ifndef ACCELERATED_INPUT_ANGLES 1197 new_node->angle =
angle;
1204 if (strcasecmp(
"cos",norm)==0) {
1207 #ifndef ACCELERATED_INPUT_ANGLES 1214 #ifndef ACCELERATED_INPUT_ANGLES 1221 #ifndef ACCELERATED_INPUT_ANGLES 1228 if (read_count == 7)
1231 #ifndef ACCELERATED_INPUT_ANGLES 1232 if (new_node->
normal == 0) {
1234 if (new_node.
normal == 0) {
1236 NAMD_die(
"ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
1238 #ifndef ACCELERATED_INPUT_ANGLES 1248 #ifndef ACCELERATED_INPUT_ANGLES 1249 new_node->
k_ub = 0.0;
1250 new_node->
r_ub = 0.0;
1252 new_node.
k_ub = 0.0;
1253 new_node.
r_ub = 0.0;
1256 #ifndef ACCELERATED_INPUT_ANGLES 1257 new_node->left = NULL;
1258 new_node->right = NULL;
1261 anglep = add_to_angle_tree(new_node, anglep);
1266 if(retval.first==
false)
1274 if ((retval.second->k != new_node.
k) ||
1275 (retval.second->theta0 != new_node.
theta0) ||
1276 (retval.second->k_ub != new_node.
k_ub) ||
1277 (retval.second->r_ub != new_node.
r_ub) ||
1278 (retval.second->normal != new_node.
normal))
1280 iout <<
"\n" <<
iWARN <<
"DUPLICATE ANGLE ENTRY FOR " 1284 <<
"\nPREVIOUS VALUES k=" 1285 << retval.second->k <<
" theta0=" 1286 << retval.second->theta0 <<
" k_ub=" 1287 << retval.second->k_ub <<
" r_ub=" 1288 << retval.second->r_ub
1289 <<
"\n USING VALUES k=" 1290 << new_node.
k <<
" theta0=" 1291 << new_node.
theta0 <<
" k_ub=" 1292 << new_node.
k_ub <<
" r_ub=" << new_node.
r_ub 1295 retval.second->k = new_node.
k;
1296 retval.second->theta0 = new_node.
theta0;
1297 retval.second->k_ub = new_node.
k_ub;
1298 retval.second->r_ub = new_node.
r_ub;
1299 retval.second->normal = new_node.
normal;
1332 #ifndef ACCELERATED_INPUT_ANGLES 1342 if (compare_code == 0)
1347 if (compare_code == 0)
1350 compare_code = strcasecmp(new_node->
atom3name,
1353 if (compare_code == 0)
1366 iout <<
"\n" <<
iWARN <<
"DUPLICATE ANGLE ENTRY FOR " 1370 <<
"\nPREVIOUS VALUES k=" 1372 << tree->
angle <<
" k_ub=" 1373 << tree->
k_ub <<
" r_ub=" 1375 <<
"\n USING VALUES k=" 1377 << new_node->
angle <<
" k_ub=" 1378 << new_node->
k_ub <<
" r_ub=" << new_node->
r_ub 1399 if (compare_code < 0)
1401 tree->left = add_to_angle_tree(new_node, tree->left);
1405 tree->right = add_to_angle_tree(new_node, tree->right);
1425 void Parameters::add_dihedral_param(
char *buf, FILE *fd)
1447 read_count=sscanf(buf,
"%*s %s %s %s %s MULTIPLE= %d %f %d %f\n",
1448 atom1name, atom2name, atom3name, atom4name, &multiplicity,
1449 &forceconstant, &periodicity, &phase_shift);
1454 read_count=sscanf(buf,
"%s %s %s %s %f %d %f\n",
1455 atom1name, atom2name, atom3name, atom4name,
1456 &forceconstant, &periodicity, &phase_shift);
1460 if ( (read_count != 4) && (read_count != 8) && (paramType ==
paraXplor) )
1464 snprintf(err_msg,
sizeof(err_msg),
1465 "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1468 else if ( (read_count != 7) && (paramType ==
paraCharmm) )
1472 snprintf(err_msg,
sizeof(err_msg),
1473 "BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
1477 if ( (read_count == 4) && (paramType ==
paraXplor) )
1480 read_count=sscanf(buf,
"%*s %*s %*s %*s %*s %f %d %f\n",
1481 &forceconstant, &periodicity, &phase_shift);
1484 if (read_count != 3)
1488 snprintf(err_msg,
sizeof(err_msg),
1489 "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
1500 snprintf(err_msg,
sizeof(err_msg),
1501 "Multiple dihedral with multiplicity of %d greater than max of %d",
1509 if (new_node == NULL)
1511 NAMD_die(
"memory allocation failed in Parameters::add_dihedral_param\n");
1519 #ifndef ACCELERATED_INPUT_DIHEDRALS 1520 strcpy(new_node->atom1name, atom1name);
1521 strcpy(new_node->atom2name, atom2name);
1522 strcpy(new_node->atom3name, atom3name);
1523 strcpy(new_node->atom4name, atom4name);
1524 new_node->atom1wild = ! strcasecmp(atom1name,
"X");
1525 new_node->atom2wild = ! strcasecmp(atom2name,
"X");
1526 new_node->atom3wild = ! strcasecmp(atom3name,
"X");
1527 new_node->atom4wild = ! strcasecmp(atom4name,
"X");
1528 new_node->multiplicity = multiplicity;
1529 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
1530 new_node->values[0].k = forceconstant;
1531 new_node->values[0].n = periodicity;
1532 new_node->values[0].delta = phase_shift;
1534 new_node->next = NULL;
1537 int cmp =strcasecmp(atom1name, atom4name);
1541 if(strcasecmp(atom2name, atom3name)>0)
1552 new_node->
key= (reverse) ?
1553 TupleString4(atom4name, atom3name, atom2name, atom1name) :
1554 TupleString4(atom1name, atom2name, atom3name, atom4name);
1556 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
1563 if (multiplicity > 1)
1565 for (i=1; i<multiplicity; i++)
1583 NAMD_die(
"EOF encoutner in middle of multiple dihedral");
1586 read_count=sscanf(buffer,
"%f %d %f\n",
1587 &forceconstant, &periodicity, &phase_shift);
1589 if (read_count != 3)
1593 snprintf(err_msg,
sizeof(err_msg),
1594 "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
1598 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
1599 #ifndef ACCELERATED_INPUT_DIHEDRALS 1600 new_node->values[i].k = forceconstant;
1601 new_node->values[i].n = periodicity;
1602 new_node->values[i].delta = phase_shift;
1615 add_to_dihedral_list(new_node);
1619 add_to_charmm_dihedral_list(new_node);
1645 void Parameters::add_to_dihedral_list(
1655 #ifndef ACCELERATED_INPUT_DIHEDRALS 1657 if (dihedralp == NULL)
1671 if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1672 (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1673 (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1674 (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1675 ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1676 (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1677 (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1678 (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
1680 if(dihedralmap==NULL)
1684 if(retval.first==
false)
1685 #endif // if duplicate 1691 #ifndef ACCELERATED_INPUT_DIHEDRALS 1692 if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
1698 #ifndef ACCELERATED_INPUT_DIHEDRALS 1699 for (i=0; i<ptr->multiplicity; i++)
1701 if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1;
break;}
1702 if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1;
break;}
1703 if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1;
break;}
1706 for (i=0; i=retval.second->multiplicity; i++)
1708 if (retval.second->values[i].k != new_node->
value.
values[i].
k) {echoWarn=1;
break;}
1709 if (retval.second->values[i].n != new_node->
value.
values[i].
n) {echoWarn=1;
break;}
1710 if (retval.second->values[i].delta != new_node->
value.
values[i].
delta) {echoWarn=1;
break;}
1716 iout <<
"\n" <<
iWARN <<
"DUPLICATE DIHEDRAL ENTRY FOR " 1717 #ifndef ACCELERATED_INPUT_DIHEDRALS 1718 << ptr->atom1name <<
"-" 1719 << ptr->atom2name <<
"-" 1720 << ptr->atom3name <<
"-" 1722 <<
"\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity <<
"\n";
1731 #ifndef ACCELERATED_INPUT_DIHEDRALS 1732 for (i=0; i<ptr->multiplicity; i++)
1734 iout <<
" k=" << ptr->values[i].k
1735 <<
" n=" << ptr->values[i].n
1736 <<
" delta=" << ptr->values[i].delta;
1738 iout <<
"\nUSING VALUES MULTIPLICITY " << new_node->multiplicity <<
"\n";
1740 for (i=0; i< retval.second->multiplicity; i++)
1742 iout <<
" k=" << retval.second->values[i].k
1743 <<
" n=" << retval.second->values[i].n
1744 <<
" delta=" << retval.second->values[i].delta;
1748 #ifndef ACCELERATED_INPUT_DIHEDRALS 1749 for (i=0; i<new_node->multiplicity; i++)
1751 iout <<
" k=" << new_node->values[i].k
1752 <<
" n=" << new_node->values[i].n
1753 <<
" delta=" << new_node->values[i].delta;
1765 #ifndef ACCELERATED_INPUT_DIHEDRALS 1766 ptr->multiplicity = new_node->multiplicity;
1768 for (i=0; i<new_node->multiplicity; i++)
1770 ptr->values[i].k = new_node->values[i].k;
1771 ptr->values[i].n = new_node->values[i].n;
1772 ptr->values[i].delta = new_node->values[i].delta;
1779 retval.second->values[i].k = new_node->
value.
values[i].
k;
1780 retval.second->values[i].n = new_node->
value.
values[i].
n;
1792 #ifndef ACCELERATED_INPUT_DIHEDRALS 1797 #ifndef ACCELERATED_INPUT_DIHEDRALS 1804 if ( new_node->atom1wild ||
1805 new_node->atom2wild ||
1806 new_node->atom3wild ||
1807 new_node->atom4wild )
1810 tail->next=new_node;
1818 new_node->next=dihedralp;
1851 void Parameters::add_to_charmm_dihedral_list(
1857 #ifndef ACCELERATED_INPUT_DIHEDRALS 1870 if (dihedralp == NULL)
1884 int same_as_last = 0;
1885 if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
1886 (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
1887 (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
1888 (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
1889 ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
1890 (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
1891 (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
1892 (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) )
1895 int same_as_last = 0;
1896 if(dihedralmap==NULL)
1899 if(retval.first==
false)
1902 #ifndef ACCELERATED_INPUT_DIHEDRALS 1909 if ( ( !strcmp(ptr->atom1name, last_dihedral.atom1name) &&
1910 !strcmp(ptr->atom2name, last_dihedral.atom2name) &&
1911 !strcmp(ptr->atom3name, last_dihedral.atom3name) &&
1912 !strcmp(ptr->atom4name, last_dihedral.atom4name)))
1915 same_as_last = (lastDihedralKey==new_node->
key) ? 1 :0;
1920 #ifndef ACCELERATED_INPUT_DIHEDRALS 1922 for (i=0; i<ptr->multiplicity; i++)
1924 if ((ptr->values[i].k == new_node->values[0].k) &&
1925 (ptr->values[i].n == new_node->values[0].n) &&
1926 (ptr->values[i].delta == new_node->values[0].delta))
1935 for (i=0; i<retval.second->multiplicity; i++)
1937 if((retval.second->values[i].k == new_node->
value.
values[0].
k) &&
1938 (retval.second->values[i].n == new_node->
value.
values[0].
n) &&
1940 {echoWarn=0;
break;}
1946 if (!same_as_last) {
1947 #ifndef ACCELERATED_INPUT_DIHEDRALS 1949 iout <<
"\n" <<
iWARN <<
"DUPLICATE DIHEDRAL ENTRY FOR " 1950 << ptr->atom1name <<
"-" 1951 << ptr->atom2name <<
"-" 1952 << ptr->atom3name <<
"-" 1954 <<
"\nPREVIOUS VALUES MULTIPLICITY: " << ptr->multiplicity <<
"\n";
1956 iout <<
"\n" <<
iWARN <<
"DUPLICATE DIHEDRAL ENTRY FOR " 1965 #ifndef ACCELERATED_INPUT_DIHEDRALS 1966 for (i=0; i<ptr->multiplicity; i++)
1968 if (!same_as_last) {
1969 iout <<
iWARN <<
"IDENTICAL PERIODICITY! " 1970 << ptr->atom1name <<
"-" 1971 << ptr->atom2name <<
"-" 1972 << ptr->atom3name <<
"-" 1973 << ptr->atom4name <<
1974 " REPLACING OLD VALUES BY: " 1975 <<
" k=" << ptr->values[i].k
1976 <<
" n=" << ptr->values[i].n
1977 <<
" delta=" << ptr->values[i].delta <<
"\n";
1979 if (ptr->values[i].n == new_node->values[0].n)
1981 ptr->values[i].k = new_node->values[0].k;
1982 ptr->values[i].delta = new_node->values[0].delta;
1983 iout <<
iWARN <<
"IDENTICAL PERIODICITY! " 1984 << ptr->atom1name <<
"-" 1985 << ptr->atom2name <<
"-" 1986 << ptr->atom3name <<
"-" 1987 << ptr->atom4name <<
1988 " REPLACING OLD VALUES BY: " 1989 <<
" k=" << ptr->values[i].k
1990 <<
" n=" << ptr->values[i].n
1991 <<
" delta=" << ptr->values[i].delta<<
"\n";
1997 for (i=0; i<retval.second->multiplicity; i++)
1999 if (!same_as_last) {
2000 iout <<
" k=" << retval.second->values[i].k
2001 <<
" n=" << retval.second->values[i].n
2002 <<
" delta=" << retval.second->values[i].delta;
2004 if (retval.second->values[i].n == new_node->
value.
values[0].
n)
2006 iout <<
iWARN <<
"IDENTICAL PERIODICITY! " 2011 " REPLACING OLD VALUES BY: ";
2012 retval.second->values[i].k = new_node->
value.
values[0].
k;
2014 iout <<
" k=" << retval.second->values[i].k
2015 <<
" n=" << retval.second->values[i].n
2016 <<
" delta=" << retval.second->values[i].delta<<
"\n";
2024 #ifndef ACCELERATED_INPUT_DIHEDRALS 2025 ptr->multiplicity += 1;
2027 retval.second->multiplicity += 1;
2031 #ifndef ACCELERATED_INPUT_DIHEDRALS 2032 multiplicity=ptr->multiplicity;
2034 multiplicity=retval.second->multiplicity;
2040 snprintf(err_msg,
sizeof(err_msg),
2041 "Multiple dihedral with multiplicity " 2042 "of %d greater than max of %d",
2046 #ifndef ACCELERATED_INPUT_DIHEDRALS 2048 iout <<
"INCREASING MULTIPLICITY TO: " << ptr->multiplicity <<
"\n";
2049 i= ptr->multiplicity - 1;
2050 ptr->values[i].k = new_node->values[0].k;
2051 ptr->values[i].n = new_node->values[0].n;
2052 ptr->values[i].delta = new_node->values[0].delta;
2055 iout <<
" k=" << ptr->values[i].k
2056 <<
" n=" << ptr->values[i].n
2057 <<
" delta=" << ptr->values[i].delta<<
"\n";
2060 iout <<
"INCREASING MULTIPLICITY TO: " << retval.second->multiplicity <<
"\n";
2061 i= retval.second->multiplicity - 1;
2062 retval.second->values[i].k = new_node->
value.
values[0].
k;
2063 retval.second->values[i].n = new_node->
value.
values[0].
n;
2066 iout <<
" k=" << retval.second->values[i].k
2067 <<
" n=" << retval.second->values[i].n
2068 <<
" delta=" << retval.second->values[i].delta<<
"\n";
2074 #ifndef ACCELERATED_INPUT_DIHEDRALS 2077 lastDihedralKey= new_node->
key;
2083 #ifdef ACCELERATED_INPUT_DIHEDRALS 2086 lastDihedralKey= new_node->
key;
2089 #ifndef ACCELERATED_INPUT_DIHEDRALS 2093 #ifndef ACCELERATED_INPUT_DIHEDRALS 2101 if ( new_node->atom1wild ||
2102 new_node->atom2wild ||
2103 new_node->atom3wild ||
2104 new_node->atom4wild )
2107 tail->next=new_node;
2116 new_node->next=dihedralp;
2142 void Parameters::add_improper_param(
char *buf, FILE *fd)
2153 #ifndef ACCELERATED_INPUT_IMPROPERS 2169 read_count=sscanf(buf,
"%*s %s %s %s %s MULTIPLE= %d %f %d %f\n",
2171 &forceconstant, &periodicity, &phase_shift);
2176 read_count=sscanf(buf,
"%s %s %s %s %f %d %f\n",
2178 &forceconstant, &periodicity, &phase_shift);
2181 orig_phase = phase_shift;
2182 if ( (read_count != 4) && (read_count != 8) && (paramType ==
paraXplor) )
2186 snprintf(err_msg,
sizeof(err_msg),
2187 "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
2190 else if ( (read_count != 7) && (paramType ==
paraCharmm) )
2194 snprintf(err_msg,
sizeof(err_msg),
2195 "BAD IMPROPER FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2199 if ( (read_count == 4) && (paramType ==
paraXplor) )
2202 read_count=sscanf(buf,
"%*s %*s %*s %*s %*s %f %d %f\n",
2203 &forceconstant, &periodicity, &phase_shift);
2206 if (read_count != 3)
2210 snprintf(err_msg,
sizeof(err_msg),
2211 "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
2222 snprintf(err_msg,
sizeof(err_msg),
2223 "Multiple improper with multiplicity of %d greater than max of %d",
2227 #ifndef ACCELERATED_INPUT_IMPROPERS 2231 if (new_node == NULL)
2233 NAMD_die(
"memory allocation failed in Parameters::add_improper_param");
2243 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
2244 new_node->
values[0].
k = forceconstant;
2245 new_node->
values[0].
n = periodicity;
2248 new_node->
next = NULL;
2266 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
2267 value.
values[0].
k = forceconstant;
2268 value.
values[0].
n = periodicity;
2295 NAMD_die(
"EOF encoutner in middle of multiple improper");
2299 read_count=sscanf(buffer,
"%f %d %f\n",
2300 &forceconstant, &periodicity, &phase_shift);
2302 if (read_count != 3)
2306 snprintf(err_msg,
sizeof(err_msg),
2307 "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
2311 if (paramType ==
paraXplor && periodicity != 0) phase_shift *= -1;
2312 #ifndef ACCELERATED_INPUT_IMPROPERS 2313 new_node->
values[i].
k = forceconstant;
2314 new_node->
values[i].
n = periodicity;
2317 value.
values[i].
k = forceconstant;
2318 value.
values[i].
n = periodicity;
2324 #ifndef ACCELERATED_INPUT_IMPROPERS 2326 add_to_improper_list(new_node);
2328 if(impropermap==NULL)
2332 if(retval.first==
false)
2339 if (retval.second->multiplicity != value.
multiplicity) {echoWarn=1;}
2343 for (i=0; i<retval.second->multiplicity; i++)
2345 if (retval.second->values[i].k != value.
values[i].
k) {echoWarn=1;
break;}
2346 if (retval.second->values[i].n != value.
values[i].
n) {echoWarn=1;
break;}
2347 if (retval.second->values[i].delta != value.
values[i].
delta) {echoWarn=1;
break;}
2353 iout <<
"\n" <<
iWARN <<
"DUPLICATE IMPROPER DIHEDRAL ENTRY FOR " 2358 <<
"\nPREVIOUS VALUES MULTIPLICITY " << retval.second->multiplicity <<
"\n";
2360 for (i=0; i<retval.second->multiplicity; i++)
2362 iout <<
" k=" << retval.second->values[i].k
2363 <<
" n=" << retval.second->values[i].n
2364 <<
" delta=" << retval.second->values[i].delta;
2367 iout <<
"\n" <<
"USING VALUES MULTIPLICITY " << value.
multiplicity <<
"\n";
2382 retval.second->values[i].k = value.
values[i].
k;
2383 retval.second->values[i].n = value.
values[i].
n;
2384 retval.second->values[i].delta = value.
values[i].
delta;
2412 void Parameters::add_to_improper_list(
struct improper_params *new_node)
2423 if (improperp == NULL)
2456 if (ptr->
values[i].
k != new_node->
values[i].
k) {echoWarn=1;
break;}
2457 if (ptr->
values[i].
n != new_node->
values[i].
n) {echoWarn=1;
break;}
2464 iout <<
"\n" <<
iWARN <<
"DUPLICATE IMPROPER DIHEDRAL ENTRY FOR " 2469 <<
"\nPREVIOUS VALUES MULTIPLICITY " << ptr->
multiplicity <<
"\n";
2478 iout <<
"\n" <<
"USING VALUES MULTIPLICITY " << new_node->
multiplicity <<
"\n";
2483 <<
" n=" << new_node->
values[i].
n 2514 if ( (strcasecmp(new_node->
atom1name,
"X") == 0) ||
2515 (strcasecmp(new_node->
atom2name,
"X") == 0) ||
2516 (strcasecmp(new_node->
atom3name,
"X") == 0) ||
2517 (strcasecmp(new_node->
atom4name,
"X") == 0) )
2520 tail->
next=new_node;
2528 new_node->
next=improperp;
2552 void Parameters::add_crossterm_param(
char *buf, FILE *fd)
2566 #ifdef ACCELERATED_INPUT_CROSSTERMS 2574 read_count=sscanf(buf,
"%s %s %s %s %s %s %s %s %d\n",
2579 if ( (read_count != 9) || dimension < 1 || dimension > 1000 )
2583 snprintf(err_msg,
sizeof(err_msg),
2584 "BAD CMAP FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
2594 #ifndef ACCELERATED_INPUT_CROSSTERMS 2604 new_node->
next = NULL;
2614 int reverseBits= (cmp21==0) ? ( (cmp22 >0 ) ? 0x1 : 0x0 ) : ((cmp21>0) ? 0x1 : 0x0);
2615 int reverseBits2= (cmp11==0) ? ( (cmp12 >0 ) ? 0x10 : 0x00 ) : ((cmp11>0) ? 0x10 : 0x00);
2618 static const short o1order1order2=0x000;
2619 static const short r1order1order2=0x100;
2620 static const short o1order1reverse2=0x001;
2621 static const short r1order1reverse2=0x101;
2622 static const short o1reverse1order2=0x010;
2623 static const short r1reverse1order2=0x110;
2624 static const short o1reverse1reverse2=0x011;
2625 static const short r1reverse1reverse2=0x111;
2626 int reverseBits3= (cmpo1==0) ? ( (cmpo2 >0 ) ? 0x100 : 0x000 ) : ((cmpo1>0) ? 0x100 : 0x000);
2627 reverseBits|=reverseBits2 | reverseBits3;
2628 switch (reverseBits)
2630 case o1order1order2:
2634 case o1order1reverse2:
2638 case o1reverse1order2:
2642 case o1reverse1reverse2:
2646 case r1order1order2:
2650 case r1order1reverse2:
2654 case r1reverse1order2:
2658 case r1reverse1reverse2:
2663 NAMD_die(
"invalid crossterm ordering value");
2670 while ( nread < nterms ) {
2674 if (ret_code == 0) {
2681 if (ret_code == 0) {
2686 if (ret_code != 0) {
2687 NAMD_die(
"EOF encoutner in middle of CMAP");
2691 read_count=sscanf(buffer,
"%lf %lf %lf %lf %lf %lf %lf %lf\n",
2692 new_node->
values + nread,
2693 new_node->
values + nread+1,
2694 new_node->
values + nread+2,
2695 new_node->
values + nread+3,
2696 new_node->
values + nread+4,
2697 new_node->
values + nread+5,
2698 new_node->
values + nread+6,
2699 new_node->
values + nread+7);
2701 nread += read_count;
2703 if (read_count == 0 || nread > nterms) {
2706 snprintf(err_msg,
sizeof(err_msg),
2707 "BAD CMAP FORMAT IN PARAMETER FILE\nLINE=*%s*\n", buffer);
2711 #ifndef ACCELERATED_INPUT_CROSSTERMS 2714 add_to_crossterm_list(new_node);
2719 NAMD_die(
"Sorry, only CMAP dimension of 24 is supported");
2723 for (
int i=0; i<N; i++) {
2724 for (
int j=0; j<N; j++) {
2729 for (
int i=0; i<N; i++) {
2730 value.
c[i][N].
d00 = value.
c[i][0].
d00;
2731 value.
c[N][i].
d00 = value.
c[0][i].
d00;
2733 value.
c[N][N].
d00 = value.
c[0][0].
d00;
2736 if(crosstermmap==NULL)
2740 if(retval.first==
false)
2746 if (retval.second->dimension != new_node->
dimension) {echoWarn=1;}
2750 int nvals = retval.second->dimension * retval.second->dimension;
2752 for (
int i=0; i<N; i++) {
2753 for (
int j=0; j<N; j++) {
2754 if(retval.second->c[i][j].d00 != new_node->
values[k]) {echoWarn=1;
break;}
2761 iout <<
"\n" <<
iWARN <<
"DUPLICATE CMAP ENTRY FOR " 2773 retval.second->dimension = new_node->
dimension;
2774 int nvals = retval.second->dimension * retval.second->dimension;
2776 for (
int i=0; i<N; i++) {
2777 for (
int j=0; j<N; j++) {
2778 retval.second->c[i][j].d00 = new_node->
values[k];
2782 for (
int i=0; i<N; i++) {
2783 value.
c[i][N].
d00 = value.
c[i][0].
d00;
2784 value.
c[N][i].
d00 = value.
c[0][i].
d00;
2786 value.
c[N][N].
d00 = value.
c[0][0].
d00;
2825 if (crosstermp == NULL)
2827 crosstermp=new_node;
2856 for (i=0; i<nvals; i++)
2858 if (ptr->
values[i] != new_node->
values[i]) {echoWarn=1;
break;}
2864 iout <<
"\n" <<
iWARN <<
"DUPLICATE CMAP ENTRY FOR " 2872 << ptr->
atom8name <<
", USING NEW VALUES\n";
2880 new_node->
values = tmpvalues;
2897 if ( (strcasecmp(new_node->
atom1name,
"X") == 0) ||
2898 (strcasecmp(new_node->
atom2name,
"X") == 0) ||
2899 (strcasecmp(new_node->
atom3name,
"X") == 0) ||
2900 (strcasecmp(new_node->
atom4name,
"X") == 0) ||
2901 (strcasecmp(new_node->
atom5name,
"X") == 0) ||
2902 (strcasecmp(new_node->
atom6name,
"X") == 0) ||
2903 (strcasecmp(new_node->
atom7name,
"X") == 0) ||
2904 (strcasecmp(new_node->
atom8name,
"X") == 0) )
2907 tail->
next=new_node;
2915 new_node->
next=crosstermp;
2916 crosstermp=new_node;
2935 void Parameters::add_vdw_param(
char *buf)
2945 #ifndef ACCELERATED_INPUT_VDW 2956 read_count=sscanf(buf,
"%*s %s %f %f %f %f\n",
atomname,
2962 read_count=sscanf(buf,
"%s %*f %f %f %*f %f %f\n",
atomname,
2967 if ((read_count != 5) && (paramType ==
paraXplor))
2971 snprintf(err_msg,
sizeof(err_msg),
2972 "BAD vdW FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
2975 else if ( ((read_count != 5) && (read_count != 3)) && (paramType ==
paraCharmm))
2979 snprintf(err_msg,
sizeof(err_msg),
2980 "BAD vdW FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
2988 sqrt26=pow(2.,(1./6.));
2991 if (read_count == 3)
3006 iout <<
iWARN <<
"Ignoring VDW parameter with positive epsilon:\n" 3007 << buf <<
"\n" <<
endi;
3010 iout <<
iWARN <<
"Ignoring VDW parameter with negative epsilon:\n" 3011 << buf <<
"\n" <<
endi;
3015 #ifndef ACCELERATED_INPUT_VDW 3019 if (new_node == NULL)
3021 NAMD_die(
"memory allocation failed in Parameters::add_vdw_param");
3031 new_node->
left = NULL;
3032 new_node->
right = NULL;
3035 vdwp=add_to_vdw_tree(new_node, vdwp);
3044 auto retval = vdwmap->insert_check(key, value);
3045 if(retval.first ==
false)
3050 if ((retval.second->sigma !=
sigma) ||
3051 (retval.second->epsilon !=
epsilon) ||
3052 (retval.second->sigma14 !=
sigma14) ||
3053 (retval.second->epsilon14 !=
epsilon14))
3056 <<
"\nPREVIOUS VALUES sigma=" << retval.second->sigma
3057 <<
" epsilon=" << retval.second->epsilon
3058 <<
" sigma14=" << retval.second->sigma14
3059 <<
" epsilon14=" << retval.second->epsilon14
3060 <<
"\n USING VALUES sigma=" <<
sigma 3066 retval.second->sigma =
sigma;
3067 retval.second->epsilon =
epsilon;
3068 retval.second->sigma14 =
sigma14;
3107 if (compare_code==0)
3118 <<
"\nPREVIOUS VALUES sigma=" << tree->
sigma 3119 <<
" epsilon=" << tree->
epsilon 3120 <<
" sigma14=" << tree->
sigma14 3122 <<
"\n USING VALUES sigma=" << new_node->
sigma 3123 <<
" epsilon=" << new_node->
epsilon 3124 <<
" sigma14=" << new_node->
sigma14 3142 if (compare_code < 0)
3144 tree->
left = add_to_vdw_tree(new_node, tree->
left);
3148 tree->
right = add_to_vdw_tree(new_node, tree->
right);
3167 void Parameters::add_table_pair_param(
char *buf)
3175 #ifndef ACCELERATED_INPUT_TABLE_PAIR 3182 read_count=sscanf(buf,
"%s %s %s\n",
atom1name,
3186 if ((read_count != 3))
3190 snprintf(err_msg,
sizeof(err_msg),
3191 "BAD TABLE PAIR FORMAT IN PARAMETER FILE\nLINE=*%s*", buf);
3194 #ifndef ACCELERATED_INPUT_TABLE_PAIR 3198 if (new_node == NULL)
3200 NAMD_die(
"memory allocation failed in Parameters::add_table_pair_param\n");
3211 snprintf(err_msg,
sizeof(err_msg),
3212 "Couldn't find table parameters for table interaction %s!\n",
3220 new_node->
next = NULL;
3224 add_to_table_pair_list(new_node);
3233 snprintf(err_msg,
sizeof(err_msg),
3234 "Couldn't find table parameters for table interaction %s!\n",
3242 bool duplicate =
false;
3244 if(ret1 != vdwmap->tupleMap.end() && ret2 != vdwmap->tupleMap.end())
3249 uint64_t k1=ret1->second;
3250 uint64_t k2=ret2->second;
3258 uint64_t key= ((k1 <<32) | k2);
3260 auto retval= tablepairmap.emplace(std::make_pair(key, value));
3261 duplicate = !retval.second;
3262 dupvalue = &(retval.first->second);
3275 if (dupvalue->
K !=
K)
3277 iout <<
iWARN <<
"DUPLICATE table PAIR ENTRY FOR " 3280 <<
"\nPREVIOUS VALUES K=" << dupvalue->
K 3281 <<
"\n USING VALUES K=" <<
K 3304 void Parameters::add_vdw_pair_param(
char *buf)
3314 #ifndef ACCELERATED_INPUT_VDW 3324 read_count=sscanf(buf,
"%*s %s %s %f %f %f %f\n",
atom1name,
3329 Real well, rmin, well14, rmin14;
3331 read_count=sscanf(buf,
"%s %s %f %f %f %f\n",
atom1name,
3332 atom2name, &well, &rmin, &well14, &rmin14);
3333 if ( read_count == 4 ) { well14 = well; rmin14 = rmin; }
3334 A = -1. * well * pow(rmin, 12.);
3335 B = -2. * well * pow(rmin, 6.);
3336 A14 = -1. * well14 * pow(rmin14, 12.);
3337 B14 = -2. * well14 * pow(rmin14, 6.);
3341 if ((read_count != 6) && (paramType ==
paraXplor))
3345 snprintf(err_msg,
sizeof(err_msg),
3346 "BAD vdW PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
3349 if ((read_count != 4) && (read_count != 6) && (paramType ==
paraCharmm))
3353 snprintf(err_msg,
sizeof(err_msg),
3354 "BAD vdW PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
3358 #ifndef ACCELERATED_INPUT_VDW 3362 if (new_node == NULL)
3364 NAMD_die(
"memory allocation failed in Parameters::add_vdw_pair_param\n");
3376 new_node->
next = NULL;
3379 add_to_vdw_pair_list(new_node);
3387 bool duplicate =
false;
3389 if(ret1 != vdwmap->tupleMap.end() && ret2 != vdwmap->tupleMap.end())
3394 uint64_t k1=ret1->second;
3395 uint64_t k2=ret2->second;
3403 uint64_t key= ((k1 <<32) | k2);
3405 auto retval= vdwpairmap.emplace(std::make_pair(key, value));
3406 duplicate = !retval.second;
3407 dupvalue = &(retval.first->second);
3420 if ((dupvalue->
A !=
A) ||
3421 (dupvalue->
B !=
B) ||
3422 (dupvalue->
A14 !=
A14) ||
3425 iout <<
iWARN <<
"DUPLICATE vdW PAIR ENTRY FOR " 3428 <<
"\nPREVIOUS VALUES A=" << dupvalue->
A 3429 <<
" B=" << dupvalue->
B 3430 <<
" A14=" << dupvalue->
A14 3431 <<
" B14" << dupvalue->
B14 3432 <<
"\n USING VALUES A=" <<
A 3461 void Parameters::add_nbthole_pair_param(
char *buf)
3468 #ifndef ACCELERATED_INPUT_NBTHOLE 3475 read_count=sscanf(buf,
"%s %s %f\n",
atom1name,
3480 if ((read_count != 3) && (paramType ==
paraCharmm))
3484 snprintf(err_msg,
sizeof(err_msg),
3485 "BAD NBTHOLE PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
3489 #ifndef ACCELERATED_INPUT_NBTHOLE 3493 if (new_node == NULL)
3495 NAMD_die(
"memory allocation failed in Parameters::nbthole_pair_param\n");
3504 new_node->
next = NULL;
3507 add_to_nbthole_pair_list(new_node);
3514 bool duplicate =
false;
3516 if(ret1 != vdwmap->tupleMap.end() && ret2 != vdwmap->tupleMap.end())
3521 uint64_t k1=ret1->second;
3522 uint64_t k2=ret2->second;
3530 uint64_t key= ((k1 <<32) | k2);
3532 auto retval= nbtholepairmap.emplace(std::make_pair(key, value));
3533 duplicate = !retval.second;
3534 dupvalue = &(retval.first->second);
3549 iout <<
iWARN <<
"DUPLICATE nbthole PAIR ENTRY FOR " 3552 <<
"\nPREVIOUS VALUES tholeij=" << dupvalue->
tholeij 3553 <<
"\n USING VALUES tholeij=" <<
tholeij 3576 void Parameters::add_hb_pair_param(
char *buf)
3588 if (sscanf(buf,
"%*s %s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
3590 snprintf(err_msg,
sizeof(err_msg),
3591 "BAD HBOND PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
3596 if (sscanf(buf,
"%s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
3598 snprintf(err_msg,
sizeof(err_msg),
3599 "BAD HBOND PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
3606 if (hbondParams.add_hbond_pair(a1n, a2n, A, B) ==
FALSE) {
3607 iout <<
"\n" <<
iWARN <<
"Duplicate HBOND parameters for types " << a1n
3608 <<
" and " << a2n <<
" found; using latest values." <<
"\n" <<
endi;
3634 if (table_pairp == NULL)
3636 table_pairp = new_node;
3649 if (compare_code == 0)
3654 if (compare_code==0)
3659 iout <<
iWARN <<
"DUPLICATE TABLE PAIR ENTRY FOR " 3677 tail->
next = new_node;
3693 void Parameters::add_to_vdw_pair_list(
struct vdw_pair_params *new_node)
3702 if (vdw_pairp == NULL)
3704 vdw_pairp = new_node;
3717 if (compare_code == 0)
3722 if (compare_code==0)
3727 if ((ptr->
A != new_node->
A) ||
3728 (ptr->
B != new_node->
B) ||
3729 (ptr->
A14 != new_node->
A14) ||
3730 (ptr->
B14 != new_node->
B14))
3732 iout <<
iWARN <<
"DUPLICATE vdW PAIR ENTRY FOR " 3735 <<
"\nPREVIOUS VALUES A=" << ptr->
A 3737 <<
" A14=" << ptr->
A14 3738 <<
" B14" << ptr->
B14 3739 <<
"\n USING VALUES A=" << new_node->
A 3740 <<
" B=" << new_node->
B 3741 <<
" A14=" << new_node->
A14 3742 <<
" B14" << new_node->
B14 3762 tail->
next = new_node;
3787 if (nbthole_pairp == NULL)
3789 nbthole_pairp = new_node;
3794 ptr = nbthole_pairp;
3796 tail->
next = new_node;
3816 AllFilesRead =
TRUE;
3821 add_bond_param(
"X DRUD 500.0 0.0\n",
FALSE);
3830 NAMD_die(
"memory allocation of bond_array failed!");
3841 NAMD_die(
"memory allocation of angle_array failed!");
3855 NAMD_die(
"memory allocation of dihedral_array failed!");
3866 NAMD_die(
"memory allocation of improper_array failed!");
3892 NAMD_die(
"memory allocation of vdw_array failed!");
3901 NAMD_die(
"memory allocation of nbthole_array failed!");
3915 index_bonds(bondp, 0);
3916 index_angles(anglep, 0);
3921 convert_nbthole_pairs();
3923 convert_vdw_pairs();
3925 convert_table_pairs();
3947 #ifndef ACCELERATED_INPUT 3953 if (tree->left != NULL)
3965 if (tree->right != NULL)
3999 #ifndef ACCELERATED_INPUT_ANGLES 4005 if (tree->left != NULL)
4023 if (tree->right != NULL)
4049 void Parameters::index_dihedrals()
4051 #ifndef ACCELERATED_INPUT_DIHEDRALS 4073 if (maxDihedralMults == NULL)
4075 NAMD_die(
"memory allocation failed in Parameters::index_dihedrals()");
4086 maxDihedralMults[index] = ptr->multiplicity;
4105 for (i=0; i<ptr->multiplicity; i++)
4124 dihedralmap->
sort();
4149 void Parameters::index_impropers()
4152 #ifndef ACCELERATED_INPUT_IMPROPERS 4174 if (maxImproperMults == NULL)
4176 NAMD_die(
"memory allocation failed in Parameters::index_impropers()");
4224 impropermap->
sort();
4253 void Parameters::index_crossterms()
4256 #ifndef ACCELERATED_INPUT_CROSSTERMS 4271 NAMD_die(
"Sorry, only CMAP dimension of 24 is supported");
4275 for (i=0; i<N; i++) {
4276 for (j=0; j<N; j++) {
4281 for (i=0; i<N; i++) {
4301 crosstermmap->
sort();
4326 #ifndef ACCELERATED_INPUT_VDW 4332 if (tree->
left != NULL)
4355 if (tree->
right != NULL)
4365 memcpy(
vdw_array, vdwmap->paramVector.data(),
sizeof(
VdwValue)*vdwmap->paramVector.size());
4366 index=vdwmap->paramVector.size();
4396 #ifdef MEM_OPT_VERSION 4403 #ifndef ACCELERATED_INPUT_VDW 4411 NAMD_die(
"Tried to assign vdw index before all parameter files were read");
4421 while (!found && (ptr!=NULL))
4423 comp_code = strcasecmp(atomtype, ptr->
atomname);
4431 else if (comp_code < 0)
4444 int64_t result = vdwmap->index(searchFor);
4455 #ifndef ACCELERATED_INPUT_VDW 4465 while (!found && (ptr!=NULL))
4470 if (windx == strlen(ptr->
atomname))
4473 comp_code = strcasecmp(atomtype, ptr->
atomname);
4477 comp_code = strncasecmp(atomtype, ptr->
atomname, windx);
4486 snprintf(errbuf,
sizeof(errbuf),
4487 "VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
4490 for(i=0; i<error_msgs.
size(); i++) {
4491 if ( strcmp(errbuf,error_msgs[i]) == 0 )
break;
4493 if ( i == error_msgs.
size() ) {
4494 char *newbuf =
new char[strlen(errbuf)+1];
4495 strcpy(newbuf,errbuf);
4496 error_msgs.
add(newbuf);
4500 else if (comp_code < 0)
4514 int64_t result = vdwmap->index(searchFor);
4529 snprintf(err_msg,
sizeof(err_msg),
4530 "DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s", atomtype);
4571 #ifndef ACCELERATED_INPUT_VDW 4574 while (!found && (ptr!=NULL))
4577 if ( (ind1 == ptr->
ind1) && (ind2 == ptr->
ind2) )
4581 else if ( (ind1 < ptr->ind1) ||
4582 ( (ind1==ptr->
ind1) && (ind2 < ptr->ind2) ) )
4607 uint64_t key= ((k1 <<32) | k2);
4608 auto ret = tablepairmap.find(key);
4609 if(ret!=tablepairmap.end())
4665 #ifndef ACCELERATED_INPUT_VDW 4668 while (!found && (ptr!=NULL))
4670 if ( (ind1 == ptr->
ind1) && (ind2 == ptr->
ind2) )
4674 else if ( (ind1 < ptr->ind1) ||
4675 ( (ind1==ptr->
ind1) && (ind2 < ptr->ind2) ) )
4703 uint64_t key= ((k1 <<32) | k2);
4704 auto ret = vdwpairmap.find(key);
4705 if(ret!=vdwpairmap.end())
4709 *A14 = ret->second.A14;
4710 *B14 = ret->second.B14;
4751 NAMD_die(
"Tried to assign bond index before all parameter files were read");
4756 if (strcasecmp(atom1, atom2) > 0)
4758 const char *tmp_name = atom1;
4762 #ifndef ACCELERATED_INPUT 4768 while (!found && (ptr!=NULL))
4770 cmp_code=strcasecmp(atom1, ptr->
atom1name);
4774 cmp_code=strcasecmp(atom2, ptr->
atom2name);
4783 else if (cmp_code < 0)
4796 int64_t result= bondmap->
index(searchFor);
4806 if ((strcmp(atom1,
"DRUD")==0 || strcmp(atom2,
"DRUD")==0)
4807 && (strcmp(atom1,
"X")!=0 && strcmp(atom2,
"X")!=0)) {
4809 char a1[8] =
"DRUD", a2[8] =
"X";
4812 else if (bond_found ==
nullptr) {
4815 snprintf(err_msg,
sizeof(err_msg),
4816 "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
4817 atom1, atom2, bond_ptr->
atom1+1, bond_ptr->
atom2+1);
4821 if (bond_found !=
nullptr) {
4822 *bond_found = bool(found);
4850 Angle *angle_ptr,
int notFoundIndex)
4860 NAMD_die(
"Tried to assign angle index before all parameter files were read");
4865 if (strcasecmp(atom1, atom3) > 0)
4867 const char *tmp_name = atom1;
4872 #ifndef ACCELERATED_INPUT_ANGLES 4878 while (!found && (ptr != NULL))
4880 comp_val = strcasecmp(atom1, ptr->
atom1name);
4885 comp_val = strcasecmp(atom2, ptr->
atom2name);
4890 comp_val = strcasecmp(atom3, ptr->
atom3name);
4900 else if (comp_val < 0)
4913 int64_t result = anglemap->
index(searchFor);
4926 snprintf(err_msg,
sizeof(err_msg),
4927 "UNABLE TO FIND ANGLE PARAMETERS FOR %s %s %s (ATOMS %i %i %i)",
4928 atom1, atom2, atom3,
4930 if ( notFoundIndex ) {
4963 const char *atom4,
Dihedral *dihedral_ptr,
4964 int multiplicity,
int notFoundIndex)
4969 #ifndef ACCELERATED_INPUT_DIHEDRALS 4976 while (!found && (ptr!=NULL))
4986 if ( ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom1)==0) ) &&
4987 ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom2)==0) ) &&
4988 ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom3)==0) ) &&
4989 ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom4)==0) ) )
4994 else if ( ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom1)==0) ) &&
4995 ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom2)==0) ) &&
4996 ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom3)==0) ) &&
4997 ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom4)==0) ) )
5009 int cmp= strcasecmp(atom1, atom4);
5013 if(strcasecmp(atom2, atom3)>0)
5024 int64_t result = dihedralmap->
index(searchFor);
5036 searchForWR[2] =
TupleString4(atom4, atom3, atom2,
"X");
5037 searchForWR[3] =
TupleString4(atom4, atom3,
"X", atom1);
5038 searchForWR[4] =
TupleString4(atom4,
"X", atom2, atom1);
5039 searchForWR[5] =
TupleString4(
"X", atom3, atom2, atom1);
5047 searchForW[2] =
TupleString4(
"X", atom2, atom3, atom4);
5048 searchForW[3] =
TupleString4(atom1,
"X", atom3, atom4);
5049 searchForW[4] =
TupleString4(atom1, atom2,
"X", atom4);
5050 searchForW[5] =
TupleString4(atom1, atom2, atom3,
"X");
5057 for(
int i=0; i<11; ++i)
5059 result = dihedralmap->
index(searchForW[i]);
5069 for(
int i=0; i<11; ++i)
5071 result = dihedralmap->
index(searchForWR[i]);
5086 #ifdef ACCELERATED_INPUT_DIHEDRALS 5089 snprintf(err_msg,
sizeof(err_msg),
5090 "UNABLE TO FIND DIHEDRAL PARAMETERS FOR %s %s %s %s " 5091 "(ATOMS %i %i %i %i)",
5092 atom4, atom3, atom2, atom1,
5093 dihedral_ptr->
atom4+1, dihedral_ptr->
atom3+1,
5094 dihedral_ptr->
atom2+1, dihedral_ptr->
atom1+1);
5098 { snprintf(err_msg,
sizeof(err_msg),
5099 "UNABLE TO FIND DIHEDRAL PARAMETERS FOR %s %s %s %s " 5100 "(ATOMS %i %i %i %i)",
5101 atom1, atom2, atom3, atom4,
5102 dihedral_ptr->
atom1+1, dihedral_ptr->
atom2+1,
5103 dihedral_ptr->
atom3+1, dihedral_ptr->
atom4+1);
5105 if ( notFoundIndex ) {
5112 #ifndef ACCELERATED_INPUT_DIHEDRALS 5113 int index = ptr->index;
5121 if (multiplicity > maxDihedralMults[index])
5125 snprintf(err_msg,
sizeof(err_msg),
5126 "Multiplicity of Paramters for dihedral bond %s %s %s %s " 5128 atom1, atom2, atom3, atom4, maxDihedralMults[index]);
5168 const char *atom4,
Improper *improper_ptr,
5173 #ifndef ACCELERATED_INPUT_IMPROPERS 5182 while (!found && (ptr!=NULL))
5192 if ( ( (strcasecmp(ptr->
atom1name, atom1)==0) ||
5193 (strcasecmp(ptr->
atom1name,
"X")==0) ) &&
5194 ( (strcasecmp(ptr->
atom2name, atom2)==0) ||
5195 (strcasecmp(ptr->
atom2name,
"X")==0) ) &&
5196 ( (strcasecmp(ptr->
atom3name, atom3)==0) ||
5197 (strcasecmp(ptr->
atom3name,
"X")==0) ) &&
5198 ( (strcasecmp(ptr->
atom4name, atom4)==0) ||
5199 (strcasecmp(ptr->
atom4name,
"X")==0) ) )
5204 else if ( ( (strcasecmp(ptr->
atom4name, atom1)==0) ||
5205 (strcasecmp(ptr->
atom4name,
"X")==0) ) &&
5206 ( (strcasecmp(ptr->
atom3name, atom2)==0) ||
5207 (strcasecmp(ptr->
atom3name,
"X")==0) ) &&
5208 ( (strcasecmp(ptr->
atom2name, atom3)==0) ||
5209 (strcasecmp(ptr->
atom2name,
"X")==0) ) &&
5210 ( (strcasecmp(ptr->
atom1name, atom4)==0) ||
5211 (strcasecmp(ptr->
atom1name,
"X")==0) ) )
5223 int cmp= strcasecmp(atom1, atom4);
5227 if(strcasecmp(atom2, atom3)>0)
5237 int64_t result = impropermap->
index(searchFor);
5249 searchForWR[2] =
TupleString4(atom4, atom3, atom2,
"X");
5250 searchForWR[3] =
TupleString4(atom4, atom3,
"X", atom1);
5251 searchForWR[4] =
TupleString4(atom4,
"X", atom2, atom1);
5252 searchForWR[5] =
TupleString4(
"X", atom3, atom2, atom1);
5262 searchForW[2] =
TupleString4(
"X", atom2, atom3, atom4);
5263 searchForW[3] =
TupleString4(atom1,
"X", atom3, atom4);
5264 searchForW[4] =
TupleString4(atom1, atom2,
"X", atom4);
5265 searchForW[5] =
TupleString4(atom1, atom2, atom3,
"X");
5273 for(
int i=0; i<13; ++i)
5275 result = impropermap->
index(searchForW[i]);
5284 for(
int i=0; i<13; ++i)
5286 result = impropermap->
index(searchForWR[i]);
5301 snprintf(err_msg,
sizeof(err_msg),
5302 "UNABLE TO FIND IMPROPER PARAMETERS FOR %s %s %s %s " 5303 "(ATOMS %i %i %i %i)",
5304 atom1, atom2, atom3, atom4,
5305 improper_ptr->
atom1+1, improper_ptr->
atom2+1,
5306 improper_ptr->
atom3+1, improper_ptr->
atom4+1);
5310 #ifndef ACCELERATED_INPUT_IMPROPERS 5324 snprintf(err_msg,
sizeof(err_msg),
5325 "Multiplicity of Paramters for improper bond %s %s %s %s " 5327 atom1, atom2, atom3, atom4, maxImproperMults[
index]);
5353 const char *atom4,
const char *atom5,
const char *atom6,
const char *atom7,
5354 const char *atom8,
Crossterm *crossterm_ptr)
5357 #ifndef ACCELERATED_INPUT_CROSSTERMS 5365 while (!found && (ptr!=NULL))
5375 if ( ( (strcasecmp(ptr->
atom1name, atom1)==0) ||
5376 (strcasecmp(ptr->
atom1name,
"X")==0) ) &&
5377 ( (strcasecmp(ptr->
atom2name, atom2)==0) ||
5378 (strcasecmp(ptr->
atom2name,
"X")==0) ) &&
5379 ( (strcasecmp(ptr->
atom3name, atom3)==0) ||
5380 (strcasecmp(ptr->
atom3name,
"X")==0) ) &&
5381 ( (strcasecmp(ptr->
atom4name, atom4)==0) ||
5382 (strcasecmp(ptr->
atom4name,
"X")==0) ) )
5387 else if ( ( (strcasecmp(ptr->
atom4name, atom1)==0) ||
5388 (strcasecmp(ptr->
atom4name,
"X")==0) ) &&
5389 ( (strcasecmp(ptr->
atom3name, atom2)==0) ||
5390 (strcasecmp(ptr->
atom3name,
"X")==0) ) &&
5391 ( (strcasecmp(ptr->
atom2name, atom3)==0) ||
5392 (strcasecmp(ptr->
atom2name,
"X")==0) ) &&
5393 ( (strcasecmp(ptr->
atom1name, atom4)==0) ||
5394 (strcasecmp(ptr->
atom1name,
"X")==0) ) )
5405 if ( ( (strcasecmp(ptr->
atom5name, atom5)==0) ||
5406 (strcasecmp(ptr->
atom5name,
"X")==0) ) &&
5407 ( (strcasecmp(ptr->
atom6name, atom6)==0) ||
5408 (strcasecmp(ptr->
atom6name,
"X")==0) ) &&
5409 ( (strcasecmp(ptr->
atom7name, atom7)==0) ||
5410 (strcasecmp(ptr->
atom7name,
"X")==0) ) &&
5411 ( (strcasecmp(ptr->
atom8name, atom8)==0) ||
5412 (strcasecmp(ptr->
atom8name,
"X")==0) ) )
5417 else if ( ( (strcasecmp(ptr->
atom8name, atom5)==0) ||
5418 (strcasecmp(ptr->
atom8name,
"X")==0) ) &&
5419 ( (strcasecmp(ptr->
atom7name, atom6)==0) ||
5420 (strcasecmp(ptr->
atom7name,
"X")==0) ) &&
5421 ( (strcasecmp(ptr->
atom6name, atom7)==0) ||
5422 (strcasecmp(ptr->
atom6name,
"X")==0) ) &&
5423 ( (strcasecmp(ptr->
atom5name, atom8)==0) ||
5424 (strcasecmp(ptr->
atom5name,
"X")==0) ) )
5437 int cmp11 = strcasecmp(atom1, atom4);
5438 int cmp12 = strcasecmp(atom2, atom3);
5439 int cmp21 = strcasecmp(atom5, atom8);
5440 int cmp22 = strcasecmp(atom6, atom7);
5443 int reverseBits= (cmp21==0) ? ( (cmp22 >0 ) ? 0x1 : 0x0 ) : ((cmp21>0) ? 0x1 : 0x0);
5444 int reverseBits2= (cmp11==0) ? ( (cmp12 >0 ) ? 0x10 : 0x00 ) : ((cmp11>0) ? 0x10 : 0x00);
5445 int cmpo1 = ( reverseBits) ? strcasecmp(atom4, atom8) : strcasecmp(atom1, atom5);
5446 int cmpo2 = (reverseBits2 ) ? strcasecmp(atom3, atom7) : strcasecmp(atom2, atom6);
5447 static const short o1order1order2=0x000;
5448 static const short r1order1order2=0x100;
5449 static const short o1order1reverse2=0x001;
5450 static const short r1order1reverse2=0x101;
5451 static const short o1reverse1order2=0x010;
5452 static const short r1reverse1order2=0x110;
5453 static const short o1reverse1reverse2=0x011;
5454 static const short r1reverse1reverse2=0x111;
5455 int reverseBits3= (cmpo1==0) ? ( (cmpo2 >0 ) ? 0x100 : 0x000 ) : ((cmpo1>0) ? 0x100 : 0x000);
5456 reverseBits|=reverseBits2 | reverseBits3;
5458 switch (reverseBits)
5460 case o1order1order2:
5462 atom5, atom6, atom7, atom8);
5464 case o1order1reverse2:
5466 atom8, atom7, atom6, atom5);
5468 case o1reverse1order2:
5470 atom5, atom6, atom7, atom8);
5472 case o1reverse1reverse2:
5474 atom8, atom7, atom6, atom5);
5476 case r1order1order2:
5478 atom1, atom2, atom3, atom4);
5480 case r1order1reverse2:
5482 atom1, atom2, atom3, atom4);
5484 case r1reverse1order2:
5486 atom4, atom3, atom2, atom1);
5488 case r1reverse1reverse2:
5490 atom4, atom3, atom2, atom1);
5493 NAMD_die(
"invalid crossterm ordering value");
5496 int64_t result = crosstermmap->
index(searchFor);
5505 switch (reverseBits)
5507 case o1order1order2 :
5508 searchForW[0] =
TupleString8(
"X", atom2, atom3, atom4, atom5, atom6, atom7,
"X");
5509 searchForW[1] =
TupleString8(atom1,
"X", atom3, atom4, atom5, atom6,
"X", atom8);
5510 searchForW[2] =
TupleString8(atom1, atom2,
"X", atom4, atom5,
"X", atom7, atom8);
5511 searchForW[3] =
TupleString8(atom1, atom2, atom3,
"X",
"X", atom6, atom7, atom8);
5512 searchForW[4] =
TupleString8(
"X", atom2, atom3,
"X",
"X", atom6, atom7,
"X");
5513 searchForW[5] =
TupleString8(atom1,
"X",
"X", atom4, atom5,
"X",
"X", atom8);
5515 case o1order1reverse2 :
5516 searchForW[0] =
TupleString8(
"X", atom2, atom3, atom4, atom8, atom7, atom6,
"X");
5517 searchForW[1] =
TupleString8(atom1,
"X", atom3, atom4, atom8, atom7,
"X", atom5);
5518 searchForW[2] =
TupleString8(atom1, atom2,
"X", atom4, atom8,
"X", atom6, atom5);
5519 searchForW[3] =
TupleString8(atom1, atom2, atom3,
"X",
"X", atom7, atom6, atom5);
5520 searchForW[4] =
TupleString8(
"X", atom2, atom3,
"X",
"X", atom7, atom6,
"X");
5521 searchForW[5] =
TupleString8(atom1,
"X",
"X", atom4, atom8,
"X",
"X", atom5);
5523 case o1reverse1order2 :
5524 searchForW[0] =
TupleString8(
"X", atom3, atom2, atom1, atom5, atom6, atom7,
"X");
5525 searchForW[1] =
TupleString8(atom4,
"X", atom2, atom1, atom5, atom6,
"X", atom8);
5526 searchForW[2] =
TupleString8(atom4, atom3,
"X", atom1, atom5,
"X", atom7, atom8);
5527 searchForW[3] =
TupleString8(atom4, atom3, atom2,
"X",
"X", atom6, atom7, atom8);
5528 searchForW[4] =
TupleString8(
"X", atom3, atom2,
"X",
"X", atom6, atom7,
"X");
5529 searchForW[5] =
TupleString8(atom4,
"X",
"X", atom1, atom5,
"X",
"X", atom8);
5531 case o1reverse1reverse2 :
5532 searchForW[0] =
TupleString8(
"X", atom3, atom2, atom1, atom8, atom7, atom6,
"X");
5533 searchForW[1] =
TupleString8(atom4,
"X", atom2, atom1, atom8, atom7,
"X", atom5);
5534 searchForW[2] =
TupleString8(atom4, atom3,
"X", atom1, atom8,
"X", atom6, atom5);
5535 searchForW[3] =
TupleString8(atom4, atom3, atom2,
"X",
"X", atom7, atom6, atom5);
5536 searchForW[4] =
TupleString8(
"X", atom3, atom2,
"X",
"X", atom7, atom6,
"X");
5537 searchForW[5] =
TupleString8(atom4,
"X",
"X", atom1, atom8,
"X",
"X", atom5);
5539 case r1reverse1reverse2:
5540 searchForW[0] =
TupleString8(
"X", atom7, atom6, atom5, atom4, atom3, atom2,
"X");
5541 searchForW[1] =
TupleString8(atom8,
"X", atom6, atom5, atom4, atom3,
"X", atom1);
5542 searchForW[2] =
TupleString8(atom8, atom7,
"X", atom5, atom4,
"X", atom2, atom1);
5543 searchForW[3] =
TupleString8(atom8, atom7, atom6,
"X",
"X", atom3, atom2, atom1);
5544 searchForW[4] =
TupleString8(
"X", atom7, atom6,
"X",
"X", atom3, atom2,
"X");
5545 searchForW[5] =
TupleString8(atom8,
"X",
"X", atom5, atom4,
"X",
"X", atom1);
5547 case r1reverse1order2:
5548 searchForW[0] =
TupleString8(
"X", atom7, atom6, atom5, atom1, atom2, atom3,
"X");
5549 searchForW[1] =
TupleString8(atom8,
"X", atom6, atom5, atom1, atom2,
"X", atom4);
5550 searchForW[2] =
TupleString8(atom8, atom7,
"X", atom5, atom1,
"X", atom3, atom4);
5551 searchForW[3] =
TupleString8(atom8, atom7, atom6,
"X",
"X", atom2, atom3, atom4);
5552 searchForW[4] =
TupleString8(
"X", atom7, atom6,
"X",
"X", atom2, atom3,
"X");
5553 searchForW[5] =
TupleString8(atom8,
"X",
"X", atom5, atom1,
"X",
"X", atom4);
5555 case r1order1reverse2:
5556 searchForW[0] =
TupleString8(
"X", atom6, atom7, atom8, atom4, atom3, atom2,
"X");
5557 searchForW[1] =
TupleString8(atom5,
"X", atom7, atom8, atom4, atom3,
"X", atom1);
5558 searchForW[2] =
TupleString8(atom5, atom6,
"X", atom8, atom4,
"X", atom2, atom1);
5559 searchForW[3] =
TupleString8(atom5, atom6, atom7,
"X",
"X", atom3, atom2, atom1);
5560 searchForW[4] =
TupleString8(
"X", atom6, atom7,
"X",
"X", atom3, atom2,
"X");
5561 searchForW[5] =
TupleString8(atom5,
"X",
"X", atom8, atom4,
"X",
"X", atom1);
5564 NAMD_die(
"invalid reverse key order for crossterms");
5566 for(
int i=0; i<7; ++i)
5568 result = crosstermmap->
index(searchForW[i]);
5584 snprintf(err_msg,
sizeof(err_msg),
5585 "UNABLE TO FIND CROSSTERM PARAMETERS FOR " 5586 "%s %s %s %s %s %s %s %s\n" 5587 "(ATOMS %i %i %i %i %i %i %i %i)",
5588 atom1, atom2, atom3, atom4, atom5, atom6, atom7, atom8,
5589 crossterm_ptr->
atom1+1, crossterm_ptr->
atom2+1,
5590 crossterm_ptr->
atom3+1, crossterm_ptr->
atom4+1,
5591 crossterm_ptr->
atom5+1, crossterm_ptr->
atom6+1,
5592 crossterm_ptr->
atom7+1, crossterm_ptr->
atom8+1);
5597 #ifndef ACCELERATED_INPUT_CROSSTERMS 5619 void Parameters::free_bond_tree(
struct bond_params *bond_ptr)
5622 #ifndef ACCELERATED_INPUT 5623 if (bond_ptr->left != NULL)
5625 free_bond_tree(bond_ptr->left);
5628 if (bond_ptr->right != NULL)
5630 free_bond_tree(bond_ptr->right);
5656 void Parameters::free_angle_tree(
struct angle_params *angle_ptr)
5659 #ifndef ACCELERATED_INPUT_ANGLES 5660 if (angle_ptr->left != NULL)
5662 free_angle_tree(angle_ptr->left);
5665 if (angle_ptr->right != NULL)
5667 free_angle_tree(angle_ptr->right);
5694 #ifndef ACCELERATED_INPUT_DIHEDRALS 5707 dihedralmap->
clear();
5729 #ifndef ACCELERATED_INPUT_IMPROPERS 5742 impropermap->
clear();
5764 #ifndef ACCELERATED_INPUT_CROSSTERMS 5777 crosstermmap->
clear();
5778 delete crosstermmap;
5799 void Parameters::free_vdw_tree(
struct vdw_params *vdw_ptr)
5802 #ifndef ACCELERATED_INPUT_VDW 5803 if (vdw_ptr->
left != NULL)
5805 free_vdw_tree(vdw_ptr->
left);
5808 if (vdw_ptr->
right != NULL)
5810 free_vdw_tree(vdw_ptr->
right);
5830 void Parameters::free_vdw_pair_list()
5832 #ifndef ACCELERATED_INPUT_VDW 5859 void Parameters::free_nbthole_pair_list()
5861 #ifndef ACCELERATED_INPUT_NBTHOLE 5875 nbthole_pairp = NULL;
5887 #ifndef ACCELERATED_INPUT_TABLE_PAIR 5888 if (table_pair_ptr->left != NULL)
5890 free_table_pair_tree(table_pair_ptr->left);
5893 if (table_pair_ptr->right != NULL)
5895 free_table_pair_tree(table_pair_ptr->right);
5898 delete table_pair_ptr;
5902 tablepairmap.clear();
5921 void Parameters::free_vdw_pair_tree(
IndexedVdwPair *vdw_pair_ptr)
5924 #ifndef ACCELERATED_INPUT_VDW 5925 if (vdw_pair_ptr->
left != NULL)
5927 free_vdw_pair_tree(vdw_pair_ptr->
left);
5930 if (vdw_pair_ptr->
right != NULL)
5932 free_vdw_pair_tree(vdw_pair_ptr->
right);
5935 delete vdw_pair_ptr;
5961 #ifndef ACCELERATED_INPUT_NBTHOLE 5962 if (nbthole_pair_ptr->
left != NULL)
5964 free_nbthole_pair_tree(nbthole_pair_ptr->
left);
5967 if (nbthole_pair_ptr->
right != NULL)
5969 free_nbthole_pair_tree(nbthole_pair_ptr->
right);
5972 delete nbthole_pair_ptr;
5975 nbtholepairmap.clear();
5993 void Parameters::traverse_bond_params(
struct bond_params *tree)
5996 #ifndef ACCELERATED_INPUT 6000 if (tree->left != NULL)
6002 traverse_bond_params(tree->left);
6007 <<
" x0=" << tree->
distance<<
"\n");
6009 if (tree->right != NULL)
6011 traverse_bond_params(tree->right);
6014 for(std::pair<TupleString2, size_t> apair : bondmap->
tupleMap)
6016 DebugM(3,
"BOND " << apair.first.getTuplePtr(0) <<
" " << apair.first.getTuplePtr(1) \
6017 <<
" index=" << apair.second <<
" k=" << bondmap->
paramVector[apair.second].k \
6018 <<
" x0=" << bondmap->
paramVector[apair.second].x0<<
"\n");
6037 void Parameters::traverse_angle_params(
struct angle_params *tree)
6040 #ifndef ACCELERATED_INPUT_ANGLES 6044 if (tree->left != NULL)
6046 traverse_angle_params(tree->left);
6051 <<
" k_ub = " << tree->
k_ub <<
" r_ub = " << tree->
r_ub \
6052 <<
" normal = " << tree->
normal \
6056 if (tree->right != NULL)
6058 traverse_angle_params(tree->right);
6061 for(std::pair<TupleString3, size_t> apair : anglemap->
tupleMap)
6063 DebugM(3,
"ANGLE " << apair.first.getTuplePtr(0) <<
" " << apair.first.getTuplePtr(1)
6064 <<
" " << apair.first.getTuplePtr(2) <<
" index = " << apair.second
6065 <<
" k = " << anglemap->
paramVector[apair.second].k
6066 <<
" theta0 = " << anglemap->
paramVector[apair.second].theta0
6067 <<
" k_ub = " << anglemap->
paramVector[apair.second].k_ub
6068 <<
" r_ub = " << anglemap->
paramVector[apair.second].r_ub
6069 <<
" normal = " << anglemap->
paramVector[apair.second].normal
6089 void Parameters::traverse_dihedral_params(
struct dihedral_params *list)
6093 #ifndef ACCELERATED_INPUT_DIHEDRALS 6094 while (list != NULL)
6096 DebugM(3,
"DIHEDRAL " << list->atom1name <<
" " \
6097 << list->atom2name <<
" " << list->atom3name \
6098 <<
" " << list->atom4name <<
" index=" \
6100 <<
" multiplicity=" << list->multiplicity);
6102 for (i=0; i<list->multiplicity; i++)
6104 iout <<
" k=" << list->values[i].k \
6105 <<
" n=" << list->values[i].n \
6106 <<
" delta=" << list->values[i].delta;
6113 for(std::pair<TupleString4, size_t> apair : dihedralmap->
tupleMap)
6115 DebugM(3,
"DIHEDRAL " << apair.first.getTuplePtr(0) <<
" " \
6116 << apair.first.getTuplePtr(1) <<
" " << apair.first.getTuplePtr(2) \
6117 <<
" " << apair.first.getTuplePtr(3) <<
" index=" \
6119 <<
" multiplicity=" << dihedralmap->
paramVector[apair.second].multiplicity);
6121 for (i=0; i< dihedralmap->
paramVector[apair.second].multiplicity; i++)
6124 <<
" n=" << dihedralmap->
paramVector[apair.second].values[i].n \
6125 <<
" delta=" << dihedralmap->
paramVector[apair.second].values[i].delta;
6148 void Parameters::traverse_improper_params(
struct improper_params *list)
6152 #ifndef ACCELERATED_INPUT_IMPROPERS 6153 while (list != NULL)
6157 <<
" " << list->
atom4name <<
" index=" \
6164 <<
" n=" << list->
values[i].
n \
6171 for(std::pair<TupleString4, size_t> apair : impropermap->
tupleMap)
6173 DebugM(3,
"Improper " << apair.first.getTuplePtr(0) <<
" " \
6174 << apair.first.getTuplePtr(1)
6175 <<
" " << apair.first.getTuplePtr(2) \
6176 <<
" " << apair.first.getTuplePtr(3)
6177 <<
" index="<< apair.second \
6178 <<
" multiplicity=" << impropermap->
paramVector[apair.second].multiplicity);
6179 for (i=0; i<impropermap->
paramVector[apair.second].multiplicity; i++)
6182 <<
" n=" << impropermap->
paramVector[apair.second].values[i].n \
6183 <<
" delta=" << impropermap->
paramVector[apair.second].values[i].delta;
6206 void Parameters::traverse_vdw_params(
struct vdw_params *tree)
6209 #ifndef ACCELERATED_INPUT_VDW 6213 if (tree->
left != NULL)
6215 traverse_vdw_params(tree->
left);
6219 <<
" sigma=" << tree->
sigma <<
" epsilon=" << \
6221 <<
" epsilon 1-4=" << tree->
epsilon14 << std::endl);
6223 if (tree->
right != NULL)
6225 traverse_vdw_params(tree->
right);
6228 for(std::pair<TupleString1, size_t> apair : vdwmap->tupleMap)
6231 DebugM(3,
"vdW " << apair.first.getTuplePtr(0) <<
" index=" << apair.second \
6232 <<
" sigma=" << vdwmap->paramVector[apair.second].sigma <<
" epsilon=" << \
6233 vdwmap->paramVector[apair.second].epsilon <<
" sigma 1-4=" << vdwmap->paramVector[apair.second].sigma14 \
6234 <<
" epsilon 1-4=" << vdwmap->paramVector[apair.second].epsilon14 << std::endl);
6253 void Parameters::traverse_vdw_pair_params(
struct vdw_pair_params *list)
6256 #ifndef ACCELERATED_INPUT_VDW 6262 <<
" B=" << list->
B <<
" A 1-4=" \
6263 << list->
A14 <<
" B 1-4=" << list->
B14 <<
"\n" 6266 traverse_vdw_pair_params(list->
next);
6268 for(
auto apair : vdwpairmap)
6270 DebugM(3,
"vdW PAIR " << apair.second.ind1 <<
" " 6271 << apair.second.ind2
6272 <<
" A=" << apair.second.A
6273 <<
" B=" << apair.second.B <<
" A 1-4=" 6274 << apair.second.A14 <<
" B 1-4=" << apair.second.B14<<
"\n");
6296 #ifndef ACCELERATED_INPUT_NBTHOLE 6304 traverse_nbthole_pair_params(list->
next);
6306 for(
auto apair : nbtholepairmap)
6309 << apair.second.ind1 <<
" " 6310 << apair.second.ind2
6311 <<
" tholeij =" << apair.second.tholeij <<
"\n");
6329 <<
"*****************************************\n" \
6332 traverse_bond_params(bondp);
6347 <<
"*****************************************\n" );
6348 traverse_angle_params(anglep);
6363 <<
"*****************************************\n" );
6365 traverse_dihedral_params(dihedralp);
6380 <<
"*****************************************\n" );
6382 traverse_improper_params(improperp);
6397 <<
"*****************************************" );
6399 traverse_vdw_params(vdwp);
6414 <<
"*****************************************" );
6416 traverse_vdw_pair_params(vdw_pairp);
6431 <<
"*****************************************" );
6433 traverse_nbthole_pair_params(nbthole_pairp);
6455 iout <<
"CROSSTERM "<< i <<
" ";
6456 for (
int j = 0; j < N; ++j) {
6457 for (
int k = 0; k < N; ++k) {
6480 iout <<
iINFO <<
"SUMMARY OF PARAMETERS:\n" 6515 free_bond_tree(bondp);
6518 free_angle_tree(anglep);
6520 if (dihedralp != NULL)
6521 free_dihedral_list(dihedralp);
6523 if (improperp != NULL)
6524 free_improper_list(improperp);
6526 if (crosstermp != NULL)
6527 free_crossterm_list(crosstermp);
6530 free_vdw_tree(vdwp);
6534 if (maxDihedralMults != NULL)
6535 delete [] maxDihedralMults;
6537 if (maxImproperMults != NULL)
6538 delete [] maxImproperMults;
6546 maxImproperMults=NULL;
6547 maxDihedralMults=NULL;
6562 Real *a1, *a2, *a3, *a4;
6583 if ( (a1 == NULL) || (a2 == NULL) )
6585 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6611 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) ||
6612 (a4 == NULL) || (i1 == NULL))
6614 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6647 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
6648 (deltavals == NULL) )
6650 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6659 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
6661 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6687 delete [] deltavals[i];
6693 delete [] deltavals;
6706 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
6707 (deltavals == NULL) )
6709 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6718 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
6720 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6746 delete [] deltavals[i];
6752 delete [] deltavals;
6773 if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
6774 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6817 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) )
6819 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6854 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
6855 (i1 == NULL) || (i2 == NULL) )
6857 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6860 vdw_pair_to_arrays(i1, i2, a1, a2, a3, a4, 0,
vdw_pair_tree);
6879 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (i1 == NULL) || (i2 == NULL) )
6881 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6910 if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
6912 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
6940 Real *a1, *a2, *a3, *a4;
6958 if ( (
bond_array == NULL) || (a1 == NULL) || (a2 == NULL) )
6960 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
6988 if ( (
angle_array == NULL) || (a1 == NULL) || (a2 == NULL) ||
6989 (a3 == NULL) || (a4 == NULL) || (i1 == NULL))
6991 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
7028 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
7031 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
7040 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
7042 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
7071 delete [] deltavals[i];
7077 delete [] deltavals;
7091 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
7094 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
7103 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
7105 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
7134 delete [] deltavals[i];
7140 delete [] deltavals;
7166 if ( (pairC6 == NULL) || (pairC12 == NULL) ) {
7167 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
7192 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
7211 if ( (
vdw_array==NULL) || (a1==NULL) || (a2==NULL) || (a3==NULL)
7212 || (a4==NULL) || (atomTypeNames==NULL))
7214 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
7249 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
7250 (i1 == NULL) || (i2 == NULL) )
7252 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
7266 if (new_node == NULL)
7268 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
7271 new_node->
ind1 = i1[i];
7272 new_node->
ind2 = i2[i];
7273 new_node->
A = a1[i];
7274 new_node->
A14 = a2[i];
7275 new_node->
B = a3[i];
7276 new_node->
B14 = a4[i];
7277 new_node->
left = NULL;
7278 new_node->
right = NULL;
7303 if ( (
nbthole_array == NULL) || (a1 == NULL) || (a2 == NULL) || (a3 == NULL)
7304 || (i1 == NULL) || (i2 == NULL) )
7306 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
7342 if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
7344 NAMD_die(
"memory allocation failed in Parameters::send_Parameters");
7355 if (new_tab_node == NULL)
7357 NAMD_die(
"memory allocation failed in Parameters::receive_Parameters");
7361 new_tab_node->
ind1 = i1[i];
7362 new_tab_node->
ind2 = i2[i];
7363 new_tab_node->
K = i3[i];
7364 #ifndef ACCELERATED_INPUT_TABLE_PAIR 7365 new_tab_node->left = NULL;
7366 new_tab_node->right = NULL;
7379 AllFilesRead =
TRUE;
7396 void Parameters::convert_vdw_pairs()
7399 #ifndef ACCELERATED_INPUT_VDW 7400 #ifdef MEM_OPT_VERSION 7405 Index index1, index2;
7417 if (new_node == NULL)
7419 NAMD_die(
"memory allocation failed in Parameters::convert_vdw_pairs");
7429 if (index1 > index2)
7431 new_node->
ind1 = index2;
7432 new_node->
ind2 = index1;
7436 new_node->
ind1 = index1;
7437 new_node->
ind2 = index2;
7440 new_node->
A = ptr->
A;
7441 new_node->
B = ptr->
B;
7442 new_node->
A14 = ptr->
A14;
7443 new_node->
B14 = ptr->
B14;
7445 new_node->
left = NULL;
7446 new_node->
right = NULL;
7477 void Parameters::convert_nbthole_pairs()
7480 #ifdef MEM_OPT_VERSION 7485 Index index1, index2;
7489 ptr = nbthole_pairp;
7497 if (new_node == NULL)
7499 NAMD_die(
"memory allocation failed in Parameters::convert_nbthole_pairs");
7509 if (index1 > index2)
7511 new_node->
ind1 = index2;
7512 new_node->
ind2 = index1;
7516 new_node->
ind1 = index1;
7517 new_node->
ind2 = index2;
7524 new_node->
left = NULL;
7525 new_node->
right = NULL;
7538 nbthole_pairp = NULL;
7554 void Parameters::convert_table_pairs()
7557 #ifndef ACCELERATED_INPUT_TABLE_PAIR 7558 #ifdef MEM_OPT_VERSION 7563 Index index1, index2;
7575 if (new_node == NULL)
7577 NAMD_die(
"memory allocation failed in Parameters::convert_table_pairs");
7589 if (index1 > index2)
7591 new_node->
ind1 = index2;
7592 new_node->
ind2 = index1;
7596 new_node->
ind1 = index1;
7597 new_node->
ind2 = index2;
7600 new_node->
K = ptr->
K;
7602 new_node->left = NULL;
7603 new_node->right = NULL;
7640 #ifndef ACCELERATED_INPUT_TABLE_PAIR 7644 if ( (new_node->
ind1 < tree->
ind1) ||
7647 tree->left = add_to_indexed_table_pairs(new_node, tree->left);
7651 tree->right = add_to_indexed_table_pairs(new_node, tree->right);
7657 uint64_t k1=new_node->
ind1;
7658 uint64_t k2=new_node->
ind2;
7659 uint64_t key= (k1 < k2) ? ((k1 <<32) | k2) : ((k2 <<32) | k1);
7661 auto retval= tablepairmap.emplace(std::make_pair(key, value));
7662 bool duplicate = !retval.second;
7670 if (dupvalue->
K != new_node->
K)
7672 iout <<
iWARN <<
"DUPLICATE table PAIR ENTRY FOR " 7673 << new_node->
ind1 <<
"-" 7675 <<
"\nPREVIOUS VALUES K=" << dupvalue->
K 7676 <<
"\n USING VALUES K=" << new_node->
K 7679 dupvalue->
K = new_node->
K;
7704 #ifndef ACCELERATED_INPUT_VDW 7708 if ( (new_node->
ind1 < tree->
ind1) ||
7711 tree->
left = add_to_indexed_vdw_pairs(new_node, tree->
left);
7715 tree->
right = add_to_indexed_vdw_pairs(new_node, tree->
right);
7721 uint64_t k1=new_node->
ind1;
7722 uint64_t k2=new_node->
ind2;
7723 uint64_t key= (k1 < k2) ? ((k1 <<32) | k2) : ((k2 <<32) | k1);
7725 auto retval= vdwpairmap.emplace(std::make_pair(key, value));
7726 bool duplicate = !retval.second;
7734 if ((dupvalue->
A != new_node->
A) ||
7735 (dupvalue->
B != new_node->
B) ||
7736 (dupvalue->
A14 != new_node->
A14) ||
7737 (dupvalue->
B14 != new_node->
B14))
7739 iout <<
iWARN <<
"DUPLICATE vdW PAIR ENTRY FOR " 7740 << new_node->
ind1 <<
"-" 7742 <<
"\nPREVIOUS VALUES A=" << dupvalue->
A 7743 <<
" B=" << dupvalue->
B 7744 <<
" A14=" << dupvalue->
A14 7745 <<
" B14" << dupvalue->
B14 7746 <<
"\n USING VALUES A=" << new_node->
A 7747 <<
" B=" << new_node->
B 7748 <<
" A14=" << new_node->
A14 7749 <<
" B14" << new_node->
B14 7752 dupvalue->
A = new_node->
A;
7753 dupvalue->
B = new_node->
B;
7754 dupvalue->
A14 = new_node->
A14;
7755 dupvalue->
B14 = new_node->
B14;
7783 if ( (new_node->
ind1 < tree->
ind1) ||
7786 tree->
left = add_to_indexed_nbthole_pairs(new_node, tree->
left);
7790 tree->
right = add_to_indexed_nbthole_pairs(new_node, tree->
right);
7818 int Parameters::vdw_pair_to_arrays(
int *ind1_array,
int *ind2_array,
7824 #ifndef ACCELERATED_INPUT_VDW 7828 ind1_array[arr_index] = tree->
ind1;
7829 ind2_array[arr_index] = tree->
ind2;
7830 A[arr_index] = tree->
A;
7831 A14[arr_index] = tree->
A14;
7832 B[arr_index] = tree->
B;
7833 B14[arr_index] = tree->
B14;
7837 arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
7838 arr_index, tree->
left);
7839 arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
7840 arr_index, tree->
right);
7845 for(
auto vdwVal : vdwpairmap)
7847 ind1_array[index] = vdwVal.second.ind1;
7848 ind2_array[index] = vdwVal.second.ind2;
7849 A[index] = vdwVal.second.A;
7850 A14[index] = vdwVal.second.A14;
7851 B[index] = vdwVal.second.B;
7852 B14[index] = vdwVal.second.B14;
7878 int Parameters::nbthole_pair_to_arrays(
int *ind1_array,
int *ind2_array,
7883 #ifndef ACCELERATED_INPUT_NBTHOLE 7887 ind1_array[arr_index] = tree->
ind1;
7888 ind2_array[arr_index] = tree->
ind2;
7889 alphai[arr_index] = tree->
alphai;
7890 alphaj[arr_index] = tree->
alphaj;
7891 tholeij[arr_index] = tree->
tholeij;
7895 arr_index = nbthole_pair_to_arrays(ind1_array, ind2_array, alphai,
7896 alphaj, tholeij, arr_index, tree->
left);
7897 arr_index = nbthole_pair_to_arrays(ind1_array, ind2_array, alphai,
7898 alphaj, tholeij, arr_index, tree->
right);
7903 for(
auto tableVal : nbtholepairmap)
7905 ind1_array[index] = tableVal.second.ind1;
7906 ind2_array[index] = tableVal.second.ind2;
7907 tholeij[index] = tableVal.second.tholeij;
7908 alphai[index] = tableVal.second.alphai;
7909 alphaj[index] = tableVal.second.alphaj;
7934 int Parameters::table_pair_to_arrays(
int *ind1_array,
int *ind2_array,
7940 #ifndef ACCELERATED_INPUT_TABLE_PAIR 7944 ind1_array[arr_index] = tree->
ind1;
7945 ind2_array[arr_index] = tree->
ind2;
7946 K[arr_index] = tree->
K;
7950 arr_index = table_pair_to_arrays(ind1_array, ind2_array,
K,
7951 arr_index, tree->left);
7952 arr_index = table_pair_to_arrays(ind1_array, ind2_array,
K,
7953 arr_index, tree->right);
7957 for(
auto tableVal : tablepairmap)
7959 ind1_array[index] = tableVal.second.ind1;
7960 ind2_array[index] = tableVal.second.ind2;
7961 K[index] = tableVal.second.K;
7982 read_parm(amber_data,vdw14);
8001 int i,j,ico,numtype,mul;
8005 NAMD_die(
"No data read from parm file yet!");
8012 NAMD_die(
"memory allocation of bond_array failed!");
8024 NAMD_die(
"memory allocation of angle_array failed!");
8046 NAMD_die(
"memory allocation of dihedral_array failed!");
8054 snprintf(err_msg,
sizeof(err_msg),
8055 "The periodicity of dihedral # %d is zero!", i+1);
8061 if (amber_data->
Pn[i] > 0)
8067 snprintf(err_msg,
sizeof(err_msg),
8068 "Multiple dihedral with multiplicity of %d greater than max of %d",
8074 NAMD_die(
"Negative periodicity in the last dihedral entry!");
8081 for (i=0; i<numtype; ++i)
8082 for (j=i; j<numtype; ++j)
8084 if (new_node == NULL)
8085 NAMD_die(
"memory allocation of vdw_pair_tree failed!");
8088 new_node->
left = new_node->
right = NULL;
8094 ico = amber_data->
Cno[numtype*i+j];
8096 { new_node->
A = amber_data->
Cn1[ico-1];
8097 new_node->
A14 = new_node->
A / vdw14;
8098 new_node->
B = amber_data->
Cn2[ico-1];
8099 new_node->
B14 = new_node->
B / vdw14;
8101 else if (amber_data->
HB12[abs(ico)-1]==0.0 && amber_data->
HB6[abs(ico)-1]==0.0)
8102 { new_node->
A = new_node->
A14 = new_node->
B = new_node->
B14 = 0.0;
8103 iout <<
iWARN <<
"Encounter 10-12 H-bond term\n";
8106 NAMD_die(
"Encounter non-zero 10-12 H-bond term!");
8119 read_parm(amber_data,vdw14);
8124 int i,j,ico,numtype,mul;
8128 NAMD_die(
"No data read from parm file yet!");
8133 NAMD_die(
"You are using a CHARMM force field in AMBER format generated by CHAMBER or ParmEd.\n" 8134 "We do not support this format. Please consider using the native format (PSF) for CHARMM force field.");
8142 NAMD_die(
"memory allocation of bond_array failed!");
8154 NAMD_die(
"memory allocation of angle_array failed!");
8176 NAMD_die(
"memory allocation of dihedral_array failed!");
8184 snprintf(err_msg,
sizeof(err_msg),
8185 "The periodicity of dihedral # %d is zero!", i+1);
8191 if (amber_data->
Pn[i] > 0)
8197 snprintf(err_msg,
sizeof(err_msg),
8198 "Multiple dihedral with multiplicity of %d greater than max of %d",
8204 NAMD_die(
"Negative periodicity in the last dihedral entry!");
8217 snprintf(err_msg,
sizeof(err_msg),
8218 "The table of %d crossterm type has a resolution(%d) " 8225 for (
int j = 0; j < N; ++j) {
8226 for (
int k = 0; k < N; ++k) {
8231 for (
int j = 0; j < N; ++j) {
8246 for (i=0; i<numtype; ++i)
8247 for (j=i; j<numtype; ++j)
8249 if (new_node == NULL)
8250 NAMD_die(
"memory allocation of vdw_pair_tree failed!");
8253 new_node->
left = new_node->
right = NULL;
8259 ico = amber_data->
Cno[numtype*i+j];
8261 { new_node->
A = amber_data->
Cn1[ico-1];
8262 new_node->
A14 = new_node->
A / vdw14;
8263 new_node->
B = amber_data->
Cn2[ico-1];
8264 new_node->
B14 = new_node->
B / vdw14;
8266 else if (amber_data->
HB12[abs(ico)-1]==0.0 && amber_data->
HB10[abs(ico)-1]==0.0)
8267 { new_node->
A = new_node->
A14 = new_node->
B = new_node->
B14 = 0.0;
8268 iout <<
iWARN <<
"Encounter 10-12 H-bond term\n";
8271 NAMD_die(
"Encounter non-zero 10-12 H-bond term!");
8321 NAMD_die(
"memory allocation of bond_array failed!");
8326 NAMD_die(
"memory allocation of dihedral_array failed!");
8331 NAMD_die(
"memory allocation of angle_array failed!");
8389 NAMD_die(
"I can't do RB dihedrals with MAX_MULTIPLICITY < 5");
8395 if(j%2 == 1) c[j] = -c[j];
8420 test2 = c[0]+1/2.*c[2]+3/8.*c[4];
8425 if(fabs(test1-test2) > 0.0001)
8426 NAMD_die(
"Internal error: failed to handle RB dihedrals");
8443 NAMD_die(
"unknown dihedral type found");
8453 for (i=0; i<numtype; i++) {
8454 for (j=i; j<numtype; j++) {
8458 if (new_node == NULL)
8459 NAMD_die(
"memory allocation of vdw_pair_tree failed!");
8462 new_node->
left = new_node->
right = NULL;
8465 &(new_node->
B14), &(new_node->
A14));
8471 if(min && ( fabs(new_node->
A) < 1.0 )) {
8475 "Adding small LJ repulsion term to some atoms.\n" <<
endi;
8487 Real *pairC6,*pairC12;
8504 const_cast<GromacsTopFile*
>(gf)->getPairLJArrays2(atom1, atom2, pairC6, pairC12);
8528 char* table_file =
simParams->tabulatedEnergiesFile;
8529 char* interp_type =
simParams->tableInterpType;
8531 enertable = fopen(table_file,
"r");
8533 if (enertable == NULL) {
8534 NAMD_die(
"ERROR: Could not open energy table file!\n");
8537 char tableline[121];
8551 iout <<
"Beginning table read\n" <<
endi;
8552 while(fgets(tableline,120,enertable)) {
8553 if (strncmp(tableline,
"#",1)==0) {
8556 readcount = sscanf(tableline,
"%i %f %f", &numtypes, &table_spacing, &maxdist);
8557 if (readcount != 3) {
8558 NAMD_die(
"ERROR: Couldn't parse table header line\n");
8563 if (maxdist < simParams->cutoff) {
8564 NAMD_die(
"Tabulated energies must at least extend to the cutoff distance\n");
8568 iout <<
"Info: Reducing max table size to match nonbond cutoff\n";
8584 int numtypesread = 0;
8586 while(fgets(tableline,120,enertable)) {
8587 if (strncmp(tableline,
"#",1)==0) {
8590 if (strncmp(tableline,
"TYPE",4)==0) {
8592 newtypename =
new char[5];
8593 int readcount = sscanf(tableline,
"%*s %s", newtypename);
8594 if (readcount != 1) {
8595 NAMD_die(
"Couldn't read interaction type from TYPE line\n");
8600 if (numtypesread == numtypes) {
8601 NAMD_die(
"Error: Number of types in table doesn't match table header\n");
8605 if (!strncasecmp(interp_type,
"linear", 6)) {
8608 snprintf(err_msg,
sizeof(err_msg),
8609 "Failed to read table energy (linear) type %s\n", newtypename);
8612 }
else if (!strncasecmp(interp_type,
"cubic", 5)) {
8615 snprintf(err_msg,
sizeof(err_msg),
8616 "Failed to read table energy (cubic) type %s\n", newtypename);
8620 NAMD_die(
"Table interpolation type must be linear or cubic\n");
8633 if (numtypesread != numtypes) {
8635 snprintf(err_msg,
sizeof(err_msg),
8636 "ERROR: Expected %i tabulated energy types but got %i\n",
8637 numtypes, numtypesread);
8643 simParams->tableSpacing = table_spacing;
8777 char tableline[120];
8794 std::vector<BigReal> dists;
8795 std::vector<BigReal> enervalues;
8796 std::vector<BigReal> forcevalues;
8808 while(fgets(tableline,120,enertable) && distbin <= (
int) (
namdnearbyint(maxdist / table_spacing) + 1)) {
8809 if (strncmp(tableline,
"#",1)==0) {
8812 if (strncmp(tableline,
"TYPE",4)==0) {
8813 fseek(enertable, -1 * (
long) strlen(tableline), SEEK_CUR);
8818 int readcount = sscanf(tableline,
"%lf %lf %lf", &readdist, &readener, &readforce);
8821 if (readcount != 3) {
8823 snprintf(err_msg,
sizeof(err_msg),
8824 "ERROR: Failed to parse table line %s!\n", tableline);
8829 if (readdist < lastdist) {
8830 NAMD_die(
"ERROR: Encountered badly ordered entries in energy table!\n");
8833 lastdist = readdist;
8834 dists.push_back(readdist);
8835 enervalues.push_back(readener);
8836 forcevalues.push_back(readforce);
8841 if (dists[0] != 0.0) {
8842 NAMD_die(
"ERROR: First data point for energy table must be at r=0\n");
8844 spacing = dists[1] - dists[0];
8845 for (i=2; i<(numentries - 1); i++) {
8847 myspacing = dists[i] - dists[i-1];
8848 if (fabs(myspacing - spacing) > 0.00001) {
8849 printf(
"Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
8850 NAMD_die(
"ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
8858 double* m =
new double[numentries*numentries];
8859 double* xe =
new double[numentries];
8860 double* xf =
new double[numentries];
8861 double* be =
new double[numentries];
8862 double* bf =
new double[numentries];
8864 be[0] = 3 * (enervalues[1] - enervalues[0]);
8865 for (i=1; i < (numentries - 1); i++) {
8867 be[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
8870 be[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
8872 bf[0] = 3 * (forcevalues[1] - forcevalues[0]);
8873 for (i=1; i < (numentries - 1); i++) {
8875 bf[i] = 3 * (forcevalues[i+1] - forcevalues[i-1]);
8878 bf[numentries - 1] = 3 * (forcevalues[numentries - 1] - forcevalues[numentries - 2]);
8880 memset(m,0,numentries*numentries*
sizeof(
double));
8884 for (i = 1; i < numentries; i++) {
8885 m[
INDEX(numentries,i-1,i)] = 1;
8886 m[
INDEX(numentries,i,i-1)] = 1;
8887 m[
INDEX(numentries,i,i)] = 4;
8889 m[
INDEX(numentries,numentries-1,numentries-1)] = 2;
8895 for (i=0; i<numentries; i++) {
8898 for (j=i+1; j<numentries; j++) {
8900 const BigReal factor = myval / baseval;
8902 for (
int k=i; k<numentries; k++) {
8904 m[
INDEX(numentries,j,k)] -= (factor * subval);
8907 be[j] -= (factor * be[i]);
8908 bf[j] -= (factor * bf[i]);
8914 for (i=numentries-1; i>=0; i--) {
8917 for (j=i+1; j<numentries; j++) {
8918 be[i] -= ( m[
INDEX(numentries,i,j)] * xe[j] );
8919 bf[i] -= ( m[
INDEX(numentries,i,j)] * xf[j] );
8922 xe[i] = be[i] / m[
INDEX(numentries,i,i)];
8923 xf[i] = bf[i] / m[
INDEX(numentries,i,i)];
8931 int entriesperseg = (int) ceil(spacing / table_spacing);
8932 int table_prefix = 2 * typeindex * (int) (
namdnearbyint(maxdist / table_spacing) + 1);
8934 for (i=0; i<numentries-1; i++) {
8937 currdist = dists[i];
8944 Ce = 3 * (enervalues[i+1] - enervalues[i]) - (2 * xe[i]) - xe[i+1];
8945 De = 2 * (enervalues[i] - enervalues[i+1]) + xe[i] + xe[i+1];
8947 Af = forcevalues[i];
8949 Cf = 3 * (forcevalues[i+1] - forcevalues[i]) - (2 * xf[i]) - xf[i+1];
8950 Df = 2 * (forcevalues[i] - forcevalues[i+1]) + xf[i] + xf[i+1];
8953 for (j=0; j<entriesperseg; j++) {
8954 const BigReal mydist = currdist + (j * table_spacing);
8955 const BigReal mydistfrac = (float) j / (entriesperseg - 1);
8957 if (distbin >= (
int)
namdnearbyint(maxdist / table_spacing))
break;
8961 energy = Ae + (Be * mydistfrac) + (Ce * mydistfrac * mydistfrac) + (De * mydistfrac * mydistfrac * mydistfrac);
8962 force = Af + (Bf * mydistfrac) + (Cf * mydistfrac * mydistfrac) + (Df * mydistfrac * mydistfrac * mydistfrac);
8965 table_ener[table_prefix + 2 * distbin] = energy;
8966 table_ener[table_prefix + 2 * distbin + 1] = force;
8969 if (currdist >= maxdist)
break;
8975 table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
8976 table_ener[table_prefix + 2 * distbin + 1] = 0.0;
8987 printf(
"Testing: %i vs %i (from %f / %f)\n", distbin, (
int) (
namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
8988 if (distbin != (
int) (
namdnearbyint(maxdist / table_spacing)))
return 1;
9014 char tableline[120];
9030 std::vector<BigReal> dists;
9031 std::vector<BigReal> enervalues;
9043 while(fgets(tableline,120,enertable) && distbin <= (
int) (
namdnearbyint(maxdist / table_spacing) + 1)) {
9044 if (strncmp(tableline,
"#",1)==0) {
9047 if (strncmp(tableline,
"TYPE",4)==0) {
9048 fseek(enertable, -1 * (
long) strlen(tableline), SEEK_CUR);
9053 int readcount = sscanf(tableline,
"%lf %lf", &readdist, &readener);
9056 if (readcount != 2) {
9058 snprintf(err_msg,
sizeof(err_msg),
9059 "ERROR: Failed to parse table line %s!\n", tableline);
9064 if (readdist < lastdist) {
9065 NAMD_die(
"ERROR: Encountered badly ordered entries in energy table!\n");
9068 lastdist = readdist;
9069 dists.push_back(readdist);
9070 enervalues.push_back(readener);
9075 if (dists[0] != 0.0) {
9076 NAMD_die(
"ERROR: First data point for energy table must be at r=0\n");
9078 spacing = dists[1] - dists[0];
9079 for (i=2; i<(numentries - 1); i++) {
9081 myspacing = dists[i] - dists[i-1];
9082 if (fabs(myspacing - spacing) > 0.00001) {
9083 printf(
"Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
9084 NAMD_die(
"ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
9091 double* m =
new double[numentries*numentries];
9092 double* x =
new double[numentries];
9093 double* b =
new double[numentries];
9095 b[0] = 3 * (enervalues[1] - enervalues[0]);
9096 for (i=1; i < (numentries - 1); i++) {
9097 printf(
"Control point %i at %f\n", i, enervalues[i]);
9098 b[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
9099 printf(
"b is %f\n", b[i]);
9101 b[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
9103 memset(m,0,numentries*numentries*
sizeof(
double));
9107 for (i = 1; i < numentries; i++) {
9108 m[
INDEX(numentries,i-1,i)] = 1;
9109 m[
INDEX(numentries,i,i-1)] = 1;
9110 m[
INDEX(numentries,i,i)] = 4;
9112 m[
INDEX(numentries,numentries-1,numentries-1)] = 2;
9116 printf(
"Solving the matrix equation: \n");
9118 for (i=0; i<numentries; i++) {
9120 for (j=0; j<numentries; j++) {
9121 printf(
" %6.3f,", m[
INDEX(numentries, i, j)]);
9123 printf(
" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
9127 for (i=0; i<numentries; i++) {
9130 for (j=i+1; j<numentries; j++) {
9132 const BigReal factor = myval / baseval;
9134 for (
int k=i; k<numentries; k++) {
9136 m[
INDEX(numentries,j,k)] -= (factor * subval);
9139 b[j] -= (factor * b[i]);
9144 printf(
" In upper diagonal form, equation is:\n");
9145 for (i=0; i<numentries; i++) {
9147 for (j=0; j<numentries; j++) {
9148 printf(
" %6.3f,", m[
INDEX(numentries, i, j)]);
9150 printf(
" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
9154 for (i=numentries-1; i>=0; i--) {
9157 for (j=i+1; j<numentries; j++) {
9158 b[i] -= ( m[
INDEX(numentries,i,j)] * x[j] );
9161 x[i] = b[i] / m[
INDEX(numentries,i,i)];
9165 printf(
" Solution vector is:\n\t(");
9166 for (i=0; i<numentries; i++) {
9167 printf(
" %6.3f ", x[i]);
9175 int entriesperseg = (int) ceil(spacing / table_spacing);
9176 int table_prefix = 2 * typeindex * (int) (
namdnearbyint(maxdist / table_spacing) + 1);
9178 for (i=0; i<numentries-1; i++) {
9180 currdist = dists[i];
9182 printf(
"Interpolating on interval %i\n", i);
9187 C = 3 * (enervalues[i+1] - enervalues[i]) - (2 * x[i]) - x[i+1];
9188 D = 2 * (enervalues[i] - enervalues[i+1]) + x[i] + x[i+1];
9190 printf(
"Coefficients for this interval: %f %f %f %f\n", A, B, C, D);
9193 for (j=0; j<entriesperseg; j++) {
9194 const BigReal mydist = currdist + (j * table_spacing);
9195 const BigReal mydistfrac = (float) j / (entriesperseg - 1);
9197 if (distbin >= (
int)
namdnearbyint(maxdist / table_spacing))
break;
9201 energy = A + (B * mydistfrac) + (C * mydistfrac * mydistfrac) + (D * mydistfrac * mydistfrac * mydistfrac);
9202 force = B + (2 * C * mydistfrac) + (3 * D * mydistfrac * mydistfrac);
9204 force *= (1.0 / spacing);
9206 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);
9207 table_ener[table_prefix + 2 * distbin] = energy;
9208 table_ener[table_prefix + 2 * distbin + 1] = force;
9211 if (currdist >= maxdist)
break;
9216 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));
9217 table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
9218 table_ener[table_prefix + 2 * distbin + 1] = 0.0;
9227 printf(
"Testing: %i vs %i (from %f / %f)\n", distbin, (
int) (
namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
9228 if (distbin != (
int) (
namdnearbyint(maxdist / table_spacing)))
return 1;
9254 char tableline[120];
9273 while(fgets(tableline,120,enertable) && distbin <= (
int) (
namdnearbyint(maxdist / table_spacing) + 1)) {
9274 printf(
"At distance %f + %f vs. %f\n", currdist, table_spacing, maxdist);
9275 if (strncmp(tableline,
"#",1)==0) {
9278 if (strncmp(tableline,
"TYPE",4)==0) {
9283 lastdist = readdist;
9284 lastener = readener;
9285 lastforce = readforce;
9286 int readcount = sscanf(tableline,
"%lf %lf %lf", &readdist, &readener, &readforce);
9287 if (distbin == -1) {
9288 if (readdist != 0.0) {
9289 NAMD_die(
"ERROR: Energy/force tables must start at d=0.0\n");
9296 if (readcount != 3) {
9298 snprintf(err_msg,
sizeof(err_msg),
9299 "ERROR: Failed to parse table line %s!\n", tableline);
9304 if (readdist < lastdist) {
9305 NAMD_die(
"ERROR: Encountered badly ordered entries in energy table!\n");
9308 currdist = lastdist;
9310 while (currdist <= readdist && distbin <= (
int) (
namdnearbyint(maxdist / table_spacing))) {
9312 int table_loc = 2 * (distbin + (typeindex * (1 + (int)
namdnearbyint(maxdist / table_spacing))));
9313 printf(
"Doing interpolation for energy between %f %f and %f %f: Dist %f\n", readener, readdist, lastener, lastdist, currdist);
9314 table_ener[table_loc] = interp_lin(readener, lastener, readdist, lastdist, currdist);
9315 table_ener[table_loc + 1] = interp_lin(readforce, lastforce, readdist, lastdist, currdist);
9316 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);
9317 currdist += table_spacing;
9323 fseek(enertable, -1 * (
long) strlen(tableline), SEEK_CUR);
9327 printf(
"Testing: %i vs %i (from %f / %f)\n", distbin, (
int) (
namdnearbyint(maxdist / table_spacing)), maxdist, table_spacing);
9328 if (distbin != (
int) (
namdnearbyint(maxdist / table_spacing)))
return 1;
9352 m = (val2 - val1) / (end2 - end1);
9354 val = ((currdist-end1) * m + val1);
struct improper_params * next
int getNumAngleParams() const
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 print_crossterm_params()
void read_charmm_parameter_file(char *)
static char ** table_types
void getDihedralParams(int num, Real *c, int *mult, int *funct) const
struct indexed_vdw_pair * right
int getNumBondParams() const
TupleString< 1 > TupleString1
std::ostream & endi(std::ostream &s)
TupleString< 2 > TupleString2
char * getTuplePtr(short index)
void crossterm_setup(CrosstermData *)
struct nbthole_pair_params * next
TupleString< 8 > TupleString8
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)
int add(const Elem &elem)
void assign_bond_index(const char *, const char *, Bond *, bool *bond_found=nullptr)
CrosstermValue * crossterm_array
IndexedVdwPair * vdw_pair_tree
int read_energy_type(FILE *, const int, BigReal *, const float, const float)
std::unordered_map< TupleString< NumStrings >, size_t, TupleStringHash< NumStrings > > tupleMap
int getNumAtomParams() const
#define INDEX(ncols, i, j)
void print_param_summary()
IndexedNbtholePair * nbthole_pair_tree
void read_parameter_file(char *)
FourBodyConsts values[MAX_MULTIPLICITY]
void getAngleParams(int num, Real *th0, Real *kth, int *funct) const
vector< int > CMAPResolution
void getVDWParams(int typea, int typeb, Real *c6, Real *c12, Real *c6pair, Real *c7) const
static Units next(Units u)
struct indexed_vdw_pair * left
int get_int_table_type(char *)
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()
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)
TupleString< 3 > TupleString3
NbtholePairValue * nbthole_array
char * atom_type_name(Index a)
std::vector< ParamValue > paramVector
FourBodyConsts values[MAX_MULTIPLICITY]
struct indexed_nbthole_pair * right
int get_table_pair_params(Index, Index, int *)
void send_Parameters(MOStream *)
void done_reading_files(Bool)
void print_nbthole_pair_params()
vector< vector< _REAL > > CMAPParameter
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)
struct table_pair_params * next
void done_reading_structure()
void getBondParams(int num, Real *b0, Real *kB, int *funct) const
struct vdw_pair_params * next
void assign_angle_index(const char *, const char *, const char *, Angle *, int)
MOStream * put(char data)
int getNumDihedralParams() const
void NAMD_remove_comment(char *str)
struct indexed_nbthole_pair * left
void print_angle_params()
int64_t const index(const TupleString< NumStrings > &findKey) const
TupleString< 4 > TupleString4
IndexedTablePair * tab_pair_tree
FourBodyConsts values[MAX_MULTIPLICITY]
struct vdw_params * right
#define MAX_ATOMTYPE_CHARS
void receive_Parameters(MIStream *)
std::pair< bool, ParamValue * > insert_check(const TupleString< NumStrings > &tKey, const ParamValue &mValue)