#include <algorithm>#include "CompressPsf.h"#include "strlib.h"#include "Molecule.h"#include "Parameters.h"#include "SimParameters.h"#include "InfoStream.h"#include "UniqueSet.h"#include "UniqueSetIter.h"Go to the source code of this file.
|
|
Referenced by Molecule::build_extra_bonds(), and getExtraBonds(). |
|
||||||||||||
|
Definition at line 2212 of file CompressPsf.C. References UniqueSet< Elem >::add(), g_mol, and Molecule::numAtoms. Referenced by buildExclusions(). 02213 {
02214 for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
02215 {
02216 vector<int> *atom1List = &eachAtomNeighbors[atom1];
02217 for(int j=0; j<atom1List->size(); j++)
02218 {
02219 int atom2 = atom1List->at(j);
02220 if(atom1<atom2)
02221 allExcls.add(Exclusion(atom1, atom2));
02222 else
02223 allExcls.add(Exclusion(atom2, atom1));
02224 }
02225 }
02226 }
|
|
||||||||||||
|
Definition at line 2228 of file CompressPsf.C. References UniqueSet< Elem >::add(), g_mol, and Molecule::numAtoms. Referenced by buildExclusions(). 02229 {
02230 for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
02231 {
02232 vector<int> *atom1List = &eachAtomNeighbors[atom1];
02233 for(int j=0; j<atom1List->size(); j++)
02234 {
02235 int atom2 = atom1List->at(j);
02236 vector<int> *atom2List = &eachAtomNeighbors[atom2];
02237 for(int k=0; k<atom2List->size(); k++)
02238 {
02239 int atom3 = atom2List->at(k);
02240 //atom1-atom2, so atom2List contains atom1 which should not be considered
02241 if(atom3 == atom1)
02242 continue;
02243 if(atom1<atom3)
02244 allExcls.add(Exclusion(atom1, atom3));
02245 else
02246 allExcls.add(Exclusion(atom3, atom1));
02247 }
02248 }
02249 }
02250 }
|
|
||||||||||||||||
|
Definition at line 2252 of file CompressPsf.C. References UniqueSet< Elem >::add(), g_mol, and Molecule::numAtoms. Referenced by buildExclusions(). 02253 {
02254 for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
02255 {
02256 vector<int> *atom1List = &eachAtomNeighbors[atom1];
02257 for(int j=0; j<atom1List->size(); j++)
02258 {
02259 int atom2 = atom1List->at(j);
02260 vector<int> *atom2List = &eachAtomNeighbors[atom2];
02261 for(int k=0; k<atom2List->size(); k++)
02262 {
02263 int atom3 = atom2List->at(k);
02264 //atom1-atom2, so atom2List contains atom1 which should not be considered
02265 if(atom3 == atom1)
02266 continue;
02267 vector<int> *atom3List = &eachAtomNeighbors[atom3];
02268 for(int l=0; l<atom3List->size(); l++)
02269 {
02270 int atom4 = atom3List->at(l);
02271 //atom1-atom2, so atom2List contains atom1 which should not be considered
02272 if(atom4 == atom2)
02273 continue;
02274 if(atom1<atom4)
02275 allExcls.add(Exclusion(atom1, atom4, modified));
02276 else
02277 allExcls.add(Exclusion(atom4, atom1, modified));
02278 }
02279 }
02280 }
02281 }
02282 }
|
|
|
Definition at line 2096 of file CompressPsf.C. References SimParameters::amberOn, atomData, BasicAtomInfo::atomSigIdx, atomSigPool, UniqueSetIter< T >::begin(), AtomSigInfo::bondSigIndices, build12Excls(), build13Excls(), build14Excls(), UniqueSet< Elem >::clear(), eachAtomExclSigs, UniqueSetIter< T >::end(), BasicAtomInfo::exclSigIdx, SimParameters::exclude, ExclSigInfo::fullExclOffset, g_mol, g_simParam, TupleSignature::isReal, lookupCstPool(), ExclSigInfo::modExclOffset, NONE, Molecule::numAtoms, TupleSignature::offset, ONEFOUR, ONETHREE, ONETWO, SimParameters::readExclusions, SCALED14, sigsOfBonds, sigsOfExclusions, and ExclSigInfo::sortExclOffset(). Referenced by compress_psf_file(). 02097 {
02098 //1. Build exclusions: mainly accomplish the function of
02099 //Molecule::build_exclusions (based on the bonds)
02100 UniqueSet<Exclusion> allExclusions;
02101
02102 int exclude_flag; //Exclusion policy
02103 exclude_flag = g_simParam->exclude;
02104 //int stripHGroupExclFlag = (simParams->splitPatch == SPLIT_PATCH_HYDROGEN);
02105
02106 //Commented now since no explicit exclusions are read
02107 // Go through the explicit exclusions and add them to the arrays
02108 //for(i=0; i<numExclusions; i++){
02109 // exclusionSet.add(exclusions[i]);
02110 //}
02111
02112 // If this is AMBER force field, and readExclusions is TRUE,
02113 // then all the exclusions were read from parm file, and we
02114 // shouldn't generate any of them.
02115 // Comment on stripHGroupExcl:
02116 // 1. Inside this function, hydrogenGroup is initialized in
02117 // build_atom_status, therefore, not available when reading psf files
02118 // 2. this function's main purpose is to reduce memory usage. Since exclusion
02119 // signatures are used, this function could be overlooked --Chao Mei
02120
02121 vector<int> *eachAtomNeighbors = new vector<int>[g_mol->numAtoms];
02122 for(int atom1=0; atom1<g_mol->numAtoms; atom1++)
02123 {
02124 AtomSigInfo *aSig = &atomSigPool[atomData[atom1].atomSigIdx];
02125 for(int j=0; j<aSig->bondSigIndices.size(); j++)
02126 {
02127 TupleSignature *tSig = &sigsOfBonds[aSig->bondSigIndices[j]];
02128 if(!tSig->isReal) continue;
02129 int atom2 = atom1+tSig->offset[0];
02130 eachAtomNeighbors[atom1].push_back(atom2);
02131 eachAtomNeighbors[atom2].push_back(atom1);
02132 }
02133 }
02134
02135 if (!g_simParam->amberOn || !g_simParam->readExclusions)
02136 { // Now calculate the bonded exlcusions based on the exclusion policy
02137 switch (exclude_flag)
02138 {
02139 case NONE:
02140 break;
02141 case ONETWO:
02142 build12Excls(allExclusions, eachAtomNeighbors);
02143 break;
02144 case ONETHREE:
02145 build12Excls(allExclusions, eachAtomNeighbors);
02146 build13Excls(allExclusions, eachAtomNeighbors);
02147 //if ( stripHGroupExclFlag ) stripHGroupExcl();
02148 break;
02149 case ONEFOUR:
02150 build12Excls(allExclusions, eachAtomNeighbors);
02151 build13Excls(allExclusions, eachAtomNeighbors);
02152 build14Excls(allExclusions, eachAtomNeighbors, 0);
02153 //if ( stripHGroupExclFlag ) stripHGroupExcl();
02154 break;
02155 case SCALED14:
02156 build12Excls(allExclusions, eachAtomNeighbors);
02157 build13Excls(allExclusions, eachAtomNeighbors);
02158 build14Excls(allExclusions, eachAtomNeighbors, 1);
02159 //if ( stripHGroupExclFlag ) stripHGroupExcl();
02160 break;
02161 }
02162 }
02163 //else if (stripHGroupExclFlag && exclude_flag!=NONE && exclude_flag!=ONETWO)
02164 // stripHGroupExcl();
02165
02166 //Commented since atomFepFlags information is not available when reading psf file
02167 //stripFepExcl();
02168
02169 for(int i=0; i<g_mol->numAtoms; i++)
02170 eachAtomNeighbors[i].clear();
02171 delete [] eachAtomNeighbors;
02172
02173 //2. Build each atom's list of exclusions
02174 UniqueSetIter<Exclusion> exclIter(allExclusions);
02175 eachAtomExclSigs = new ExclSigInfo[g_mol->numAtoms];
02176 for(exclIter=exclIter.begin(); exclIter!=exclIter.end(); exclIter++)
02177 {
02178 int atom1 = exclIter->atom1;
02179 int atom2 = exclIter->atom2;
02180 int offset21 = atom2-atom1;
02181 if(exclIter->modified)
02182 {
02183 eachAtomExclSigs[atom1].modExclOffset.push_back(offset21);
02184 eachAtomExclSigs[atom2].modExclOffset.push_back(-offset21);
02185 }
02186 else
02187 {
02188 eachAtomExclSigs[atom1].fullExclOffset.push_back(offset21);
02189 eachAtomExclSigs[atom2].fullExclOffset.push_back(-offset21);
02190 }
02191 }
02192 allExclusions.clear();
02193
02194 //3. Build up exclusion signatures and determine each atom's
02195 //exclusion signature index
02196 for(int i=0; i<g_mol->numAtoms; i++)
02197 {
02198 eachAtomExclSigs[i].sortExclOffset();
02199 int poolIndex = lookupCstPool(sigsOfExclusions, eachAtomExclSigs[i]);
02200 if(poolIndex==-1)
02201 {
02202 poolIndex = sigsOfExclusions.size();
02203 sigsOfExclusions.push_back(eachAtomExclSigs[i]);
02204 }
02205 atomData[i].exclSigIdx = poolIndex;
02206 }
02207 delete [] eachAtomExclSigs;
02208 eachAtomExclSigs = NULL;
02209 printf("Exclusion signatures: %d\n", (int)sigsOfExclusions.size());
02210 }
|
|
|
Definition at line 302 of file CompressPsf.C. References atomNamePool, atomTypePool, chargePool, eachClusterSize, massPool, resNamePool, segNamePool, sigsOfAngles, sigsOfBonds, sigsOfDihedrals, sigsOfExclusions, and sigsOfImpropers. 00303 {
00304 segNamePool.clear();
00305 resNamePool.clear();
00306 atomNamePool.clear();
00307 atomTypePool.clear();
00308 chargePool.clear();
00309 massPool.clear();
00310 sigsOfBonds.clear();
00311 sigsOfAngles.clear();
00312 sigsOfDihedrals.clear();
00313 sigsOfImpropers.clear();
00314 sigsOfExclusions.clear();
00315
00316 eachClusterSize.clear();
00317 }
|
|
||||||||||||||||||||||||
|
Definition at line 319 of file CompressPsf.C. References buildExclusions(), g_cfgList, g_mol, g_param, g_simParam, integrateAllAtomSigs(), outputPsfFile(), and readPsfFile(). Referenced by Molecule::Molecule(). 00320 {
00321
00322 g_mol = mol;
00323 g_param = param;
00324 g_simParam = simParam; //used for building exclusions
00325 g_cfgList = cfgList; //used for integrating extra bonds
00326
00327 //read psf files
00328 readPsfFile(psfFileName);
00329
00330 integrateAllAtomSigs();
00331
00332 buildExclusions();
00333
00334
00335 char *outFileName = new char[strlen(psfFileName)+10];
00336 sprintf(outFileName, "%s.inter", psfFileName);
00337 FILE *ofp = fopen(outFileName, "w");
00338 delete [] outFileName;
00339 //output compressed psf file
00340 outputPsfFile(ofp);
00341 fclose(ofp);
00342 }
|
|
|
Definition at line 1900 of file CompressPsf.C. References bond::atom1, bond::atom2, Bond, g_mol, NAMD_die(), NAMD_read_int(), Molecule::numAcceptors, and Molecule::numAtoms. Referenced by readPsfFile(). 01901 {
01902 int d[2]; // temporary storage of atom index
01903 register int j; // Loop counter
01904 int num_read=0; // Number of bonds read so far
01905 int num_no_ante=0; // number of pairs with no antecedent
01906
01907 int numAcceptors = g_mol->numAcceptors;
01908
01909 /* Allocate the array to hold the bonds */
01910 Bond *acceptors=new Bond[numAcceptors];
01911
01912 if (acceptors == NULL)
01913 {
01914 NAMD_die("memory allocations failed in Molecule::read_acceptors");
01915 }
01916
01917 /* Loop through and read in all the acceptors */
01918 while (num_read < numAcceptors)
01919 {
01920 /* Loop and read in the two atom indexes */
01921 for (j=0; j<2; j++)
01922 {
01923 /* Read the atom number from the file. */
01924 /* Subtract 1 to convert the index from the */
01925 /* 1 to NumAtoms used in the file to the */
01926 /* 0 to NumAtoms-1 that we need */
01927 d[j]=NAMD_read_int(fd, "ACCEPTORS")-1;
01928
01929 /* Check to make sure the index isn't too big */
01930 if (d[j] >= g_mol->numAtoms)
01931 {
01932 char err_msg[128];
01933
01934 sprintf(err_msg, "ACCEPTOR INDEX %d GREATER THAN NATOM %d IN DONOR # %d IN PSF FILE", d[j]+1, g_mol->numAtoms, num_read+1);
01935 NAMD_die(err_msg);
01936 }
01937
01938 /* Check if there is an antecedent given */
01939 if (d[j] < 0)
01940 num_no_ante++;
01941 }
01942
01943 /* Assign the atom indexes to the array element */
01944 Bond *b = &(acceptors[num_read]);
01945 b->atom1=d[0];
01946 b->atom2=d[1];
01947
01948 num_read++;
01949 }
01950
01951 delete [] acceptors;
01952 }
|
|
|
Definition at line 1437 of file CompressPsf.C. References ANGLE, Angle, angle::angle_type, AtomSigInfo::angleSigIndices, Parameters::assign_angle_index(), angle::atom1, angle::atom2, angle::atom3, atomData, BasicAtomInfo::atomTypeIdx, atomTypePool, eachAtomSigs, extraAngles, g_mol, g_param, Parameters::get_angle_params(), iout, iWARN(), lookupCstPool(), NAMD_die(), NAMD_read_int(), Molecule::numAngles, Molecule::numAtoms, Real, and sigsOfAngles. Referenced by readPsfFile(). 01438 {
01439
01440 int atom_nums[3]; // Atom numbers for the three atoms
01441 char atom1name[11]; // Atom type for atom 1
01442 char atom2name[11]; // Atom type for atom 2
01443 char atom3name[11]; // Atom type for atom 3
01444 register int j; // Loop counter
01445 int num_read=0; // Number of angles read so far
01446
01447 int numAngles = g_mol->numAngles;
01448 int origNumAngles = numAngles; // Number of angles in file
01449 /* Alloc the array of angles */
01450 Angle *angles=new Angle[numAngles + extraAngles.size()];
01451
01452 if (angles == NULL)
01453 {
01454 NAMD_die("memory allocation failed in Molecule::read_angles");
01455 }
01456
01457 /* Loop through and read all the angles */
01458 while (num_read < numAngles)
01459 {
01460 /* Loop through the 3 atom indexes in the current angle*/
01461 for (j=0; j<3; j++)
01462 {
01463 /* Read the atom number from the file. */
01464 /* Subtract 1 to convert the index from the */
01465 /* 1 to NumAtoms used in the file to the */
01466 /* 0 to NumAtoms-1 that we need */
01467 atom_nums[j]=NAMD_read_int(fd, "ANGLES")-1;
01468
01469 /* Check to make sure the atom index doesn't */
01470 /* exceed the Number of Atoms */
01471 if (atom_nums[j] >= g_mol->numAtoms)
01472 {
01473 char err_msg[128];
01474
01475 sprintf(err_msg, "ANGLES INDEX %d GREATER THAN NATOM %d IN ANGLES # %d IN PSF FILE", atom_nums[j]+1, g_mol->numAtoms, num_read+1);
01476 NAMD_die(err_msg);
01477 }
01478 }
01479
01480 /* Place the bond name that is alphabetically first */
01481 /* in the atom1name. This is OK since the order of */
01482 /* atom1 and atom3 are interchangable. And to search */
01483 /* the tree of angle parameters, we need the order */
01484 /* to be predictable. */
01485
01486 const char *atom1Type = atomTypePool[atomData[atom_nums[0]].atomTypeIdx].c_str();
01487 const char *atom2Type = atomTypePool[atomData[atom_nums[1]].atomTypeIdx].c_str();
01488 const char *atom3Type = atomTypePool[atomData[atom_nums[2]].atomTypeIdx].c_str();
01489
01490 if (strcasecmp(atom1Type,atom2Type) < 0)
01491 {
01492 strcpy(atom1name, atom1Type);
01493 strcpy(atom2name, atom2Type);
01494 strcpy(atom3name, atom3Type);
01495 }
01496 else
01497 {
01498 strcpy(atom1name, atom1Type);
01499 strcpy(atom2name, atom2Type);
01500 strcpy(atom3name, atom3Type);
01501 }
01502
01503 /* Get the constant values for this bond from the */
01504 /* parameter object */
01505 g_param->assign_angle_index(atom1name, atom2name,
01506 atom3name, &(angles[num_read]));
01507
01508 /* Assign the three atom indices */
01509 angles[num_read].atom1=atom_nums[0];
01510 angles[num_read].atom2=atom_nums[1];
01511 angles[num_read].atom3=atom_nums[2];
01512
01513 /* Make sure this isn't a fake angle meant for shake in x-plor. */
01514 Real k, t0, k_ub, r_ub;
01515 g_param->get_angle_params(&k,&t0,&k_ub,&r_ub,angles[num_read].angle_type);
01516 if ( k == 0. && k_ub == 0. )
01517 --numAngles; // fake angle
01518 else
01519 ++num_read; // real angle
01520 }
01521
01522 /* Tell user about our subterfuge */
01523 if ( numAngles != origNumAngles )
01524 {
01525 iout << iWARN << "Ignored " << origNumAngles - numAngles <<
01526 " angles with zero force constants.\n" << endi;
01527 }
01528
01529 //copy extra angles to angles structure
01530 for(int i=0; i<extraAngles.size(); i++)
01531 angles[numAngles+i] = extraAngles[i];
01532 numAngles += extraAngles.size();
01533 extraAngles.clear();
01534
01535 g_mol->numAngles = numAngles;
01536
01537 //create angles' tupleSignature
01538 for(int i=0; i<numAngles; i++)
01539 {
01540 Angle *tuple = angles+i;
01541 TupleSignature oneSig(2,ANGLE,tuple->angle_type);
01542 int offset[2];
01543 offset[0] = tuple->atom2 - tuple->atom1;
01544 offset[1] = tuple->atom3 - tuple->atom1;
01545 oneSig.setOffsets(offset);
01546
01547 int poolIndex = lookupCstPool(sigsOfAngles, oneSig);
01548 if(poolIndex == -1)
01549 {
01550 sigsOfAngles.push_back(oneSig);
01551 poolIndex = (short)sigsOfAngles.size()-1;
01552 }
01553 eachAtomSigs[tuple->atom1].angleSigIndices.push_back(poolIndex);
01554 }
01555
01556 delete [] angles;
01557
01558 }
|
|
|
Definition at line 870 of file CompressPsf.C. References atomData, BasicAtomInfo::atomNameIdx, atomNamePool, BasicAtomInfo::atomTypeIdx, atomTypePool, charge, BasicAtomInfo::chargeIdx, chargePool, g_mol, lookupCstPool(), BasicAtomInfo::massIdx, massPool, NAMD_blank_string(), NAMD_die(), NAMD_read_line(), Molecule::numAtoms, Real, BasicAtomInfo::resID, BasicAtomInfo::resNameIdx, resNamePool, BasicAtomInfo::segNameIdx, and segNamePool. Referenced by readPsfFile(). 00871 {
00872 char buffer[512]; //Buffer for reading from file
00873 int atomID=0;
00874 int lastAtomID=0; //to ensure atoms are in order
00875 char segmentName[11];
00876 int residueID;
00877 char residueName[11];
00878 char atomName[11];
00879 char atomType[11];
00880 Real charge;
00881 Real mass;
00882 int fieldsCnt; //Number of fields read by sscanf
00883
00884 int numAtoms = g_mol->numAtoms;
00885
00886 //1. parse atom data to build constant pool (atom name, mass, charge etc.)
00887 atomData = new BasicAtomInfo[numAtoms];
00888 while(atomID < numAtoms)
00889 {
00890 NAMD_read_line(ifp, buffer);
00891 //If it is blank or a comment, skip it
00892 if(NAMD_blank_string(buffer) || buffer[0]=='!')
00893 continue;
00894
00895 //Fields are arranged as follows:
00896 //atomID, segName, resID, resName, atomName, atomType, charge, mass
00897 fieldsCnt = sscanf(buffer, "%d %s %d %s %s %s %f %f",
00898 &atomID, segmentName, &residueID, residueName,
00899 atomName, atomType, &charge, &mass);
00900
00901 //check to make sure we found what we were expecting
00902 if(fieldsCnt!=8)
00903 {
00904 char errMsg[128];
00905 sprintf(errMsg, "BAD ATOM LINE FORMAT IN PSF FILE FOR ATOM %d (%s)", lastAtomID+1, buffer);
00906 NAMD_die(errMsg);
00907 }
00908
00909 // check if this is in XPLOR format
00910 int atomTypeNum;
00911 if(sscanf(atomType, "%d", &atomTypeNum)>0)
00912 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.");
00913
00914 //make sure the atoms were in sequence
00915 if(atomID!=lastAtomID+1)
00916 {
00917 char errMsg[128];
00918 sprintf(errMsg, "ATOM NUMBERS OUT OF ORDER AT ATOM #%d IN FSF FILE",
00919 lastAtomID+1);
00920 NAMD_die(errMsg);
00921 }
00922
00923 lastAtomID++;
00924
00925 //building constant pool
00926 int poolIndex;
00927 string fieldName;
00928 fieldName.assign(segmentName);
00929 poolIndex = lookupCstPool(segNamePool, fieldName);
00930 if(poolIndex==-1)
00931 {
00932 segNamePool.push_back(fieldName);
00933 poolIndex = segNamePool.size()-1;
00934 }
00935 atomData[atomID-1].segNameIdx = poolIndex;
00936
00937 //poolIndex = lookupCstPool(resIDPool, residueID);
00938 //if(poolIndex==-1)
00939 // resIDPool.push_back(residueID);
00940 atomData[atomID-1].resID = residueID;
00941
00942 fieldName.assign(residueName);
00943 poolIndex = lookupCstPool(resNamePool, fieldName);
00944 if(poolIndex==-1)
00945 {
00946 resNamePool.push_back(fieldName);
00947 poolIndex = resNamePool.size()-1;
00948 }
00949 atomData[atomID-1].resNameIdx = poolIndex;
00950
00951 fieldName.assign(atomName);
00952 poolIndex = lookupCstPool(atomNamePool, fieldName);
00953 if(poolIndex==-1)
00954 {
00955 atomNamePool.push_back(fieldName);
00956 poolIndex = atomNamePool.size()-1;
00957 }
00958 atomData[atomID-1].atomNameIdx = poolIndex;
00959
00960 fieldName.assign(atomType);
00961 poolIndex = lookupCstPool(atomTypePool, fieldName);
00962 if(poolIndex==-1)
00963 {
00964 atomTypePool.push_back(fieldName);
00965 poolIndex = atomTypePool.size()-1;
00966 }
00967 atomData[atomID-1].atomTypeIdx = poolIndex;
00968
00969 poolIndex = lookupCstPool(chargePool, charge);
00970 if(poolIndex==-1)
00971 {
00972 chargePool.push_back(charge);
00973 poolIndex = chargePool.size()-1;
00974 }
00975 atomData[atomID-1].chargeIdx = poolIndex;
00976
00977 poolIndex = lookupCstPool(massPool, mass);
00978 if(poolIndex==-1)
00979 {
00980 massPool.push_back(mass);
00981 poolIndex = massPool.size()-1;
00982 }
00983 atomData[atomID-1].massIdx = poolIndex;
00984 }
00985 }
|
|
|
Definition at line 1182 of file CompressPsf.C. References Parameters::assign_bond_index(), bond::atom1, bond::atom2, atomData, BasicAtomInfo::atomTypeIdx, atomTypePool, BOND, Bond, bond::bond_type, AtomSigInfo::bondSigIndices, eachAtomClusterID, eachAtomSigs, eachClusterSize, extraBonds, g_mol, g_param, Parameters::get_bond_params(), iout, iWARN(), lookupCstPool(), NAMD_die(), NAMD_read_int(), Molecule::numAtoms, Molecule::numBonds, TupleSignature::offset, Real, and sigsOfBonds. Referenced by readPsfFile(). 01183 {
01184
01185 //first read bonds
01186 int atom_nums[2]; // Atom indexes for the bonded atoms
01187 char atom1name[11]; // Atom type for atom #1
01188 char atom2name[11]; // Atom type for atom #2
01189 register int j; // Loop counter
01190 int num_read=0; // Number of bonds read so far
01191
01192 int numBonds = g_mol->numBonds;
01193 int origNumBonds = numBonds; // number of bonds in file header
01194
01195 /* Allocate the array to hold the bonds */
01196 Bond *bonds=new Bond[numBonds + extraBonds.size()];
01197
01198 if (bonds == NULL)
01199 {
01200 NAMD_die("memory allocations failed in Molecule::read_bonds");
01201 }
01202
01203 /* Loop through and read in all the bonds */
01204 while (num_read < numBonds)
01205 {
01206 /* Loop and read in the two atom indexes */
01207 for (j=0; j<2; j++)
01208 {
01209 /* Read the atom number from the file. */
01210 /* Subtract 1 to convert the index from the */
01211 /* 1 to NumAtoms used in the file to the */
01212 /* 0 to NumAtoms-1 that we need */
01213 atom_nums[j]=NAMD_read_int(fd, "BONDS")-1;
01214
01215 /* Check to make sure the index isn't too big */
01216 if (atom_nums[j] >= g_mol->numAtoms)
01217 {
01218 char err_msg[128];
01219
01220 sprintf(err_msg, "BOND INDEX %d GREATER THAN NATOM %d IN BOND # %d IN PSF FILE", atom_nums[j]+1, g_mol->numAtoms, num_read+1);
01221 NAMD_die(err_msg);
01222 }
01223 }
01224
01225 if(atom_nums[0] == atom_nums[1])
01226 { //bond to itself
01227 char err_msg[128];
01228 sprintf(err_msg, "ATOM %d is bonded to itself!", atom_nums[0]+1);
01229 NAMD_die(err_msg);
01230 }
01231
01232 /* Get the atom type for the two atoms. When we query */
01233 /* the parameter object, we need to send the atom type */
01234 /* that is alphabetically first as atom 1. */
01235 const char *atom1Type = atomTypePool[atomData[atom_nums[0]].atomTypeIdx].c_str();
01236 const char *atom2Type = atomTypePool[atomData[atom_nums[1]].atomTypeIdx].c_str();
01237 if (strcasecmp(atom1Type,atom2Type) < 0)
01238 {
01239 strcpy(atom1name, atom1Type);
01240 strcpy(atom2name, atom2Type);
01241 }
01242 else
01243 {
01244 strcpy(atom2name, atom1Type);
01245 strcpy(atom1name, atom2Type);
01246 }
01247
01248 /* Query the parameter object for the constants for */
01249 /* this bond */
01250 Bond *b = &(bonds[num_read]);
01251 g_param->assign_bond_index(atom1name, atom2name, b);
01252
01253 /* Assign the atom indexes to the array element */
01254 b->atom1=atom_nums[0];
01255 b->atom2=atom_nums[1];
01256
01257 /* Make sure this isn't a fake bond meant for shake in x-plor. */
01258 Real k, x0;
01259 g_param->get_bond_params(&k,&x0,b->bond_type);
01260 if ( k == 0. )
01261 --numBonds; // fake bond
01262 else
01263 ++num_read; // real bond
01264 }
01265
01266 /* Tell user about our subterfuge */
01267 if ( numBonds != origNumBonds )
01268 {
01269 iout << iWARN << "Ignored " << origNumBonds - numBonds <<
01270 " bonds with zero force constants.\n" << endi;
01271 iout << iWARN <<
01272 "Will get H-H distance in rigid H2O from H-O-H angle.\n" << endi;
01273 }
01274
01275 //copy extra bonds to bonds structure
01276 int numRealBonds = numBonds;
01277 for(int i=0; i<extraBonds.size(); i++)
01278 bonds[numBonds+i] = extraBonds[i];
01279 numBonds += extraBonds.size();
01280 extraBonds.clear();
01281
01282 g_mol->numBonds = numBonds;
01283
01284 //then creating bond's tupleSignature
01285 for(int i=0; i<numBonds; i++)
01286 {
01287 Bond *b = bonds+i;
01288 TupleSignature oneSig(1,BOND,b->bond_type);
01289 oneSig.offset[0] = b->atom2 - b->atom1;
01290 oneSig.isReal = (i<numRealBonds );
01291
01292 int poolIndex = lookupCstPool(sigsOfBonds, oneSig);
01293 int newSig=0;
01294 if(poolIndex == -1)
01295 {
01296 sigsOfBonds.push_back(oneSig);
01297 poolIndex = (short)sigsOfBonds.size()-1;
01298 newSig=1;
01299 }
01300
01301 if(!newSig)
01302 {//check duplicate bonds in the form of (a, b) && (a, b);
01303 int dupIdx = lookupCstPool(eachAtomSigs[b->atom1].bondSigIndices, (short)poolIndex);
01304 if(dupIdx!=-1)
01305 {
01306 char err_msg[128];
01307 sprintf(err_msg, "Duplicate bond %d-%d!", b->atom1+1, b->atom2+1);
01308 NAMD_die(err_msg);
01309 }
01310 }
01311 eachAtomSigs[b->atom1].bondSigIndices.push_back(poolIndex);
01312 }
01313
01314 //check duplicate bonds in the form of (a, b) && (b, a)
01315 for(int i=0; i<numBonds; i++)
01316 {
01317 Bond *b=bonds+i;
01318 int atom2 = b->atom2;
01319 int thisOffset = atom2 - b->atom1;
01320 for(int j=0; j<eachAtomSigs[atom2].bondSigIndices.size(); j++)
01321 {
01322 short atom2BondId = eachAtomSigs[atom2].bondSigIndices[j];
01323 TupleSignature *secSig = &(sigsOfBonds[atom2BondId]);
01324 if(thisOffset== -(secSig->offset[0]))
01325 {
01326 char err_msg[128];
01327 sprintf(err_msg, "Duplicate bond %d-%d because two atoms are just reversed!", b->atom1+1, atom2+1);
01328 NAMD_die(err_msg);
01329 }
01330 }
01331 }
01332
01333 //building clusters for this simulation system in two steps
01334 //1. create a list for each atom where each atom in the list is bonded with that atom
01335 vector<int> *atomListOfBonded = new vector<int>[g_mol->numAtoms];
01336
01337 for(int i=0; i<numRealBonds; i++)
01338 {
01339 Bond *b=bonds+i;
01340 int atom1 = b->atom1;
01341 int atom2 = b->atom2;
01342 atomListOfBonded[atom1].push_back(atom2);
01343 atomListOfBonded[atom2].push_back(atom1);
01344 }
01345
01346 delete [] bonds;
01347
01348 //2. using breadth-first-search to build the clusters. Here, we avoid recursive call
01349 // because the depth of calls may be of thousands which will blow up the stack, and
01350 //recursive call is slower than the stack-based BFS.
01351 //Considering such structure
01352 //1->1245; 7->1243; 1243->1245
01353 eachAtomClusterID = new int[g_mol->numAtoms];
01354 for(int i=0; i<g_mol->numAtoms; i++)
01355 eachAtomClusterID[i] = -1;
01356
01357 for(int i=0; i<g_mol->numAtoms; i++)
01358 {
01359 int curClusterID=eachAtomClusterID[i];
01360 if(curClusterID==-1)
01361 {
01362 curClusterID=i;
01363 }
01364
01365 deque<int> toVisitAtoms;
01366 toVisitAtoms.push_back(i);
01367 while(!toVisitAtoms.empty())
01368 {
01369 int visAtomID = toVisitAtoms.front();
01370 toVisitAtoms.pop_front();
01371 eachAtomClusterID[visAtomID] = curClusterID;
01372 for(int j=0; j<atomListOfBonded[visAtomID].size(); j++)
01373 {
01374 int otherAtom = atomListOfBonded[visAtomID][j];
01375 if(eachAtomClusterID[otherAtom]!=curClusterID)
01376 toVisitAtoms.push_back(otherAtom);
01377 }
01378 }
01379 }
01380
01381 //Now the clusterID of each atom should be in the non-decreasing order, otherwise
01382 //report a psf format error. Well, it's not a real problem, it can be fixed the problem
01383 //later. Here, such assumption is made for the ease of programming
01384 int curClusterID;
01385 int prevClusterID=eachAtomClusterID[0];
01386 int curClusterSize=1;
01387 for(int i=1; i<g_mol->numAtoms; i++)
01388 {
01389 curClusterID = eachAtomClusterID[i];
01390 if(curClusterID > prevClusterID)
01391 {
01392 eachClusterSize.push_back(curClusterSize);
01393 curClusterSize=1;
01394 }
01395 else if(curClusterID == prevClusterID)
01396 {
01397 curClusterSize++;
01398 }
01399 else
01400 { //psf format error
01401 char errMsg[256];
01402 sprintf(errMsg, "Cluster ID of each atom should be in increasing order (error for atom %d)!\n", i+1);
01403 NAMD_die(errMsg);
01404 }
01405 prevClusterID = curClusterID;
01406 }
01407 eachClusterSize.push_back(curClusterSize); //record the last sequence of cluster
01408
01409 //Now iterate over the cluster size again to filter out the repeating cluster size.
01410 //Since the cluster id of atoms is monotonically increasing, the size of the cluster can be used
01411 //as this cluster's signature.
01412 //After this, eachAtomClusterID will store the cluster signature. Only the atom whose id is "clustersize"
01413 //will store the size. Others will be -1
01414
01415 int aid=0;
01416 for(int clusterIdx=0; clusterIdx<eachClusterSize.size(); clusterIdx++)
01417 {
01418 int curSize = eachClusterSize[clusterIdx];
01419 eachAtomClusterID[aid] = curSize;
01420 for(int i=aid+1; i<aid+curSize; i++)
01421 eachAtomClusterID[i] = -1;
01422 aid += curSize;
01423 }
01424
01425 //check whether cluster is built correctly
01426 /*printf("num clusters: %d\n", eachClusterSize.size());
01427 FILE *checkFile = fopen("cluster.opt", "w");
01428 for(int i=0; i<g_mol->numAtoms; i++) fprintf(checkFile, "%d\n", eachAtomClusterID[i]);
01429 fclose(checkFile);*/
01430
01431 eachClusterSize.clear();
01432 for(int i=0; i<g_mol->numAtoms; i++)
01433 atomListOfBonded[i].clear();
01434 delete [] atomListOfBonded;
01435 }
|
|
|
Definition at line 1963 of file CompressPsf.C. References Parameters::assign_crossterm_index(), crossterm::atom1, crossterm::atom2, crossterm::atom3, crossterm::atom4, crossterm::atom5, crossterm::atom6, crossterm::atom7, crossterm::atom8, atomData, BasicAtomInfo::atomTypeIdx, atomTypePool, Bool, CROSSTERM, Crossterm, crossterm::crossterm_type, AtomSigInfo::crosstermSigIndices, eachAtomSigs, g_mol, g_param, iout, iWARN(), lookupCstPool(), NAMD_die(), NAMD_read_int(), Molecule::numAtoms, Molecule::numCrossterms, and sigsOfCrossterms. Referenced by readPsfFile(). 01964 {
01965 int atom_nums[8]; // Atom indexes for the 4 atoms
01966 int last_atom_nums[8]; // Atom indexes from previous bond
01967 char atom1name[11]; // Atom type for atom 1
01968 char atom2name[11]; // Atom type for atom 2
01969 char atom3name[11]; // Atom type for atom 3
01970 char atom4name[11]; // Atom type for atom 4
01971 char atom5name[11]; // Atom type for atom 5
01972 char atom6name[11]; // Atom type for atom 6
01973 char atom7name[11]; // Atom type for atom 7
01974 char atom8name[11]; // Atom type for atom 8
01975 int j; // Loop counter
01976 int num_read=0; // Number of items read so far
01977 Bool duplicate_bond; // Is this a duplicate of the last bond
01978
01979 // Initialize the array used to look for duplicate crossterm
01980 // entries. Set them all to -1 so we know nothing will match
01981 for (j=0; j<8; j++)
01982 last_atom_nums[j] = -1;
01983
01984 int numCrossterms = g_mol->numCrossterms;
01985 int numAtoms = g_mol->numAtoms;
01986
01987 /* Allocate the array to hold the cross-terms */
01988 Crossterm* crossterms=new Crossterm[numCrossterms];
01989
01990 if (crossterms == NULL)
01991 {
01992 NAMD_die("memory allocation failed in getCrosstermData when compressing psf file");
01993 }
01994
01995 /* Loop through and read all the cross-terms */
01996 while (num_read < numCrossterms)
01997 {
01998 duplicate_bond = TRUE;
01999
02000 /* Loop through the 8 indexes for this cross-term */
02001 for (j=0; j<8; j++)
02002 {
02003 /* Read the atom number from the file. */
02004 /* Subtract 1 to convert the index from the */
02005 /* 1 to NumAtoms used in the file to the */
02006 /* 0 to NumAtoms-1 that we need */
02007 atom_nums[j]=NAMD_read_int(fd, "CROSS-TERMS")-1;
02008
02009 /* Check to make sure the index isn't too big */
02010 if (atom_nums[j] >= numAtoms)
02011 {
02012 char err_msg[128];
02013
02014 sprintf(err_msg, "CROSS-TERM INDEX %d GREATER THAN NATOM %d IN CROSS-TERMS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
02015 NAMD_die(err_msg);
02016 }
02017
02018 if (atom_nums[j] != last_atom_nums[j]){
02019 duplicate_bond = FALSE;
02020 }
02021
02022 last_atom_nums[j] = atom_nums[j];
02023 }
02024
02025 /* Get the atom types so we can look up the parameters */
02026 const char *atom1Type = atomTypePool[atomData[atom_nums[0]].atomTypeIdx].c_str();
02027 const char *atom2Type = atomTypePool[atomData[atom_nums[1]].atomTypeIdx].c_str();
02028 const char *atom3Type = atomTypePool[atomData[atom_nums[2]].atomTypeIdx].c_str();
02029 const char *atom4Type = atomTypePool[atomData[atom_nums[3]].atomTypeIdx].c_str();
02030 const char *atom5Type = atomTypePool[atomData[atom_nums[4]].atomTypeIdx].c_str();
02031 const char *atom6Type = atomTypePool[atomData[atom_nums[5]].atomTypeIdx].c_str();
02032 const char *atom7Type = atomTypePool[atomData[atom_nums[6]].atomTypeIdx].c_str();
02033 const char *atom8Type = atomTypePool[atomData[atom_nums[7]].atomTypeIdx].c_str();
02034 strcpy(atom1name, atom1Type);
02035 strcpy(atom2name, atom2Type);
02036 strcpy(atom3name, atom3Type);
02037 strcpy(atom4name, atom4Type);
02038 strcpy(atom5name, atom5Type);
02039 strcpy(atom6name, atom6Type);
02040 strcpy(atom7name, atom7Type);
02041 strcpy(atom8name, atom8Type);
02042
02043 // Check to see if this is a duplicate term
02044 if (duplicate_bond)
02045 {
02046 iout << iWARN << "Duplicate cross-term detected.\n" << endi;
02047 }
02048
02049 /* Look up the constants for this bond */
02050 g_param->assign_crossterm_index(atom1name, atom2name,
02051 atom3name, atom4name, atom5name, atom6name,
02052 atom7name, atom8name, &(crossterms[num_read]));
02053
02054 /* Assign the atom indexes */
02055 crossterms[num_read].atom1=atom_nums[0];
02056 crossterms[num_read].atom2=atom_nums[1];
02057 crossterms[num_read].atom3=atom_nums[2];
02058 crossterms[num_read].atom4=atom_nums[3];
02059 crossterms[num_read].atom5=atom_nums[4];
02060 crossterms[num_read].atom6=atom_nums[5];
02061 crossterms[num_read].atom7=atom_nums[6];
02062 crossterms[num_read].atom8=atom_nums[7];
02063
02064 num_read++;
02065 }
02066
02067 numCrossterms = num_read;
02068
02069 //create crossterm's tupleSignature
02070 for(int i=0; i<numCrossterms; i++)
02071 {
02072 Crossterm *tuple = crossterms+i;
02073 TupleSignature oneSig(7, CROSSTERM, tuple->crossterm_type);
02074 int offset[7];
02075 offset[0] = tuple->atom2 - tuple->atom1;
02076 offset[1] = tuple->atom3 - tuple->atom1;
02077 offset[2] = tuple->atom4 - tuple->atom1;
02078 offset[3] = tuple->atom5 - tuple->atom1;
02079 offset[4] = tuple->atom6 - tuple->atom1;
02080 offset[5] = tuple->atom7 - tuple->atom1;
02081 offset[6] = tuple->atom8 - tuple->atom1;
02082 oneSig.setOffsets(offset);
02083
02084 int poolIndex = lookupCstPool(sigsOfCrossterms, oneSig);
02085 if(poolIndex == -1)
02086 {
02087 sigsOfCrossterms.push_back(oneSig);
02088 poolIndex = (short)sigsOfCrossterms.size()-1;
02089 }
02090 eachAtomSigs[tuple->atom1].crosstermSigIndices.push_back(poolIndex);
02091 }
02092
02093 delete[] crossterms;
02094 }
|
|
|
Definition at line 1560 of file CompressPsf.C. References Parameters::assign_dihedral_index(), dihedral::atom1, dihedral::atom2, dihedral::atom3, dihedral::atom4, atomData, BasicAtomInfo::atomTypeIdx, atomTypePool, Bool, DIHEDRAL, Dihedral, dihedral::dihedral_type, AtomSigInfo::dihedralSigIndices, eachAtomSigs, extraDihedrals, g_mol, g_param, iout, iWARN(), lookupCstPool(), NAMD_die(), NAMD_read_int(), Molecule::numAtoms, Molecule::numDihedrals, Molecule::numMultipleDihedrals, and sigsOfDihedrals. Referenced by readPsfFile(). 01561 {
01562
01563 int atom_nums[4]; // The 4 atom indexes
01564 int last_atom_nums[4]; // Atom numbers from previous bond
01565 char atom1name[11]; // Atom type for atom 1
01566 char atom2name[11]; // Atom type for atom 2
01567 char atom3name[11]; // Atom type for atom 3
01568 char atom4name[11]; // Atom type for atom 4
01569 register int j; // loop counter
01570 int num_read=0; // number of dihedrals read so far
01571 int multiplicity=1; // multiplicity of the current bond
01572 Bool duplicate_bond; // Is this a duplicate of the last bond
01573 int num_unique=0; // Number of unique dihedral bonds
01574
01575 // Initialize the array used to check for duplicate dihedrals
01576 for (j=0; j<4; j++)
01577 last_atom_nums[j] = -1;
01578
01579 int numDihedrals = g_mol->numDihedrals;
01580 int numAtoms = g_mol->numAtoms;
01581
01582 /* Allocate an array to hold the Dihedrals */
01583 Dihedral *dihedrals = new Dihedral[numDihedrals + extraDihedrals.size()];
01584
01585 if (dihedrals == NULL)
01586 {
01587 NAMD_die("memory allocation failed in Molecule::read_dihedrals");
01588 }
01589
01590 /* Loop through and read all the dihedrals */
01591 while (num_read < numDihedrals)
01592 {
01593 duplicate_bond = TRUE;
01594
01595 /* Loop through and read the 4 indexes for this bond */
01596 for (j=0; j<4; j++)
01597 {
01598 /* Read the atom number from the file. */
01599 /* Subtract 1 to convert the index from the */
01600 /* 1 to NumAtoms used in the file to the */
01601 /* 0 to NumAtoms-1 that we need */
01602 atom_nums[j]=NAMD_read_int(fd, "DIHEDRALS")-1;
01603
01604 /* Check for an atom index that is too large */
01605 if (atom_nums[j] >= numAtoms)
01606 {
01607 char err_msg[128];
01608
01609 sprintf(err_msg, "DIHEDRALS INDEX %d GREATER THAN NATOM %d IN DIHEDRALS # %d IN PSF FILE", atom_nums[j]+1, numAtoms, num_read+1);
01610 NAMD_die(err_msg);
01611 }
01612
01613 // Check to see if this atom matches the last bond
01614 if (atom_nums[j] != last_atom_nums[j])
01615 {
01616 duplicate_bond = FALSE;
01617 }
01618
01619 last_atom_nums[j] = atom_nums[j];
01620 }
01621
01622 /* Get the atom types for the 4 atoms so we can look */
01623 /* up the constants in the parameter object */
01624 const char *atom1Type = atomTypePool[atomData[atom_nums[0]].atomTypeIdx].c_str();
01625 const char *atom2Type = atomTypePool[atomData[atom_nums[1]].atomTypeIdx].c_str();
01626 const char *atom3Type = atomTypePool[atomData[atom_nums[2]].atomTypeIdx].c_str();
01627 const char *atom4Type = atomTypePool[atomData[atom_nums[3]].atomTypeIdx].c_str();
01628 strcpy(atom1name, atom1Type);
01629 strcpy(atom2name, atom2Type);
01630 strcpy(atom3name, atom3Type);
01631 strcpy(atom4name, atom4Type);
01632
01633 // Check to see if this is really a new bond or just
01634 // a repeat of the last one
01635 if (duplicate_bond)
01636 {
01637 // This is a duplicate, so increase the multiplicity
01638 multiplicity++;
01639
01640 if (multiplicity == 2)
01641 {
01642 g_mol->numMultipleDihedrals++;
01643 }
01644 }
01645 else
01646 {
01647 multiplicity=1;
01648 num_unique++;
01649 }
01650
01651 /* Get the constants for this dihedral bond */
01652 g_param->assign_dihedral_index(atom1name, atom2name,
01653 atom3name, atom4name, &(dihedrals[num_unique-1]),
01654 multiplicity);
01655
01656 /* Assign the atom indexes */
01657 dihedrals[num_unique-1].atom1=atom_nums[0];
01658 dihedrals[num_unique-1].atom2=atom_nums[1];
01659 dihedrals[num_unique-1].atom3=atom_nums[2];
01660 dihedrals[num_unique-1].atom4=atom_nums[3];
01661
01662 num_read++;
01663 }
01664
01665 //copy extra dihedrals to dihedrals structure
01666 if(extraDihedrals.size()>0){
01667 iout << iWARN << "It's your responsibility to ensure there is no duplicates among these extra dihedrals\n" << endi;
01668 }
01669 for(int i=0; i<extraDihedrals.size(); i++)
01670 dihedrals[num_unique+i] = extraDihedrals[i];
01671 num_unique += extraDihedrals.size();
01672 extraDihedrals.clear();
01673
01674 g_mol->numDihedrals = num_unique;
01675
01676 //create dihedrals' tupleSignature
01677 for(int i=0; i<num_unique; i++)
01678 {
01679 Dihedral *tuple = dihedrals+i;
01680 TupleSignature oneSig(3,DIHEDRAL,tuple->dihedral_type);
01681 int offset[3];
01682 offset[0] = tuple->atom2 - tuple->atom1;
01683 offset[1] = tuple->atom3 - tuple->atom1;
01684 offset[2] = tuple->atom4 - tuple->atom1;
01685 oneSig.setOffsets(offset);
01686
01687 int poolIndex = lookupCstPool(sigsOfDihedrals, oneSig);
01688 if(poolIndex == -1)
01689 {
01690 sigsOfDihedrals.push_back(oneSig);
01691 poolIndex = (short)sigsOfDihedrals.size()-1;
01692 }
01693 eachAtomSigs[tuple->atom1].dihedralSigIndices.push_back(poolIndex);
01694 |