00001
00008
00009
00010
00011
00012
00013
00014
00015 #include <stdio.h>
00016 #include <string.h>
00017 #include <stdlib.h>
00018 #include <ctype.h>
00019
00020 #include "InfoStream.h"
00021 #include "Molecule.h"
00022 #include "strlib.h"
00023 #include "MStream.h"
00024 #include "Communicate.h"
00025
00026 #include "ObjectArena.h"
00027 #include "Parameters.h"
00028 #include "PDB.h"
00029 #include "SimParameters.h"
00030 #include "Hydrogen.h"
00031 #include "UniqueSetIter.h"
00032 #include "charm++.h"
00033
00034 #include "GridForceGrid.h"
00035
00036
00037 #define MIN_DEBUG_LEVEL 3
00038
00039 #include "Debug.h"
00040
00041 #include "CompressPsf.h"
00042 #include <deque>
00043 #include <algorithm>
00044 using namespace std;
00045
00046 #ifdef MEM_OPT_VERSION
00047 template int lookupCstPool<AtomSignature>(const vector<AtomSignature>&, const AtomSignature&);
00048 template int lookupCstPool<ExclusionSignature>(const vector<ExclusionSignature>&, const ExclusionSignature&);
00049 #endif
00050
00051 class ResidueLookupElem
00052 {
00053 public:
00054 char mySegid[11];
00055 ResidueLookupElem *next;
00056 int firstResid;
00057 int lastResid;
00058 ResizeArray<int> atomIndex;
00059
00060 ResidueLookupElem(void) { next = 0; firstResid = -1; lastResid = -1; }
00061 ~ResidueLookupElem(void) { delete next; }
00062 int lookup(const char *segid, int resid, int *begin, int *end) const;
00063 ResidueLookupElem* append(const char *segid, int resid, int aid);
00064 };
00065
00066 int ResidueLookupElem::lookup(
00067 const char *segid, int resid, int *begin, int *end) const {
00068 const ResidueLookupElem *elem = this;
00069 int rval = -1;
00070 while ( elem && strcasecmp(elem->mySegid,segid) ) elem = elem->next;
00071 if ( elem && (resid >= elem->firstResid) && (resid <= elem->lastResid) ) {
00072 *begin = elem->atomIndex[resid - elem->firstResid];
00073 *end = elem->atomIndex[resid - elem->firstResid + 1];
00074 rval = 0;
00075 }
00076 return rval;
00077 }
00078
00079 ResidueLookupElem* ResidueLookupElem::append(
00080 const char *segid, int resid, int aid) {
00081 ResidueLookupElem *rval = this;
00082 if ( firstResid == -1 ) {
00083 strcpy(mySegid,segid);
00084 firstResid = resid;
00085 lastResid = resid;
00086 atomIndex.add(aid);
00087 atomIndex.add(aid+1);
00088 } else if ( ! strcasecmp(mySegid,segid) ) {
00089 if ( resid == lastResid ) {
00090 atomIndex[lastResid - firstResid + 1] = aid + 1;
00091 } else if ( resid < lastResid ) {
00092
00093 iout << iWARN << "Residue " << resid <<
00094 " out of order in segment " << segid <<
00095 ", lookup for additional residues in this segment disabled.\n" << endi;
00096 rval = next = new ResidueLookupElem;
00097 next->append(segid,resid,aid);
00098 } else {
00099 for ( ; lastResid < resid; ++lastResid ) atomIndex.add(aid);
00100 atomIndex[lastResid - firstResid + 1] = aid + 1;
00101 }
00102 } else {
00103 rval = next = new ResidueLookupElem;
00104 next->append(segid,resid,aid);
00105 }
00106 return rval;
00107 }
00108
00109
00110
00111 int Molecule::get_atom_from_name(
00112 const char *segid, int resid, const char *aname) const {
00113
00114 if (atomNames == NULL || resLookup == NULL)
00115 {
00116 NAMD_die("Tried to find atom from name on node other than node 0");
00117 }
00118
00119 int i = 0;
00120 int end = 0;
00121 if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
00122 for ( ; i < end; ++i ) {
00123 #ifdef MEM_OPT_VERSION
00124 Index idx = atomNames[i].atomnameIdx;
00125 if(!strcasecmp(aname, atomNamePool[idx])) return i;
00126 #else
00127 if ( ! strcasecmp(aname,atomNames[i].atomname) ) return i;
00128 #endif
00129 }
00130 return -1;
00131 }
00132
00133
00134 int Molecule::get_residue_size(
00135 const char *segid, int resid) const {
00136
00137 if (atomNames == NULL || resLookup == NULL)
00138 {
00139 NAMD_die("Tried to find atom from name on node other than node 0");
00140 }
00141 int i = 0;
00142 int end = 0;
00143 if ( resLookup->lookup(segid,resid,&i,&end) ) return 0;
00144 return ( end - i );
00145 }
00146
00147
00148 int Molecule::get_atom_from_index_in_residue(
00149 const char *segid, int resid, int index) const {
00150
00151 if (atomNames == NULL || resLookup == NULL)
00152 {
00153 NAMD_die("Tried to find atom from name on node other than node 0");
00154 }
00155 int i = 0;
00156 int end = 0;
00157 if ( resLookup->lookup(segid,resid,&i,&end) ) return -1;
00158 if ( index >= 0 && index < ( end - i ) ) return ( index + i );
00159 return -1;
00160 }
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 void Molecule::initialize(SimParameters *simParams, Parameters *param)
00174 {
00175 if ( sizeof(int32) != 4 ) { NAMD_bug("sizeof(int32) != 4"); }
00176 this->simParams = simParams;
00177 this->params = param;
00178
00179 atoms=NULL;
00180 atomNames=NULL;
00181 resLookup=NULL;
00182 if ( simParams->globalForcesOn ) {
00183 resLookup = new ResidueLookupElem;
00184 }
00185
00186 #ifdef MEM_OPT_VERSION
00187 eachAtomSig = NULL;
00188 atomSigPoolSize = 0;
00189 atomSigPool = NULL;
00190 massPoolSize = 0;
00191 atomMassPool = NULL;
00192 eachAtomMass = NULL;
00193 chargePoolSize = 0;
00194 atomChargePool = NULL;
00195 eachAtomCharge = NULL;
00196 #else
00197 bonds=NULL;
00198 angles=NULL;
00199 dihedrals=NULL;
00200 impropers=NULL;
00201 crossterms=NULL;
00202 #endif
00203
00204 donors=NULL;
00205 acceptors=NULL;
00206
00207
00208 #ifndef MEM_OPT_VERSION
00209 tmpArena=NULL;
00210 exclusions=NULL;
00211 bondsWithAtom=NULL;
00212 bondsByAtom=NULL;
00213 anglesByAtom=NULL;
00214 dihedralsByAtom=NULL;
00215 impropersByAtom=NULL;
00216 crosstermsByAtom=NULL;
00217 #endif
00218
00219 #ifdef MEM_OPT_VERSION
00220 exclSigPool = NULL;
00221 exclChkSigPool = NULL;
00222 exclSigPoolSize = 0;
00223 eachAtomExclSig = NULL;
00224 #else
00225 exclusionsByAtom=NULL;
00226 fullExclusionsByAtom=NULL;
00227 modExclusionsByAtom=NULL;
00228 all_exclusions=NULL;
00229 #endif
00230
00231 langevinParams=NULL;
00232 fixedAtomFlags=NULL;
00233
00234 #ifdef MEM_OPT_VERSION
00235 clusterSigs=NULL;
00236 #else
00237 cluster=NULL;
00238 clusterSize=NULL;
00239 #endif
00240
00241 exPressureAtomFlags=NULL;
00242 rigidBondLengths=NULL;
00243 consIndexes=NULL;
00244 consParams=NULL;
00245
00246 gridfrcIndexes=NULL;
00247 gridfrcParams=NULL;
00248 gridfrcGrid=NULL;
00249
00250 stirIndexes=NULL;
00251 stirParams=NULL;
00252 movDragIndexes=NULL;
00253 movDragParams=NULL;
00254 rotDragIndexes=NULL;
00255 rotDragParams=NULL;
00256 consTorqueIndexes=NULL;
00257 consTorqueParams=NULL;
00258 consForceIndexes=NULL;
00259 consForce=NULL;
00260
00261 fepAtomFlags=NULL;
00262
00263
00264 nameArena = new ObjectArena<char>;
00265
00266
00267 #ifndef MEM_OPT_VERSION
00268 arena = new ObjectArena<int32>;
00269 exclArena = new ObjectArena<char>;
00270 #endif
00271
00272
00273
00274 numAtoms=0;
00275 numRealBonds=0;
00276 numBonds=0;
00277 numAngles=0;
00278 numDihedrals=0;
00279 numImpropers=0;
00280 numCrossterms=0;
00281 numDonors=0;
00282 numAcceptors=0;
00283 numExclusions=0;
00284 numLP=0;
00285 numConstraints=0;
00286 numStirredAtoms=0;
00287 numMovDrag=0;
00288 numRotDrag=0;
00289 numConsTorque=0;
00290 numConsForce=0;
00291 numFixedAtoms=0;
00292 numFixedGroups=0;
00293 numExPressureAtoms=0;
00294 numRigidBonds=0;
00295 numFixedRigidBonds=0;
00296 numMultipleDihedrals=0;
00297 numMultipleImpropers=0;
00298 numCalcBonds=0;
00299 numCalcAngles=0;
00300 numCalcDihedrals=0;
00301 numCalcImpropers=0;
00302 numCalcCrossterms=0;
00303 numCalcExclusions=0;
00304
00305
00306 numFepInitial = 0;
00307 numFepFinal = 0;
00308
00309
00310 }
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 Molecule::Molecule(SimParameters *simParams, Parameters *param)
00323 {
00324 initialize(simParams,param);
00325 }
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 Molecule::Molecule(SimParameters *simParams, Parameters *param, char *filename, ConfigList *cfgList)
00336 {
00337 initialize(simParams,param);
00338
00339 if(simParams->useCompressedPsf)
00340 read_compressed_psf_file(filename, param);
00341 else if(simParams->genCompressedPsf){
00342 compress_psf_file(this, filename, param, simParams, cfgList);
00343 }
00344 else
00345 read_psf_file(filename, param);
00346 }
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360 Molecule::~Molecule()
00361 {
00362
00363
00364 if (atoms != NULL)
00365 delete [] atoms;
00366
00367 if (atomNames != NULL)
00368 {
00369
00370 delete [] atomNames;
00371 }
00372 delete nameArena;
00373
00374 if (resLookup != NULL)
00375 delete resLookup;
00376
00377 #ifdef MEM_OPT_VERSION
00378 if(eachAtomSig) delete [] eachAtomSig;
00379 if(atomSigPool) delete [] atomSigPool;
00380 #else
00381 if (bonds != NULL)
00382 delete [] bonds;
00383
00384 if (angles != NULL)
00385 delete [] angles;
00386
00387 if (dihedrals != NULL)
00388 delete [] dihedrals;
00389
00390 if (impropers != NULL)
00391 delete [] impropers;
00392
00393 if (crossterms != NULL)
00394 delete [] crossterms;
00395
00396 if (exclusions != NULL)
00397 delete [] exclusions;
00398 #endif
00399
00400 if (donors != NULL)
00401 delete [] donors;
00402
00403 if (acceptors != NULL)
00404 delete [] acceptors;
00405
00406 #ifdef MEM_OPT_VERSION
00407 if(exclSigPool) delete [] exclSigPool;
00408 if(exclChkSigPool) delete [] exclChkSigPool;
00409 if(eachAtomExclSig) delete [] eachAtomExclSig;
00410 #else
00411 if (bondsByAtom != NULL)
00412 delete [] bondsByAtom;
00413
00414 if (anglesByAtom != NULL)
00415 delete [] anglesByAtom;
00416
00417 if (dihedralsByAtom != NULL)
00418 delete [] dihedralsByAtom;
00419
00420 if (impropersByAtom != NULL)
00421 delete [] impropersByAtom;
00422
00423 if (crosstermsByAtom != NULL)
00424 delete [] crosstermsByAtom;
00425
00426 if (exclusionsByAtom != NULL)
00427 delete [] exclusionsByAtom;
00428
00429 if (fullExclusionsByAtom != NULL)
00430 delete [] fullExclusionsByAtom;
00431
00432 if (modExclusionsByAtom != NULL)
00433 delete [] modExclusionsByAtom;
00434
00435 if (all_exclusions != NULL)
00436 delete [] all_exclusions;
00437 #endif
00438
00439
00440 if (fixedAtomFlags != NULL)
00441 delete [] fixedAtomFlags;
00442
00443 if (stirIndexes != NULL)
00444 delete [] stirIndexes;
00445
00446
00447 #ifdef MEM_OPT_VERSION
00448 if(clusterSigs != NULL){
00449 delete [] clusterSigs;
00450 }
00451 #else
00452 if (cluster != NULL)
00453 delete [] cluster;
00454 if (clusterSize != NULL)
00455 delete [] clusterSize;
00456 #endif
00457
00458 if (exPressureAtomFlags != NULL)
00459 delete [] exPressureAtomFlags;
00460
00461 if (rigidBondLengths != NULL)
00462 delete [] rigidBondLengths;
00463
00464
00465 if (fepAtomFlags != NULL)
00466 delete [] fepAtomFlags;
00467
00468
00469
00470 #ifndef MEM_OPT_VERSION
00471 delete arena;
00472 delete exclArena;
00473 #endif
00474 }
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494 void Molecule::read_psf_file(char *fname, Parameters *params)
00495
00496 {
00497 #ifdef MEM_OPT_VERSION
00498 return;
00499 #else
00500 char err_msg[512];
00501 char buffer[512];
00502 int i;
00503 int NumTitle;
00504 FILE *psf_file;
00505 int ret_code;
00506
00507
00508 if ( (psf_file = Fopen(fname, "r")) == NULL)
00509 {
00510 sprintf(err_msg, "UNABLE TO OPEN .psf FILE %s", fname);
00511 NAMD_die(err_msg);
00512 }
00513
00514
00515 ret_code = NAMD_read_line(psf_file, buffer);
00516
00517 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00518 {
00519 ret_code = NAMD_read_line(psf_file, buffer);
00520 }
00521
00522
00523
00524 if (ret_code!=0)
00525 {
00526 sprintf(err_msg, "EMPTY .psf FILE %s", fname);
00527 NAMD_die(err_msg);
00528 }
00529
00530
00531
00532 if (!NAMD_find_word(buffer, "psf"))
00533 {
00534 sprintf(err_msg, "UNABLE TO FIND \"PSF\" STRING IN PSF FILE %s",
00535 fname);
00536 NAMD_die(err_msg);
00537 }
00538
00539
00540 ret_code = NAMD_read_line(psf_file, buffer);
00541
00542 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00543 {
00544 ret_code = NAMD_read_line(psf_file, buffer);
00545 }
00546
00547
00548
00549
00550 if (ret_code!=0)
00551 {
00552 sprintf(err_msg, "MISSING EVERYTHING BUT PSF FROM %s", fname);
00553 NAMD_die(err_msg);
00554 }
00555
00556
00557
00558 if (!NAMD_find_word(buffer, "NTITLE"))
00559 {
00560 sprintf(err_msg,"CAN NOT FIND \"NTITLE\" STRING IN PSF FILE %s",
00561 fname);
00562 NAMD_die(err_msg);
00563 }
00564
00565 sscanf(buffer, "%d", &NumTitle);
00566
00567
00568
00569 i=0;
00570
00571 while ( ((ret_code=NAMD_read_line(psf_file, buffer)) == 0) &&
00572 (i<NumTitle) )
00573 {
00574 if (!NAMD_blank_string(buffer))
00575 i++;
00576 }
00577
00578
00579 if (ret_code!=0)
00580 {
00581 sprintf(err_msg, "FOUND EOF INSTEAD OF NATOM IN PSF FILE %s",
00582 fname);
00583 NAMD_die(err_msg);
00584 }
00585
00586 while (NAMD_blank_string(buffer))
00587 {
00588 NAMD_read_line(psf_file, buffer);
00589 }
00590
00591
00592 if (!NAMD_find_word(buffer, "NATOM"))
00593 {
00594 sprintf(err_msg, "DIDN'T FIND \"NATOM\" IN PSF FILE %s",
00595 fname);
00596 NAMD_die(err_msg);
00597 }
00598
00599
00600 sscanf(buffer, "%d", &numAtoms);
00601
00602 read_atoms(psf_file, params);
00603
00604
00605 ret_code = NAMD_read_line(psf_file, buffer);
00606
00607 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00608 {
00609 ret_code = NAMD_read_line(psf_file, buffer);
00610 }
00611
00612
00613 if (ret_code != 0)
00614 {
00615 NAMD_die("EOF ENCOUNTERED LOOKING FOR NBONDS IN PSF");
00616 }
00617
00618
00619 if (!NAMD_find_word(buffer, "NBOND"))
00620 {
00621 NAMD_die("DID NOT FIND NBOND AFTER ATOM LIST IN PSF");
00622 }
00623
00624
00625 sscanf(buffer, "%d", &numBonds);
00626
00627 if (numBonds)
00628 read_bonds(psf_file, params);
00629
00630
00631 ret_code = NAMD_read_line(psf_file, buffer);
00632
00633 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00634 {
00635 ret_code = NAMD_read_line(psf_file, buffer);
00636 }
00637
00638
00639 if (ret_code != 0)
00640 {
00641 NAMD_die("EOF ENCOUNTERED LOOKING FOR NTHETA IN PSF");
00642 }
00643
00644
00645 if (!NAMD_find_word(buffer, "NTHETA"))
00646 {
00647 NAMD_die("DID NOT FIND NTHETA AFTER BOND LIST IN PSF");
00648 }
00649
00650
00651 sscanf(buffer, "%d", &numAngles);
00652
00653 if (numAngles)
00654 read_angles(psf_file, params);
00655
00656
00657 ret_code = NAMD_read_line(psf_file, buffer);
00658
00659 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00660 {
00661 ret_code = NAMD_read_line(psf_file, buffer);
00662 }
00663
00664
00665 if (ret_code != 0)
00666 {
00667 NAMD_die("EOF ENCOUNTERED LOOKING FOR NPHI IN PSF");
00668 }
00669
00670
00671 if (!NAMD_find_word(buffer, "NPHI"))
00672 {
00673 NAMD_die("DID NOT FIND NPHI AFTER ANGLE LIST IN PSF");
00674 }
00675
00676
00677 sscanf(buffer, "%d", &numDihedrals);
00678
00679 if (numDihedrals)
00680 read_dihedrals(psf_file, params);
00681
00682
00683 ret_code = NAMD_read_line(psf_file, buffer);
00684
00685 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00686 {
00687 ret_code = NAMD_read_line(psf_file, buffer);
00688 }
00689
00690
00691 if (ret_code != 0)
00692 {
00693 NAMD_die("EOF ENCOUNTERED LOOKING FOR NIMPHI IN PSF");
00694 }
00695
00696
00697 if (!NAMD_find_word(buffer, "NIMPHI"))
00698 {
00699 NAMD_die("DID NOT FIND NIMPHI AFTER ATOM LIST IN PSF");
00700 }
00701
00702
00703 sscanf(buffer, "%d", &numImpropers);
00704
00705 if (numImpropers)
00706 read_impropers(psf_file, params);
00707
00708
00709 ret_code = NAMD_read_line(psf_file, buffer);
00710
00711 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00712 {
00713 ret_code = NAMD_read_line(psf_file, buffer);
00714 }
00715
00716
00717 if (ret_code != 0)
00718 {
00719 NAMD_die("EOF ENCOUNTERED LOOKING FOR NDON IN PSF");
00720 }
00721
00722
00723 if (!NAMD_find_word(buffer, "NDON"))
00724 {
00725 NAMD_die("DID NOT FIND NDON AFTER ATOM LIST IN PSF");
00726 }
00727
00728
00729 sscanf(buffer, "%d", &numDonors);
00730
00731 if (numDonors)
00732 read_donors(psf_file);
00733
00734
00735 ret_code = NAMD_read_line(psf_file, buffer);
00736
00737 while ( (ret_code==0) && (NAMD_blank_string(buffer)) )
00738 {
00739 ret_code = NAMD_read_line(psf_file, buffer);
00740 }
00741
00742
00743 if (ret_code != 0)
00744 {
00745 NAMD_die("EOF ENCOUNTERED LOOKING FOR NACC IN PSF");
00746 }
00747
00748
00749 if (!NAMD_find_word(buffer, "NACC"))
00750 {
00751 NAMD_die("DID NOT FIND NACC AFTER ATOM LIST IN PSF");
00752 }
00753
00754
00755 sscanf(buffer, "%d", &numAcceptors);
00756
00757 if (numAcceptors)
00758 read_acceptors(psf_file);
00759
00760
00761 while (!NAMD_find_word(buffer, "NNB"))
00762 {
00763 ret_code = NAMD_read_line(psf_file, buffer);
00764
00765 if (ret_code != 0)
00766 {
00767 NAMD_die("EOF ENCOUNTERED LOOKING FOR NNB IN PSF FILE");
00768 }
00769 }
00770
00771
00772 sscanf(buffer, "%d", &numExclusions);
00773
00774 if (numExclusions)
00775 read_exclusions(psf_file);
00776
00777
00778 int crossterms_present = 1;
00779 while (!NAMD_find_word(buffer, "NCRTERM"))
00780 {
00781 ret_code = NAMD_read_line(psf_file, buffer);
00782
00783 if (ret_code != 0)
00784 {
00785
00786 crossterms_present = 0;
00787 break;
00788 }
00789 }
00790
00791 if ( crossterms_present) {
00792
00793
00794 sscanf(buffer, "%d", &numCrossterms);
00795
00796 if (numCrossterms)
00797 read_crossterms(psf_file, params);
00798
00799 }
00800
00801
00802 Fclose(psf_file);
00803
00804
00805 numRealBonds = numBonds;
00806 build_atom_status();
00807
00808 return;
00809 #endif
00810 }
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829 void Molecule::read_compressed_psf_file(char *fname, Parameters *params){
00830 #ifndef MEM_OPT_VERSION
00831 return;
00832 #else
00833 FILE *psf_file;
00834 int ret_code;
00835 char buffer[512];
00836
00837
00838 if ( (psf_file = Fopen(fname, "r")) == NULL)
00839 {
00840 char err_msg[512];
00841 sprintf(err_msg, "UNABLE TO OPEN THE COMPRESSED .psf FILE %s", fname);
00842 NAMD_die(err_msg);
00843 }
00844
00845 char strBuf[12];
00846
00847
00848
00849 NAMD_read_line(psf_file, buffer);
00850 if(!NAMD_find_word(buffer, "NSEGMENTNAMES"))
00851 NAMD_die("UNABLE TO FIND NSEGMENTNAMES");
00852 sscanf(buffer, "%d", &segNamePoolSize);
00853 if(segNamePoolSize!=0)
00854 segNamePool = new char *[segNamePoolSize];
00855 for(int i=0; i<segNamePoolSize; i++){
00856 NAMD_read_line(psf_file, buffer);
00857 sscanf(buffer, "%s", strBuf);
00858 segNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
00859 strcpy(segNamePool[i], strBuf);
00860 }
00861
00862
00863 NAMD_read_line(psf_file, buffer);
00864 if(!NAMD_find_word(buffer, "NRESIDUENAMES"))
00865 NAMD_die("UNABLE TO FIND NRESIDUENAMES");
00866 sscanf(buffer, "%d", &resNamePoolSize);
00867 if(resNamePoolSize!=0)
00868 resNamePool = new char *[resNamePoolSize];
00869 for(int i=0; i<resNamePoolSize; i++){
00870 NAMD_read_line(psf_file, buffer);
00871 sscanf(buffer, "%s", strBuf);
00872 resNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
00873 strcpy(resNamePool[i], strBuf);
00874 }
00875
00876
00877 NAMD_read_line(psf_file, buffer);
00878 if(!NAMD_find_word(buffer, "NATOMNAMES"))
00879 NAMD_die("UNABLE TO FIND NATOMNAMES");
00880 sscanf(buffer, "%d", &atomNamePoolSize);
00881 if(atomNamePoolSize!=0)
00882 atomNamePool = new char *[atomNamePoolSize];
00883 for(int i=0; i<atomNamePoolSize; i++){
00884 NAMD_read_line(psf_file, buffer);
00885 sscanf(buffer, "%s", strBuf);
00886 atomNamePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
00887 strcpy(atomNamePool[i], strBuf);
00888 }
00889
00890
00891 NAMD_read_line(psf_file, buffer);
00892 if(!NAMD_find_word(buffer, "NATOMTYPES"))
00893 NAMD_die("UNABLE TO FIND NATOMTYPES");
00894 sscanf(buffer, "%d", &atomTypePoolSize);
00895 if(atomTypePoolSize!=0)
00896 atomTypePool = new char *[atomTypePoolSize];
00897 for(int i=0; i<atomTypePoolSize; i++){
00898 NAMD_read_line(psf_file, buffer);
00899 sscanf(buffer, "%s", strBuf);
00900 atomTypePool[i] = nameArena->getNewArray(strlen(strBuf)+1);
00901 strcpy(atomTypePool[i], strBuf);
00902 }
00903
00904
00905 NAMD_read_line(psf_file, buffer);
00906 if(!NAMD_find_word(buffer, "NCHARGES"))
00907 NAMD_die("UNABLE TO FIND NCHARGES");
00908 sscanf(buffer, "%d", &chargePoolSize);
00909 if(chargePoolSize!=0)
00910 atomChargePool = new Real[chargePoolSize];
00911 for(int i=0; i<chargePoolSize; i++){
00912 NAMD_read_line(psf_file, buffer);
00913 sscanf(buffer, "%f", atomChargePool+i);
00914 }
00915
00916
00917 NAMD_read_line(psf_file, buffer);
00918 if(!NAMD_find_word(buffer, "NMASSES"))
00919 NAMD_die("UNABLE TO FIND NMASSES");
00920 sscanf(buffer, "%d", &massPoolSize);
00921 if(massPoolSize!=0)
00922 atomMassPool = new Real[massPoolSize];
00923 for(int i=0; i<massPoolSize; i++){
00924 NAMD_read_line(psf_file, buffer);
00925 sscanf(buffer, "%f", atomMassPool+i);
00926 }
00927
00928
00929 NAMD_read_line(psf_file, buffer);
00930 if(!NAMD_find_word(buffer, "ATOMSIGS"))
00931 NAMD_die("UNABLE TO FIND ATOMSIGS");
00932 sscanf(buffer, "%d", &atomSigPoolSize);
00933 atomSigPool = new AtomSignature[atomSigPoolSize];
00934 int typeCnt;
00935 int tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
00936 int tisReal;
00937 int ttype;
00938 for(int i=0; i<atomSigPoolSize; i++){
00939
00940 NAMD_read_line(psf_file, buffer);
00941 if(!NAMD_find_word(buffer, "NBONDSIGS"))
00942 NAMD_die("UNABLE TO FIND NBONDSIGS");
00943 sscanf(buffer, "%d", &typeCnt);
00944 if(typeCnt!=0){
00945 atomSigPool[i].bondCnt = typeCnt;
00946 atomSigPool[i].bondSigs = new TupleSignature[typeCnt];
00947 }
00948 for(int j=0; j<typeCnt; j++){
00949 NAMD_read_line(psf_file, buffer);
00950 sscanf(buffer, "%d | %d | %d", &tmp1, &ttype, &tisReal);
00951 TupleSignature oneSig(1, BOND, (Index)ttype, (char)tisReal);
00952 oneSig.offset[0] = tmp1;
00953 atomSigPool[i].bondSigs[j]=oneSig;
00954 if(tisReal) numRealBonds++;
00955 }
00956
00957
00958 NAMD_read_line(psf_file, buffer);
00959 if(!NAMD_find_word(buffer, "NTHETASIGS"))
00960 NAMD_die("UNABLE TO FIND NTHETASIGS");
00961 sscanf(buffer, "%d", &typeCnt);
00962 if(typeCnt!=0){
00963 atomSigPool[i].angleCnt = typeCnt;
00964 atomSigPool[i].angleSigs = new TupleSignature[typeCnt];
00965 }
00966 for(int j=0; j<typeCnt; j++){
00967 NAMD_read_line(psf_file, buffer);
00968 sscanf(buffer, "%d %d | %d | %d", &tmp1, &tmp2, &ttype, &tisReal);
00969 TupleSignature oneSig(2,ANGLE,(Index)ttype, (char)tisReal);
00970 oneSig.offset[0] = tmp1;
00971 oneSig.offset[1] = tmp2;
00972 atomSigPool[i].angleSigs[j] = oneSig;
00973 }
00974
00975 NAMD_read_line(psf_file, buffer);
00976 if(!NAMD_find_word(buffer, "NPHISIGS"))
00977 NAMD_die("UNABLE TO FIND NPHISIGS");
00978 sscanf(buffer, "%d", &typeCnt);
00979 if(typeCnt!=0){
00980 atomSigPool[i].dihedralCnt = typeCnt;
00981 atomSigPool[i].dihedralSigs = new TupleSignature[typeCnt];
00982 }
00983 for(int j=0; j<typeCnt; j++){
00984 NAMD_read_line(psf_file, buffer);
00985 sscanf(buffer, "%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
00986 TupleSignature oneSig(3,DIHEDRAL,(Index)ttype, (char)tisReal);
00987 oneSig.offset[0] = tmp1;
00988 oneSig.offset[1] = tmp2;
00989 oneSig.offset[2] = tmp3;
00990 atomSigPool[i].dihedralSigs[j] = oneSig;
00991 }
00992
00993 NAMD_read_line(psf_file, buffer);
00994 if(!NAMD_find_word(buffer, "NIMPHISIGS"))
00995 NAMD_die("UNABLE TO FIND NIMPHISIGS");
00996 sscanf(buffer, "%d", &typeCnt);
00997 if(typeCnt!=0){
00998 atomSigPool[i].improperCnt = typeCnt;
00999 atomSigPool[i].improperSigs = new TupleSignature[typeCnt];
01000 }
01001 for(int j=0; j<typeCnt; j++){
01002 NAMD_read_line(psf_file, buffer);
01003 sscanf(buffer, "%d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &ttype, &tisReal);
01004 TupleSignature oneSig(3,IMPROPER,(Index)ttype, (char)tisReal);
01005 oneSig.offset[0] = tmp1;
01006 oneSig.offset[1] = tmp2;
01007 oneSig.offset[2] = tmp3;
01008 atomSigPool[i].improperSigs[j] = oneSig;
01009 }
01010
01011 NAMD_read_line(psf_file, buffer);
01012 if(!NAMD_find_word(buffer, "NCRTERMSIGS"))
01013 NAMD_die("UNABLE TO FIND NCRTERMSIGS");
01014 sscanf(buffer, "%d", &typeCnt);
01015 if(typeCnt!=0){
01016 atomSigPool[i].crosstermCnt = typeCnt;
01017 atomSigPool[i].crosstermSigs = new TupleSignature[typeCnt];
01018 }
01019 for(int j=0; j<typeCnt; j++){
01020 NAMD_read_line(psf_file, buffer);
01021 sscanf(buffer, "%d %d %d %d %d %d %d | %d | %d", &tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6, &tmp7, &ttype, &tisReal);
01022 TupleSignature oneSig(7,CROSSTERM,(Index)ttype, (char)tisReal);
01023 oneSig.offset[0] = tmp1;
01024 oneSig.offset[1] = tmp2;
01025 oneSig.offset[2] = tmp3;
01026 oneSig.offset[3] = tmp4;
01027 oneSig.offset[4] = tmp5;
01028 oneSig.offset[5] = tmp6;
01029 oneSig.offset[6] = tmp7;
01030 atomSigPool[i].crosstermSigs[j] = oneSig;
01031 }
01032 }
01033
01034
01035 NAMD_read_line(psf_file, buffer);
01036 if(!NAMD_find_word(buffer, "NEXCLSIGS")){
01037 NAMD_die("UNABLE TO FIND NEXCLSIGS");
01038 }
01039 sscanf(buffer, "%d", &exclSigPoolSize);
01040 if(exclSigPoolSize>0) exclSigPool = new ExclusionSignature[exclSigPoolSize];
01041 vector<int> fullExcls;
01042 vector<int> modExcls;
01043 for(int i=0; i<exclSigPoolSize; i++){
01044 int fullExclCnt = NAMD_read_int(psf_file, buffer);
01045 for(int j=0; j<fullExclCnt; j++)
01046 fullExcls.push_back(NAMD_read_int(psf_file, buffer));
01047 int modExclCnt = NAMD_read_int(psf_file, buffer);
01048 for(int j=0; j<modExclCnt; j++)
01049 modExcls.push_back(NAMD_read_int(psf_file, buffer));
01050
01051
01052
01053 exclSigPool[i].setOffsets(fullExcls, modExcls);
01054
01055 fullExcls.clear();
01056 modExcls.clear();
01057 }
01058
01059
01060
01061
01062 NAMD_read_line(psf_file, buffer);
01063 if(!NAMD_find_word(buffer, "NATOM"))
01064 NAMD_die("UNABLE TO FIND NATOM");
01065 sscanf(buffer, "%d", &numAtoms);
01066
01067 if(numAtoms>0){
01068 atoms = new AtomCstInfo[numAtoms];
01069 atomNames = new AtomNameIdx[numAtoms];
01070 eachAtomMass = new Index[numAtoms];
01071 eachAtomCharge = new Index[numAtoms];
01072 eachAtomSig = new Index[numAtoms];
01073 eachAtomExclSig = new Index[numAtoms];
01074
01075 clusterSigs = new int[numAtoms];
01076
01077 hydrogenGroup.resize(0);
01078 ResidueLookupElem *tmpResLookup = resLookup;
01079
01080 int residue_number;
01081 char *segment_name;
01082 int read_count;
01083 int idx[10];
01084
01085
01086 for(int i=0; i<numAtoms; i++){
01087 NAMD_read_line(psf_file, buffer);
01088 read_count = sscanf(buffer, "%d %d %d %d %d %d %d %d %d %d",
01089 idx, idx+1, idx+2, idx+3, idx+4,
01090 idx+5, idx+6, idx+7, idx+8, idx+9);
01091 if(read_count!=10){
01092 char err_msg[128];
01093 sprintf(err_msg, "BAD ATOM LINE FORMAT IN COMPRESSED PSF FILE IN ATOM LINE %d\nLINE=%s",i+1, buffer);
01094 NAMD_die(err_msg);
01095 }
01096
01097 segment_name = segNamePool[idx[0]];
01098 residue_number = idx[1];
01099 atomNames[i].resnameIdx = (Index)idx[2];
01100 atomNames[i].atomnameIdx = (Index)idx[3];
01101 atomNames[i].atomtypeIdx = (Index)idx[4];
01102 eachAtomCharge[i] = (Index)idx[5];
01103 eachAtomMass[i] = (Index)idx[6];
01104 eachAtomSig[i] = (Index)idx[7];
01105 eachAtomExclSig[i] = (Index)idx[8];
01106
01107 clusterSigs[i] = idx[9];
01108
01109
01110
01111
01112 if(tmpResLookup) tmpResLookup =
01113 tmpResLookup->append(segment_name, residue_number, i);
01114
01115 Real thisAtomMass = atomMassPool[eachAtomMass[i]];
01116 if (thisAtomMass < 0.1) {
01117 atoms[i].status = LonepairAtom;
01118 } else if(thisAtomMass <= 3.5){
01119 atoms[i].status = HydrogenAtom;
01120 }else if(atomNamePool[atomNames[i].atomnameIdx][0]=='O' &&
01121 (thisAtomMass >= 14.0) && (thisAtomMass <= 18.0)){
01122 atoms[i].status = OxygenAtom;
01123 }
01124
01125
01126 params->assign_vdw_index(atomTypePool[atomNames[i].atomtypeIdx],
01127 &(atoms[i]));
01128 }
01129
01130
01131 }
01132
01133 numBonds=0;
01134 numAngles=0;
01135 numDihedrals=0;
01136 numImpropers=0;
01137 numCrossterms=0;
01138 numTotalExclusions=0;
01139 for(int i=0; i<numAtoms; i++){
01140 AtomSignature *thisSig = &atomSigPool[eachAtomSig[i]];
01141 numBonds += thisSig->bondCnt;
01142 numAngles += thisSig->angleCnt;
01143 numDihedrals += thisSig->dihedralCnt;
01144 numImpropers += thisSig->improperCnt;
01145 numCrossterms += thisSig->crosstermCnt;
01146
01147 ExclusionSignature *exclSig = &exclSigPool[eachAtomExclSig[i]];
01148 numTotalExclusions += (exclSig->fullExclCnt + exclSig->modExclCnt);
01149 }
01150
01151 numTotalExclusions /= 2;
01152
01153
01154 if(simParams->extraBondsOn){
01155 int numExtraParams=0;
01156
01157
01158 NAMD_read_line(psf_file, buffer);
01159 if(!NAMD_find_word(buffer, "NEXTRABONDPARAMS")){
01160 NAMD_die("UNABLE TO FIND NEXTRABONDPARAMS");
01161 }
01162 sscanf(buffer, "%d", &numExtraParams);
01163 if(numExtraParams>0){
01164 BondValue *newParams = new BondValue[params->NumBondParams+numExtraParams];
01165 memcpy(newParams, params->bond_array, params->NumBondParams*sizeof(BondValue));
01166 delete [] params->bond_array;
01167
01168 int curNumPs = params->NumBondParams;
01169 for(int i=0; i<numExtraParams; i++){
01170 Real k, x0;
01171 NAMD_read_line(psf_file, buffer);
01172 sscanf(buffer, "%f %f", &k, &x0);
01173 newParams[curNumPs+i].k = k;
01174 newParams[curNumPs+i].x0 = x0;
01175 }
01176 params->bond_array = newParams;
01177 params->NumBondParams += numExtraParams;
01178 }
01179
01180
01181 NAMD_read_line(psf_file, buffer);
01182 if(!NAMD_find_word(buffer, "NEXTRAANGLEPARAMS")){
01183 NAMD_die("UNABLE TO FIND NEXTRAANGLEPARAMS");
01184 }
01185 sscanf(buffer, "%d", &numExtraParams);
01186 if(numExtraParams>0){
01187 NAMD_die("CURRENTLY NOT SUPPORT EXTRA ANGLES");
01188 }
01189
01190
01191 NAMD_read_line(psf_file, buffer);
01192 if(!NAMD_find_word(buffer, "NEXTRADIHEDRALPARAMS")){
01193 NAMD_die("UNABLE TO FIND NEXTRADIHEDRALPARAMS");
01194 }
01195 sscanf(buffer, "%d", &numExtraParams);
01196 if(numExtraParams>0){
01197 NAMD_die("CURRENTLY NOT SUPPORT EXTRA DIHEDRALS");
01198 }
01199
01200
01201 NAMD_read_line(psf_file, buffer);
01202 if(!NAMD_find_word(buffer, "NEXTRAIMPROPERPARAMS")){
01203 NAMD_die("UNABLE TO FIND NEXTRAIMPROPERPARAMS");
01204 }
01205 sscanf(buffer, "%d", &numExtraParams);
01206 if(numExtraParams>0){
01207 NAMD_die("CURRENTLY NOT SUPPORT EXTRA IMPROPERS");
01208 }
01209 }
01210
01211
01212 NAMD_read_line(psf_file, buffer);
01213 if(!NAMD_find_word(buffer, "DIHEDRALPARAMARRAY"))
01214 NAMD_die("UNABLE TO FIND DIHEDRALPARAMARRAY");
01215 for(int i=0; i<params->NumDihedralParams; i++){
01216 params->dihedral_array[i].multiplicity = NAMD_read_int(psf_file, buffer);
01217 }
01218
01219 NAMD_read_line(psf_file, buffer);
01220 NAMD_read_line(psf_file, buffer);
01221 if(!NAMD_find_word(buffer, "IMPROPERPARAMARRAY"))
01222 NAMD_die("UNABLE TO FIND IMPROPERPARAMARRAY");
01223 for(int i=0; i<params->NumImproperParams; i++){
01224 params->improper_array[i].multiplicity = NAMD_read_int(psf_file, buffer);
01225 }
01226
01227 Fclose(psf_file);
01228
01229
01230
01231 build_atom_status();
01232 #endif
01233 }
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253 void Molecule::read_atoms(FILE *fd, Parameters *params)
01254
01255 {
01256 #ifdef MEM_OPT_VERSION
01257 return;
01258 #else
01259 char buffer[512];
01260 int atom_number=0;
01261 int last_atom_number=0;
01262
01263 char segment_name[11];
01264 char residue_number[11];
01265 char residue_name[11];
01266 char atom_name[11];
01267 char atom_type[11];
01268 Real charge;
01269 Real mass;
01270 int read_count;
01271
01272
01273 atoms = new Atom[numAtoms];
01274 atomNames = new AtomNameInfo[numAtoms];
01275 hydrogenGroup.resize(0);
01276
01277 if (atoms == NULL || atomNames == NULL )
01278 {
01279 NAMD_die("memory allocation failed in Molecule::read_atoms");
01280 }
01281
01282 ResidueLookupElem *tmpResLookup = resLookup;
01283
01284
01285 while (atom_number < numAtoms)
01286 {
01287
01288 NAMD_read_line(fd, buffer);
01289
01290
01291 if ( (NAMD_blank_string(buffer)) || (buffer[0] == '!') )
01292 continue;
01293
01294
01295 read_count=sscanf(buffer, "%d %s %s %s %s %s %f %f",
01296 &atom_number, segment_name, residue_number,
01297 residue_name, atom_name, atom_type, &charge, &mass);
01298
01299
01300 if (read_count != 8)
01301 {
01302 char err_msg[128];
01303
01304 sprintf(err_msg, "BAD ATOM LINE FORMAT IN PSF FILE IN ATOM LINE %d\nLINE=%s",
01305 last_atom_number+1, buffer);
01306 NAMD_die(err_msg);
01307 }
01308
01309
01310 int atom_type_num;
01311 if ( sscanf(atom_type, "%d", &atom_type_num) > 0 )
01312 {
01313 NAMD_die("Structure (psf) file is either in CHARMM format (with numbers for atoms types, the X-PLOR format using names is required) or the segment name field is empty.");
01314 }
01315
01316
01317 if (atom_number != last_atom_number+1)
01318 {
01319 char err_msg[128];
01320
01321 sprintf(err_msg, "ATOM NUMBERS OUT OF ORDER AT ATOM #%d OF PSF FILE",
01322 last_atom_number+1);
01323 NAMD_die(err_msg);
01324 }
01325
01326 last_atom_number++;
01327
01328
01329
01330
01331 int reslength = strlen(residue_name)+1;
01332 int namelength = strlen(atom_name)+1;
01333 int typelength = strlen(atom_type)+1;
01334
01335 atomNames[atom_number-1].resname = nameArena->getNewArray(reslength);
01336 atomNames[atom_number-1].atomname = nameArena->getNewArray(namelength);
01337 atomNames[atom_number-1].atomtype = nameArena->getNewArray(typelength);
01338
01339 if (atomNames[atom_number-1].resname == NULL)
01340 {
01341 NAMD_die("memory allocation failed in Molecule::read_atoms");
01342 }
01343
01344
01345 strcpy(atomNames[atom_number-1].resname, residue_name);
01346 strcpy(atomNames[atom_number-1].atomname, atom_name);
01347 strcpy(atomNames[atom_number-1].atomtype, atom_type);
01348 atoms[atom_number-1].mass = mass;
01349 atoms[atom_number-1].charge = charge;
01350 atoms[atom_number-1].status = UnknownAtom;
01351
01352
01353 if ( tmpResLookup ) tmpResLookup =
01354 tmpResLookup->append(segment_name, atoi(residue_number), atom_number-1);
01355
01356
01357 if (atoms[atom_number-1].mass < 0.1) {
01358 atoms[atom_number-1].status |= LonepairAtom;
01359 } else if (atoms[atom_number-1].mass <=3.5) {
01360 atoms[atom_number-1].status |= HydrogenAtom;
01361 } else if ((atomNames[atom_number-1].atomname[0] == 'O') &&
01362 (atoms[atom_number-1].mass >= 14.0) &&
01363 (atoms[atom_number-1].mass <= 18.0)) {
01364 atoms[atom_number-1].status |= OxygenAtom;
01365 }
01366
01367
01368 params->assign_vdw_index(atomNames[atom_number-1].atomtype,
01369 &(atoms[atom_number-1]));
01370 }
01371
01372 return;
01373 #endif
01374 }
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389 void Molecule::read_bonds(FILE *fd, Parameters *params)
01390
01391 {
01392 #ifdef MEM_OPT_VERSION
01393 return;
01394 #else
01395 int atom_nums[2];
01396 char atom1name[11];
01397 char atom2name[11];
01398 register int j;
01399 int num_read=0;
01400 int origNumBonds = numBonds;
01401
01402
01403 bonds=new Bond[numBonds];
01404
01405 if (bonds == NULL)
01406 {
01407 NAMD_die("memory allocations failed in Molecule::read_bonds");
01408 }
01409
01410
01411 while (num_read < numBonds)
01412 {
01413
01414 for (j=0; j<2; j++)
01415 {
01416
01417
01418
01419
01420 atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
01421
01422
01423 if (atom_nums[j] >= numAtoms)
01424 {
01425 char err_msg[128];
01426
01427 sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01428 NAMD_die(err_msg);
01429 }
01430 }
01431
01432
01433
01434
01435 if (strcasecmp(atomNames[atom_nums[0]].atomtype,
01436 atomNames[atom_nums[1]].atomtype) < 0)
01437 {
01438 strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
01439 strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
01440 }
01441 else
01442 {
01443 strcpy(atom2name, atomNames[atom_nums[0]].atomtype);
01444 strcpy(atom1name, atomNames[atom_nums[1]].atomtype);
01445 }
01446
01447
01448
01449 Bond *b = &(bonds[num_read]);
01450 params->assign_bond_index(atom1name, atom2name, b);
01451
01452
01453 b->atom1=atom_nums[0];
01454 b->atom2=atom_nums[1];
01455
01456
01457 Real k, x0;
01458 params->get_bond_params(&k,&x0,b->bond_type);
01459 if ( k == 0. ) --numBonds;
01460 else ++num_read;
01461 }
01462
01463
01464 if ( numBonds != origNumBonds ) {
01465 iout << iWARN << "Ignored " << origNumBonds - numBonds <<
01466 " bonds with zero force constants.\n" << endi;
01467 iout << iWARN <<
01468 "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
01469 }
01470
01471 return;
01472 #endif
01473 }
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492 void Molecule::read_angles(FILE *fd, Parameters *params)
01493
01494 {
01495 #ifdef MEM_OPT_VERSION
01496 return;
01497 #else
01498 int atom_nums[3];
01499 char atom1name[11];
01500 char atom2name[11];
01501 char atom3name[11];
01502 register int j;
01503 int num_read=0;
01504 int origNumAngles = numAngles;
01505
01506 angles=new Angle[numAngles];
01507
01508 if (angles == NULL)
01509 {
01510 NAMD_die("memory allocation failed in Molecule::read_angles");
01511 }
01512
01513
01514 while (num_read < numAngles)
01515 {
01516
01517 for (j=0; j<3; j++)
01518 {
01519
01520
01521
01522
01523 atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
01524
01525
01526
01527 if (atom_nums[j] >= numAtoms)
01528 {
01529 char err_msg[128];
01530
01531 sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01532 NAMD_die(err_msg);
01533 }
01534 }
01535
01536
01537
01538
01539
01540
01541 if (strcasecmp(atomNames[atom_nums[0]].atomtype,
01542 atomNames[atom_nums[2]].atomtype) < 0)
01543 {
01544 strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
01545 strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
01546 strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
01547 }
01548 else
01549 {
01550 strcpy(atom1name, atomNames[atom_nums[2]].atomtype);
01551 strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
01552 strcpy(atom3name, atomNames[atom_nums[0]].atomtype);
01553 }
01554
01555
01556
01557 params->assign_angle_index(atom1name, atom2name,
01558 atom3name, &(angles[num_read]));
01559
01560
01561 angles[num_read].atom1=atom_nums[0];
01562 angles[num_read].atom2=atom_nums[1];
01563 angles[num_read].atom3=atom_nums[2];
01564
01565
01566 Real k, t0, k_ub, r_ub;
01567 params->get_angle_params(&k,&t0,&k_ub,&r_ub,angles[num_read].angle_type);
01568 if ( k == 0. && k_ub == 0. ) --numAngles;
01569 else ++num_read;
01570 }
01571
01572
01573 if ( numAngles != origNumAngles ) {
01574 iout << iWARN << "Ignored " << origNumAngles - numAngles <<
01575 " angles with zero force constants.\n" << endi;
01576 }
01577
01578 return;
01579 #endif
01580 }
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598 void Molecule::read_dihedrals(FILE *fd, Parameters *params)
01599 {
01600 #ifdef MEM_OPT_VERSION
01601 return;
01602 #else
01603 int atom_nums[4];
01604 int last_atom_nums[4];
01605 char atom1name[11];
01606 char atom2name[11];
01607 char atom3name[11];
01608 char atom4name[11];
01609 register int j;
01610 int num_read=0;
01611 int multiplicity=1;
01612 Bool duplicate_bond;
01613 int num_unique=0;
01614
01615
01616 for (j=0; j<4; j++)
01617 last_atom_nums[j] = -1;
01618
01619
01620 dihedrals = new Dihedral[numDihedrals];
01621
01622 if (dihedrals == NULL)
01623 {
01624 NAMD_die("memory allocation failed in Molecule::read_dihedrals");
01625 }
01626
01627
01628 while (num_read < numDihedrals)
01629 {
01630 duplicate_bond = TRUE;
01631
01632
01633 for (j=0; j<4; j++)
01634 {
01635
01636
01637
01638
01639 atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
01640
01641
01642 if (atom_nums[j] >= numAtoms)
01643 {
01644 char err_msg[128];
01645
01646 sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01647 NAMD_die(err_msg);
01648 }
01649
01650
01651 if (atom_nums[j] != last_atom_nums[j])
01652 {
01653 duplicate_bond = FALSE;
01654 }
01655
01656 last_atom_nums[j] = atom_nums[j];
01657 }
01658
01659
01660
01661 strcpy(atom1name, atomNames[atom_nums[0]].atomtype);
01662 strcpy(atom2name, atomNames[atom_nums[1]].atomtype);
01663 strcpy(atom3name, atomNames[atom_nums[2]].atomtype);
01664 strcpy(atom4name, atomNames[atom_nums[3]].atomtype);
01665
01666
01667
01668 if (duplicate_bond)
01669 {
01670
01671 multiplicity++;
01672
01673 if (multiplicity == 2)
01674 {
01675 numMultipleDihedrals++;
01676 }
01677 }
01678 else
01679 {
01680 multiplicity=1;
01681 num_unique++;
01682 }
01683
01684
01685 params->assign_dihedral_index(atom1name, atom2name,
01686 atom3name, atom4name, &(dihedrals[num_unique-1]),
01687 multiplicity);
01688
01689
01690 dihedrals[num_unique-1].atom1=atom_nums[0];
01691 dihedrals[num_unique-1].atom2=atom_nums[1];
01692 dihedrals[num_unique-1].atom3=atom_nums[2];
01693 dihedrals[num_unique-1].atom4=atom_nums[3];
01694
01695 num_read++;
01696 }
01697
01698 numDihedrals = num_unique;
01699
01700 return;
01701 #endif
01702 }
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720 void Molecule::read_impropers(FILE *fd, Parameters *params)
01721 {
01722 #ifdef MEM_OPT_VERSION
01723 return;
01724 #else
01725 int atom_nums[4];
01726 int last_atom_nums[4];
01727 char atom1name[11];
01728 char atom2name[11];
01729 char atom3name[11];
01730 char atom4name[11];
01731 register int j;
01732 int num_read=0;
01733 int multiplicity=1;
01734 Bool duplicate_bond;
01735 int num_unique=0;
01736
01737
01738
01739 for (j=0; j<4; j++)
01740 last_atom_nums[j] = -1;
01741
01742
01743 impropers=new Improper[numImpropers];
01744
01745 if (impropers == NULL)
01746 {
01747 NAMD_die("memory allocation failed in Molecule::read_impropers");
01748 }
01749
01750