00001
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include <stdio.h>
00016 #include <stdlib.h>
00017 #include <string.h>
00018 #include <vector>
00019 #ifndef WIN32
00020 #include <strings.h>
00021 #endif
00022 #include "InfoStream.h"
00023 #include <charm++.h>
00024 #include "Parameters.h"
00025 #include "Communicate.h"
00026 #include "ConfigList.h"
00027
00028 #include "SimParameters.h"
00029
00030
00031 #define MIN_DEBUG_LEVEL 3
00032
00033 #include "Debug.h"
00034
00035 #define INDEX(ncols,i,j) ((i)*ncols + (j))
00036
00037 #define ENABLETABLES
00038
00039 static char** table_types;
00040
00041
00042
00043
00044
00045 struct bond_params
00046 {
00047 char atom1name[11];
00048 char atom2name[11];
00049 Real forceconstant;
00050 Real distance;
00051 Index index;
00052 struct bond_params *left;
00053 struct bond_params *right;
00054 };
00055
00056
00057
00058
00059
00060 struct angle_params
00061 {
00062 char atom1name[11];
00063 char atom2name[11];
00064 char atom3name[11];
00065 Real forceconstant;
00066 int normal;
00067 Real angle;
00068 Real k_ub;
00069 Real r_ub;
00070 Index index;
00071 struct angle_params *left;
00072 struct angle_params *right;
00073 };
00074
00075
00076
00077
00078
00079
00080
00081 struct dihedral_params
00082 {
00083 char atom1name[11];
00084 char atom2name[11];
00085 char atom3name[11];
00086 char atom4name[11];
00087 char atom1wild;
00088 char atom2wild;
00089 char atom3wild;
00090 char atom4wild;
00091 int multiplicity;
00092 FourBodyConsts values[MAX_MULTIPLICITY];
00093 Index index;
00094 struct dihedral_params *next;
00095 dihedral_params() { memset(this, 0, sizeof(dihedral_params)); }
00096 };
00097
00098
00099
00100
00101
00102
00103
00104 struct improper_params
00105 {
00106 char atom1name[11];
00107 char atom2name[11];
00108 char atom3name[11];
00109 char atom4name[11];
00110 int multiplicity;
00111 FourBodyConsts values[MAX_MULTIPLICITY];
00112 Index index;
00113 struct improper_params *next;
00114 };
00115
00116 struct crossterm_params
00117 {
00118 crossterm_params(int dim) : dimension(dim) {
00119 values = new double[dimension*dimension];
00120 }
00121 ~crossterm_params() {
00122 delete [] values;
00123 }
00124 char atom1name[11];
00125 char atom2name[11];
00126 char atom3name[11];
00127 char atom4name[11];
00128 char atom5name[11];
00129 char atom6name[11];
00130 char atom7name[11];
00131 char atom8name[11];
00132 int dimension;
00133 double *values;
00134 Index index;
00135 struct crossterm_params *next;
00136 };
00137
00138
00139
00140
00141 struct vdw_params
00142 {
00143 char atomname[11];
00144 Real sigma;
00145 Real epsilon;
00146 Real sigma14;
00147 Real epsilon14;
00148 Index index;
00149 struct vdw_params *left;
00150 struct vdw_params *right;
00151 };
00152
00153
00154
00155
00156 struct vdw_pair_params
00157 {
00158 char atom1name[11];
00159 char atom2name[11];
00160 Real A;
00161 Real B;
00162 Real A14;
00163 Real B14;
00164 struct vdw_pair_params *next;
00165 };
00166
00167 struct table_pair_params
00168 {
00169 char atom1name[11];
00170 char atom2name[11];
00171 int K;
00172 struct table_pair_params *next;
00173 };
00174
00175 struct nbthole_pair_params
00176 {
00177 char atom1name[11];
00178 char atom2name[11];
00179 Real alphai;
00180 Real alphaj;
00181 Real tholeij;
00182 Index index;
00183 struct nbthole_pair_params *next;
00184 };
00185
00186 Parameters::Parameters() {
00187 initialize();
00188 }
00189
00190 void Parameters::initialize() {
00191
00192 paramType = -1;
00193
00194
00195 atomTypeNames=NULL;
00196 bondp=NULL;
00197 anglep=NULL;
00198 improperp=NULL;
00199 dihedralp=NULL;
00200 crosstermp=NULL;
00201 vdwp=NULL;
00202 vdw_pairp=NULL;
00203 nbthole_pairp=NULL;
00204 table_pairp=NULL;
00205 bond_array=NULL;
00206 angle_array=NULL;
00207 dihedral_array=NULL;
00208 improper_array=NULL;
00209 crossterm_array=NULL;
00210 vdw_array=NULL;
00211 vdw_pair_tree=NULL;
00212 nbthole_pair_tree=NULL;
00213 tab_pair_tree=NULL;
00214 maxDihedralMults=NULL;
00215 maxImproperMults=NULL;
00216 table_ener = NULL;
00217
00218
00219 NumBondParams=0;
00220 NumAngleParams=0;
00221 NumDihedralParams=0;
00222 NumImproperParams=0;
00223 NumCrosstermParams=0;
00224 NumVdwParams=0;
00225 NumVdwPairParams=0;
00226 NumNbtholePairParams=0;
00227 NumTablePairParams=0;
00228 NumCosAngles=0;
00229 numenerentries=0;
00230 }
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 Parameters::Parameters(SimParameters *simParams, StringList *f)
00244 {
00245 initialize();
00246
00248 if (simParams->paraTypeXplorOn)
00249 {
00250 paramType = paraXplor;
00251 }
00252 else if (simParams->paraTypeCharmmOn)
00253 {
00254 paramType = paraCharmm;
00255 }
00256
00257
00258 if (simParams->cosAngles) {
00259 cosAngles = true;
00260 } else {
00261 cosAngles = false;
00262 }
00263
00264 if (simParams->tabulatedEnergies) {
00265 CkPrintf("Working on tables\n");
00266 read_ener_table(simParams);
00267 }
00268
00269
00270
00271
00272
00273 AllFilesRead = FALSE;
00274
00275 if (NULL != f)
00276 {
00277 do
00278 {
00279
00280 if (paramType == paraXplor)
00281 {
00282 read_parameter_file(f->data);
00283 }
00284 else if (paramType == paraCharmm)
00285 {
00286 read_charmm_parameter_file(f->data);
00287 }
00288
00289 f = f->next;
00290 } while ( f != NULL );
00291
00292 done_reading_files();
00293 }
00294
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 Parameters::~Parameters()
00308
00309 {
00310 if (atomTypeNames)
00311 delete [] atomTypeNames;
00312
00313 if (bondp != NULL)
00314 free_bond_tree(bondp);
00315
00316 if (anglep != NULL)
00317 free_angle_tree(anglep);
00318
00319 if (dihedralp != NULL)
00320 free_dihedral_list(dihedralp);
00321
00322 if (improperp != NULL)
00323 free_improper_list(improperp);
00324
00325 if (crosstermp != NULL)
00326 free_crossterm_list(crosstermp);
00327
00328 if (vdwp != NULL)
00329 free_vdw_tree(vdwp);
00330
00331 if (vdw_pairp != NULL)
00332 free_vdw_pair_list();
00333
00334 if (nbthole_pairp != NULL)
00335 free_nbthole_pair_list();
00336
00337 if (bond_array != NULL)
00338 delete [] bond_array;
00339
00340 if (angle_array != NULL)
00341 delete [] angle_array;
00342
00343 if (dihedral_array != NULL)
00344 delete [] dihedral_array;
00345
00346 if (improper_array != NULL)
00347 delete [] improper_array;
00348
00349 if (crossterm_array != NULL)
00350 delete [] crossterm_array;
00351
00352 if (vdw_array != NULL)
00353 delete [] vdw_array;
00354
00355 if (tab_pair_tree != NULL)
00356 free_table_pair_tree(tab_pair_tree);
00357
00358 if (vdw_pair_tree != NULL)
00359 free_vdw_pair_tree(vdw_pair_tree);
00360
00361 if (nbthole_pair_tree != NULL)
00362 free_nbthole_pair_tree(nbthole_pair_tree);
00363
00364 if (maxDihedralMults != NULL)
00365 delete [] maxDihedralMults;
00366
00367 if (maxImproperMults != NULL)
00368 delete [] maxImproperMults;
00369
00370 for( int i = 0; i < error_msgs.size(); ++i ) {
00371 delete [] error_msgs[i];
00372 }
00373 error_msgs.resize(0);
00374 }
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393 void Parameters::read_parameter_file(char *fname)
00394
00395 {
00396 char buffer[512];
00397 char first_word[512];
00398 FILE *pfile;
00399
00400
00401
00402 if (AllFilesRead)
00403 {
00404 NAMD_die("Tried to read another parameter file after being told that all files were read!");
00405 }
00406
00407
00408 if ( (pfile = Fopen(fname, "r")) == NULL)
00409 {
00410 char err_msg[256];
00411
00412 sprintf(err_msg, "UNABLE TO OPEN XPLOR PARAMETER FILE %s\n", fname);
00413 NAMD_die(err_msg);
00414 }
00415
00416
00417 while (NAMD_read_line(pfile, buffer) != -1)
00418 {
00419
00420 NAMD_find_first_word(buffer, first_word);
00421
00422
00423
00424
00425
00426
00427
00428 if ((buffer[0] != '!') &&
00429 !NAMD_blank_string(buffer) &&
00430 (strncasecmp(first_word, "REMARK", 6) != 0) &&
00431 (strcasecmp(first_word, "set")!=0) &&
00432 (strncasecmp(first_word, "AEXP", 4) != 0) &&
00433 (strncasecmp(first_word, "REXP", 4) != 0) &&
00434 (strncasecmp(first_word, "HAEX", 4) != 0) &&
00435 (strncasecmp(first_word, "AAEX", 4) != 0) &&
00436 (strncasecmp(first_word, "NBOND", 5) != 0) &&
00437 (strncasecmp(first_word, "CUTNB", 5) != 0) &&
00438 (strncasecmp(first_word, "END", 3) != 0) &&
00439 (strncasecmp(first_word, "CTONN", 5) != 0) &&
00440 (strncasecmp(first_word, "EPS", 3) != 0) &&
00441 (strncasecmp(first_word, "VSWI", 4) != 0) &&
00442 (strncasecmp(first_word, "NBXM", 4) != 0) &&
00443 (strncasecmp(first_word, "INHI", 4) != 0) )
00444 {
00445
00446
00447 if (strncasecmp(first_word, "bond", 4)==0)
00448 {
00449 add_bond_param(buffer);
00450 NumBondParams++;
00451 }
00452 else if (strncasecmp(first_word, "angl", 4)==0)
00453 {
00454 add_angle_param(buffer);
00455 NumAngleParams++;
00456 }
00457 else if (strncasecmp(first_word, "dihe", 4)==0)
00458 {
00459 add_dihedral_param(buffer, pfile);
00460 NumDihedralParams++;
00461 }
00462 else if (strncasecmp(first_word, "impr", 4)==0)
00463 {
00464 add_improper_param(buffer, pfile);
00465 NumImproperParams++;
00466 }
00467 else if (strncasecmp(first_word, "nonb", 4)==0)
00468 {
00469 add_vdw_param(buffer);
00470 NumVdwParams++;
00471 }
00472 else if (strncasecmp(first_word, "nbfi", 4)==0)
00473 {
00474 add_vdw_pair_param(buffer);
00475 NumVdwPairParams++;
00476 }
00477 else if (strncasecmp(first_word, "nbta", 4)==0)
00478 {
00479
00480 if (table_ener == NULL) {
00481 continue;
00482 }
00483
00484 add_table_pair_param(buffer);
00485 NumTablePairParams++;
00486 }
00487 else if (strncasecmp(first_word, "hbon", 4)==0)
00488 {
00489 add_hb_pair_param(buffer);
00490 }
00491 else
00492 {
00493
00494
00495 char err_msg[512];
00496
00497 sprintf(err_msg, "UNKNOWN PARAMETER IN XPLOR PARAMETER FILE %s\nLINE=*%s*",
00498 fname, buffer);
00499 NAMD_die(err_msg);
00500 }
00501 }
00502 }
00503
00504
00505 Fclose(pfile);
00506
00507 return;
00508 }
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528 void Parameters::read_charmm_parameter_file(char *fname)
00529
00530 {
00531 int par_type=0;
00532
00533 int skipline;
00534 int skipall = 0;
00535 char buffer[512];
00536 char first_word[512];
00537 FILE *pfile;
00538
00539
00540
00541 if (AllFilesRead)
00542 {
00543 NAMD_die("Tried to read another parameter file after being told that all files were read!");
00544 }
00545
00546
00547 if ( (pfile = fopen(fname, "r")) == NULL)
00548 {
00549 char err_msg[256];
00550
00551 sprintf(err_msg, "UNABLE TO OPEN CHARMM PARAMETER FILE %s\n", fname);
00552 NAMD_die(err_msg);
00553 }
00554
00555
00556 while (NAMD_read_line(pfile, buffer) != -1)
00557 {
00558
00559 NAMD_find_first_word(buffer, first_word);
00560 skipline=0;
00561
00562
00563
00564
00565 if (!NAMD_blank_string(buffer) &&
00566 (strncmp(first_word, "!", 1) != 0) &&
00567 (strncmp(first_word, "*", 1) != 0) &&
00568 (strncasecmp(first_word, "END", 3) != 0))
00569 {
00570 if ( skipall ) {
00571 iout << iWARN << "SKIPPING PART OF PARAMETER FILE AFTER RETURN STATEMENT\n" << endi;
00572 break;
00573 }
00574
00575 if (strncasecmp(first_word, "bond", 4)==0)
00576 {
00577 par_type=1; skipline=1;
00578 }
00579 else if (strncasecmp(first_word, "angl", 4)==0)
00580 {
00581 par_type=2; skipline=1;
00582 }
00583 else if (strncasecmp(first_word, "thet", 4)==0)
00584 {
00585 par_type=2; skipline=1;
00586 }
00587 else if (strncasecmp(first_word, "dihe", 4)==0)
00588 {
00589 par_type=3; skipline=1;
00590 }
00591 else if (strncasecmp(first_word, "phi", 3)==0)
00592 {
00593 par_type=3; skipline=1;
00594 }
00595 else if (strncasecmp(first_word, "impr", 4)==0)
00596 {
00597 par_type=4; skipline=1;
00598 }
00599 else if (strncasecmp(first_word, "imph", 4)==0)
00600 {
00601 par_type=4; skipline=1;
00602 }
00603 else if (strncasecmp(first_word, "nbon", 4)==0)
00604 {
00605 par_type=5; skipline=1;
00606 }
00607 else if (strncasecmp(first_word, "nonb", 4)==0)
00608 {
00609 par_type=5; skipline=1;
00610 }
00611 else if (strncasecmp(first_word, "nbfi", 4)==0)
00612 {
00613 par_type=6; skipline=1;
00614 }
00615 else if (strncasecmp(first_word, "hbon", 4)==0)
00616 {
00617 par_type=7; skipline=1;
00618 }
00619 else if (strncasecmp(first_word, "cmap", 4)==0)
00620 {
00621 par_type=8; skipline=1;
00622 }
00623 else if (strncasecmp(first_word, "nbta", 4)==0)
00624 {
00625 par_type=9; skipline=1;
00626 }
00627 else if (strncasecmp(first_word, "thol", 4)==0)
00628 {
00629 par_type=10; skipline=1;
00630 }
00631 else if (strncasecmp(first_word, "atom", 4)==0)
00632 {
00633 par_type=11; skipline=1;
00634 }
00635 else if (strncasecmp(first_word, "read", 4)==0)
00636 {
00637 skip_stream_read(buffer, pfile); skipline=1;
00638 }
00639 else if (strncasecmp(first_word, "return", 4)==0)
00640 {
00641 skipall=1; skipline=1;
00642 }
00643 else if ((strncasecmp(first_word, "nbxm", 4) == 0) ||
00644 (strncasecmp(first_word, "grou", 4) == 0) ||
00645 (strncasecmp(first_word, "cdie", 4) == 0) ||
00646 (strncasecmp(first_word, "shif", 4) == 0) ||
00647 (strncasecmp(first_word, "vgro", 4) == 0) ||
00648 (strncasecmp(first_word, "vdis", 4) == 0) ||
00649 (strncasecmp(first_word, "vswi", 4) == 0) ||
00650 (strncasecmp(first_word, "cutn", 4) == 0) ||
00651 (strncasecmp(first_word, "ctof", 4) == 0) ||
00652 (strncasecmp(first_word, "cton", 4) == 0) ||
00653 (strncasecmp(first_word, "eps", 3) == 0) ||
00654 (strncasecmp(first_word, "e14f", 4) == 0) ||
00655 (strncasecmp(first_word, "wmin", 4) == 0) ||
00656 (strncasecmp(first_word, "aexp", 4) == 0) ||
00657 (strncasecmp(first_word, "rexp", 4) == 0) ||
00658 (strncasecmp(first_word, "haex", 4) == 0) ||
00659 (strncasecmp(first_word, "aaex", 4) == 0) ||
00660 (strncasecmp(first_word, "noac", 4) == 0) ||
00661 (strncasecmp(first_word, "hbno", 4) == 0) ||
00662 (strncasecmp(first_word, "cuth", 4) == 0) ||
00663 (strncasecmp(first_word, "ctof", 4) == 0) ||
00664 (strncasecmp(first_word, "cton", 4) == 0) ||
00665 (strncasecmp(first_word, "cuth", 4) == 0) ||
00666 (strncasecmp(first_word, "ctof", 4) == 0) ||
00667 (strncasecmp(first_word, "cton", 4) == 0) )
00668 {
00669 if ((par_type != 5) && (par_type != 6) && (par_type != 7) && (par_type != 9))
00670 {
00671 char err_msg[512];
00672
00673 sprintf(err_msg, "ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00674 NAMD_die(err_msg);
00675 }
00676 else
00677 {
00678 skipline = 1;
00679 }
00680 }
00681 else if (par_type == 0)
00682 {
00683
00684
00685 char err_msg[512];
00686
00687 sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00688 NAMD_die(err_msg);
00689 }
00690 }
00691 else
00692 {
00693 skipline=1;
00694 }
00695
00696 if ( (par_type != 0) && (!skipline) )
00697 {
00698
00699
00700
00701 if (par_type == 1)
00702 {
00703 add_bond_param(buffer);
00704 NumBondParams++;
00705 }
00706 else if (par_type == 2)
00707 {
00708 add_angle_param(buffer);
00709 NumAngleParams++;
00710 }
00711 else if (par_type == 3)
00712 {
00713 add_dihedral_param(buffer, pfile);
00714 NumDihedralParams++;
00715 }
00716 else if (par_type == 4)
00717 {
00718 add_improper_param(buffer, pfile);
00719 NumImproperParams++;
00720 }
00721 else if (par_type == 5)
00722 {
00723 add_vdw_param(buffer);
00724 NumVdwParams++;
00725 }
00726 else if (par_type == 6)
00727 {
00728 add_vdw_pair_param(buffer);
00729 NumVdwPairParams++;
00730 }
00731 else if (par_type == 7)
00732 {
00733 add_hb_pair_param(buffer);
00734 }
00735 else if (par_type == 8)
00736 {
00737 add_crossterm_param(buffer, pfile);
00738 NumCrosstermParams++;
00739 }
00740 else if (par_type == 9)
00741 {
00742
00743 if (table_ener == NULL) {
00744 continue;
00745 }
00746
00747 add_table_pair_param(buffer);
00748 NumTablePairParams++;
00749 }
00750 else if (par_type == 10)
00751 {
00752 add_nbthole_pair_param(buffer);
00753 NumNbtholePairParams++;
00754 }
00755 else if (par_type == 11)
00756 {
00757 if ( strncasecmp(first_word, "mass", 4) != 0 ) {
00758 char err_msg[512];
00759 sprintf(err_msg, "UNKNOWN PARAMETER IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00760 NAMD_die(err_msg);
00761 }
00762 }
00763 else
00764 {
00765
00766
00767
00768 char err_msg[512];
00769
00770 sprintf(err_msg, "INTERNAL ERROR IN CHARMM PARAMETER FILE %s\nLINE=*%s*",fname, buffer);
00771 NAMD_die(err_msg);
00772 }
00773 }
00774 }
00775
00776
00777 fclose(pfile);
00778
00779 return;
00780 }
00781
00782
00783
00784
00785 void Parameters::skip_stream_read(char *buf, FILE *fd) {
00786
00787 char buffer[513];
00788 char first_word[513];
00789 char s1[128];
00790 char s2[128];
00791 int rval = sscanf(buf, "%s %s", s1, s2);
00792 if (rval != 2) {
00793 char err_msg[512];
00794 sprintf(err_msg, "BAD FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
00795 NAMD_die(err_msg);
00796 }
00797 if ( ! strncasecmp(s2,"PARA",4) ) return;
00798
00799 iout << iINFO << "SKIPPING " << s2 << " SECTION IN STREAM FILE\n" << endi;
00800
00801 while (NAMD_read_line(fd, buffer) != -1)
00802 {
00803
00804 NAMD_find_first_word(buffer, first_word);
00805 if (!NAMD_blank_string(buffer) &&
00806 (strncmp(first_word, "!", 1) != 0) &&
00807 (strncmp(first_word, "*", 1) != 0) &&
00808 (strncasecmp(first_word, "END", 3) == 0) ) return;
00809 }
00810
00811 }
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827 void Parameters::add_bond_param(char *buf)
00828
00829 {
00830 char atom1name[11];
00831 char atom2name[11];
00832 Real forceconstant;
00833 Real distance;
00834 int read_count;
00835 struct bond_params *new_node;
00836
00837
00838
00839 if (paramType == paraXplor)
00840 {
00841
00842 read_count=sscanf(buf, "%*s %s %s %f %f\n", atom1name, atom2name,
00843 &forceconstant, &distance);
00844 }
00845 else if (paramType == paraCharmm)
00846 {
00847
00848 read_count=sscanf(buf, "%s %s %f %f\n", atom1name, atom2name,
00849 &forceconstant, &distance);
00850 }
00851
00852
00853
00854 if (read_count != 4)
00855 {
00856 char err_msg[512];
00857
00858 if (paramType == paraXplor)
00859 {
00860 sprintf(err_msg, "BAD BOND FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
00861 }
00862 else if (paramType == paraCharmm)
00863 {
00864 sprintf(err_msg, "BAD BOND FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
00865 }
00866 NAMD_die(err_msg);
00867 }
00868
00869
00870 new_node = new bond_params;
00871
00872 if (new_node == NULL)
00873 {
00874 NAMD_die("memory allocation failed in Parameters::add_bond_param\n");
00875 }
00876
00877
00878
00879
00880
00881 if (strcasecmp(atom1name, atom2name) < 0)
00882 {
00883 strcpy(new_node->atom1name, atom1name);
00884 strcpy(new_node->atom2name, atom2name);
00885 }
00886 else
00887 {
00888 strcpy(new_node->atom2name, atom1name);
00889 strcpy(new_node->atom1name, atom2name);
00890 }
00891
00892
00893 new_node->forceconstant = forceconstant;
00894 new_node->distance = distance;
00895
00896
00897 new_node->left = NULL;
00898 new_node->right = NULL;
00899
00900
00901
00902 bondp=add_to_bond_tree(new_node, bondp);
00903
00904 return;
00905 }
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926 struct bond_params *Parameters::add_to_bond_tree(struct bond_params *new_node,
00927 struct bond_params *tree)
00928
00929 {
00930 int compare_code;
00931
00932
00933
00934 if (tree == NULL)
00935 return(new_node);
00936
00937
00938
00939 compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
00940
00941
00942 if (compare_code == 0)
00943 {
00944
00945 compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
00946
00947
00948 if (compare_code == 0)
00949 {
00950
00951
00952
00953
00954
00955 if ((tree->forceconstant != new_node->forceconstant) ||
00956 (tree->distance != new_node->distance))
00957 {
00958 iout << "\n" << iWARN << "DUPLICATE BOND ENTRY FOR "
00959 << new_node->atom1name << "-"
00960 << new_node->atom2name
00961 << "\nPREVIOUS VALUES k=" << tree->forceconstant
00962 << " x0=" << tree->distance
00963 << "\n USING VALUES k=" << new_node->forceconstant
00964 << " x0=" << new_node->distance
00965 << "\n" << endi;
00966
00967 tree->forceconstant=new_node->forceconstant;
00968 tree->distance=new_node->distance;
00969 }
00970
00971
00972 delete new_node;
00973
00974 return(tree);
00975 }
00976 }
00977
00978
00979
00980
00981 if (compare_code < 0)
00982 {
00983 tree->left = add_to_bond_tree(new_node, tree->left);
00984 }
00985 else
00986 {
00987 tree->right = add_to_bond_tree(new_node, tree->right);
00988 }
00989
00990 return(tree);
00991 }
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007 void Parameters::add_angle_param(char *buf)
01008
01009 {
01010 char atom1name[11];
01011 char atom2name[11];
01012 char atom3name[11];
01013 char norm[4]="xxx";
01014 Real forceconstant;
01015 Real angle;
01016 Real k_ub;
01017 Real r_ub;
01018 int read_count;
01019 struct angle_params *new_node;
01020
01021
01022
01023 if (paramType == paraXplor)
01024 {
01025
01026 read_count=sscanf(buf, "%*s %s %s %s %f %f UB %f %f\n",
01027 atom1name, atom2name, atom3name, &forceconstant, &angle,
01028 &k_ub, &r_ub);
01029 }
01030 else if ((paramType == paraCharmm) && cosAngles) {
01031 read_count=sscanf(buf, "%s %s %s %f %f %3s %f %f\n",
01032 atom1name, atom2name, atom3name, &forceconstant, &angle, norm,
01033 &k_ub, &r_ub);
01034
01035
01036 }
01037 else if (paramType == paraCharmm)
01038 {
01039
01040 read_count=sscanf(buf, "%s %s %s %f %f %f %f\n",
01041 atom1name, atom2name, atom3name, &forceconstant, &angle,
01042 &k_ub, &r_ub);
01043
01044
01045 }
01046
01047
01048
01049 if ( (read_count != 5) && (read_count != 7) && !(cosAngles && read_count == 6))
01050 {
01051 char err_msg[512];
01052
01053 if (paramType == paraXplor)
01054 {
01055 sprintf(err_msg, "BAD ANGLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01056 }
01057 else if (paramType == paraCharmm)
01058 {
01059 sprintf(err_msg, "BAD ANGLE FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
01060 }
01061 NAMD_die(err_msg);
01062 }
01063
01064
01065 new_node = new angle_params;
01066
01067 if (new_node == NULL)
01068 {
01069 NAMD_die("memory allocation failed in Parameters::add_angle_param");
01070 }
01071
01072
01073
01074
01075 if (strcasecmp(atom1name, atom3name) < 0)
01076 {
01077 strcpy(new_node->atom1name, atom1name);
01078 strcpy(new_node->atom2name, atom2name);
01079 strcpy(new_node->atom3name, atom3name);
01080 }
01081 else
01082 {
01083 strcpy(new_node->atom3name, atom1name);
01084 strcpy(new_node->atom2name, atom2name);
01085 strcpy(new_node->atom1name, atom3name);
01086 }
01087
01088
01089 new_node->forceconstant = forceconstant;
01090 new_node->angle = angle;
01091
01092 if (cosAngles) {
01093 if (strcasecmp("cos",norm)==0) {
01094
01095 NumCosAngles++;
01096 new_node->normal = 0;
01097 } else {
01098
01099 new_node->normal = 1;
01100 }
01101 } else {
01102 new_node->normal = 1;
01103 }
01104
01105 if (read_count == 7)
01106 {
01107
01108 if (new_node->normal == 0) {
01109 NAMD_die("ERROR: Urey-Bradley angles can't be used with cosine-based terms\n");
01110 }
01111 new_node->k_ub = k_ub;
01112 new_node->r_ub = r_ub;
01113 }
01114 else
01115 {
01116 new_node->k_ub = 0.0;
01117 new_node->r_ub = 0.0;
01118 }
01119
01120 new_node->left = NULL;
01121 new_node->right = NULL;
01122
01123
01124 anglep = add_to_angle_tree(new_node, anglep);
01125
01126 return;
01127 }
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150 struct angle_params *Parameters::add_to_angle_tree(struct angle_params *new_node,
01151 struct angle_params *tree)
01152
01153 {
01154 int compare_code;
01155
01156
01157 if (tree == NULL)
01158 return(new_node);
01159
01160
01161 compare_code = strcasecmp(new_node->atom1name, tree->atom1name);
01162
01163 if (compare_code == 0)
01164 {
01165
01166 compare_code = strcasecmp(new_node->atom2name, tree->atom2name);
01167
01168 if (compare_code == 0)
01169 {
01170
01171 compare_code = strcasecmp(new_node->atom3name,
01172 tree->atom3name);
01173
01174 if (compare_code == 0)
01175 {
01176
01177
01178
01179
01180
01181
01182 if ((tree->forceconstant != new_node->forceconstant) ||
01183 (tree->angle != new_node->angle) ||
01184 (tree->k_ub != new_node->k_ub) ||
01185 (tree->r_ub != new_node->r_ub) || (tree->normal != new_node->normal))
01186 {
01187 iout << "\n" << iWARN << "DUPLICATE ANGLE ENTRY FOR "
01188 << new_node->atom1name << "-"
01189 << new_node->atom2name << "-"
01190 << new_node->atom3name
01191 << "\nPREVIOUS VALUES k="
01192 << tree->forceconstant << " theta0="
01193 << tree->angle << " k_ub="
01194 << tree->k_ub << " r_ub="
01195 << tree->r_ub
01196 << "\n USING VALUES k="
01197 << new_node->forceconstant << " theta0="
01198 << new_node->angle << " k_ub="
01199 << new_node->k_ub << " r_ub=" << new_node->r_ub
01200 << "\n" << endi;
01201
01202 tree->forceconstant=new_node->forceconstant;
01203 tree->angle=new_node->angle;
01204 tree->k_ub=new_node->k_ub;
01205 tree->r_ub=new_node->r_ub;
01206 tree->normal=new_node->normal;
01207 }
01208
01209
01210 delete new_node;
01211
01212 return(tree);
01213 }
01214 }
01215 }
01216
01217
01218
01219
01220 if (compare_code < 0)
01221 {
01222 tree->left = add_to_angle_tree(new_node, tree->left);
01223 }
01224 else
01225 {
01226 tree->right = add_to_angle_tree(new_node, tree->right);
01227 }
01228
01229 return(tree);
01230 }
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 void Parameters::add_dihedral_param(char *buf, FILE *fd)
01247
01248 {
01249 char atom1name[11];
01250 char atom2name[11];
01251 char atom3name[11];
01252 char atom4name[11];
01253 Real forceconstant;
01254 int periodicity;
01255 Real phase_shift;
01256 int read_count;
01257 struct dihedral_params *new_node;
01258 int multiplicity;
01259 int i;
01260 char buffer[513];
01261 int ret_code;
01262
01263
01264
01265 if (paramType == paraXplor)
01266 {
01267
01268 read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n",
01269 atom1name, atom2name, atom3name, atom4name, &multiplicity,
01270 &forceconstant, &periodicity, &phase_shift);
01271 }
01272 else if (paramType == paraCharmm)
01273 {
01274
01275 read_count=sscanf(buf, "%s %s %s %s %f %d %f\n",
01276 atom1name, atom2name, atom3name, atom4name,
01277 &forceconstant, &periodicity, &phase_shift);
01278 multiplicity=1;
01279 }
01280
01281 if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
01282 {
01283 char err_msg[512];
01284
01285 sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01286 NAMD_die(err_msg);
01287 }
01288 else if ( (read_count != 7) && (paramType == paraCharmm) )
01289 {
01290 char err_msg[512];
01291
01292 sprintf(err_msg, "BAD DIHEDRAL FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
01293 NAMD_die(err_msg);
01294 }
01295
01296 if ( (read_count == 4) && (paramType == paraXplor) )
01297
01298 {
01299 read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n",
01300 &forceconstant, &periodicity, &phase_shift);
01301
01302
01303 if (read_count != 3)
01304 {
01305 char err_msg[512];
01306
01307 sprintf(err_msg, "BAD DIHEDRAL FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01308 NAMD_die(err_msg);
01309 }
01310
01311 multiplicity = 1;
01312 }
01313
01314 if (multiplicity > MAX_MULTIPLICITY)
01315 {
01316 char err_msg[181];
01317
01318 sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
01319 multiplicity, MAX_MULTIPLICITY);
01320 NAMD_die(err_msg);
01321 }
01322
01323
01324 new_node = new dihedral_params;
01325
01326 if (new_node == NULL)
01327 {
01328 NAMD_die("memory allocation failed in Parameters::add_dihedral_param\n");
01329 }
01330
01331
01332
01333
01334
01335 strcpy(new_node->atom1name, atom1name);
01336 strcpy(new_node->atom2name, atom2name);
01337 strcpy(new_node->atom3name, atom3name);
01338 strcpy(new_node->atom4name, atom4name);
01339 new_node->atom1wild = ! strcasecmp(atom1name, "X");
01340 new_node->atom2wild = ! strcasecmp(atom2name, "X");
01341 new_node->atom3wild = ! strcasecmp(atom3name, "X");
01342 new_node->atom4wild = ! strcasecmp(atom4name, "X");
01343 new_node->multiplicity = multiplicity;
01344 if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01345 new_node->values[0].k = forceconstant;
01346 new_node->values[0].n = periodicity;
01347 new_node->values[0].delta = phase_shift;
01348
01349 new_node->next = NULL;
01350
01351
01352 if (multiplicity > 1)
01353 {
01354 for (i=1; i<multiplicity; i++)
01355 {
01356 ret_code = NAMD_read_line(fd, buffer);
01357
01358
01359 if (ret_code == 0)
01360 {
01361 NAMD_remove_comment(buffer);
01362 }
01363
01364
01365 while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
01366 {
01367 ret_code = NAMD_read_line(fd, buffer);
01368 }
01369
01370 if (ret_code != 0)
01371 {
01372 NAMD_die("EOF encoutner in middle of multiple dihedral");
01373 }
01374
01375 read_count=sscanf(buffer, "%f %d %f\n",
01376 &forceconstant, &periodicity, &phase_shift);
01377
01378 if (read_count != 3)
01379 {
01380 char err_msg[512];
01381
01382 sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
01383 NAMD_die(err_msg);
01384 }
01385
01386 if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01387 new_node->values[i].k = forceconstant;
01388 new_node->values[i].n = periodicity;
01389 new_node->values[i].delta = phase_shift;
01390 }
01391 }
01392
01393
01394
01395 if (paramType == paraXplor)
01396 {
01397 add_to_dihedral_list(new_node);
01398 }
01399 else if (paramType == paraCharmm)
01400 {
01401 add_to_charmm_dihedral_list(new_node);
01402 }
01403
01404
01405 return;
01406 }
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427 void Parameters::add_to_dihedral_list(
01428 struct dihedral_params *new_node)
01429
01430 {
01431 static struct dihedral_params *ptr;
01432 static struct dihedral_params *tail;
01433
01434
01435
01436 int i;
01437
01438
01439 if (dihedralp == NULL)
01440 {
01441 dihedralp=new_node;
01442 tail=new_node;
01443
01444 return;
01445 }
01446
01447
01448 ptr=dihedralp;
01449
01450 while (ptr != NULL)
01451 {
01452 if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01453 (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01454 (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01455 (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01456 ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01457 (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01458 (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01459 (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
01460 {
01461
01462
01463
01464 int echoWarn=0;
01465
01466 if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
01467
01468 if (!echoWarn)
01469 {
01470 for (i=0; i<ptr->multiplicity; i++)
01471 {
01472 if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
01473 if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
01474 if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
01475 }
01476 }
01477
01478 if (echoWarn)
01479 {
01480 iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
01481 << ptr->atom1name << "-"
01482 << ptr->atom2name << "-"
01483 << ptr->atom3name << "-"
01484 << ptr->atom4name
01485 << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
01486
01487 for (i=0; i<ptr->multiplicity; i++)
01488 {
01489 iout << " k=" << ptr->values[i].k
01490 << " n=" << ptr->values[i].n
01491 << " delta=" << ptr->values[i].delta;
01492 }
01493
01494 iout << "\nUSING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
01495
01496 for (i=0; i<new_node->multiplicity; i++)
01497 {
01498 iout << " k=" << new_node->values[i].k
01499 << " n=" << new_node->values[i].n
01500 << " delta=" << new_node->values[i].delta;
01501 }
01502
01503 iout << endi;
01504
01505 ptr->multiplicity = new_node->multiplicity;
01506
01507 for (i=0; i<new_node->multiplicity; i++)
01508 {
01509 ptr->values[i].k = new_node->values[i].k;
01510 ptr->values[i].n = new_node->values[i].n;
01511 ptr->values[i].delta = new_node->values[i].delta;
01512 }
01513
01514 }
01515
01516
01517 delete new_node;
01518
01519 return;
01520 }
01521
01522 ptr=ptr->next;
01523 }
01524
01525
01526
01527
01528
01529
01530
01531 if ( new_node->atom1wild ||
01532 new_node->atom2wild ||
01533 new_node->atom3wild ||
01534 new_node->atom4wild )
01535 {
01536
01537 tail->next=new_node;
01538 tail=new_node;
01539
01540 return;
01541 }
01542 else
01543 {
01544
01545 new_node->next=dihedralp;
01546 dihedralp=new_node;
01547
01548 return;
01549 }
01550
01551 }
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576 void Parameters::add_to_charmm_dihedral_list(
01577 struct dihedral_params *new_node)
01578
01579 {
01580 static struct dihedral_params *ptr;
01581 static struct dihedral_params *tail;
01582
01583
01584
01585 int i;
01586 int replace;
01587
01588
01589
01590 static struct dihedral_params last_dihedral;
01591
01592
01593 if (dihedralp == NULL)
01594 {
01595 dihedralp=new_node;
01596 tail=new_node;
01597 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01598
01599 return;
01600 }
01601
01602
01603 ptr=dihedralp;
01604
01605 while (ptr != NULL)
01606 {
01607 int same_as_last = 0;
01608 if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01609 (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01610 (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01611 (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01612 ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01613 (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01614 (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01615 (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) )
01616 )
01617 {
01618
01619
01620
01621
01622
01623
01624 if ( ( !strcmp(ptr->atom1name, last_dihedral.atom1name) &&
01625 !strcmp(ptr->atom2name, last_dihedral.atom2name) &&
01626 !strcmp(ptr->atom3name, last_dihedral.atom3name) &&
01627 !strcmp(ptr->atom4name, last_dihedral.atom4name)))
01628 same_as_last = 1;
01629
01630
01631
01632 int echoWarn=1;
01633
01634
01635 for (i=0; i<ptr->multiplicity; i++)
01636 {
01637 if ((ptr->values[i].k == new_node->values[0].k) &&
01638 (ptr->values[i].n == new_node->values[0].n) &&
01639 (ptr->values[i].delta == new_node->values[0].delta))
01640 {
01641
01642 echoWarn=0;
01643 break;
01644 }
01645
01646 }
01647
01648 if (echoWarn)
01649 {
01650 if (!same_as_last) {
01651 iout << "\n" << iWARN << "DUPLICATE DIHEDRAL ENTRY FOR "
01652 << ptr->atom1name << "-"
01653 << ptr->atom2name << "-"
01654 << ptr->atom3name << "-"
01655 << ptr->atom4name
01656 << "\nPREVIOUS VALUES MULTIPLICITY: " << ptr->multiplicity << "\n";
01657 }
01658 replace=0;
01659
01660 for (i=0; i<ptr->multiplicity; i++)
01661 {
01662 if (!same_as_last) {
01663 iout << " k=" << ptr->values[i].k
01664 << " n=" << ptr->values[i].n
01665 << " delta=" << ptr->values[i].delta << "\n";
01666 }
01667 if (ptr->values[i].n == new_node->values[0].n)
01668 {
01669 iout << iWARN << "IDENTICAL PERIODICITY! REPLACING OLD VALUES BY: \n";
01670 ptr->values[i].k = new_node->values[0].k;
01671 ptr->values[i].delta = new_node->values[0].delta;
01672 iout << " k=" << ptr->values[i].k
01673 << " n=" << ptr->values[i].n
01674 << " delta=" << ptr->values[i].delta<< "\n";
01675 replace=1;
01676 break;
01677 }
01678 }
01679
01680 if (!replace)
01681 {
01682 ptr->multiplicity += 1;
01683
01684 if (ptr->multiplicity > MAX_MULTIPLICITY)
01685 {
01686 char err_msg[181];
01687
01688 sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
01689 ptr->multiplicity, MAX_MULTIPLICITY);
01690 NAMD_die(err_msg);
01691 }
01692 if (!same_as_last)
01693 iout << "INCREASING MULTIPLICITY TO: " << ptr->multiplicity << "\n";
01694
01695 i= ptr->multiplicity - 1;
01696 ptr->values[i].k = new_node->values[0].k;
01697 ptr->values[i].n = new_node->values[0].n;
01698 ptr->values[i].delta = new_node->values[0].delta;
01699
01700 if (!same_as_last)
01701 iout << " k=" << ptr->values[i].k
01702 << " n=" << ptr->values[i].n
01703 << " delta=" << ptr->values[i].delta<< "\n";
01704 }
01705
01706 iout << endi;
01707 }
01708
01709
01710 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01711 delete new_node;
01712
01713 return;
01714 }
01715
01716 ptr=ptr->next;
01717 }
01718
01719
01720
01721
01722
01723
01724
01725
01726 if ( new_node->atom1wild ||
01727 new_node->atom2wild ||
01728 new_node->atom3wild ||
01729 new_node->atom4wild )
01730 {
01731
01732 tail->next=new_node;
01733 tail=new_node;
01734
01735 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01736 return;
01737 }
01738 else
01739 {
01740
01741 new_node->next=dihedralp;
01742 dihedralp=new_node;
01743
01744 memcpy(&last_dihedral, new_node, sizeof(dihedral_params));
01745 return;
01746 }
01747
01748 }
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765 void Parameters::add_improper_param(char *buf, FILE *fd)
01766
01767 {
01768 char atom1name[11];
01769 char atom2name[11];
01770 char atom3name[11];
01771 char atom4name[11];
01772 Real forceconstant;
01773 int periodicity;
01774 Real phase_shift;
01775 int read_count;
01776 struct improper_params *new_node;
01777 int multiplicity;
01778 int i;
01779 char buffer[513];
01780 int ret_code;
01781
01782
01783
01784 if (paramType == paraXplor)
01785 {
01786
01787 read_count=sscanf(buf, "%*s %s %s %s %s MULTIPLE= %d %f %d %f\n",
01788 atom1name, atom2name, atom3name, atom4name, &multiplicity,
01789 &forceconstant, &periodicity, &phase_shift);
01790 }
01791 else if (paramType == paraCharmm)
01792 {
01793
01794 read_count=sscanf(buf, "%s %s %s %s %f %d %f\n",
01795 atom1name, atom2name, atom3name, atom4name,
01796 &forceconstant, &periodicity, &phase_shift);
01797 multiplicity=1;
01798 }
01799
01800 if ( (read_count != 4) && (read_count != 8) && (paramType == paraXplor) )
01801 {
01802 char err_msg[512];
01803
01804 sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
01805 NAMD_die(err_msg);
01806 }
01807 else if ( (read_count != 7) && (paramType == paraCharmm) )
01808 {
01809 char err_msg[512];
01810
01811 sprintf(err_msg, "BAD IMPROPER FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
01812 NAMD_die(err_msg);
01813 }
01814
01815 if ( (read_count == 4) && (paramType == paraXplor) )
01816
01817 {
01818 read_count=sscanf(buf, "%*s %*s %*s %*s %*s %f %d %f\n",
01819 &forceconstant, &periodicity, &phase_shift);
01820
01821
01822 if (read_count != 3)
01823 {
01824 char err_msg[512];
01825
01826 sprintf(err_msg, "BAD IMPROPER FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
01827 NAMD_die(err_msg);
01828 }
01829
01830 multiplicity = 1;
01831 }
01832
01833 if (multiplicity > MAX_MULTIPLICITY)
01834 {
01835 char err_msg[181];
01836
01837 sprintf(err_msg, "Multiple improper with multiplicity of %d greater than max of %d",
01838 multiplicity, MAX_MULTIPLICITY);
01839 NAMD_die(err_msg);
01840 }
01841
01842
01843 new_node = new improper_params;
01844
01845 if (new_node == NULL)
01846 {
01847 NAMD_die("memory allocation failed in Parameters::add_improper_param");
01848 }
01849
01850
01851
01852 strcpy(new_node->atom1name, atom1name);
01853 strcpy(new_node->atom2name, atom2name);
01854 strcpy(new_node->atom3name, atom3name);
01855 strcpy(new_node->atom4name, atom4name);
01856 new_node->multiplicity = multiplicity;
01857 if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01858 new_node->values[0].k = forceconstant;
01859 new_node->values[0].n = periodicity;
01860 new_node->values[0].delta = phase_shift;
01861
01862 new_node->next = NULL;
01863
01864
01865 if (multiplicity > 1)
01866 {
01867
01868 for (i=1; i<multiplicity; i++)
01869 {
01870 ret_code = NAMD_read_line(fd, buffer);
01871
01872
01873 if (ret_code == 0)
01874 {
01875 NAMD_remove_comment(buffer);
01876 }
01877
01878
01879 while ( (ret_code == 0) && (NAMD_blank_string(buffer)) )
01880 {
01881 ret_code = NAMD_read_line(fd, buffer);
01882 }
01883
01884 if (ret_code != 0)
01885 {
01886 NAMD_die("EOF encoutner in middle of multiple improper");
01887 }
01888
01889
01890 read_count=sscanf(buffer, "%f %d %f\n",
01891 &forceconstant, &periodicity, &phase_shift);
01892
01893 if (read_count != 3)
01894 {
01895 char err_msg[512];
01896
01897 sprintf(err_msg, "BAD MULTIPLE FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buffer);
01898 NAMD_die(err_msg);
01899 }
01900
01901 if (paramType == paraXplor && periodicity != 0) phase_shift *= -1;
01902 new_node->values[i].k = forceconstant;
01903 new_node->values[i].n = periodicity;
01904 new_node->values[i].delta = phase_shift;
01905 }
01906 }
01907
01908
01909 add_to_improper_list(new_node);
01910
01911 return;
01912 }
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933 void Parameters::add_to_improper_list(struct improper_params *new_node)
01934
01935 {
01936 int i;
01937 static struct improper_params *ptr;
01938 static struct improper_params *tail;
01939
01940
01941
01942
01943
01944 if (improperp == NULL)
01945 {
01946 improperp=new_node;
01947 tail=new_node;
01948
01949 return;
01950 }
01951
01952
01953 ptr=improperp;
01954
01955 while (ptr != NULL)
01956 {
01957 if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
01958 (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
01959 (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
01960 (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) ) ||
01961 ( (strcasecmp(new_node->atom4name, ptr->atom1name) == 0) &&
01962 (strcasecmp(new_node->atom3name, ptr->atom2name) == 0) &&
01963 (strcasecmp(new_node->atom2name, ptr->atom3name) == 0) &&
01964 (strcasecmp(new_node->atom1name, ptr->atom4name) == 0) ) )
01965 {
01966
01967
01968
01969 int echoWarn=0;
01970
01971 if (ptr->multiplicity != new_node->multiplicity) {echoWarn=1;}
01972
01973 if (!echoWarn)
01974 {
01975 for (i=0; i<ptr->multiplicity; i++)
01976 {
01977 if (ptr->values[i].k != new_node->values[i].k) {echoWarn=1; break;}
01978 if (ptr->values[i].n != new_node->values[i].n) {echoWarn=1; break;}
01979 if (ptr->values[i].delta != new_node->values[i].delta) {echoWarn=1; break;}
01980 }
01981 }
01982
01983 if (echoWarn)
01984 {
01985 iout << "\n" << iWARN << "DUPLICATE IMPROPER DIHEDRAL ENTRY FOR "
01986 << ptr->atom1name << "-"
01987 << ptr->atom2name << "-"
01988 << ptr->atom3name << "-"
01989 << ptr->atom4name
01990 << "\nPREVIOUS VALUES MULTIPLICITY " << ptr->multiplicity << "\n";
01991
01992 for (i=0; i<ptr->multiplicity; i++)
01993 {
01994 iout << " k=" << ptr->values[i].k
01995 << " n=" << ptr->values[i].n
01996 << " delta=" << ptr->values[i].delta;
01997 }
01998
01999 iout << "\n" << "USING VALUES MULTIPLICITY " << new_node->multiplicity << "\n";
02000
02001 for (i=0; i<new_node->multiplicity; i++)
02002 {
02003 iout << " k=" << new_node->values[i].k
02004 << " n=" << new_node->values[i].n
02005 << " delta=" << new_node->values[i].delta;
02006 }
02007
02008 iout << endi;
02009
02010 ptr->multiplicity = new_node->multiplicity;
02011
02012 for (i=0; i<new_node->multiplicity; i++)
02013 {
02014 ptr->values[i].k = new_node->values[i].k;
02015 ptr->values[i].n = new_node->values[i].n;
02016 ptr->values[i].delta = new_node->values[i].delta;
02017 }
02018 }
02019
02020
02021 delete new_node;
02022
02023 return;
02024 }
02025
02026 ptr=ptr->next;
02027 }
02028
02029
02030
02031
02032
02033
02034
02035 if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
02036 (strcasecmp(new_node->atom2name, "X") == 0) ||
02037 (strcasecmp(new_node->atom3name, "X") == 0) ||
02038 (strcasecmp(new_node->atom4name, "X") == 0) )
02039 {
02040
02041 tail->next=new_node;
02042 tail=new_node;
02043
02044 return;
02045 }
02046 else
02047 {
02048
02049 new_node->next=improperp;
02050 improperp=new_node;
02051
02052 return;
02053 }
02054 }
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070 void Parameters::add_crossterm_param(char *buf, FILE *fd)
02071
02072 {
02073 char atom1name[11];
02074 char atom2name[11];
02075 char atom3name[11];
02076 char atom4name[11];
02077 char atom5name[11];
02078 char atom6name[11];
02079 char atom7name[11];
02080 char atom8name[11];
02081 int dimension;
02082 int read_count;
02083 struct crossterm_params *new_node;
02084 char buffer[513];
02085 int ret_code;
02086
02087
02088 read_count=sscanf(buf, "%s %s %s %s %s %s %s %s %d\n",
02089 atom1name, atom2name, atom3name, atom4name,
02090 atom5name, atom6name, atom7name, atom8name,
02091 &dimension);
02092
02093 if ( (read_count != 9) || dimension < 1 || dimension > 1000 )
02094 {
02095 char err_msg[512];
02096
02097 sprintf(err_msg, "BAD CMAP FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
02098 NAMD_die(err_msg);
02099 }
02100
02101
02102 new_node = new crossterm_params(dimension);
02103
02104
02105
02106 strcpy(new_node->atom1name, atom1name);
02107 strcpy(new_node->atom2name, atom2name);
02108 strcpy(new_node->atom3name, atom3name);
02109 strcpy(new_node->atom4name, atom4name);
02110 strcpy(new_node->atom5name, atom5name);
02111 strcpy(new_node->atom6name, atom6name);
02112 strcpy(new_node->atom7name, atom7name);
02113 strcpy(new_node->atom8name, atom8name);
02114
02115 new_node->next = NULL;
02116
02117 int nterms = dimension * dimension;
02118 int nread = 0;
02119
02120
02121 while ( nread < nterms ) {
02122 ret_code = NAMD_read_line(fd, buffer);
02123
02124
02125 if (ret_code == 0) {
02126 NAMD_remove_comment(buffer);
02127 }
02128
02129
02130 while ( (ret_code == 0) && (NAMD_blank_string(buffer)) ) {
02131 ret_code = NAMD_read_line(fd, buffer);
02132 if (ret_code == 0) {
02133 NAMD_remove_comment(buffer);
02134 }
02135 }
02136
02137 if (ret_code != 0) {
02138 NAMD_die("EOF encoutner in middle of CMAP");
02139 }
02140
02141
02142 read_count=sscanf(buffer, "%lf %lf %lf %lf %lf %lf %lf %lf\n",
02143 new_node->values + nread,
02144 new_node->values + nread+1,
02145 new_node->values + nread+2,
02146 new_node->values + nread+3,
02147 new_node->values + nread+4,
02148 new_node->values + nread+5,
02149 new_node->values + nread+6,
02150 new_node->values + nread+7);
02151
02152 nread += read_count;
02153
02154 if (read_count == 0 || nread > nterms) {
02155 char err_msg[512];
02156
02157 sprintf(err_msg, "BAD CMAP FORMAT IN PARAMETER FILE\nLINE=*%s*\n", buffer);
02158 NAMD_die(err_msg);
02159 }
02160 }
02161
02162
02163 add_to_crossterm_list(new_node);
02164
02165 return;
02166 }
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187 void Parameters::add_to_crossterm_list(struct crossterm_params *new_node)
02188
02189 {
02190 int i;
02191 static struct crossterm_params *ptr;
02192 static struct crossterm_params *tail;
02193
02194
02195
02196
02197
02198 if (crosstermp == NULL)
02199 {
02200 crosstermp=new_node;
02201 tail=new_node;
02202
02203 return;
02204 }
02205
02206
02207 ptr=crosstermp;
02208
02209 while (ptr != NULL)
02210 {
02211 if ( ( (strcasecmp(new_node->atom1name, ptr->atom1name) == 0) &&
02212 (strcasecmp(new_node->atom2name, ptr->atom2name) == 0) &&
02213 (strcasecmp(new_node->atom3name, ptr->atom3name) == 0) &&
02214 (strcasecmp(new_node->atom4name, ptr->atom4name) == 0) &&
02215 (strcasecmp(new_node->atom5name, ptr->atom5name) == 0) &&
02216 (strcasecmp(new_node->atom6name, ptr->atom6name) == 0) &&
02217 (strcasecmp(new_node->atom7name, ptr->atom7name) == 0) &&
02218 (strcasecmp(new_node->atom8name, ptr->atom8name) == 0) ) )
02219 {
02220
02221
02222 int echoWarn=0;
02223
02224 if (ptr->dimension != new_node->dimension) {echoWarn=1;}
02225
02226 if (!echoWarn)
02227 {
02228 int nvals = ptr->dimension * ptr->dimension;
02229 for (i=0; i<nvals; i++)
02230 {
02231 if (ptr->values[i] != new_node->values[i]) {echoWarn=1; break;}
02232 }
02233 }
02234
02235 if (echoWarn)
02236 {
02237 iout << "\n" << iWARN << "DUPLICATE CMAP ENTRY FOR "
02238 << ptr->atom1name << "-"
02239 << ptr->atom2name << "-"
02240 << ptr->atom3name << "-"
02241 << ptr->atom4name << " "
02242 << ptr->atom5name << "-"
02243 << ptr->atom6name << "-"
02244 << ptr->atom7name << "-"
02245 << ptr->atom8name << ", USING NEW VALUES\n";
02246
02247 iout << endi;
02248
02249 ptr->dimension = new_node->dimension;
02250
02251 BigReal *tmpvalues = ptr->values;
02252 ptr->values = new_node->values;
02253 new_node->values = tmpvalues;
02254 }
02255
02256 delete new_node;
02257
02258 return;
02259 }
02260
02261 ptr=ptr->next;
02262 }
02263
02264
02265
02266
02267
02268
02269
02270 if ( (strcasecmp(new_node->atom1name, "X") == 0) ||
02271 (strcasecmp(new_node->atom2name, "X") == 0) ||
02272 (strcasecmp(new_node->atom3name, "X") == 0) ||
02273 (strcasecmp(new_node->atom4name, "X") == 0) ||
02274 (strcasecmp(new_node->atom5name, "X") == 0) ||
02275 (strcasecmp(new_node->atom6name, "X") == 0) ||
02276 (strcasecmp(new_node->atom7name, "X") == 0) ||
02277 (strcasecmp(new_node->atom8name, "X") == 0) )
02278 {
02279
02280 tail->next=new_node;
02281 tail=new_node;
02282
02283 return;
02284 }
02285 else
02286 {
02287
02288 new_node->next=crosstermp;
02289 crosstermp=new_node;
02290
02291 return;
02292 }
02293 }
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308 void Parameters::add_vdw_param(char *buf)
02309
02310 {
02311 char atomname[11];
02312 Real sigma;
02313 Real epsilon;
02314 Real sigma14;
02315 Real epsilon14;
02316 Real sqrt26;
02317 int read_count;
02318 struct vdw_params *new_node;
02319
02320
02321
02322 if (paramType == paraXplor)
02323 {
02324
02325 read_count=sscanf(buf, "%*s %s %f %f %f %f\n", atomname,
02326 &epsilon, &sigma, &epsilon14, &sigma14);
02327 }
02328 else if (paramType == paraCharmm)
02329 {
02330
02331 read_count=sscanf(buf, "%s %*f %f %f %*f %f %f\n", atomname,
02332 &epsilon, &sigma, &epsilon14, &sigma14);
02333 }
02334
02335
02336 if ((read_count != 5) && (paramType == paraXplor))
02337 {
02338 char err_msg[512];
02339
02340 sprintf(err_msg, "BAD vdW FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*\n", buf);
02341 NAMD_die(err_msg);
02342 }
02343 else if ( ((read_count != 5) && (read_count != 3)) && (paramType == paraCharmm))
02344 {
02345 char err_msg[512];
02346
02347 sprintf(err_msg, "BAD vdW FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*\n", buf);
02348 NAMD_die(err_msg);
02349 }
02350
02351 if (paramType == paraCharmm)
02352 {
02353
02354 epsilon*=-1.;
02355 sqrt26=pow(2.,(1./6.));
02356 sigma=2.*sigma/sqrt26;
02357
02358 if (read_count == 3)
02359 {
02360 epsilon14=epsilon;
02361 sigma14=sigma;
02362 }
02363 else
02364 {
02365 epsilon14*=-1.;
02366 sigma14=2.*sigma14/sqrt26;
02367 }
02368 }
02369
02370
02371 if ( epsilon < 0. || epsilon14 < 0. ) {
02372 iout << iWARN << "Ignoring VDW parameter with negative epsilon:\n"
02373 << buf << "\n" << endi;
02374 return;
02375 }
02376
02377
02378 new_node = new vdw_params;
02379
02380 if (new_node == NULL)
02381 {
02382 NAMD_die("memory allocation failed in Parameters::add_vdw_param");
02383 }
02384
02385
02386 strcpy(new_node->atomname, atomname);
02387 new_node->sigma = sigma;
02388 new_node->sigma14 = sigma14;
02389 new_node->epsilon = epsilon;
02390 new_node->epsilon14 = epsilon14;
02391
02392 new_node->left = NULL;
02393 new_node->right = NULL;
02394
02395
02396 vdwp=add_to_vdw_tree(new_node, vdwp);
02397
02398 return;
02399 }
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418 struct vdw_params *Parameters::add_to_vdw_tree(struct vdw_params *new_node,
02419 struct vdw_params *tree)
02420
02421 {
02422 int compare_code;
02423
02424
02425 if (tree == NULL)
02426 return(new_node);
02427
02428 compare_code = strcasecmp(new_node->atomname, tree->atomname);
02429
02430
02431 if (compare_code==0)
02432 {
02433
02434
02435
02436 if ((tree->sigma != new_node->sigma) ||
02437 (tree->epsilon != new_node->epsilon) ||
02438 (tree->sigma14 != new_node->sigma14) ||
02439 (tree->epsilon14 != new_node->epsilon14))
02440 {
02441 iout << iWARN << "DUPLICATE vdW ENTRY FOR " << tree->atomname
02442 << "\nPREVIOUS VALUES sigma=" << tree->sigma
02443 << " epsilon=" << tree->epsilon
02444 << " sigma14=" << tree->sigma14
02445 << " epsilon14=" << tree->epsilon14
02446 << "\n USING VALUES sigma=" << new_node->sigma
02447 << " epsilon=" << new_node->epsilon
02448 << " sigma14=" << new_node->sigma14
02449 << " epsilon14=" << new_node->epsilon14
02450 << "\n" << endi;
02451
02452 tree->sigma=new_node->sigma;
02453 tree->epsilon=new_node->epsilon;
02454 tree->sigma14=new_node->sigma14;
02455 tree->epsilon14=new_node->epsilon14;
02456 }
02457
02458 delete new_node;
02459
02460 return(tree);
02461 }
02462
02463
02464
02465
02466 if (compare_code < 0)
02467 {
02468 tree->left = add_to_vdw_tree(new_node, tree->left);
02469 }
02470 else
02471 {
02472 tree->right = add_to_vdw_tree(new_node, tree->right);
02473 }
02474
02475 return(tree);
02476 }
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491 void Parameters::add_table_pair_param(char *buf)
02492
02493 {
02494 char atom1name[11];
02495 char atom2name[11];
02496 char tabletype[11];
02497 int K;
02498 int read_count;
02499 struct table_pair_params *new_node;
02500
02501
02502 read_count=sscanf(buf, "%s %s %s\n", atom1name,
02503 atom2name, tabletype);
02504
02505
02506 if ((read_count != 3))
02507 {
02508 char err_msg[512];
02509
02510 sprintf(err_msg, "BAD TABLE PAIR FORMAT IN PARAMETER FILE\nLINE=*%s*", buf);
02511 NAMD_die(err_msg);
02512 }
02513
02514
02515 new_node = new table_pair_params;
02516
02517 if (new_node == NULL)
02518 {
02519 NAMD_die("memory allocation failed in Parameters::add_table_pair_param\n");
02520 }
02521
02522 strcpy(new_node->atom1name, atom1name);
02523 strcpy(new_node->atom2name, atom2name);
02524
02525
02526 K = get_int_table_type(tabletype);
02527
02528 if (K < 0) {
02529 char err_msg[512];
02530 sprintf(err_msg, "Couldn't find table parameters for table interaction %s!\n", tabletype);
02531 NAMD_die(err_msg);
02532 }
02533
02534
02535 new_node->K = K;
02536
02537 new_node->next = NULL;
02538
02539
02540
02541 add_to_table_pair_list(new_node);
02542
02543 return;
02544 }
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559 void Parameters::add_vdw_pair_param(char *buf)
02560
02561 {
02562 char atom1name[11];
02563 char atom2name[11];
02564 Real A;
02565 Real B;
02566 Real A14;
02567 Real B14;
02568 Real sqrt26;
02569 Real expo;
02570 int read_count;
02571 struct vdw_pair_params *new_node;
02572
02573
02574 if (paramType == paraXplor)
02575 {
02576
02577 read_count=sscanf(buf, "%*s %s %s %f %f %f %f\n", atom1name,
02578 atom2name, &A, &B, &A14, &B14);
02579 }
02580 else if (paramType == paraCharmm)
02581 {
02582
02583
02584 read_count=sscanf(buf, "%s %s %f %f\n", atom1name,
02585 atom2name, &A, &B);
02586
02587
02588
02589
02590
02591
02592 A14=-A;
02593 sqrt26=pow(2.,(1./6.));
02594 B14=B/sqrt26;
02595 expo=12.;
02596 A=pow(B14,expo);
02597 A=A*4.*A14;
02598 expo=6.;
02599 B=pow(B14,expo);
02600 B=B*4.*A14;
02601 A14=A;
02602 B14=B;
02603 }
02604
02605
02606 if ((read_count != 6) && (paramType == paraXplor))
02607 {
02608 char err_msg[512];
02609
02610 sprintf(err_msg, "BAD vdW PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
02611 NAMD_die(err_msg);
02612 }
02613 if ((read_count != 4) && (paramType == paraCharmm))
02614 {
02615 char err_msg[512];
02616
02617 sprintf(err_msg, "BAD vdW PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
02618 NAMD_die(err_msg);
02619 }
02620
02621
02622
02623 new_node = new vdw_pair_params;
02624
02625 if (new_node == NULL)
02626 {
02627 NAMD_die("memory allocation failed in Parameters::add_vdw_pair_param\n");
02628 }
02629
02630 strcpy(new_node->atom1name, atom1name);
02631 strcpy(new_node->atom2name, atom2name);
02632
02633
02634 new_node->A = A;
02635 new_node->A14 = A14;
02636 new_node->B = B;
02637 new_node->B14 = B14;
02638
02639 new_node->next = NULL;
02640
02641
02642 add_to_vdw_pair_list(new_node);
02643
02644 return;
02645 }
02646
02647
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660 void Parameters::add_nbthole_pair_param(char *buf)
02661
02662 {
02663 char atom1name[11];
02664 char atom2name[11];
02665 Real tholeij;
02666 int read_count;
02667 struct nbthole_pair_params *new_node;
02668
02669
02670 if (paramType == paraCharmm)
02671 {
02672
02673 read_count=sscanf(buf, "%s %s %f\n", atom1name,
02674 atom2name, &tholeij);
02675 }
02676
02677
02678 if ((read_count != 3) && (paramType == paraCharmm))
02679 {
02680 char err_msg[512];
02681
02682 sprintf(err_msg, "BAD NBTHOLE PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
02683 NAMD_die(err_msg);
02684 }
02685
02686
02687
02688 new_node = new nbthole_pair_params;
02689
02690 if (new_node == NULL)
02691 {
02692 NAMD_die("memory allocation failed in Parameters::nbthole_pair_param\n");
02693 }
02694
02695 strcpy(new_node->atom1name, atom1name);
02696 strcpy(new_node->atom2name, atom2name);
02697
02698
02699 new_node->tholeij = tholeij;
02700
02701 new_node->next = NULL;
02702
02703
02704 add_to_nbthole_pair_list(new_node);
02705
02706 return;
02707 }
02708
02709
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722 void Parameters::add_hb_pair_param(char *buf)
02723
02724 {
02725 #if 0
02726 char a1n[11];
02727 char a2n[11];
02728 Real A, B;
02729
02730
02732
02733 if (paramType == paraXplor) {
02734 if (sscanf(buf, "%*s %s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
02735 char err_msg[512];
02736 sprintf(err_msg, "BAD HBOND PAIR FORMAT IN XPLOR PARAMETER FILE\nLINE=*%s*", buf);
02737 NAMD_die(err_msg);
02738 }
02739 }
02740 else if (paramType == paraCharmm) {
02741 if (sscanf(buf, "%s %s %f %f\n", a1n, a2n, &A, &B) != 4) {
02742 char err_msg[512];
02743 sprintf(err_msg, "BAD HBOND PAIR FORMAT IN CHARMM PARAMETER FILE\nLINE=*%s*", buf);
02744 NAMD_die(err_msg);
02745 }
02746 }
02747
02748
02749
02750 if (hbondParams.add_hbond_pair(a1n, a2n, A, B) == FALSE) {
02751 iout << "\n" << iWARN << "Duplicate HBOND parameters for types " << a1n
02752 << " and " << a2n << " found; using latest values." << "\n" << endi;
02753 }
02754 #endif
02755 }
02756
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766
02767
02768
02769 void Parameters::add_to_table_pair_list(struct table_pair_params *new_node)
02770
02771 {
02772 static struct table_pair_params *tail=NULL;
02773 struct table_pair_params *ptr;
02774 int compare_code;
02775
02776
02777
02778 if (table_pairp == NULL)
02779 {
02780 table_pairp = new_node;
02781 tail = new_node;
02782 return;
02783 }
02784
02785 ptr = table_pairp;
02786
02787
02788 while (ptr!=NULL)
02789 {
02790
02791 compare_code = strncasecmp(new_node->atom1name, ptr->atom1name, 5);
02792
02793 if (compare_code == 0)
02794 {
02795
02796 compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
02797
02798 if (compare_code==0)
02799 {
02800
02801
02802
02803 iout << iWARN << "DUPLICATE TABLE PAIR ENTRY FOR "
02804 << new_node->atom1name << "-"
02805 << new_node->atom2name
02806 << "\n" << endi;
02807
02808 ptr->K=new_node->K;
02809
02810 delete new_node;
02811
02812 return;
02813 }
02814 }
02815
02816 ptr = ptr->next;
02817 }
02818
02819
02820
02821 tail->next = new_node;
02822 tail = new_node;
02823 }
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835
02836
02837 void Parameters::add_to_vdw_pair_list(struct vdw_pair_params *new_node)
02838
02839 {
02840 static struct vdw_pair_params *tail=NULL;
02841 struct vdw_pair_params *ptr;
02842 int compare_code;
02843
02844
02845
02846 if (vdw_pairp == NULL)
02847 {
02848 vdw_pairp = new_node;
02849 tail = new_node;
02850 return;
02851 }
02852
02853 ptr = vdw_pairp;
02854
02855
02856 while (ptr!=NULL)
02857 {
02858
02859 compare_code = strcasecmp(new_node->atom1name, ptr->atom1name);
02860
02861 if (compare_code == 0)
02862 {
02863
02864 compare_code = strcasecmp(new_node->atom2name, ptr->atom2name);
02865
02866 if (compare_code==0)
02867 {
02868
02869
02870
02871 if ((ptr->A != new_node->A) ||
02872 (ptr->B != new_node->B) ||
02873 (ptr->A14 != new_node->A14) ||
02874 (ptr->B14 != new_node->B14))
02875 {
02876 iout << iWARN << "DUPLICATE vdW PAIR ENTRY FOR "
02877 << new_node->atom1name << "-"
02878 << new_node->atom2name
02879 << "\nPREVIOUS VALUES A=" << ptr->A
02880 << " B=" << ptr->B
02881 << " A14=" << ptr->A14
02882 << " B14" << ptr->B14
02883 << "\n USING VALUES A=" << new_node->A
02884 << " B=" << new_node->B
02885 << " A14=" << new_node->A14
02886 << " B14" << new_node->B14
02887 << "\n" << endi;
02888
02889 ptr->A=new_node->A;
02890 ptr->B=new_node->B;
02891 ptr->A14=new_node->A14;
02892 ptr->B14=new_node->B14;
02893 }
02894
02895 delete new_node;
02896
02897 return;
02898 }
02899 }
02900
02901 ptr = ptr->next;
02902 }
02903
02904
02905
02906 tail->next = new_node;
02907 tail = new_node;
02908 }
02909
02910
02911
02912
02913
02914
02915
02916
02917
02918
02919
02920
02921
02922 void Parameters::add_to_nbthole_pair_list(struct nbthole_pair_params *new_node)
02923
02924 {
02925 static struct nbthole_pair_params *tail=NULL;
02926 struct nbthole_pair_params *ptr;
02927 int compare_code;
02928
02929
02930
02931 if (nbthole_pairp == NULL)
02932 {
02933 nbthole_pairp = new_node;
02934 tail = new_node;
02935 return;
02936 }
02937
02938 ptr = nbthole_pairp;
02939
02940 tail->next = new_node;
02941 tail = new_node;
02942 }
02943
02944
02945
02946
02947
02948
02949
02950
02951
02952
02953
02954
02955
02956
02957 void Parameters::done_reading_files()
02958
02959 {
02960 AllFilesRead = TRUE;
02961
02962
02963 if (NumBondParams)
02964 {
02965 bond_array = new BondValue[NumBondParams];
02966
02967 if (bond_array == NULL)
02968 {
02969 NAMD_die("memory allocation of bond_array failed!");
02970 }
02971 }
02972
02973 if (NumAngleParams)
02974 {
02975 angle_array = new AngleValue[NumAngleParams];
02976
02977 if (angle_array == NULL)
02978 {
02979 NAMD_die("memory allocation of angle_array failed!");
02980 }
02981 }
02982
02983 if (NumDihedralParams)
02984 {
02985 dihedral_array = new DihedralValue[NumDihedralParams];
02986
02987 if (dihedral_array == NULL)
02988 {
02989 NAMD_die("memory allocation of dihedral_array failed!");
02990 }
02991 memset(dihedral_array, 0, NumDihedralParams*sizeof(DihedralValue));
02992 }
02993
02994 if (NumImproperParams)
02995 {
02996 improper_array = new ImproperValue[NumImproperParams];
02997
02998 if (improper_array == NULL)
02999 {
03000 NAMD_die("memory allocation of improper_array failed!");
03001 }
03002 memset(improper_array, 0, NumImproperParams*sizeof(ImproperValue));
03003 }
03004
03005 if (NumCrosstermParams)
03006 {
03007 crossterm_array = new CrosstermValue[NumCrosstermParams];
03008 memset(crossterm_array, 0, NumCrosstermParams*sizeof(CrosstermValue));
03009 }
03010
03011 if (NumVdwParams)
03012 {
03013 atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
03014 vdw_array = new VdwValue[NumVdwParams];
03015
03016 if (vdw_array == NULL)
03017 {
03018 NAMD_die("memory allocation of vdw_array failed!");
03019 }
03020 }
03021 if (NumNbtholePairParams)
03022 {
03023 nbthole_array = new NbtholePairValue[NumNbtholePairParams];
03024
03025 if(nbthole_array == NULL)
03026 {
03027 NAMD_die("memory allocation of nbthole_array failed!");
03028 }
03029 }
03030
03031
03032
03033 index_bonds(bondp, 0);
03034 index_angles(anglep, 0);
03035 NumVdwParamsAssigned = index_vdw(vdwp, 0);
03036 index_dihedrals();
03037 index_impropers();
03038 index_crossterms();
03039 convert_nbthole_pairs();
03040
03041 convert_vdw_pairs();
03042 convert_table_pairs();
03043 }
03044
03045
03046
03047
03048
03049
03050
03051
03052
03053
03054
03055
03056
03057
03058
03059
03060
03061 Index Parameters::index_bonds(struct bond_params *tree, Index index)
03062
03063 {
03064
03065 if (tree==NULL)
03066 return(index);
03067
03068
03069 if (tree->left != NULL)
03070 {
03071 index=index_bonds(tree->left, index);
03072 }
03073
03074
03075 tree->index = index;
03076 bond_array[index].k = tree->forceconstant;
03077 bond_array[index].x0 = tree->distance;
03078 index++;
03079
03080
03081 if (tree->right != NULL)
03082 {
03083 index=index_bonds(tree->right, index);
03084 }
03085
03086 return(index);
03087 }
03088
03089
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103
03104
03105 Index Parameters::index_angles(struct angle_params *tree, Index index)
03106
03107 {
03108
03109 if (tree==NULL)
03110 return(index);
03111
03112
03113 if (tree->left != NULL)
03114 {
03115 index=index_angles(tree->left, index);
03116 }
03117
03118
03119 tree->index = index;
03120
03121 angle_array[index].k = tree->forceconstant;
03122 angle_array[index].k_ub = tree->k_ub;
03123 angle_array[index].r_ub = tree->r_ub;
03124 angle_array[index].normal = tree->normal;
03125
03126
03127 angle_array[index].theta0 = (tree->angle*PI)/180.0;
03128 index++;
03129
03130
03131 if (tree->right != NULL)
03132 {
03133 index=index_angles(tree->right, index);
03134 }
03135
03136 return(index);
03137 }
03138
03139
03140
03141
03142
03143
03144
03145
03146
03147
03148
03149
03150 void Parameters::index_dihedrals()
03151
03152 {
03153 struct dihedral_params *ptr;
03154 Index index=0;
03155 int i;
03156
03157
03158
03159
03160
03161
03162
03163
03164
03165
03166
03167
03168
03169
03170
03171
03172
03173 maxDihedralMults = new int[NumDihedralParams];
03174
03175 if (maxDihedralMults == NULL)
03176 {
03177 NAMD_die("memory allocation failed in Parameters::index_dihedrals()");
03178 }
03179
03180
03181 ptr = dihedralp;
03182
03183 while (ptr != NULL)
03184 {
03185
03186
03187
03188 maxDihedralMults[index] = ptr->multiplicity;
03189
03190
03191
03192 if (paramType == paraXplor)
03193 {
03194
03195
03196 dihedral_array[index].multiplicity = -1;
03197 }
03198 else if (paramType == paraCharmm)
03199 {
03200
03201
03202
03203 dihedral_array[index].multiplicity = ptr->multiplicity;
03204 }
03205
03206
03207 for (i=0; i<ptr->multiplicity; i++)
03208 {
03209 dihedral_array[index].values[i].k = ptr->values[i].k;
03210 dihedral_array[index].values[i].n = ptr->values[i].n;
03211
03212
03213 dihedral_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
03214 }
03215
03216 ptr->index = index;
03217
03218 index++;
03219 ptr=ptr->next;
03220 }
03221 }
03222
03223
03224
03225
03226
03227
03228
03229
03230
03231
03232
03233
03234 void Parameters::index_impropers()
03235
03236 {
03237 struct improper_params *ptr;
03238 Index index=0;
03239 int i;
03240
03241
03242
03243
03244
03245
03246
03247
03248
03249
03250
03251
03252
03253
03254
03255
03256 maxImproperMults = new int[NumImproperParams];
03257
03258 if (maxImproperMults == NULL)
03259 {
03260 NAMD_die("memory allocation failed in Parameters::index_impropers()");
03261 }
03262
03263
03264 ptr = improperp;
03265
03266 while (ptr != NULL)
03267 {
03268
03269
03270
03271 maxImproperMults[index] = ptr->multiplicity;
03272
03273
03274
03275 improper_array[index].multiplicity = -1;
03276
03277 for (i=0; i<ptr->multiplicity; i++)
03278 {
03279 improper_array[index].values[i].k = ptr->values[i].k;
03280 improper_array[index].values[i].n = ptr->values[i].n;
03281
03282
03283 improper_array[index].values[i].delta = ptr->values[i].delta*PI/180.0;
03284 }
03285
03286 ptr->index=index;
03287
03288 index++;
03289 ptr=ptr->next;
03290 }
03291 }
03292
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305 void crossterm_setup(CrosstermData *);
03306
03307 void Parameters::index_crossterms()
03308
03309 {
03310 struct crossterm_params *ptr;
03311 Index index=0;
03312 int i,j,k;
03313
03314
03315 ptr = crosstermp;
03316
03317 while (ptr != NULL)
03318 {
03319
03320
03321 int N = CrosstermValue::dim - 1;
03322
03323 if ( ptr->dimension != N ) {
03324 NAMD_die("Sorry, only CMAP dimension of 24 is supported");
03325 }
03326
03327 k = 0;
03328 for (i=0; i<N; i++) {
03329 for (j=0; j<N; j++) {
03330 crossterm_array[index].c[i][j].d00 = ptr->values[k];
03331 ++k;
03332 }
03333 }
03334 for (i=0; i<N; i++) {
03335 crossterm_array[index].c[i][N].d00 =
03336 crossterm_array[index].c[i][0].d00;
03337 crossterm_array[index].c[N][i].d00 =
03338 crossterm_array[index].c[0][i].d00;
03339 }
03340 crossterm_array[index].c[N][N].d00 =
03341 crossterm_array[index].c[0][0].d00;
03342
03343 crossterm_setup(&crossterm_array[index].c[0][0]);
03344
03345 ptr->index=index;
03346
03347 index++;
03348 ptr=ptr->next;
03349 }
03350 }
03351
03352
03353
03354
03355
03356
03357
03358
03359
03360
03361
03362
03363
03364
03365
03366
03367
03368 Index Parameters::index_vdw(struct vdw_params *tree, Index index)
03369
03370 {
03371
03372 if (tree==NULL)
03373 return(index);
03374
03375
03376 if (tree->left != NULL)
03377 {
03378 index=index_vdw(tree->left, index);
03379 }
03380
03381
03382 tree->index = index;
03383
03384 vdw_array[index].sigma = tree->sigma;
03385 vdw_array[index].epsilon = tree->epsilon;
03386 vdw_array[index].sigma14 = tree->sigma14;
03387 vdw_array[index].epsilon14 = tree->epsilon14;
03388
03389 char *nameloc = atom_type_name(index);
03390 strncpy(nameloc, tree->atomname, MAX_ATOMTYPE_CHARS);
03391 nameloc[MAX_ATOMTYPE_CHARS] = '\0';
03392
03393
03394
03395
03396 index++;
03397
03398
03399 if (tree->right != NULL)
03400 {
03401 index=index_vdw(tree->right, index);
03402 }
03403
03404 return(index);
03405 }
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416
03417
03418
03419
03420
03421
03422
03423
03424
03425
03426 #ifdef MEM_OPT_VERSION
03427 void Parameters::assign_vdw_index(char *atomtype, AtomCstInfo *atom_ptr)
03428 #else
03429 void Parameters::assign_vdw_index(char *atomtype, Atom *atom_ptr)
03430 #endif
03431 {
03432 struct vdw_params *ptr;
03433 int found=0;
03434 int comp_code;
03435
03436
03437 if (!AllFilesRead)
03438 {
03439 NAMD_die("Tried to assign vdw index before all parameter files were read");
03440 }
03441
03442
03443 ptr=vdwp;
03444
03445
03446
03447
03448
03449 while (!found && (ptr!=NULL))
03450 {
03451 comp_code = strcasecmp(atomtype, ptr->atomname);
03452
03453 if (comp_code == 0)
03454 {
03455
03456 atom_ptr->vdw_type=ptr->index;
03457 found=1;
03458 }
03459 else if (comp_code < 0)
03460 {
03461
03462 ptr=ptr->left;
03463 }
03464 else
03465 {
03466
03467 ptr=ptr->right;
03468 }
03469 }
03470
03471
03472 if (!found)
03473 {
03474
03475
03476
03477 size_t windx;
03478
03479
03480 ptr=vdwp;
03481
03482 while (!found && (ptr!=NULL))
03483 {
03484
03485
03486 windx= strcspn(ptr->atomname,"*");
03487 if (windx == strlen(ptr->atomname))
03488 {
03489
03490 comp_code = strcasecmp(atomtype, ptr->atomname);
03491 }
03492 else
03493 {
03494 comp_code = strncasecmp(atomtype, ptr->atomname, windx);
03495 }
03496
03497 if (comp_code == 0)
03498 {
03499
03500 atom_ptr->vdw_type=ptr->index;
03501 found=1;
03502 char errbuf[100];
03503 sprintf(errbuf,"VDW TYPE NAME %s MATCHES PARAMETER TYPE NAME %s",
03504 atomtype, ptr->atomname);
03505 int i;
03506 for(i=0; i<error_msgs.size(); i++) {
03507 if ( strcmp(errbuf,error_msgs[i]) == 0 ) break;
03508 }
03509 if ( i == error_msgs.size() ) {
03510 char *newbuf = new char[strlen(errbuf)+1];
03511 strcpy(newbuf,errbuf);
03512 error_msgs.add(newbuf);
03513 iout << iWARN << newbuf << "\n" << endi;
03514 }
03515 }
03516 else if (comp_code < 0)
03517 {
03518
03519 ptr=ptr->left;
03520 }
03521 else
03522 {
03523
03524 ptr=ptr->right;
03525 }
03526
03527 }
03528
03529 }
03530
03531
03532
03533 if (!found)
03534 {
03535 char err_msg[100];
03536
03537 sprintf(err_msg, "DIDN'T FIND vdW PARAMETER FOR ATOM TYPE %s",
03538 atomtype);
03539 NAMD_die(err_msg);
03540 }
03541
03542 return;
03543 }
03544
03545
03546
03547
03548
03549
03550
03551
03552
03553
03554
03555
03556
03557
03558
03559
03560
03561
03562
03563
03564 int Parameters::get_table_pair_params(Index ind1, Index ind2, int *K) {
03565 IndexedTablePair *ptr;
03566 Index temp;
03567 int found=FALSE;
03568
03569 ptr=tab_pair_tree;
03570
03571
03572
03573 if (ind1 > ind2)
03574 {
03575 temp = ind1;
03576 ind1 = ind2;
03577 ind2 = temp;
03578 }
03579
03580
03581
03582 while (!found && (ptr!=NULL))
03583 {
03584
03585 if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
03586 {
03587 found = TRUE;
03588 }
03589 else if ( (ind1 < ptr->ind1) ||
03590 ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
03591 {
03592
03593 ptr=ptr->left;
03594 }
03595 else
03596 {
03597
03598 ptr=ptr->right;
03599 }
03600 }
03601
03602
03603 if (found)
03604 {
03605 *K = ptr->K;
03606 return(TRUE);
03607 }
03608 else
03609 {
03610 return(FALSE);
03611 }
03612 }
03613
03614
03615
03616
03617
03618
03619
03620
03621
03622
03623
03624
03625
03626
03627
03628
03629
03630
03631
03632
03633
03634
03635
03636
03637
03638
03639 int Parameters::get_vdw_pair_params(Index ind1, Index ind2, Real *A,
03640 Real *B, Real *A14, Real *B14)
03641
03642 {
03643 IndexedVdwPair *ptr;
03644 Index temp;
03645
03646 int found=FALSE;
03647
03648 ptr=vdw_pair_tree;
03649
03650
03651
03652 if (ind1 > ind2)
03653 {
03654 temp = ind1;
03655 ind1 = ind2;
03656 ind2 = temp;
03657 }
03658
03659
03660
03661 while (!found && (ptr!=NULL))
03662 {
03663 if ( (ind1 == ptr->ind1) && (ind2 == ptr->ind2) )
03664 {
03665 found = TRUE;
03666 }
03667 else if ( (ind1 < ptr->ind1) ||
03668 ( (ind1==ptr->ind1) && (ind2 < ptr->ind2) ) )
03669 {
03670
03671 ptr=ptr->left;
03672 }
03673 else
03674 {
03675
03676 ptr=ptr->right;
03677 }
03678 }
03679
03680
03681 if (found)
03682 {
03683 *A = ptr->A;
03684 *B = ptr->B;
03685 *A14 = ptr->A14;
03686 *B14 = ptr->B14;
03687
03688 return(TRUE);
03689 }
03690 else
03691 {
03692 return(FALSE);
03693 }
03694 }
03695
03696
03697
03698
03699
03700
03701
03702
03703
03704
03705
03706
03707
03708
03709
03710
03711
03712
03713
03714
03715
03716
03717 void Parameters::assign_bond_index(char *atom1, char *atom2, Bond *bond_ptr)
03718
03719 {
03720 struct bond_params *ptr;
03721 int found=0;
03722 int cmp_code;
03723 char tmp_name[15];
03724
03725
03726 if (!AllFilesRead)
03727 {
03728 NAMD_die("Tried to assign bond index before all parameter files were read");
03729 }
03730
03731
03732
03733 if (strcasecmp(atom1, atom2) > 0)
03734 {
03735 strcpy(tmp_name, atom1);
03736 strcpy(atom1, atom2);
03737 strcpy(atom2, tmp_name);
03738 }
03739
03740
03741 ptr=bondp;
03742
03743
03744
03745 while (!found && (ptr!=NULL))
03746 {
03747 cmp_code=strcasecmp(atom1, ptr->atom1name);
03748
03749 if (cmp_code == 0)
03750 {
03751 cmp_code=strcasecmp(atom2, ptr->atom2name);
03752 }
03753
03754 if (cmp_code == 0)
03755 {
03756
03757 found=1;
03758 bond_ptr->bond_type = ptr->index;
03759 }
03760 else if (cmp_code < 0)
03761 {
03762
03763 ptr=ptr->left;
03764 }
03765 else
03766 {
03767
03768 ptr=ptr->right;
03769 }
03770 }
03771
03772
03773 if (!found)
03774 {
03775 if ((strcmp(atom1, "DRUD")==0 || strcmp(atom2, "DRUD")==0)
03776 && (strcmp(atom1, "X")!=0 && strcmp(atom2, "X")!=0)) {
03777
03778 char a1[8] = "DRUD", a2[8] = "X";
03779 return assign_bond_index(a1, a2, bond_ptr);
03780 }
03781 else {
03782 char err_msg[128];
03783
03784 sprintf(err_msg, "UNABLE TO FIND BOND PARAMETERS FOR %s %s (ATOMS %i %i)",
03785 atom1, atom2, bond_ptr->atom1+1, bond_ptr->atom2+1);
03786 NAMD_die(err_msg);
03787 }
03788 }
03789
03790 return;
03791 }
03792
03793
03794
03795
03796
03797
03798
03799
03800
03801
03802
03803
03804
03805
03806
03807
03808
03809
03810
03811
03812
03813
03814 void Parameters::assign_angle_index(char *atom1, char *atom2, char*atom3,
03815 Angle *angle_ptr, int notFoundIndex)
03816
03817 {
03818 struct angle_params *ptr;
03819 int comp_val;
03820 int found=0;
03821 char tmp_name[15];
03822
03823
03824 if (!AllFilesRead)
03825 {
03826 NAMD_die("Tried to assign angle index before all parameter files were read");
03827 }
03828
03829
03830
03831 if (strcasecmp(atom1, atom3) > 0)
03832 {
03833 strcpy(tmp_name, atom1);
03834 strcpy(atom1, atom3);
03835 strcpy(atom3, tmp_name);
03836 }
03837
03838
03839 ptr=anglep;
03840
03841
03842
03843 while (!found && (ptr != NULL))
03844 {
03845 comp_val = strcasecmp(atom1, ptr->atom1name);
03846
03847 if (comp_val == 0)
03848 {
03849
03850 comp_val = strcasecmp(atom2, ptr->atom2name);
03851
03852 if (comp_val == 0)
03853 {
03854
03855 comp_val = strcasecmp(atom3, ptr->atom3name);
03856 }
03857 }
03858
03859 if (comp_val == 0)
03860 {
03861
03862 found = 1;
03863 angle_ptr->angle_type = ptr->index;
03864 }
03865 else if (comp_val < 0)
03866 {
03867
03868 ptr=ptr->left;
03869 }
03870 else
03871 {
03872
03873 ptr=ptr->right;
03874 }
03875 }
03876
03877
03878 if (!found)
03879 {
03880 char err_msg[128];
03881
03882 sprintf(err_msg, "UNABLE TO FIND ANGLE PARAMETERS FOR %s %s %s (ATOMS %i %i %i)",
03883 atom1, atom2, atom3, angle_ptr->atom1+1, angle_ptr->atom2+1, angle_ptr->atom3+1);
03884
03885 if ( notFoundIndex ) {
03886 angle_ptr->angle_type = notFoundIndex;
03887 iout << iWARN << err_msg << "\n" << endi;
03888 return;
03889 } else NAMD_die(err_msg);
03890 }
03891
03892 return;
03893 }
03894
03895
03896
03897
03898
03899
03900
03901
03902
03903
03904
03905
03906
03907
03908
03909
03910
03911
03912
03913
03914
03915
03916
03917 void Parameters::assign_dihedral_index(char *atom1, char *atom2, char *atom3,
03918 char *atom4, Dihedral *dihedral_ptr,
03919 int multiplicity, int notFoundIndex)
03920
03921 {
03922 struct dihedral_params *ptr;
03923 int found=0;
03924
03925
03926 ptr=dihedralp;
03927
03928
03929
03930 while (!found && (ptr!=NULL))
03931 {
03932
03933
03934
03935
03936
03937
03938
03939
03940 if ( ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom1)==0) ) &&
03941 ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom2)==0) ) &&
03942 ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom3)==0) ) &&
03943 ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom4)==0) ) )
03944 {
03945
03946 found=1;
03947 }
03948 else if ( ( ptr->atom4wild || (strcasecmp(ptr->atom4name, atom1)==0) ) &&
03949 ( ptr->atom3wild || (strcasecmp(ptr->atom3name, atom2)==0) ) &&
03950 ( ptr->atom2wild || (strcasecmp(ptr->atom2name, atom3)==0) ) &&
03951 ( ptr->atom1wild || (strcasecmp(ptr->atom1name, atom4)==0) ) )
03952 {
03953
03954 found=1;
03955 }
03956 else
03957 {
03958
03959 ptr=ptr->next;
03960 }
03961 }
03962
03963
03964 if (!found)
03965 {
03966 char err_msg[128];
03967
03968 sprintf(err_msg, "UNABLE TO FIND DIHEDRAL PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
03969 atom1, atom2, atom3, atom4, dihedral_ptr->atom1+1, dihedral_ptr->atom2+1,
03970 dihedral_ptr->atom3+1, dihedral_ptr->atom4+1);
03971
03972 if ( notFoundIndex ) {
03973 dihedral_ptr->dihedral_type = notFoundIndex;
03974 iout << iWARN << err_msg << "\n" << endi;
03975 return;
03976 } else NAMD_die(err_msg);
03977 }
03978
03979
03980
03981
03982 if (multiplicity > maxDihedralMults[ptr->index])
03983 {
03984 char err_msg[257];
03985
03986 sprintf(err_msg, "Multiplicity of Paramters for diehedral bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxDihedralMults[ptr->index]);
03987 NAMD_die(err_msg);
03988 }
03989
03990
03991
03992 if (multiplicity > dihedral_array[ptr->index].multiplicity)
03993 {
03994 dihedral_array[ptr->index].multiplicity = multiplicity;
03995 }
03996
03997 dihedral_ptr->dihedral_type = ptr->index;
03998
03999 return;
04000 }
04001
04002
04003
04004
04005
04006
04007
04008
04009
04010
04011
04012
04013
04014
04015
04016
04017
04018
04019
04020
04021
04022
04023
04024 void Parameters::assign_improper_index(char *atom1, char *atom2, char *atom3,
04025 char *atom4, Improper *improper_ptr,
04026 int multiplicity)
04027
04028 {
04029 struct improper_params *ptr;
04030 int found=0;
04031
04032
04033 ptr=improperp;
04034
04035
04036
04037 while (!found && (ptr!=NULL))
04038 {
04039
04040
04041
04042
04043
04044
04045
04046
04047 if ( ( (strcasecmp(ptr->atom1name, atom1)==0) ||
04048 (strcasecmp(ptr->atom1name, "X")==0) ) &&
04049 ( (strcasecmp(ptr->atom2name, atom2)==0) ||
04050 (strcasecmp(ptr->atom2name, "X")==0) ) &&
04051 ( (strcasecmp(ptr->atom3name, atom3)==0) ||
04052 (strcasecmp(ptr->atom3name, "X")==0) ) &&
04053 ( (strcasecmp(ptr->atom4name, atom4)==0) ||
04054 (strcasecmp(ptr->atom4name, "X")==0) ) )
04055 {
04056
04057 found=1;
04058 }
04059 else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) ||
04060 (strcasecmp(ptr->atom4name, "X")==0) ) &&
04061 ( (strcasecmp(ptr->atom3name, atom2)==0) ||
04062 (strcasecmp(ptr->atom3name, "X")==0) ) &&
04063 ( (strcasecmp(ptr->atom2name, atom3)==0) ||
04064 (strcasecmp(ptr->atom2name, "X")==0) ) &&
04065 ( (strcasecmp(ptr->atom1name, atom4)==0) ||
04066 (strcasecmp(ptr->atom1name, "X")==0) ) )
04067 {
04068
04069 found=1;
04070 }
04071 else
04072 {
04073
04074 ptr=ptr->next;
04075 }
04076 }
04077
04078
04079 if (!found)
04080 {
04081 char err_msg[128];
04082
04083 sprintf(err_msg, "UNABLE TO FIND IMPROPER PARAMETERS FOR %s %s %s %s (ATOMS %i %i %i %i)",
04084 atom1, atom2, atom3, atom4, improper_ptr->atom1+1, improper_ptr->atom2+1,
04085 improper_ptr->atom3+1, improper_ptr->atom4+1);
04086
04087 NAMD_die(err_msg);
04088 }
04089
04090
04091
04092
04093 if (multiplicity > maxImproperMults[ptr->index])
04094 {
04095 char err_msg[257];
04096
04097 sprintf(err_msg, "Multiplicity of Paramters for improper bond %s %s %s %s of %d exceeded", atom1, atom2, atom3, atom4, maxImproperMults[ptr->index]);
04098 NAMD_die(err_msg);
04099 }
04100
04101
04102
04103 if (multiplicity > improper_array[ptr->index].multiplicity)
04104 {
04105 improper_array[ptr->index].multiplicity = multiplicity;
04106 }
04107
04108
04109 improper_ptr->improper_type = ptr->index;
04110
04111 return;
04112 }
04113
04114
04115
04116
04117
04118
04119
04120
04121 void Parameters::assign_crossterm_index(char *atom1, char *atom2, char *atom3,
04122 char *atom4, char *atom5, char *atom6, char *atom7,
04123 char *atom8, Crossterm *crossterm_ptr)
04124 {
04125 struct crossterm_params *ptr;
04126 int found=0;
04127
04128
04129 ptr=crosstermp;
04130
04131
04132
04133 while (!found && (ptr!=NULL))
04134 {
04135
04136
04137
04138
04139
04140
04141
04142
04143 if ( ( (strcasecmp(ptr->atom1name, atom1)==0) ||
04144 (strcasecmp(ptr->atom1name, "X")==0) ) &&
04145 ( (strcasecmp(ptr->atom2name, atom2)==0) ||
04146 (strcasecmp(ptr->atom2name, "X")==0) ) &&
04147 ( (strcasecmp(ptr->atom3name, atom3)==0) ||
04148 (strcasecmp(ptr->atom3name, "X")==0) ) &&
04149 ( (strcasecmp(ptr->atom4name, atom4)==0) ||
04150 (strcasecmp(ptr->atom4name, "X")==0) ) )
04151 {
04152
04153 found=1;
04154 }
04155 else if ( ( (strcasecmp(ptr->atom4name, atom1)==0) ||
04156 (strcasecmp(ptr->atom4name, "X")==0) ) &&
04157 ( (strcasecmp(ptr->atom3name, atom2)==0) ||
04158 (strcasecmp(ptr->atom3name, "X")==0) ) &&
04159 ( (strcasecmp(ptr->atom2name, atom3)==0) ||
04160 (strcasecmp(ptr->atom2name, "X")==0) ) &&
04161 ( (strcasecmp(ptr->atom1name, atom4)==0) ||
04162 (strcasecmp(ptr->atom1name, "X")==0) ) )
04163 {
04164
04165 found=1;
04166 }
04167 if ( ! found ) {
04168
04169 ptr=ptr->next;
04170 continue;
04171 }
04172 found = 0;
04173 if ( ( (strcasecmp(ptr->atom5name, atom5)==0) ||
04174 (strcasecmp(ptr->atom5name, "X")==0) ) &&
04175 ( (strcasecmp(ptr->atom6name, atom6)==0) ||
04176 (strcasecmp(ptr->atom6name, "X")==0) ) &&
04177 ( (strcasecmp(ptr->atom7name, atom7)==0) ||
04178 (strcasecmp(ptr->atom7name, "X")==0) ) &&
04179 ( (strcasecmp(ptr->atom8name, atom8)==0) ||
04180 (strcasecmp(ptr->atom8name, "X")==0) ) )
04181 {
04182
04183 found=1;
04184 }
04185 else if ( ( (strcasecmp(ptr->atom8name, atom5)==0) ||
04186 (strcasecmp(ptr->atom8name, "X")==0) ) &&
04187 ( (strcasecmp(ptr->atom7name, atom6)==0) ||
04188 (strcasecmp(ptr->atom7name, "X")==0) ) &&
04189 ( (strcasecmp(ptr->atom6name, atom7)==0) ||
04190 (strcasecmp(ptr->atom6name, "X")==0) ) &&
04191 ( (strcasecmp(ptr->atom5name, atom8)==0) ||
04192 (strcasecmp(ptr->atom5name, "X")==0) ) )
04193 {
04194
04195 found=1;
04196 }
04197 if ( ! found ) {
04198
04199 ptr=ptr->next;
04200 }
04201 }
04202
04203
04204 if (!found)
04205 {
04206 char err_msg[128];
04207
04208 sprintf(err_msg, "UNABLE TO FIND CROSSTERM PARAMETERS FOR %s %s %s %s %s %s %s %s\n"
04209 "(ATOMS %i %i %i %i %i %i %i %i)",
04210 atom1, atom2, atom3, atom4, atom5, atom6, atom7, atom8,
04211 crossterm_ptr->atom1+1, crossterm_ptr->atom1+2, crossterm_ptr->atom1+3, crossterm_ptr->atom4+1,
04212 crossterm_ptr->atom1+5, crossterm_ptr->atom1+6, crossterm_ptr->atom1+7, crossterm_ptr->atom8+1);
04213
04214 NAMD_die(err_msg);
04215 }
04216
04217
04218 crossterm_ptr->crossterm_type = ptr->index;
04219
04220 return;
04221 }
04222
04223
04224
04225
04226
04227
04228
04229
04230
04231
04232
04233
04234
04235
04236
04237
04238 void Parameters::free_bond_tree(struct bond_params *bond_ptr)
04239
04240 {
04241 if (bond_ptr->left != NULL)
04242 {
04243 free_bond_tree(bond_ptr->left);
04244 }
04245
04246 if (bond_ptr->right != NULL)
04247 {
04248 free_bond_tree(bond_ptr->right);
04249 }
04250
04251 delete bond_ptr;
04252
04253 return;
04254 }
04255
04256
04257
04258
04259
04260
04261
04262
04263
04264
04265
04266
04267
04268
04269
04270
04271 void Parameters::free_angle_tree(struct angle_params *angle_ptr)
04272
04273 {
04274 if (angle_ptr->left != NULL)
04275 {
04276 free_angle_tree(angle_ptr->left);
04277 }
04278
04279 if (angle_ptr->right != NULL)
04280 {
04281 free_angle_tree(angle_ptr->right);
04282 }
04283
04284 delete angle_ptr;
04285
04286 return;
04287 }
04288
04289
04290
04291
04292
04293
04294
04295
04296
04297
04298
04299
04300
04301
04302 void Parameters::free_dihedral_list(struct dihedral_params *dih_ptr)
04303
04304 {
04305 struct dihedral_params *ptr;
04306 struct dihedral_params *next;
04307
04308 ptr=dih_ptr;
04309
04310 while (ptr != NULL)
04311 {
04312 next=ptr->next;
04313 delete ptr;
04314 ptr=next;
04315 }
04316
04317 return;
04318 }
04319
04320
04321
04322
04323
04324
04325
04326
04327
04328
04329
04330
04331
04332
04333 void Parameters::free_improper_list(struct improper_params *imp_ptr)
04334
04335 {
04336 struct improper_params *ptr;
04337 struct improper_params *next;
04338
04339 ptr=imp_ptr;
04340
04341 while (ptr != NULL)
04342 {
04343 next=ptr->next;
04344 delete ptr;
04345 ptr=next;
04346 }
04347
04348 return;
04349 }
04350
04351
04352
04353
04354
04355
04356
04357
04358
04359
04360
04361
04362
04363
04364 void Parameters::free_crossterm_list(struct crossterm_params *imp_ptr)
04365
04366 {
04367 struct crossterm_params *ptr;
04368 struct crossterm_params *next;
04369
04370 ptr=imp_ptr;
04371
04372 while (ptr != NULL)
04373 {
04374 next=ptr->next;
04375 delete ptr;
04376 ptr=next;
04377 }
04378
04379 return;
04380 }
04381
04382
04383
04384
04385
04386
04387
04388
04389
04390
04391
04392
04393
04394
04395
04396
04397 void Parameters::free_vdw_tree(struct vdw_params *vdw_ptr)
04398
04399 {
04400 if (vdw_ptr->left != NULL)
04401 {
04402 free_vdw_tree(vdw_ptr->left);
04403 }
04404
04405 if (vdw_ptr->right != NULL)
04406 {
04407 free_vdw_tree(vdw_ptr->right);
04408 }
04409
04410 delete vdw_ptr;
04411
04412 return;
04413 }
04414
04415
04416
04417
04418
04419
04420
04421
04422
04423
04424 void Parameters::free_vdw_pair_list()
04425 {
04426 struct vdw_pair_params *ptr, *next;
04427
04428 ptr=vdw_pairp;
04429
04430 while (ptr != NULL)
04431 {
04432 next = ptr->next;
04433
04434 delete ptr;
04435
04436 ptr = next;
04437 }
04438
04439 vdw_pairp = NULL;
04440 }
04441
04442
04443
04444
04445
04446
04447
04448
04449
04450
04451 void Parameters::free_nbthole_pair_list()
04452 {
04453 struct nbthole_pair_params *ptr, *next;
04454
04455 ptr=nbthole_pairp;
04456
04457 while (ptr != NULL)
04458 {
04459 next = ptr->next;
04460
04461 delete ptr;
04462
04463 ptr = next;
04464 }
04465
04466 nbthole_pairp = NULL;
04467 }
04468
04469
04470
04471
04472
04473
04474
04475
04476 void Parameters::free_table_pair_tree(IndexedTablePair *table_pair_ptr) {
04477 if (table_pair_ptr->left != NULL)
04478 {
04479 free_table_pair_tree(table_pair_ptr->left);
04480 }
04481
04482 if (table_pair_ptr->right != NULL)
04483 {
04484 free_table_pair_tree(table_pair_ptr->right);
04485 }
04486
04487 delete table_pair_ptr;
04488
04489 return;
04490 }
04491
04492
04493
04494
04495
04496
04497
04498
04499
04500
04501
04502
04503
04504
04505
04506
04507 void Parameters::free_vdw_pair_tree(IndexedVdwPair *vdw_pair_ptr)
04508
04509 {
04510 if (vdw_pair_ptr->left != NULL)
04511 {
04512 free_vdw_pair_tree(vdw_pair_ptr->left);
04513 }
04514
04515 if (vdw_pair_ptr->right != NULL)
04516 {
04517 free_vdw_pair_tree(vdw_pair_ptr->right);
04518 }
04519
04520 delete vdw_pair_ptr;
04521
04522 return;
04523 }
04524
04525
04526
04527
04528
04529
04530
04531
04532
04533
04534
04535
04536
04537
04538
04539
04540 void Parameters::free_nbthole_pair_tree(IndexedNbtholePair *nbthole_pair_ptr)
04541
04542 {
04543 if (nbthole_pair_ptr->left != NULL)
04544 {
04545 free_nbthole_pair_tree(nbthole_pair_ptr->left);
04546 }
04547
04548 if (nbthole_pair_ptr->right != NULL)
04549 {
04550 free_nbthole_pair_tree(nbthole_pair_ptr->right);
04551 }
04552
04553 delete nbthole_pair_ptr;
04554
04555 return;
04556 }
04557
04558
04559
04560
04561
04562
04563
04564
04565
04566
04567
04568
04569
04570
04571
04572 void Parameters::traverse_bond_params(struct bond_params *tree)
04573
04574 {
04575 if (tree==NULL)
04576 return;
04577
04578 if (tree->left != NULL)
04579 {
04580 traverse_bond_params(tree->left);
04581 }
04582
04583 DebugM(3,"BOND " << tree->atom1name << " " << tree->atom2name \
04584 << " index=" << tree->index << " k=" << tree->forceconstant \
04585 << " x0=" << tree->distance);
04586
04587 if (tree->right != NULL)
04588 {
04589 traverse_bond_params(tree->right);
04590 }
04591 }
04592
04593
04594
04595
04596
04597
04598
04599
04600
04601
04602
04603
04604
04605
04606
04607 void Parameters::traverse_angle_params(struct angle_params *tree)
04608
04609 {
04610 if (tree==NULL)
04611 return;
04612
04613 if (tree->left != NULL)
04614 {
04615 traverse_angle_params(tree->left);
04616 }
04617 DebugM(3,"ANGLE " << tree->atom1name << " " << tree->atom2name \
04618 << " " << tree->atom3name << " index=" << tree->index \
04619 << " k=" << tree->forceconstant << " theta0=" << tree->angle \
04620 );
04621
04622 if (tree->right != NULL)
04623 {
04624 traverse_angle_params(tree->right);
04625 }
04626 }
04627
04628
04629
04630
04631
04632
04633
04634
04635
04636
04637
04638
04639
04640
04641
04642 void Parameters::traverse_dihedral_params(struct dihedral_params *list)
04643
04644 {
04645 int i;
04646
04647 while (list != NULL)
04648 {
04649 DebugM(3,"DIHEDRAL " << list->atom1name << " " \
04650 << list->atom2name << " " << list->atom3name \
04651 << " " << list->atom4name << " index=" \
04652 << list->index \
04653 << " multiplicity=" << list->multiplicity << "\n");
04654
04655 for (i=0; i<list->multiplicity; i++)
04656 {
04657 DebugM(3,"k=" << list->values[i].k \
04658 << " n=" << list->values[i].n \
04659 << " delta=" << list->values[i].delta);
04660 }
04661
04662 list=list->next;
04663 }
04664 }
04665
04666
04667
04668
04669
04670
04671
04672
04673
04674
04675
04676
04677
04678
04679
04680 void Parameters::traverse_improper_params(struct improper_params *list)
04681
04682 {
04683 int i;
04684
04685 while (list != NULL)
04686 {
04687 DebugM(3,"Improper " << list->atom1name << " " \
04688 << list->atom2name << " " << list->atom3name \
04689 << " " << list->atom4name << " index=" \
04690 << list->index \
04691 << " multiplicity=" << list->multiplicity << "\n");
04692
04693 for (i=0; i<list->multiplicity; i++)
04694 {
04695 DebugM(3,"k=" << list->values[i].k \
04696 << " n=" << list->values[i].n \
04697 << " delta=" << list->values[i].delta);
04698 }
04699
04700 list=list->next;
04701 }
04702 }
04703
04704
04705
04706
04707
04708
04709
04710
04711
04712
04713
04714
04715
04716
04717
04718
04719 void Parameters::traverse_vdw_params(struct vdw_params *tree)
04720
04721 {
04722 if (tree==NULL)
04723 return;
04724
04725 if (tree->left != NULL)
04726 {
04727 traverse_vdw_params(tree->left);
04728 }
04729
04730 DebugM(3,"vdW " << tree->atomname << " index=" << tree->index \
04731 << " sigma=" << tree->sigma << " epsilon=" << \
04732 tree->epsilon << " sigma 1-4=" << tree->sigma14 \
04733 << " epsilon 1-4=" << tree->epsilon14 << std::endl);
04734
04735 if (tree->right != NULL)
04736 {
04737 traverse_vdw_params(tree->right);
04738 }
04739 }
04740
04741
04742
04743
04744
04745
04746
04747
04748
04749
04750
04751
04752
04753
04754 void Parameters::traverse_vdw_pair_params(struct vdw_pair_params *list)
04755
04756 {
04757 if (list==NULL)
04758 return;
04759
04760 DebugM(3,"vdW PAIR " << list->atom1name << " " \
04761 << list->atom2name << " A=" << list->A \
04762 << " B=" << list->B << " A 1-4=" \
04763 << list->A14 << " B 1-4=" << list->B14 \
04764 );
04765
04766 traverse_vdw_pair_params(list->next);
04767 }
04768
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779
04780
04781 void Parameters::traverse_nbthole_pair_params(struct nbthole_pair_params *list)
04782
04783 {
04784 if (list==NULL)
04785 return;
04786
04787 DebugM(3,"NBTHOLE PAIR " << list->atom1name << " " \
04788 << list->atom2name << " tholeij =" << list->tholeij \
04789 );
04790
04791 traverse_nbthole_pair_params(list->next);
04792 }
04793
04794
04795
04796
04797
04798
04799
04800
04801
04802
04803
04804 void Parameters::print_bond_params()
04805 {
04806 DebugM(3,NumBondParams << " BOND PARAMETERS\n" \
04807 << "*****************************************" \
04808 );
04809
04810 traverse_bond_params(bondp);
04811 }
04812
04813
04814
04815
04816
04817
04818
04819
04820
04821
04822 void Parameters::print_angle_params()
04823 {
04824 DebugM(3,NumAngleParams << " ANGLE PARAMETERS\n"
04825 << "*****************************************" );
04826 traverse_angle_params(anglep);
04827 }
04828
04829
04830
04831
04832
04833
04834
04835
04836
04837
04838 void Parameters::print_dihedral_params()
04839 {
04840 DebugM(3,NumDihedralParams << " DIHEDRAL PARAMETERS\n" \
04841 << "*****************************************" );
04842
04843 traverse_dihedral_params(dihedralp);
04844 }
04845
04846
04847
04848
04849
04850
04851
04852
04853
04854
04855 void Parameters::print_improper_params()
04856 {
04857 DebugM(3,NumImproperParams << " IMPROPER PARAMETERS\n" \
04858 << "*****************************************" );
04859
04860 traverse_improper_params(improperp);
04861 }
04862
04863
04864
04865
04866
04867
04868
04869
04870
04871
04872 void Parameters::print_vdw_params()
04873 {
04874 DebugM(3,NumVdwParams << " vdW PARAMETERS\n" \
04875 << "*****************************************" );
04876
04877 traverse_vdw_params(vdwp);
04878 }
04879
04880
04881
04882
04883
04884
04885
04886
04887
04888
04889 void Parameters::print_vdw_pair_params()
04890 {
04891 DebugM(3,NumVdwPairParams << " vdW PAIR PARAMETERS\n" \
04892 << "*****************************************" );
04893
04894 traverse_vdw_pair_params(vdw_pairp);
04895 }
04896
04897
04898
04899
04900
04901
04902
04903
04904
04905
04906 void Parameters::print_nbthole_pair_params()
04907 {
04908 DebugM(3,NumNbtholePairParams << " NBTHOLE PAIR PARAMETERS\n" \
04909 << "*****************************************" );
04910
04911 traverse_nbthole_pair_params(nbthole_pairp);
04912 }
04913
04914
04915
04916
04917
04918
04919
04920
04921
04922
04923 void Parameters::print_param_summary()
04924 {
04925 iout << iINFO << "SUMMARY OF PARAMETERS:\n"
04926 << iINFO << NumBondParams << " BONDS\n"
04927 << iINFO << NumAngleParams << " ANGLES\n" << endi;
04928 if (cosAngles) {
04929 iout << iINFO << " " << (NumAngleParams - NumCosAngles) << " HARMONIC\n"
04930 << iINFO << " " << NumCosAngles << " COSINE-BASED\n" << endi;
04931 }
04932 iout << iINFO << NumDihedralParams << " DIHEDRAL\n"
04933 << iINFO << NumImproperParams << " IMPROPER\n"
04934 << iINFO << NumCrosstermParams << " CROSSTERM\n"
04935 << iINFO << NumVdwParams << " VDW\n"
04936 << iINFO << NumVdwPairParams << " VDW_PAIRS\n"
04937 << iINFO << NumNbtholePairParams << " NBTHOLE_PAIRS\n" << endi;
04938 }
04939
04940
04941
04942
04943
04944
04945
04946
04947
04948
04949
04950
04951
04952
04953
04954
04955
04956 void Parameters::done_reading_structure()
04957
04958 {
04959 if (bondp != NULL)
04960 free_bond_tree(bondp);
04961
04962 if (anglep != NULL)
04963 free_angle_tree(anglep);
04964
04965 if (dihedralp != NULL)
04966 free_dihedral_list(dihedralp);
04967
04968 if (improperp != NULL)
04969 free_improper_list(improperp);
04970
04971 if (crosstermp != NULL)
04972 free_crossterm_list(crosstermp);
04973
04974 if (vdwp != NULL)
04975 free_vdw_tree(vdwp);
04976
04977
04978
04979 if (maxDihedralMults != NULL)
04980 delete [] maxDihedralMults;
04981
04982 if (maxImproperMults != NULL)
04983 delete [] maxImproperMults;
04984
04985 bondp=NULL;
04986 anglep=NULL;
04987 dihedralp=NULL;
04988 improperp=NULL;
04989 crosstermp=NULL;
04990 vdwp=NULL;
04991 maxImproperMults=NULL;
04992 maxDihedralMults=NULL;
04993 }
04994
04995
04996
04997
04998
04999
05000
05001
05002
05003
05004
05005 void Parameters::send_Parameters(MOStream *msg)
05006 {
05007 Real *a1, *a2, *a3, *a4;
05008 int *i1, *i2, *i3;
05009 int i, j;
05010 Real **kvals;
05011 int **nvals;
05012 Real **deltavals;
05013
05014
05015
05016
05017
05018
05019
05020 msg->put(NumBondParams);
05021
05022 if (NumBondParams)
05023 {
05024 a1 = new Real[NumBondParams];
05025 a2 = new Real[NumBondParams];
05026
05027 if ( (a1 == NULL) || (a2 == NULL) )
05028 {
05029 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05030 }
05031
05032 for (i=0; i<NumBondParams; i++)
05033 {
05034 a1[i] = bond_array[i].k;
05035 a2[i] = bond_array[i].x0;
05036 }
05037
05038 msg->put(NumBondParams, a1)->put(NumBondParams, a2);
05039
05040 delete [] a1;
05041 delete [] a2;
05042 }
05043
05044
05045 msg->put(NumAngleParams);
05046
05047 if (NumAngleParams)
05048 {
05049 a1 = new Real[NumAngleParams];
05050 a2 = new Real[NumAngleParams];
05051 a3 = new Real[NumAngleParams];
05052 a4 = new Real[NumAngleParams];
05053 i1 = new int[NumAngleParams];
05054
05055 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) ||
05056 (a4 == NULL) || (i1 == NULL))
05057 {
05058 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05059 }
05060
05061 for (i=0; i<NumAngleParams; i++)
05062 {
05063 a1[i] = angle_array[i].k;
05064 a2[i] = angle_array[i].theta0;
05065 a3[i] = angle_array[i].k_ub;
05066 a4[i] = angle_array[i].r_ub;
05067 i1[i] = angle_array[i].normal;
05068 }
05069
05070 msg->put(NumAngleParams, a1)->put(NumAngleParams, a2);
05071 msg->put(NumAngleParams, a3)->put(NumAngleParams, a4);
05072 msg->put(NumAngleParams, i1);
05073
05074 delete [] a1;
05075 delete [] a2;
05076 delete [] a3;
05077 delete [] a4;
05078 delete [] i1;
05079 }
05080
05081
05082 msg->put(NumDihedralParams);
05083
05084 if (NumDihedralParams)
05085 {
05086 i1 = new int[NumDihedralParams];
05087 kvals = new Real *[MAX_MULTIPLICITY];
05088 nvals = new int *[MAX_MULTIPLICITY];
05089 deltavals = new Real *[MAX_MULTIPLICITY];
05090
05091 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
05092 (deltavals == NULL) )
05093 {
05094 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05095 }
05096
05097 for (i=0; i<MAX_MULTIPLICITY; i++)
05098 {
05099 kvals[i] = new Real[NumDihedralParams];
05100 nvals[i] = new int[NumDihedralParams];
05101 deltavals[i] = new Real[NumDihedralParams];
05102
05103 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05104 {
05105 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05106 }
05107 }
05108
05109 for (i=0; i<NumDihedralParams; i++)
05110 {
05111 i1[i] = dihedral_array[i].multiplicity;
05112
05113 for (j=0; j<MAX_MULTIPLICITY; j++)
05114 {
05115 kvals[j][i] = dihedral_array[i].values[j].k;
05116 nvals[j][i] = dihedral_array[i].values[j].n;
05117 deltavals[j][i] = dihedral_array[i].values[j].delta;
05118 }
05119 }
05120
05121 msg->put(NumDihedralParams, i1);
05122
05123 for (i=0; i<MAX_MULTIPLICITY; i++)
05124 {
05125 msg->put(NumDihedralParams, kvals[i]);
05126 msg->put(NumDihedralParams, nvals[i]);
05127 msg->put(NumDihedralParams, deltavals[i]);
05128
05129 delete [] kvals[i];
05130 delete [] nvals[i];
05131 delete [] deltavals[i];
05132 }
05133
05134 delete [] i1;
05135 delete [] kvals;
05136 delete [] nvals;
05137 delete [] deltavals;
05138 }
05139
05140
05141 msg->put(NumImproperParams);
05142
05143 if (NumImproperParams)
05144 {
05145 i1 = new int[NumImproperParams];
05146 kvals = new Real *[MAX_MULTIPLICITY];
05147 nvals = new int *[MAX_MULTIPLICITY];
05148 deltavals = new Real *[MAX_MULTIPLICITY];
05149
05150 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
05151 (deltavals == NULL) )
05152 {
05153 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05154 }
05155
05156 for (i=0; i<MAX_MULTIPLICITY; i++)
05157 {
05158 kvals[i] = new Real[NumImproperParams];
05159 nvals[i] = new int[NumImproperParams];
05160 deltavals[i] = new Real[NumImproperParams];
05161
05162 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05163 {
05164 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05165 }
05166 }
05167
05168 for (i=0; i<NumImproperParams; i++)
05169 {
05170 i1[i] = improper_array[i].multiplicity;
05171
05172 for (j=0; j<MAX_MULTIPLICITY; j++)
05173 {
05174 kvals[j][i] = improper_array[i].values[j].k;
05175 nvals[j][i] = improper_array[i].values[j].n;
05176 deltavals[j][i] = improper_array[i].values[j].delta;
05177 }
05178 }
05179
05180 msg->put(NumImproperParams, i1);
05181
05182 for (i=0; i<MAX_MULTIPLICITY; i++)
05183 {
05184 msg->put(NumImproperParams, kvals[i]);
05185 msg->put(NumImproperParams, nvals[i]);
05186 msg->put(NumImproperParams, deltavals[i]);
05187
05188 delete [] kvals[i];
05189 delete [] nvals[i];
05190 delete [] deltavals[i];
05191 }
05192
05193 delete [] i1;
05194 delete [] kvals;
05195 delete [] nvals;
05196 delete [] deltavals;
05197 }
05198
05199
05200 msg->put(NumCrosstermParams);
05201
05202 if (NumCrosstermParams)
05203 {
05204 for (i=0; i<NumCrosstermParams; ++i) {
05205 int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
05206 msg->put(nvals,&crossterm_array[i].c[0][0].d00);
05207 }
05208 }
05209
05210
05211 msg->put(numenerentries);
05212
05213 if (numenerentries) {
05214
05215
05216
05217
05218
05219
05220
05221
05222 msg->put(numenerentries, table_ener);
05223 }
05224
05225
05226 msg->put(NumVdwParams);
05227 msg->put(NumVdwParamsAssigned);
05228
05229 if (NumVdwParams)
05230 {
05231 a1 = new Real[NumVdwParams];
05232 a2 = new Real[NumVdwParams];
05233 a3 = new Real[NumVdwParams];
05234 a4 = new Real[NumVdwParams];
05235
05236 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) )
05237 {
05238 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05239 }
05240
05241 for (i=0; i<NumVdwParams; i++)
05242 {
05243 a1[i] = vdw_array[i].sigma;
05244 a2[i] = vdw_array[i].epsilon;
05245 a3[i] = vdw_array[i].sigma14;
05246 a4[i] = vdw_array[i].epsilon14;
05247 }
05248
05249 msg->put(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
05250 msg->put(NumVdwParams, a1);
05251 msg->put(NumVdwParams, a2);
05252 msg->put(NumVdwParams, a3);
05253 msg->put(NumVdwParams, a4);
05254
05255 delete [] a1;
05256 delete [] a2;
05257 delete [] a3;
05258 delete [] a4;
05259 }
05260
05261
05262 msg->put(NumVdwPairParams);
05263
05264 if (NumVdwPairParams)
05265 {
05266 a1 = new Real[NumVdwPairParams];
05267 a2 = new Real[NumVdwPairParams];
05268 a3 = new Real[NumVdwPairParams];
05269 a4 = new Real[NumVdwPairParams];
05270 i1 = new int[NumVdwPairParams];
05271 i2 = new int[NumVdwPairParams];
05272
05273 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
05274 (i1 == NULL) || (i2 == NULL) )
05275 {
05276 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05277 }
05278
05279 vdw_pair_to_arrays(i1, i2, a1, a2, a3, a4, 0, vdw_pair_tree);
05280
05281 msg->put(NumVdwPairParams, i1)->put(NumVdwPairParams, i2);
05282 msg->put(NumVdwPairParams, a1);
05283 msg->put(NumVdwPairParams, a2)->put(NumVdwPairParams, a3);
05284 msg->put(NumVdwPairParams, a4);
05285 }
05286
05287
05288 msg->put(NumNbtholePairParams);
05289
05290 if (NumNbtholePairParams)
05291 {
05292 a1 = new Real[NumNbtholePairParams];
05293 a2 = new Real[NumNbtholePairParams];
05294 a3 = new Real[NumNbtholePairParams];
05295 i1 = new int[NumNbtholePairParams];
05296 i2 = new int[NumNbtholePairParams];
05297
05298 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05299 {
05300 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05301 }
05302
05303 nbthole_pair_to_arrays(i1, i2, a1, a2, a3, 0, nbthole_pair_tree);
05304
05305 for (i=0; i<NumNbtholePairParams; i++)
05306 {
05307 nbthole_array[i].ind1 = i1[i];
05308 nbthole_array[i].ind2 = i2[i];
05309 nbthole_array[i].alphai = a1[i];
05310 nbthole_array[i].alphaj = a2[i];
05311 nbthole_array[i].tholeij = a3[i];
05312 }
05313
05314 msg->put(NumNbtholePairParams, i1)->put(NumNbtholePairParams, i2);
05315 msg->put(NumNbtholePairParams, a1);
05316 msg->put(NumNbtholePairParams, a2)->put(NumNbtholePairParams, a3);
05317 }
05318
05319
05320
05321 msg->put(NumTablePairParams);
05322
05323 if (NumTablePairParams)
05324 {
05325 i1 = new int[NumTablePairParams];
05326 i2 = new int[NumTablePairParams];
05327 i3 = new int[NumTablePairParams];
05328
05329 if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05330 {
05331 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05332 }
05333
05334 table_pair_to_arrays(i1, i2, i3, 0, tab_pair_tree);
05335
05336 msg->put(NumTablePairParams, i1)->put(NumTablePairParams, i2);
05337 msg->put(NumTablePairParams, i3);
05338 }
05339
05340
05341
05342 msg->end();
05343 delete msg;
05344 }
05345
05346
05347
05348
05349
05350
05351
05352
05353
05354
05355 void Parameters::receive_Parameters(MIStream *msg)
05356
05357 {
05358 int i, j;
05359 Real *a1, *a2, *a3, *a4;
05360 int *i1, *i2, *i3;
05361 IndexedVdwPair *new_node;
05362 IndexedTablePair *new_tab_node;
05363 Real **kvals;
05364 int **nvals;
05365 Real **deltavals;
05366
05367
05368 msg->get(NumBondParams);
05369
05370 if (NumBondParams)
05371 {
05372 bond_array = new BondValue[NumBondParams];
05373 a1 = new Real[NumBondParams];
05374 a2 = new Real[NumBondParams];
05375
05376 if ( (bond_array == NULL) || (a1 == NULL) || (a2 == NULL) )
05377 {
05378 NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05379 }
05380
05381 msg->get(NumBondParams, a1);
05382 msg->get(NumBondParams, a2);
05383
05384 for (i=0; i<NumBondParams; i++)
05385 {
05386 bond_array[i].k = a1[i];
05387 bond_array[i].x0 = a2[i];
05388 }
05389
05390 delete [] a1;
05391 delete [] a2;
05392 }
05393
05394
05395 msg->get(NumAngleParams);
05396
05397 if (NumAngleParams)
05398 {
05399 angle_array = new AngleValue[NumAngleParams];
05400 a1 = new Real[NumAngleParams];
05401 a2 = new Real[NumAngleParams];
05402 a3 = new Real[NumAngleParams];
05403 a4 = new Real[NumAngleParams];
05404 i1 = new int[NumAngleParams];
05405
05406 if ( (angle_array == NULL) || (a1 == NULL) || (a2 == NULL) ||
05407 (a3 == NULL) || (a4 == NULL) || (i1 == NULL))
05408 {
05409 NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05410 }
05411
05412 msg->get(NumAngleParams, a1);
05413 msg->get(NumAngleParams, a2);
05414 msg->get(NumAngleParams, a3);
05415 msg->get(NumAngleParams, a4);
05416 msg->get(NumAngleParams, i1);
05417
05418 for (i=0; i<NumAngleParams; i++)
05419 {
05420 angle_array[i].k = a1[i];
05421 angle_array[i].theta0 = a2[i];
05422 angle_array[i].k_ub = a3[i];
05423 angle_array[i].r_ub = a4[i];
05424 angle_array[i].normal = i1[i];
05425 }
05426
05427 delete [] a1;
05428 delete [] a2;
05429 delete [] a3;
05430 delete [] a4;
05431 delete [] i1;
05432 }
05433
05434
05435 msg->get(NumDihedralParams);
05436
05437 if (NumDihedralParams)
05438 {
05439 dihedral_array = new DihedralValue[NumDihedralParams];
05440
05441 i1 = new int[NumDihedralParams];
05442 kvals = new Real *[MAX_MULTIPLICITY];
05443 nvals = new int *[MAX_MULTIPLICITY];
05444 deltavals = new Real *[MAX_MULTIPLICITY];
05445
05446 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
05447 (deltavals == NULL) || (dihedral_array == NULL) )
05448 {
05449 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05450 }
05451
05452 for (i=0; i<MAX_MULTIPLICITY; i++)
05453 {
05454 kvals[i] = new Real[NumDihedralParams];
05455 nvals[i] = new int[NumDihedralParams];
05456 deltavals[i] = new Real[NumDihedralParams];
05457
05458 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05459 {
05460 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05461 }
05462 }
05463
05464 msg->get(NumDihedralParams, i1);
05465
05466 for (i=0; i<MAX_MULTIPLICITY; i++)
05467 {
05468 msg->get(NumDihedralParams, kvals[i]);
05469 msg->get(NumDihedralParams, nvals[i]);
05470 msg->get(NumDihedralParams, deltavals[i]);
05471 }
05472
05473 for (i=0; i<NumDihedralParams; i++)
05474 {
05475 dihedral_array[i].multiplicity = i1[i];
05476
05477 for (j=0; j<MAX_MULTIPLICITY; j++)
05478 {
05479 dihedral_array[i].values[j].k = kvals[j][i];
05480 dihedral_array[i].values[j].n = nvals[j][i];
05481 dihedral_array[i].values[j].delta = deltavals[j][i];
05482 }
05483 }
05484
05485 for (i=0; i<MAX_MULTIPLICITY; i++)
05486 {
05487 delete [] kvals[i];
05488 delete [] nvals[i];
05489 delete [] deltavals[i];
05490 }
05491
05492 delete [] i1;
05493 delete [] kvals;
05494 delete [] nvals;
05495 delete [] deltavals;
05496 }
05497
05498
05499 msg->get(NumImproperParams);
05500
05501 if (NumImproperParams)
05502 {
05503 improper_array = new ImproperValue[NumImproperParams];
05504 i1 = new int[NumImproperParams];
05505 kvals = new Real *[MAX_MULTIPLICITY];
05506 nvals = new int *[MAX_MULTIPLICITY];
05507 deltavals = new Real *[MAX_MULTIPLICITY];
05508
05509 if ( (i1 == NULL) || (kvals == NULL) || (nvals == NULL) ||
05510 (deltavals == NULL) || (improper_array==NULL) )
05511 {
05512 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05513 }
05514
05515 for (i=0; i<MAX_MULTIPLICITY; i++)
05516 {
05517 kvals[i] = new Real[NumImproperParams];
05518 nvals[i] = new int[NumImproperParams];
05519 deltavals[i] = new Real[NumImproperParams];
05520
05521 if ( (kvals[i] == NULL) || (nvals[i] == NULL) || (deltavals[i] == NULL) )
05522 {
05523 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05524 }
05525 }
05526
05527 msg->get(NumImproperParams,i1);
05528
05529 for (i=0; i<MAX_MULTIPLICITY; i++)
05530 {
05531 msg->get(NumImproperParams,kvals[i]);
05532 msg->get(NumImproperParams,nvals[i]);
05533 msg->get(NumImproperParams,deltavals[i]);
05534 }
05535
05536 for (i=0; i<NumImproperParams; i++)
05537 {
05538 improper_array[i].multiplicity = i1[i];
05539
05540 for (j=0; j<MAX_MULTIPLICITY; j++)
05541 {
05542 improper_array[i].values[j].k = kvals[j][i];
05543 improper_array[i].values[j].n = nvals[j][i];
05544 improper_array[i].values[j].delta = deltavals[j][i];
05545 }
05546 }
05547
05548 for (i=0; i<MAX_MULTIPLICITY; i++)
05549 {
05550 delete [] kvals[i];
05551 delete [] nvals[i];
05552 delete [] deltavals[i];
05553 }
05554
05555 delete [] i1;
05556 delete [] kvals;
05557 delete [] nvals;
05558 delete [] deltavals;
05559 }
05560
05561
05562 msg->get(NumCrosstermParams);
05563
05564 if (NumCrosstermParams)
05565 {
05566 crossterm_array = new CrosstermValue[NumCrosstermParams];
05567
05568 for (i=0; i<NumCrosstermParams; ++i) {
05569 int nvals = CrosstermValue::dim * CrosstermValue::dim * 2 * 2;
05570 msg->get(nvals,&crossterm_array[i].c[0][0].d00);
05571 }
05572 }
05573
05574
05575 msg->get(numenerentries);
05576 if (numenerentries > 0) {
05577
05578
05579 table_ener = new BigReal[numenerentries];
05580
05581 if (table_ener==NULL)
05582 {
05583 NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05584 }
05585
05586 msg->get(numenerentries, table_ener);
05587 }
05588
05589
05590 msg->get(NumVdwParams);
05591 msg->get(NumVdwParamsAssigned);
05592
05593 if (NumVdwParams)
05594 {
05595 atomTypeNames = new char[NumVdwParams*(MAX_ATOMTYPE_CHARS+1)];
05596 vdw_array = new VdwValue[NumVdwParams];
05597 a1 = new Real[NumVdwParams];
05598 a2 = new Real[NumVdwParams];
05599 a3 = new Real[NumVdwParams];
05600 a4 = new Real[NumVdwParams];
05601
05602 if ( (vdw_array==NULL) || (a1==NULL) || (a2==NULL) || (a3==NULL)
05603 || (a4==NULL) || (atomTypeNames==NULL))
05604 {
05605 NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05606 }
05607
05608 msg->get(NumVdwParams * (MAX_ATOMTYPE_CHARS+1), atomTypeNames);
05609 msg->get(NumVdwParams, a1);
05610 msg->get(NumVdwParams, a2);
05611 msg->get(NumVdwParams, a3);
05612 msg->get(NumVdwParams, a4);
05613
05614 for (i=0; i<NumVdwParams; i++)
05615 {
05616 vdw_array[i].sigma = a1[i];
05617 vdw_array[i].epsilon = a2[i];
05618 vdw_array[i].sigma14 = a3[i];
05619 vdw_array[i].epsilon14 = a4[i];
05620 }
05621
05622 delete [] a1;
05623 delete [] a2;
05624 delete [] a3;
05625 delete [] a4;
05626 }
05627
05628
05629 msg->get(NumVdwPairParams);
05630
05631 if (NumVdwPairParams)
05632 {
05633 a1 = new Real[NumVdwPairParams];
05634 a2 = new Real[NumVdwPairParams];
05635 a3 = new Real[NumVdwPairParams];
05636 a4 = new Real[NumVdwPairParams];
05637 i1 = new int[NumVdwPairParams];
05638 i2 = new int[NumVdwPairParams];
05639
05640 if ( (a1 == NULL) || (a2 == NULL) || (a3 == NULL) || (a4 == NULL) ||
05641 (i1 == NULL) || (i2 == NULL) )
05642 {
05643 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05644 }
05645
05646 msg->get(NumVdwPairParams, i1);
05647 msg->get(NumVdwPairParams, i2);
05648 msg->get(NumVdwPairParams, a1);
05649 msg->get(NumVdwPairParams, a2);
05650 msg->get(NumVdwPairParams, a3);
05651 msg->get(NumVdwPairParams, a4);
05652
05653 for (i=0; i<NumVdwPairParams; i++)
05654 {
05655 new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
05656
05657 if (new_node == NULL)
05658 {
05659 NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05660 }
05661
05662 new_node->ind1 = i1[i];
05663 new_node->ind2 = i2[i];
05664 new_node->A = a1[i];
05665 new_node->A14 = a2[i];
05666 new_node->B = a3[i];
05667 new_node->B14 = a4[i];
05668 new_node->left = NULL;
05669 new_node->right = NULL;
05670
05671 vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
05672 }
05673
05674 delete [] i1;
05675 delete [] i2;
05676 delete [] a1;
05677 delete [] a2;
05678 delete [] a3;
05679 delete [] a4;
05680 }
05681
05682
05683 msg->get(NumNbtholePairParams);
05684
05685 if (NumNbtholePairParams)
05686 {
05687 nbthole_array = new NbtholePairValue[NumNbtholePairParams];
05688 a1 = new Real[NumNbtholePairParams];
05689 a2 = new Real[NumNbtholePairParams];
05690 a3 = new Real[NumNbtholePairParams];
05691 i1 = new int[NumNbtholePairParams];
05692 i2 = new int[NumNbtholePairParams];
05693
05694 if ( (nbthole_array == NULL) || (a1 == NULL) || (a2 == NULL) || (a3 == NULL)
05695 || (i1 == NULL) || (i2 == NULL) )
05696 {
05697 NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05698 }
05699
05700 msg->get(NumNbtholePairParams, i1);
05701 msg->get(NumNbtholePairParams, i2);
05702 msg->get(NumNbtholePairParams, a1);
05703 msg->get(NumNbtholePairParams, a2);
05704 msg->get(NumNbtholePairParams, a3);
05705
05706 for (i=0; i<NumNbtholePairParams; i++)
05707 {
05708
05709 nbthole_array[i].ind1 = i1[i];
05710 nbthole_array[i].ind2 = i2[i];
05711 nbthole_array[i].alphai = a1[i];
05712 nbthole_array[i].alphaj = a2[i];
05713 nbthole_array[i].tholeij = a3[i];
05714
05715 }
05716
05717 delete [] i1;
05718 delete [] i2;
05719 delete [] a1;
05720 delete [] a2;
05721 delete [] a3;
05722 }
05723
05724
05725 msg->get(NumTablePairParams);
05726
05727 if (NumTablePairParams)
05728 {
05729 i1 = new int[NumTablePairParams];
05730 i2 = new int[NumTablePairParams];
05731 i3 = new int[NumTablePairParams];
05732
05733 if ( (i3 == NULL) || (i1 == NULL) || (i2 == NULL) )
05734 {
05735 NAMD_die("memory allocation failed in Parameters::send_Parameters");
05736 }
05737
05738 msg->get(NumTablePairParams, i1);
05739 msg->get(NumTablePairParams, i2);
05740 msg->get(NumTablePairParams, i3);
05741
05742 for (i=0; i<NumTablePairParams; i++)
05743 {
05744 new_tab_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
05745
05746 if (new_tab_node == NULL)
05747 {
05748 NAMD_die("memory allocation failed in Parameters::receive_Parameters");
05749 }
05750
05751
05752 new_tab_node->ind1 = i1[i];
05753 new_tab_node->ind2 = i2[i];
05754 new_tab_node->K = i3[i];
05755 new_tab_node->left = NULL;
05756 new_tab_node->right = NULL;
05757
05758 tab_pair_tree = add_to_indexed_table_pairs(new_tab_node, tab_pair_tree);
05759 }
05760
05761 delete [] i1;
05762 delete [] i2;
05763 delete [] i3;
05764 }
05765
05766
05767
05768
05769 AllFilesRead = TRUE;
05770
05771 delete msg;
05772 }
05773
05774
05775
05776
05777
05778
05779
05780
05781
05782
05783
05784
05785
05786 void Parameters::convert_vdw_pairs()
05787
05788 {
05789 #ifdef MEM_OPT_VERSION
05790 AtomCstInfo atom_struct;
05791 #else
05792 Atom atom_struct;
05793 #endif
05794 Index index1, index2;
05795 IndexedVdwPair *new_node;
05796 struct vdw_pair_params *ptr, *next;
05797
05798 ptr = vdw_pairp;
05799
05800
05801
05802 while (ptr != NULL)
05803 {
05804 new_node = (IndexedVdwPair *) malloc(sizeof(IndexedVdwPair));
05805
05806 if (new_node == NULL)
05807 {
05808 NAMD_die("memory allocation failed in Parameters::convert_vdw_pairs");
05809 }
05810
05811
05812
05813 assign_vdw_index(ptr->atom1name, &atom_struct);
05814 index1 = atom_struct.vdw_type;
05815 assign_vdw_index(ptr->atom2name, &atom_struct);
05816 index2 = atom_struct.vdw_type;
05817
05818 if (index1 > index2)
05819 {
05820 new_node->ind1 = index2;
05821 new_node->ind2 = index1;
05822 }
05823 else
05824 {
05825 new_node->ind1 = index1;
05826 new_node->ind2 = index2;
05827 }
05828
05829 new_node->A = ptr->A;
05830 new_node->B = ptr->B;
05831 new_node->A14 = ptr->A14;
05832 new_node->B14 = ptr->B14;
05833
05834 new_node->left = NULL;
05835 new_node->right = NULL;
05836
05837
05838 vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
05839
05840
05841 next = ptr->next;
05842
05843 delete ptr;
05844
05845 ptr = next;
05846 }
05847
05848 vdw_pairp = NULL;
05849
05850 }
05851
05852
05853
05854
05855
05856
05857
05858
05859
05860
05861
05862
05863
05864 void Parameters::convert_nbthole_pairs()
05865
05866 {
05867 #ifdef MEM_OPT_VERSION
05868 AtomCstInfo atom_struct;
05869 #else
05870 Atom atom_struct;
05871 #endif
05872 Index index1, index2;
05873 IndexedNbtholePair *new_node;
05874 struct nbthole_pair_params *ptr, *next;
05875
05876 ptr = nbthole_pairp;
05877
05878
05879
05880 while (ptr != NULL)
05881 {
05882 new_node = (IndexedNbtholePair *) malloc(sizeof(IndexedNbtholePair));
05883
05884 if (new_node == NULL)
05885 {
05886 NAMD_die("memory allocation failed in Parameters::convert_nbthole_pairs");
05887 }
05888
05889
05890
05891 assign_vdw_index(ptr->atom1name, &atom_struct);
05892 index1 = atom_struct.vdw_type;
05893 assign_vdw_index(ptr->atom2name, &atom_struct);
05894 index2 = atom_struct.vdw_type;
05895
05896 if (index1 > index2)
05897 {
05898 new_node->ind1 = index2;
05899 new_node->ind2 = index1;
05900 }
05901 else
05902 {
05903 new_node->ind1 = index1;
05904 new_node->ind2 = index2;
05905 }
05906
05907 new_node->alphai = ptr->alphai;
05908 new_node->alphaj = ptr->alphaj;
05909 new_node->tholeij = ptr->tholeij;
05910
05911 new_node->left = NULL;
05912 new_node->right = NULL;
05913
05914
05915 nbthole_pair_tree = add_to_indexed_nbthole_pairs(new_node, nbthole_pair_tree);
05916
05917
05918 next = ptr->next;
05919
05920 delete ptr;
05921
05922 ptr = next;
05923 }
05924
05925 nbthole_pairp = NULL;
05926
05927 }
05928
05929
05930
05931
05932
05933
05934
05935
05936
05937
05938
05939
05940
05941 void Parameters::convert_table_pairs()
05942
05943 {
05944 #ifdef MEM_OPT_VERSION
05945 AtomCstInfo atom_struct;
05946 #else
05947 Atom atom_struct;
05948 #endif
05949 Index index1, index2;
05950 IndexedTablePair *new_node;
05951 struct table_pair_params *ptr, *next;
05952
05953 ptr = table_pairp;
05954
05955
05956
05957 while (ptr != NULL)
05958 {
05959 new_node = (IndexedTablePair *) malloc(sizeof(IndexedTablePair));
05960
05961 if (new_node == NULL)
05962 {
05963 NAMD_die("memory allocation failed in Parameters::convert_table_pairs");
05964 }
05965
05966
05967
05968 assign_vdw_index(ptr->atom1name, &atom_struct);
05969 index1 = atom_struct.vdw_type;
05970 assign_vdw_index(ptr->atom2name, &atom_struct);
05971 index2 = atom_struct.vdw_type;
05972
05973
05974
05975 if (index1 > index2)
05976 {
05977 new_node->ind1 = index2;
05978 new_node->ind2 = index1;
05979 }
05980 else
05981 {
05982 new_node->ind1 = index1;
05983 new_node->ind2 = index2;
05984 }
05985
05986 new_node->K = ptr->K;
05987
05988 new_node->left = NULL;
05989 new_node->right = NULL;
05990
05991
05992 tab_pair_tree = add_to_indexed_table_pairs(new_node, tab_pair_tree);
05993
05994
05995 next = ptr->next;
05996
05997 delete ptr;
05998
05999 ptr = next;
06000 }
06001
06002 table_pairp = NULL;
06003
06004 }
06005
06006
06007
06008
06009
06010
06011
06012
06013
06014
06015
06016
06017
06018
06019
06020 IndexedTablePair *Parameters::add_to_indexed_table_pairs(IndexedTablePair *new_node,
06021 IndexedTablePair *tree)
06022
06023 {
06024 if (tree == NULL)
06025 return(new_node);
06026
06027 if ( (new_node->ind1 < tree->ind1) ||
06028 ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
06029 {
06030 tree->left = add_to_indexed_table_pairs(new_node, tree->left);
06031 }
06032 else
06033 {
06034 tree->right = add_to_indexed_table_pairs(new_node, tree->right);
06035 }
06036
06037 return(tree);
06038 }
06039
06040
06041
06042
06043
06044
06045
06046
06047
06048
06049
06050
06051
06052
06053
06054 IndexedVdwPair *Parameters::add_to_indexed_vdw_pairs(IndexedVdwPair *new_node,
06055 IndexedVdwPair *tree)
06056
06057 {
06058 if (tree == NULL)
06059 return(new_node);
06060
06061 if ( (new_node->ind1 < tree->ind1) ||
06062 ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
06063 {
06064 tree->left = add_to_indexed_vdw_pairs(new_node, tree->left);
06065 }
06066 else
06067 {
06068 tree->right = add_to_indexed_vdw_pairs(new_node, tree->right);
06069 }
06070
06071 return(tree);
06072 }
06073
06074
06075
06076
06077
06078
06079
06080
06081
06082
06083
06084
06085
06086
06087
06088 IndexedNbtholePair *Parameters::add_to_indexed_nbthole_pairs(IndexedNbtholePair *new_node,
06089 IndexedNbtholePair *tree)
06090
06091 {
06092 if (tree == NULL)
06093 return(new_node);
06094
06095 if ( (new_node->ind1 < tree->ind1) ||
06096 ((new_node->ind1 == tree->ind1) && (new_node->ind2 < tree->ind2)) )
06097 {
06098 tree->left = add_to_indexed_nbthole_pairs(new_node, tree->left);
06099 }
06100 else
06101 {
06102 tree->right = add_to_indexed_nbthole_pairs(new_node, tree->right);
06103 }
06104
06105 return(tree);
06106 }
06107
06108
06109
06110
06111
06112
06113
06114
06115
06116
06117
06118
06119
06120
06121
06122
06123
06124
06125
06126
06127
06128
06129
06130 int Parameters::vdw_pair_to_arrays(int *ind1_array, int *ind2_array,
06131 Real *A, Real *A14,
06132 Real *B, Real *B14,
06133 int arr_index, IndexedVdwPair *tree)
06134
06135 {
06136 if (tree == NULL)
06137 return(arr_index);
06138
06139 ind1_array[arr_index] = tree->ind1;
06140 ind2_array[arr_index] = tree->ind2;
06141 A[arr_index] = tree->A;
06142 A14[arr_index] = tree->A14;
06143 B[arr_index] = tree->B;
06144 B14[arr_index] = tree->B14;
06145
06146 arr_index++;
06147
06148 arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
06149 arr_index, tree->left);
06150 arr_index = vdw_pair_to_arrays(ind1_array, ind2_array, A, A14, B, B14,
06151 arr_index, tree->right);
06152
06153 return(arr_index);
06154 }
06155
06156
06157
06158
06159
06160
06161
06162
06163
06164
06165
06166
06167
06168
06169
06170
06171
06172
06173
06174
06175 int Parameters::nbthole_pair_to_arrays(int *ind1_array, int *ind2_array,
06176 Real *alphai, Real *alphaj, Real *tholeij,
06177 int arr_index, IndexedNbtholePair *tree)
06178
06179 {
06180 if (tree == NULL)
06181 return(arr_index);
06182
06183 ind1_array[arr_index] = tree->ind1;
06184 ind2_array[arr_index] = tree->ind2;
06185 alphai[arr_index] = tree->alphai;
06186 alphaj[arr_index] = tree->alphaj;
06187 tholeij[arr_index] = tree->tholeij;
06188
06189 arr_index++;
06190
06191 arr_index = nbthole_pair_to_arrays(ind1_array, ind2_array, alphai,
06192 alphaj, tholeij, arr_index, tree->left);
06193 arr_index = nbthole_pair_to_arrays(ind1_array, ind2_array, alphai,
06194 alphaj, tholeij, arr_index, tree->right);
06195
06196 return(arr_index);
06197 }
06198
06199
06200
06201
06202
06203
06204
06205
06206
06207
06208
06209
06210
06211
06212
06213
06214
06215
06216 int Parameters::table_pair_to_arrays(int *ind1_array, int *ind2_array,
06217 int *K,
06218 int arr_index, IndexedTablePair *tree)
06219
06220 {
06221 if (tree == NULL)
06222 return(arr_index);
06223
06224 ind1_array[arr_index] = tree->ind1;
06225 ind2_array[arr_index] = tree->ind2;
06226 K[arr_index] = tree->K;
06227
06228 arr_index++;
06229
06230 arr_index = table_pair_to_arrays(ind1_array, ind2_array, K,
06231 arr_index, tree->left);
06232 arr_index = table_pair_to_arrays(ind1_array, ind2_array, K,
06233 arr_index, tree->right);
06234
06235 return(arr_index);
06236 }
06237
06238
06239
06240
06241
06242
06243
06244
06245
06246
06247 Parameters::Parameters(Ambertoppar *amber_data, BigReal vdw14)
06248 {
06249 initialize();
06250
06251
06252 read_parm(amber_data,vdw14);
06253 }
06254
06255
06256
06257
06258
06259
06260
06261
06262
06263
06264
06265
06266
06267
06268
06269 void Parameters::read_parm(Ambertoppar *amber_data, BigReal vdw14)
06270 {
06271 int i,j,ico,numtype,mul;
06272 IndexedVdwPair *new_node;
06273
06274 if (!amber_data->data_read)
06275 NAMD_die("No data read from parm file yet!");
06276
06277
06278 NumBondParams = amber_data->Numbnd;
06279 if (NumBondParams)
06280 { bond_array = new BondValue[NumBondParams];
06281 if (bond_array == NULL)
06282 NAMD_die("memory allocation of bond_array failed!");
06283 }
06284 for (i=0; i<NumBondParams; ++i)
06285 { bond_array[i].k = amber_data->Rk[i];
06286 bond_array[i].x0 = amber_data->Req[i];
06287 }
06288
06289
06290 NumAngleParams = amber_data->Numang;
06291 if (NumAngleParams)
06292 { angle_array = new AngleValue[NumAngleParams];
06293 if (angle_array == NULL)
06294 NAMD_die("memory allocation of angle_array failed!");
06295 }
06296 for (i=0; i<NumAngleParams; ++i)
06297 { angle_array[i].k = amber_data->Tk[i];
06298 angle_array[i].theta0 = amber_data->Teq[i];
06299
06300 angle_array[i].k_ub = angle_array[i].r_ub = 0;
06301
06302 angle_array[i].normal = 1;
06303 }
06304
06305
06306
06307
06308
06309
06310
06311
06312 NumDihedralParams = amber_data->Nptra;
06313 if (NumDihedralParams)
06314 { dihedral_array = new DihedralValue[amber_data->Nptra];
06315 if (dihedral_array == NULL)
06316 NAMD_die("memory allocation of dihedral_array failed!");
06317 }
06318 mul = 0;
06319 for (i=0; i<NumDihedralParams; ++i)
06320 { dihedral_array[i-mul].values[mul].k = amber_data->Pk[i];
06321 dihedral_array[i-mul].values[mul].n = int(fabs(amber_data->Pn[i])+0.5);
06322 if (dihedral_array[i-mul].values[mul].n == 0)
06323 { char err_msg[128];
06324 sprintf(err_msg, "The periodicity of dihedral # %d is zero!", i+1);
06325 NAMD_die(err_msg);
06326 }
06327 dihedral_array[i-mul].values[mul].delta = amber_data->Phase[i];
06328
06329
06330 if (amber_data->Pn[i] > 0)
06331 { dihedral_array[i-mul].multiplicity = mul+1;
06332 mul = 0;
06333 }
06334 else if (++mul >= MAX_MULTIPLICITY)
06335 { char err_msg[181];
06336 sprintf(err_msg, "Multiple dihedral with multiplicity of %d greater than max of %d",
06337 mul+1, MAX_MULTIPLICITY);
06338 NAMD_die(err_msg);
06339 }
06340 }
06341 if (mul > 0)
06342 NAMD_die("Negative periodicity in the last dihedral entry!");
06343
06344
06345
06346
06347 NumVdwParamsAssigned = numtype = amber_data->Ntypes;
06348 NumVdwPairParams = numtype * (numtype+1) / 2;
06349 for (i=0; i<numtype; ++i)
06350 for (j=i; j<numtype; ++j)
06351 { new_node = new IndexedVdwPair;
06352 if (new_node == NULL)
06353 NAMD_die("memory allocation of vdw_pair_tree failed!");
06354 new_node->ind1 = i;
06355 new_node->ind2 = j;
06356 new_node->left = new_node->right = NULL;
06357
06358
06359
06360
06361
06362 ico = amber_data->Cno[numtype*i+j];
06363 if (ico>0)
06364 { new_node->A = amber_data->Cn1[ico-1];
06365 new_node->A14 = new_node->A / vdw14;
06366 new_node->B = amber_data->Cn2[ico-1];
06367 new_node->B14 = new_node->B / vdw14;
06368 }
06369 else if (amber_data->HB12[abs(ico)-1]==0.0 && amber_data->HB6[abs(ico)-1]==0.0)
06370 { new_node->A = new_node->A14 = new_node->B = new_node->B14 = 0.0;
06371 iout << iWARN << "Encounter 10-12 H-bond term\n";
06372 }
06373 else
06374 NAMD_die("Encounter non-zero 10-12 H-bond term!");
06375
06376 vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
06377 }
06378 }
06379
06380
06381
06382
06383
06384
06385
06386
06387
06388
06389 Parameters::Parameters(const GromacsTopFile *gf, Bool min)
06390 {
06391 initialize();
06392
06393
06394 read_parm(gf,min);
06395 }
06396
06397
06398
06399
06400
06401
06402
06403
06404
06405
06406
06407
06408
06409
06410
06411 void Parameters::read_parm(const GromacsTopFile *gf, Bool min)
06412 {
06413 int numtype;
06414 IndexedVdwPair *new_node;
06415 int i,j,funct;
06416 Real test1,test2;
06417
06418
06419 NumBondParams = gf->getNumBondParams();
06420 NumAngleParams = gf->getNumAngleParams();
06421 NumDihedralParams = gf->getNumDihedralParams();
06422 if (NumBondParams) {
06423 bond_array = new BondValue[NumBondParams];
06424 if (bond_array == NULL)
06425 NAMD_die("memory allocation of bond_array failed!");
06426 }
06427 if (NumDihedralParams) {
06428 dihedral_array = new DihedralValue[NumDihedralParams];
06429 if (dihedral_array == NULL)
06430 NAMD_die("memory allocation of dihedral_array failed!");
06431 }
06432 if (NumAngleParams) {
06433 angle_array = new AngleValue[NumAngleParams];
06434 if (angle_array == NULL)
06435 NAMD_die("memory allocation of angle_array failed!");
06436 }
06437
06438
06439
06440
06441 for (i=0;i<NumBondParams;i++) {
06442 Real x0,k;
06443 gf->getBondParams(i,
06444 &x0,
06445 &k,
06446 &funct);
06447 bond_array[i].x0 = x0;
06448 bond_array[i].k = k;
06449 }
06450
06451
06452
06453
06454 for (i=0;i<NumAngleParams;i++) {
06455 Real theta0,k;
06456 gf->getAngleParams(i,
06457 &theta0,
06458 &k,
06459 &funct);
06460 angle_array[i].theta0 = theta0*PI/180;
06461 angle_array[i].k = k;
06462
06463 angle_array[i].k_ub = angle_array[i].r_ub = 0;
06464 angle_array[i].normal = 1;
06465 }
06466
06467
06468
06469 for (i=0; i<NumDihedralParams; ++i) {
06470 Real c[6];
06471 int num_periods;
06472 int funct;
06473
06474 gf->getDihedralParams(i,c,&num_periods,&funct);
06475
06476 switch(funct) {
06477 case 1:
06478 dihedral_array[i].values[0].delta = c[0]*PI/180;
06479 dihedral_array[i].values[0].k = c[1];
06480 dihedral_array[i].values[0].n = num_periods;
06481 dihedral_array[i].multiplicity = 1;
06482 break;
06483 case 2:
06484 dihedral_array[i].values[0].delta = c[0]*PI/180;
06485 dihedral_array[i].values[0].k = c[1];
06486 dihedral_array[i].values[0].n = 0;
06487 dihedral_array[i].multiplicity = 1;
06488 break;
06489 case 3:
06490
06491
06492 if(MAX_MULTIPLICITY < 5)
06493 NAMD_die("I can't do RB dihedrals with MAX_MULTIPLICITY < 5");
06494 dihedral_array[i].multiplicity = 5;
06495
06496
06497
06498 for(j=0;j<6;j++) {
06499 if(j%2 == 1) c[j] = -c[j];
06500 }
06501
06502
06503
06504 for(j=0;j<5;j++) dihedral_array[i].values[j].delta = 0;
06505 for(j=0;j<5;j++) dihedral_array[i].values[j].n = j+1;
06506
06507
06508
06509
06510 dihedral_array[i].values[0].k = 1*c[1] + 3/4.*c[3] + 10/16.*c[5];
06511 dihedral_array[i].values[1].k = 1/2.*c[2] + 4/8.*c[4] ;
06512 dihedral_array[i].values[2].k = 1/4.*c[3] + 5/16.*c[5];
06513 dihedral_array[i].values[3].k = 1/8.*c[4] ;
06514 dihedral_array[i].values[4].k = 1/16.*c[5];
06515
06516
06517
06518 test1 = 0;
06519 for(j=5;j>=0;j--) {
06520 test1 *= cos(0.5);
06521 test1 += c[j];
06522 }
06523
06524 test2 = c[0]+1/2.*c[2]+3/8.*c[4];
06525 for(j=0;j<5;j++) {
06526 test2 += dihedral_array[i].values[j].k * cos((j+1)*0.5);
06527 }
06528
06529 if(fabs(test1-test2) > 0.0001)
06530 NAMD_die("Internal error: failed to handle RB dihedrals");
06531
06532
06533
06534
06535
06536
06537
06538
06539
06540
06541
06542
06543
06544
06545 break;
06546 default:
06547 NAMD_die("unknown dihedral type found");
06548 }
06549 }
06550
06551
06552
06553 Bool warned=false;
06554
06555 NumVdwParamsAssigned = numtype = gf->getNumAtomParams();
06556 NumVdwPairParams = numtype * (numtype+1) / 2;
06557 for (i=0; i<numtype; i++) {
06558 for (j=i; j<numtype; j++) {
06559
06560
06561 new_node = new IndexedVdwPair;
06562 if (new_node == NULL)
06563 NAMD_die("memory allocation of vdw_pair_tree failed!");
06564 new_node->ind1 = i;
06565 new_node->ind2 = j;
06566 new_node->left = new_node->right = NULL;
06567
06568 gf->getVDWParams(i,j, &(new_node->B), &(new_node->A),
06569 &(new_node->B14), &(new_node->A14));
06570
06571
06572
06573
06574
06575 if(min && ( fabs(new_node->A) < 1.0 )) {
06576 new_node->A = 1.0;
06577 if(!warned) {
06578 iout << iWARN <<
06579 "Adding small LJ repulsion term to some atoms.\n" << endi;
06580 warned=true;
06581 }
06582 }
06583
06584 vdw_pair_tree = add_to_indexed_vdw_pairs(new_node, vdw_pair_tree);
06585 }
06586 }
06587 }
06588
06589
06590
06591
06592
06593
06594
06595
06596
06597
06598
06599
06600
06601 void Parameters::read_ener_table(SimParameters *simParams) {
06602 char* table_file = simParams->tabulatedEnergiesFile;
06603 char* interp_type = simParams->tableInterpType;
06604 FILE* enertable;
06605 enertable = fopen(table_file, "r");
06606
06607 if (enertable == NULL) {
06608 NAMD_die("ERROR: Could not open energy table file!\n");
06609 }
06610
06611 char tableline[121];
06612 char* newtypename;
06613 int numtypes;
06614 int atom1;
06615 int atom2;
06616 int distbin;
06617 int readcount;
06618 Real dist;
06619 Real energy;
06620 Real force;
06621 Real table_spacing;
06622 Real maxdist;
06623
06624
06625 iout << "Beginning table read\n" << endi;
06626 while(fgets(tableline,120,enertable)) {
06627 if (strncmp(tableline,"#",1)==0) {
06628 continue;
06629 }
06630 readcount = sscanf(tableline, "%i %f %f", &numtypes, &table_spacing, &maxdist);
06631 if (readcount != 3) {
06632 NAMD_die("ERROR: Couldn't parse table header line\n");
06633 }
06634 break;
06635 }
06636
06637 if (maxdist < simParams->cutoff) {
06638 NAMD_die("Tabulated energies must at least extend to the cutoff distance\n");
06639 }
06640
06641 if (maxdist > simParams->cutoff) {
06642 iout << "Info: Reducing max table size to match nonbond cutoff\n";
06643 maxdist = ceil(simParams->cutoff);
06644 }
06645
06646
06647 numenerentries = 2 * numtypes * int (mynearbyint(maxdist/table_spacing) + 1);
06648
06649
06650 fprintf(stdout, "Table has %i entries\n",numenerentries);
06651
06652 if ( table_ener ) delete [] table_ener;
06653 table_ener = new BigReal[numenerentries];
06654 if ( table_types ) delete [] table_types;
06655 table_types = new char*[numtypes];
06656
06657
06658 int numtypesread = 0;
06659
06660 while(fgets(tableline,120,enertable)) {
06661 if (strncmp(tableline,"#",1)==0) {
06662 continue;
06663 }
06664 if (strncmp(tableline,"TYPE",4)==0) {
06665
06666 newtypename = new char[5];
06667 int readcount = sscanf(tableline, "%*s %s", newtypename);
06668 if (readcount != 1) {
06669 NAMD_die("Couldn't read interaction type from TYPE line\n");
06670 }
06671
06672 table_types[numtypesread] = newtypename;
06673
06674 if (numtypesread == numtypes) {
06675 NAMD_die("Error: Number of types in table doesn't match table header\n");
06676 }
06677
06678
06679 if (!strncasecmp(interp_type, "linear", 6)) {
06680 if (read_energy_type(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
06681 char err_msg[512];
06682 sprintf(err_msg, "Failed to read table energy (linear) type %s\n", newtypename);
06683 NAMD_die(err_msg);
06684 }
06685 } else if (!strncasecmp(interp_type, "cubic", 5)) {
06686 if (read_energy_type_bothcubspline(enertable, numtypesread, table_ener, table_spacing, maxdist) != 0) {
06687 char err_msg[512];
06688 sprintf(err_msg, "Failed to read table energy (cubic) type %s\n", newtypename);
06689 NAMD_die(err_msg);
06690 }
06691 } else {
06692 NAMD_die("Table interpolation type must be linear or cubic\n");
06693 }
06694
06695 numtypesread++;
06696 continue;
06697 }
06698
06699
06700
06701 }
06702
06703
06704 if (numtypesread != numtypes) {
06705 char err_msg[512];
06706 sprintf(err_msg, "ERROR: Expected %i tabulated energy types but got %i\n", numtypes, numtypesread);
06707 NAMD_die(err_msg);
06708 }
06709
06710
06711 simParams->tableNumTypes = numtypes;
06712 simParams->tableSpacing = table_spacing;
06713 simParams->tableMaxDist = maxdist;
06714 tablenumtypes = numtypes;
06715
06716
06717
06718
06719
06720
06721
06722
06723
06724
06725
06726
06727
06728
06729
06730
06731
06732
06733
06734
06735
06736
06737
06738
06739
06740
06741
06742
06743
06744
06745
06746
06747
06748
06749
06750
06751
06752
06753
06754
06755
06756
06757
06758
06759
06760
06761
06762
06763
06764
06765
06766
06767
06768
06769
06770
06771
06772
06773
06774
06775
06776
06777
06778
06779
06780
06781
06782
06783
06784
06785
06786
06787
06788
06789
06790
06791
06792
06793
06794
06795
06796
06797
06798
06799
06800
06801
06802
06803
06804
06805
06806
06807
06808
06809
06810
06811
06812
06813
06814
06815
06816
06817
06818
06819
06820
06821 }
06822
06823
06824
06825
06826
06827
06828
06829
06830
06831
06832
06833
06834
06835
06836
06837
06838
06839
06840
06841
06842
06843
06844 int Parameters::read_energy_type_bothcubspline(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
06845
06846 char tableline[120];
06847 int i,j;
06848
06849
06850 BigReal readdist;
06851 BigReal readener;
06852 BigReal readforce;
06853 BigReal spacing;
06854
06855 BigReal lastdist;
06856
06857
06858
06859
06860
06861
06862
06863 std::vector<BigReal> dists;
06864 std::vector<BigReal> enervalues;
06865 std::vector<BigReal> forcevalues;
06866 int numentries = 0;
06867
06868
06869
06870 BigReal currdist;
06871 int distbin;
06872 currdist = 0.0;
06873 lastdist = -1.0;
06874 distbin = 0;
06875
06876
06877 while(fgets(tableline,120,enertable) && distbin <= (int) (mynearbyint(maxdist / table_spacing) + 1)) {
06878 if (strncmp(tableline,"#",1)==0) {
06879 continue;
06880 }
06881 if (strncmp(tableline,"TYPE",4)==0) {
06882 fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR);
06883 break;
06884 }
06885
06886
06887 int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
06888
06889
06890 if (readcount != 3) {
06891 char err_msg[512];
06892 sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
06893 NAMD_die(err_msg);
06894 }
06895
06896
06897 if (readdist < lastdist) {
06898 NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
06899 }
06900
06901 lastdist = readdist;
06902 dists.push_back(readdist);
06903 enervalues.push_back(readener);
06904 forcevalues.push_back(readforce);
06905 numentries++;
06906 }
06907
06908
06909 if (dists[0] != 0.0) {
06910 NAMD_die("ERROR: First data point for energy table must be at r=0\n");
06911 }
06912 spacing = dists[1] - dists[0];
06913 for (i=2; i<(numentries - 1); i++) {
06914 BigReal myspacing;
06915 myspacing = dists[i] - dists[i-1];
06916 if (fabs(myspacing - spacing) > 0.00001) {
06917 printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
06918 NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
06919 }
06920 }
06921
06922
06923
06924
06925
06926 double* m = new double[numentries*numentries];
06927 double* xe = new double[numentries];
06928 double* xf = new double[numentries];
06929 double* be = new double[numentries];
06930 double* bf = new double[numentries];
06931
06932 be[0] = 3 * (enervalues[1] - enervalues[0]);
06933 for (i=1; i < (numentries - 1); i++) {
06934
06935 be[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
06936
06937 }
06938 be[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
06939
06940 bf[0] = 3 * (forcevalues[1] - forcevalues[0]);
06941 for (i=1; i < (numentries - 1); i++) {
06942
06943 bf[i] = 3 * (forcevalues[i+1] - forcevalues[i-1]);
06944
06945 }
06946 bf[numentries - 1] = 3 * (forcevalues[numentries - 1] - forcevalues[numentries - 2]);
06947
06948 memset(m,0,numentries*numentries*sizeof(double));
06949
06950
06951 m[0] = 2;
06952 for (i = 1; i < numentries; i++) {
06953 m[INDEX(numentries,i-1,i)] = 1;
06954 m[INDEX(numentries,i,i-1)] = 1;
06955 m[INDEX(numentries,i,i)] = 4;
06956 }
06957 m[INDEX(numentries,numentries-1,numentries-1)] = 2;
06958
06959
06960
06961
06962
06963 for (i=0; i<numentries; i++) {
06964
06965 const BigReal baseval = m[INDEX(numentries,i,i)];
06966 for (j=i+1; j<numentries; j++) {
06967 const BigReal myval = m[INDEX(numentries,j,i)];
06968 const BigReal factor = myval / baseval;
06969
06970 for (int k=i; k<numentries; k++) {
06971 const BigReal subval = m[INDEX(numentries,i,k)];
06972 m[INDEX(numentries,j,k)] -= (factor * subval);
06973 }
06974
06975 be[j] -= (factor * be[i]);
06976 bf[j] -= (factor * bf[i]);
06977
06978 }
06979 }
06980
06981
06982 for (i=numentries-1; i>=0; i--) {
06983
06984
06985 for (j=i+1; j<numentries; j++) {
06986 be[i] -= ( m[INDEX(numentries,i,j)] * xe[j] );
06987 bf[i] -= ( m[INDEX(numentries,i,j)] * xf[j] );
06988 }
06989
06990 xe[i] = be[i] / m[INDEX(numentries,i,i)];
06991 xf[i] = bf[i] / m[INDEX(numentries,i,i)];
06992
06993 }
06994
06995
06996
06997
06998 distbin = 0;
06999 int entriesperseg = (int) ceil(spacing / table_spacing);
07000 int table_prefix = 2 * typeindex * (int) (mynearbyint(maxdist / table_spacing) + 1);
07001
07002 for (i=0; i<numentries-1; i++) {
07003 BigReal Ae,Be,Ce,De;
07004 BigReal Af,Bf,Cf,Df;
07005 currdist = dists[i];
07006
07007
07008
07009
07010 Ae = enervalues[i];
07011 Be = xe[i];
07012 Ce = 3 * (enervalues[i+1] - enervalues[i]) - (2 * xe[i]) - xe[i+1];
07013 De = 2 * (enervalues[i] - enervalues[i+1]) + xe[i] + xe[i+1];
07014
07015 Af = forcevalues[i];
07016 Bf = xf[i];
07017 Cf = 3 * (forcevalues[i+1] - forcevalues[i]) - (2 * xf[i]) - xf[i+1];
07018 Df = 2 * (forcevalues[i] - forcevalues[i+1]) + xf[i] + xf[i+1];
07019
07020
07021 for (j=0; j<entriesperseg; j++) {
07022 const BigReal mydist = currdist + (j * table_spacing);
07023 const BigReal mydistfrac = (float) j / (entriesperseg - 1);
07024 distbin = (int) mynearbyint(mydist / table_spacing);
07025 if (distbin >= (int) mynearbyint(maxdist / table_spacing)) break;
07026 BigReal energy;
07027 BigReal force;
07028
07029 energy = Ae + (Be * mydistfrac) + (Ce * mydistfrac * mydistfrac) + (De * mydistfrac * mydistfrac * mydistfrac);
07030 force = Af + (Bf * mydistfrac) + (Cf * mydistfrac * mydistfrac) + (Df * mydistfrac * mydistfrac * mydistfrac);
07031
07032
07033 table_ener[table_prefix + 2 * distbin] = energy;
07034 table_ener[table_prefix + 2 * distbin + 1] = force;
07035 distbin++;
07036 }
07037 if (currdist >= maxdist) break;
07038 }
07039
07040
07041 distbin = (int) mynearbyint(maxdist / table_spacing);
07042
07043 table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
07044 table_ener[table_prefix + 2 * distbin + 1] = 0.0;
07045 distbin++;
07046
07047
07048
07049 delete m;
07050 delete xe;
07051 delete xf;
07052 delete be;
07053 delete bf;
07054 distbin--;
07055 printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (mynearbyint(maxdist / table_spacing)), maxdist, table_spacing);
07056 if (distbin != (int) (mynearbyint(maxdist / table_spacing))) return 1;
07057 return 0;
07058 }
07059
07060
07061
07062
07063
07064
07065
07066
07067
07068
07069
07070
07071
07072
07073
07074
07075
07076
07077
07078
07079
07080 int Parameters::read_energy_type_cubspline(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
07081
07082 char tableline[120];
07083 int i,j;
07084
07085
07086 BigReal readdist;
07087 BigReal readener;
07088 BigReal spacing;
07089
07090 BigReal lastdist;
07091
07092
07093
07094
07095
07096
07097
07098 std::vector<BigReal> dists;
07099 std::vector<BigReal> enervalues;
07100 int numentries = 0;
07101
07102
07103
07104 BigReal currdist;
07105 int distbin;
07106 currdist = 0.0;
07107 lastdist = -1.0;
07108 distbin = 0;
07109
07110
07111 while(fgets(tableline,120,enertable) && distbin <= (int) (mynearbyint(maxdist / table_spacing) + 1)) {
07112 if (strncmp(tableline,"#",1)==0) {
07113 continue;
07114 }
07115 if (strncmp(tableline,"TYPE",4)==0) {
07116 fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR);
07117 break;
07118 }
07119
07120
07121 int readcount = sscanf(tableline, "%lf %lf", &readdist, &readener);
07122
07123
07124 if (readcount != 2) {
07125 char err_msg[512];
07126 sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
07127 NAMD_die(err_msg);
07128 }
07129
07130
07131 if (readdist < lastdist) {
07132 NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
07133 }
07134
07135 lastdist = readdist;
07136 dists.push_back(readdist);
07137 enervalues.push_back(readener);
07138 numentries++;
07139 }
07140
07141
07142 if (dists[0] != 0.0) {
07143 NAMD_die("ERROR: First data point for energy table must be at r=0\n");
07144 }
07145 spacing = dists[1] - dists[0];
07146 for (i=2; i<(numentries - 1); i++) {
07147 BigReal myspacing;
07148 myspacing = dists[i] - dists[i-1];
07149 if (fabs(myspacing - spacing) > 0.00001) {
07150 printf("Bad spacing in table: %f should be %f (between distances %f and %f)\n", myspacing, spacing, dists[i-1], dists[i]);
07151 NAMD_die("ERROR: Nonuniform table encountered on cubic interpolation. Use a uniformly spaced table or switch to linear interpolation.\n");
07152 }
07153 }
07154
07155
07156
07157
07158 double* m = new double[numentries*numentries];
07159 double* x = new double[numentries];
07160 double* b = new double[numentries];
07161
07162 b[0] = 3 * (enervalues[1] - enervalues[0]);
07163 for (i=1; i < (numentries - 1); i++) {
07164 printf("Control point %i at %f\n", i, enervalues[i]);
07165 b[i] = 3 * (enervalues[i+1] - enervalues[i-1]);
07166 printf("b is %f\n", b[i]);
07167 }
07168 b[numentries - 1] = 3 * (enervalues[numentries - 1] - enervalues[numentries - 2]);
07169
07170 memset(m,0,numentries*numentries*sizeof(double));
07171
07172
07173 m[0] = 2;
07174 for (i = 1; i < numentries; i++) {
07175 m[INDEX(numentries,i-1,i)] = 1;
07176 m[INDEX(numentries,i,i-1)] = 1;
07177 m[INDEX(numentries,i,i)] = 4;
07178 }
07179 m[INDEX(numentries,numentries-1,numentries-1)] = 2;
07180
07181
07182
07183 printf("Solving the matrix equation: \n");
07184
07185 for (i=0; i<numentries; i++) {
07186 printf(" ( ");
07187 for (j=0; j<numentries; j++) {
07188 printf(" %6.3f,", m[INDEX(numentries, i, j)]);
07189 }
07190 printf(" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
07191 }
07192
07193
07194 for (i=0; i<numentries; i++) {
07195
07196 const BigReal baseval = m[INDEX(numentries,i,i)];
07197 for (j=i+1; j<numentries; j++) {
07198 const BigReal myval = m[INDEX(numentries,j,i)];
07199 const BigReal factor = myval / baseval;
07200
07201 for (int k=i; k<numentries; k++) {
07202 const BigReal subval = m[INDEX(numentries,i,k)];
07203 m[INDEX(numentries,j,k)] -= (factor * subval);
07204 }
07205
07206 b[j] -= (factor * b[i]);
07207
07208 }
07209 }
07210
07211 printf(" In upper diagonal form, equation is:\n");
07212 for (i=0; i<numentries; i++) {
07213 printf(" ( ");
07214 for (j=0; j<numentries; j++) {
07215 printf(" %6.3f,", m[INDEX(numentries, i, j)]);
07216 }
07217 printf(" ) ( D%-3i ) = ( %6.3f )\n", i, b[i]);
07218 }
07219
07220
07221 for (i=numentries-1; i>=0; i--) {
07222
07223
07224 for (j=i+1; j<numentries; j++) {
07225 b[i] -= ( m[INDEX(numentries,i,j)] * x[j] );
07226 }
07227
07228 x[i] = b[i] / m[INDEX(numentries,i,i)];
07229
07230 }
07231
07232 printf(" Solution vector is:\n\t(");
07233 for (i=0; i<numentries; i++) {
07234 printf(" %6.3f ", x[i]);
07235 }
07236 printf(" ) \n");
07237
07238
07239
07240
07241 distbin = 0;
07242 int entriesperseg = (int) ceil(spacing / table_spacing);
07243 int table_prefix = 2 * typeindex * (int) (mynearbyint(maxdist / table_spacing) + 1);
07244
07245 for (i=0; i<numentries-1; i++) {
07246 BigReal A,B,C,D;
07247 currdist = dists[i];
07248
07249 printf("Interpolating on interval %i\n", i);
07250
07251
07252 A = enervalues[i];
07253 B = x[i];
07254 C = 3 * (enervalues[i+1] - enervalues[i]) - (2 * x[i]) - x[i+1];
07255 D = 2 * (enervalues[i] - enervalues[i+1]) + x[i] + x[i+1];
07256
07257 printf("Coefficients for this interval: %f %f %f %f\n", A, B, C, D);
07258
07259
07260 for (j=0; j<entriesperseg; j++) {
07261 const BigReal mydist = currdist + (j * table_spacing);
07262 const BigReal mydistfrac = (float) j / (entriesperseg - 1);
07263 distbin = (int) mynearbyint(mydist / table_spacing);
07264 if (distbin >= (int) mynearbyint(maxdist / table_spacing)) break;
07265 BigReal energy;
07266 BigReal force;
07267
07268 energy = A + (B * mydistfrac) + (C * mydistfrac * mydistfrac) + (D * mydistfrac * mydistfrac * mydistfrac);
07269 force = B + (2 * C * mydistfrac) + (3 * D * mydistfrac * mydistfrac);
07270
07271 force *= (1.0 / spacing);
07272
07273 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);
07274 table_ener[table_prefix + 2 * distbin] = energy;
07275 table_ener[table_prefix + 2 * distbin + 1] = force;
07276 distbin++;
07277 }
07278 if (currdist >= maxdist) break;
07279 }
07280
07281
07282 distbin = (int) mynearbyint(maxdist / table_spacing);
07283 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));
07284 table_ener[table_prefix + 2 * distbin] = enervalues[numentries - 1];
07285 table_ener[table_prefix + 2 * distbin + 1] = 0.0;
07286 distbin++;
07287
07288
07289
07290 delete m;
07291 delete x;
07292 delete b;
07293 distbin--;
07294 printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (mynearbyint(maxdist / table_spacing)), maxdist, table_spacing);
07295 if (distbin != (int) (mynearbyint(maxdist / table_spacing))) return 1;
07296 return 0;
07297 }
07298
07299
07300
07301
07302
07303
07304
07305
07306
07307
07308
07309
07310
07311
07312
07313
07314
07315
07316
07317
07318
07319 int Parameters::read_energy_type(FILE* enertable, const int typeindex, BigReal* table_ener, const float table_spacing, const float maxdist) {
07320
07321 char tableline[120];
07322
07323
07324 BigReal readdist;
07325 BigReal readener;
07326 BigReal readforce;
07327 BigReal lastdist;
07328 BigReal lastener;
07329 BigReal lastforce;
07330 readdist = -1.0;
07331 readener = 0.0;
07332 readforce = 0.0;
07333
07334
07335 float currdist;
07336 int distbin;
07337 currdist = -1.0;
07338 distbin = -1;
07339
07340 while(fgets(tableline,120,enertable) && distbin <= (int) (mynearbyint(maxdist / table_spacing) + 1)) {
07341 printf("At distance %f + %f vs. %f\n", currdist, table_spacing, maxdist);
07342 if (strncmp(tableline,"#",1)==0) {
07343 continue;
07344 }
07345 if (strncmp(tableline,"TYPE",4)==0) {
07346 break;
07347 }
07348
07349
07350 lastdist = readdist;
07351 lastener = readener;
07352 lastforce = readforce;
07353 int readcount = sscanf(tableline, "%lf %lf %lf", &readdist, &readener, &readforce);
07354 if (distbin == -1) {
07355 if (readdist != 0.0) {
07356 NAMD_die("ERROR: Energy/force tables must start at d=0.0\n");
07357 } else {
07358 distbin = 0;
07359 continue;
07360 }
07361 }
07362
07363 if (readcount != 3) {
07364 char err_msg[512];
07365 sprintf(err_msg, "ERROR: Failed to parse table line %s!\n", tableline);
07366 NAMD_die(err_msg);
07367 }
07368
07369
07370 if (readdist < lastdist) {
07371 NAMD_die("ERROR: Encountered badly ordered entries in energy table!\n");
07372 }
07373
07374 currdist = lastdist;
07375
07376 while (currdist <= readdist && distbin <= (int) (mynearbyint(maxdist / table_spacing))) {
07377 distbin = (int) (mynearbyint(currdist / table_spacing));
07378 int table_loc = 2 * (distbin + (typeindex * (1 + (int) mynearbyint(maxdist / table_spacing))));
07379 printf("Doing interpolation for energy between %f %f and %f %f: Dist %f\n", readener, readdist, lastener, lastdist, currdist);
07380 table_ener[table_loc] = interp_lin(readener, lastener, readdist, lastdist, currdist);
07381 table_ener[table_loc + 1] = interp_lin(readforce, lastforce, readdist, lastdist, currdist);
07382 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);
07383 currdist += table_spacing;
07384 distbin++;
07385 }
07386 }
07387
07388
07389 fseek(enertable, -1 * (long) strlen(tableline), SEEK_CUR);
07390
07391
07392 distbin--;
07393 printf("Testing: %i vs %i (from %f / %f)\n", distbin, (int) (mynearbyint(maxdist / table_spacing)), maxdist, table_spacing);
07394 if (distbin != (int) (mynearbyint(maxdist / table_spacing))) return 1;
07395 return 0;
07396 }
07397
07398
07399
07400
07401
07402
07403
07404
07405
07406
07407
07408
07409
07410
07411
07412
07413 BigReal Parameters::interp_lin(BigReal val1, BigReal val2, BigReal end1, BigReal end2, BigReal currdist) {
07414
07415 BigReal m;
07416 BigReal val;
07417
07418 m = (val2 - val1) / (end2 - end1);
07419
07420 val = ((currdist-end1) * m + val1);
07421 return val;
07422 }
07423
07424
07425
07426
07427
07428
07429
07430
07431
07432
07433
07434
07435
07436
07437 int Parameters::get_int_table_type(char* tabletype) {
07438 for (int i=0; i<tablenumtypes; i++) {
07439 if (!strncmp(tabletype, table_types[i], 5)) {
07440 return i;
07441 }
07442 }
07443
07444 return -1;
07445 }
07446
07447