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);
  9382 #ifndef ACCELERATED_INPUT_VDW  9385   std::remove_reference<decltype((*vdwmap))>::type::KeyType k;
  9386   const bool key_found = vdwmap->get_key_by_index(a, k);
  9388     return std::string(k.getCatKey());
  9390     return std::string();
 
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
 
std::string atom_type_name_str(Index a) const
 
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
 
int64_t const index(const KeyType &findKey) const
 
IndexedVdwPair * vdw_pair_tree
 
int read_energy_type(FILE *, const int, BigReal *, const float, const float)
 
std::unordered_map< KeyType, size_t, TupleStringHash< NumStrings > > tupleMap
 
int getNumAtomParams() const
 
std::pair< bool, ParamValue * > insert_check(const KeyType &tKey, const ParamValue &mValue)
 
#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 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()
 
TupleString< 4 > TupleString4
 
IndexedTablePair * tab_pair_tree
 
FourBodyConsts values[MAX_MULTIPLICITY]
 
void done_reading_files(const bool addDrudeBond)
 
struct vdw_params * right
 
#define MAX_ATOMTYPE_CHARS
 
void receive_Parameters(MIStream *)